MFEM  v3.1 Finite element discretization library
fe_coll.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
10 // Software Foundation) version 2.1 dated February 1999.
11
12 #include "fem.hpp"
13 #include <cstdlib>
14 #include <cstring>
15 #include <cstdio>
16 #ifdef _WIN32
17 #define snprintf _snprintf_s
18 #endif
19
20 namespace mfem
21 {
22
23 using namespace std;
24
25 int FiniteElementCollection::HasFaceDofs(int GeomType) const
26 {
27  switch (GeomType)
28  {
29  case Geometry::TETRAHEDRON: return DofForGeometry (Geometry::TRIANGLE);
30  case Geometry::CUBE: return DofForGeometry (Geometry::SQUARE);
31  default:
32  mfem_error ("FiniteElementCollection::HasFaceDofs:"
33  " unknown geometry type.");
34  }
35  return 0;
36 }
37
39 {
40  MFEM_ABORT("this method is not implemented in this derived class!");
41  return NULL;
42 }
43
45 {
46  FiniteElementCollection *fec = NULL;
47
48  if (!strcmp(name, "Linear"))
49  {
50  fec = new LinearFECollection;
51  }
53  {
55  }
57  {
59  }
60  else if (!strcmp(name, "Cubic"))
61  {
62  fec = new CubicFECollection;
63  }
64  else if (!strcmp(name, "Const3D"))
65  {
66  fec = new Const3DFECollection;
67  }
68  else if (!strcmp(name, "Const2D"))
69  {
70  fec = new Const2DFECollection;
71  }
72  else if (!strcmp(name, "LinearDiscont2D"))
73  {
75  }
76  else if (!strcmp(name, "GaussLinearDiscont2D"))
77  {
79  }
81  {
83  }
85  {
87  }
89  {
91  }
93  {
95  }
96  else if (!strcmp(name, "CubicDiscont2D"))
97  {
99  }
100  else if (!strcmp(name, "LinearDiscont3D"))
101  {
102  fec = new LinearDiscont3DFECollection;
103  }
105  {
107  }
108  else if (!strcmp(name, "LinearNonConf3D"))
109  {
110  fec = new LinearNonConf3DFECollection;
111  }
112  else if (!strcmp(name, "CrouzeixRaviart"))
113  {
114  fec = new CrouzeixRaviartFECollection;
115  }
116  else if (!strcmp(name, "ND1_3D"))
117  {
118  fec = new ND1_3DFECollection;
119  }
120  else if (!strcmp(name, "RT0_2D"))
121  {
122  fec = new RT0_2DFECollection;
123  }
124  else if (!strcmp(name, "RT1_2D"))
125  {
126  fec = new RT1_2DFECollection;
127  }
128  else if (!strcmp(name, "RT2_2D"))
129  {
130  fec = new RT2_2DFECollection;
131  }
132  else if (!strcmp(name, "RT0_3D"))
133  {
134  fec = new RT0_3DFECollection;
135  }
136  else if (!strcmp(name, "RT1_3D"))
137  {
138  fec = new RT1_3DFECollection;
139  }
140  else if (!strncmp(name, "H1_Trace_", 9))
141  {
142  fec = new H1_Trace_FECollection(atoi(name + 13), atoi(name + 9));
143  }
144  else if (!strncmp(name, "H1_", 3))
145  {
146  fec = new H1_FECollection(atoi(name + 7), atoi(name + 3));
147  }
148  else if (!strncmp(name, "H1Pos_Trace_", 12))
149  {
150  fec = new H1_Trace_FECollection(atoi(name + 16), atoi(name + 12), 1);
151  }
152  else if (!strncmp(name, "H1Pos_", 6))
153  {
154  fec = new H1Pos_FECollection(atoi(name + 10), atoi(name + 6));
155  }
156  else if (!strncmp(name, "L2_T", 4))
157  fec = new L2_FECollection(atoi(name + 10), atoi(name + 6),
158  atoi(name + 4));
159  else if (!strncmp(name, "L2_", 3))
160  {
161  fec = new L2_FECollection(atoi(name + 7), atoi(name + 3));
162  }
163  else if (!strncmp(name, "RT_Trace_", 9))
164  {
165  fec = new RT_Trace_FECollection(atoi(name + 13), atoi(name + 9));
166  }
167  else if (!strncmp(name, "RT_ValTrace_", 12))
168  {
169  fec = new RT_Trace_FECollection(atoi(name + 16), atoi(name + 12),
171  }
172  else if (!strncmp(name, "DG_Iface_", 9))
173  {
174  fec = new DG_Interface_FECollection(atoi(name + 13), atoi(name + 9));
175  }
176  else if (!strncmp(name, "DG_IntIface_", 12))
177  {
178  fec = new DG_Interface_FECollection(atoi(name + 16), atoi(name + 12),
180  }
181  else if (!strncmp(name, "RT_", 3))
182  {
183  fec = new RT_FECollection(atoi(name + 7), atoi(name + 3));
184  }
185  else if (!strncmp(name, "ND_Trace_", 9))
186  {
187  fec = new ND_Trace_FECollection(atoi(name + 13), atoi(name + 9));
188  }
189  else if (!strncmp(name, "ND_", 3))
190  {
191  fec = new ND_FECollection(atoi(name + 7), atoi(name + 3));
192  }
193  else if (!strncmp(name, "Local_", 6))
194  {
195  fec = new Local_FECollection(name + 6);
196  }
197  else if (!strncmp(name, "NURBS", 5))
198  {
199  fec = new NURBSFECollection(atoi(name + 5));
200  }
201  else
202  mfem_error ("FiniteElementCollection::New : "
203  "Unknown FiniteElementCollection!");
204
205  return fec;
206 }
207
208 const FiniteElement *
210 {
211  switch (GeomType)
212  {
213  case Geometry::POINT: return &PointFE;
214  case Geometry::SEGMENT: return &SegmentFE;
215  case Geometry::TRIANGLE: return &TriangleFE;
217  case Geometry::TETRAHEDRON: return &TetrahedronFE;
218  case Geometry::CUBE: return &ParallelepipedFE;
219  default:
220  mfem_error ("LinearFECollection: unknown geometry type.");
221  }
222  return &SegmentFE; // Make some compilers happy
223 }
224
225 int LinearFECollection::DofForGeometry(int GeomType) const
226 {
227  switch (GeomType)
228  {
229  case Geometry::POINT: return 1;
230  case Geometry::SEGMENT: return 0;
231  case Geometry::TRIANGLE: return 0;
232  case Geometry::SQUARE: return 0;
233  case Geometry::TETRAHEDRON: return 0;
234  case Geometry::CUBE: return 0;
235  default:
236  mfem_error ("LinearFECollection: unknown geometry type.");
237  }
238  return 0; // Make some compilers happy
239 }
240
241 int * LinearFECollection::DofOrderForOrientation(int GeomType, int Or) const
242 {
243  return NULL;
244 }
245
246
247 const FiniteElement *
249 {
250  switch (GeomType)
251  {
252  case Geometry::POINT: return &PointFE;
253  case Geometry::SEGMENT: return &SegmentFE;
254  case Geometry::TRIANGLE: return &TriangleFE;
256  case Geometry::TETRAHEDRON: return &TetrahedronFE;
257  case Geometry::CUBE: return &ParallelepipedFE;
258  default:
259  mfem_error ("QuadraticFECollection: unknown geometry type.");
260  }
261  return &SegmentFE; // Make some compilers happy
262 }
263
265 {
266  switch (GeomType)
267  {
268  case Geometry::POINT: return 1;
269  case Geometry::SEGMENT: return 1;
270  case Geometry::TRIANGLE: return 0;
271  case Geometry::SQUARE: return 1;
272  case Geometry::TETRAHEDRON: return 0;
273  case Geometry::CUBE: return 1;
274  default:
275  mfem_error ("QuadraticFECollection: unknown geometry type.");
276  }
277  return 0; // Make some compilers happy
278 }
279
280 int * QuadraticFECollection::DofOrderForOrientation(int GeomType, int Or) const
281 {
282  static int indexes[] = { 0 };
283
284  return indexes;
285 }
286
287
288 const FiniteElement *
290 {
291  switch (GeomType)
292  {
293  case Geometry::SEGMENT: return &SegmentFE;
295  default:
296  mfem_error ("QuadraticPosFECollection: unknown geometry type.");
297  }
298  return NULL; // Make some compilers happy
299 }
300
302 {
303  switch (GeomType)
304  {
305  case Geometry::POINT: return 1;
306  case Geometry::SEGMENT: return 1;
307  case Geometry::SQUARE: return 1;
308  default:
309  mfem_error ("QuadraticPosFECollection: unknown geometry type.");
310  }
311  return 0; // Make some compilers happy
312 }
313
315 const
316 {
317  static int indexes[] = { 0 };
318
319  return indexes;
320 }
321
322
323 const FiniteElement *
325 {
326  switch (GeomType)
327  {
328  case Geometry::POINT: return &PointFE;
329  case Geometry::SEGMENT: return &SegmentFE;
330  case Geometry::TRIANGLE: return &TriangleFE;
332  case Geometry::TETRAHEDRON: return &TetrahedronFE;
333  case Geometry::CUBE: return &ParallelepipedFE;
334  default:
335  mfem_error ("CubicFECollection: unknown geometry type.");
336  }
337  return &SegmentFE; // Make some compilers happy
338 }
339
340 int CubicFECollection::DofForGeometry(int GeomType) const
341 {
342  switch (GeomType)
343  {
344  case Geometry::POINT: return 1;
345  case Geometry::SEGMENT: return 2;
346  case Geometry::TRIANGLE: return 1;
347  case Geometry::SQUARE: return 4;
348  case Geometry::TETRAHEDRON: return 0;
349  case Geometry::CUBE: return 8;
350  default:
351  mfem_error ("CubicFECollection: unknown geometry type.");
352  }
353  return 0; // Make some compilers happy
354 }
355
356 int * CubicFECollection::DofOrderForOrientation(int GeomType, int Or) const
357 {
358  if (GeomType == Geometry::SEGMENT)
359  {
360  static int ind_pos[] = { 0, 1 };
361  static int ind_neg[] = { 1, 0 };
362
363  if (Or < 0)
364  {
365  return ind_neg;
366  }
367  return ind_pos;
368  }
369  else if (GeomType == Geometry::TRIANGLE)
370  {
371  static int indexes[] = { 0 };
372
373  return indexes;
374  }
375  else if (GeomType == Geometry::SQUARE)
376  {
377  static int sq_ind[8][4] = {{0, 1, 2, 3}, {0, 2, 1, 3},
378  {2, 0, 3, 1}, {1, 0, 3, 2},
379  {3, 2, 1, 0}, {3, 1, 2, 0},
380  {1, 3, 0, 2}, {2, 3, 0, 1}
381  };
382  return sq_ind[Or];
383  }
384
385  return NULL;
386 }
387
388
389 const FiniteElement *
391 {
392  switch (GeomType)
393  {
394  case Geometry::SEGMENT: return &SegmentFE;
395  case Geometry::TRIANGLE: return &TriangleFE;
397  default:
398  mfem_error ("CrouzeixRaviartFECollection: unknown geometry type.");
399  }
400  return &SegmentFE; // Make some compilers happy
401 }
402
404 {
405  switch (GeomType)
406  {
407  case Geometry::POINT: return 0;
408  case Geometry::SEGMENT: return 1;
409  case Geometry::TRIANGLE: return 0;
410  case Geometry::SQUARE: return 0;
411  default:
412  mfem_error ("CrouzeixRaviartFECollection: unknown geometry type.");
413  }
414  return 0; // Make some compilers happy
415 }
416
418 const
419 {
420  static int indexes[] = { 0 };
421
422  return indexes;
423 }
424
425
426 const FiniteElement *
428 {
429  switch (GeomType)
430  {
431  case Geometry::SEGMENT: return &SegmentFE;
432  case Geometry::TRIANGLE: return &TriangleFE;
434  default:
435  mfem_error ("RT0_2DFECollection: unknown geometry type.");
436  }
437  return &SegmentFE; // Make some compilers happy
438 }
439
440 int RT0_2DFECollection::DofForGeometry(int GeomType) const
441 {
442  switch (GeomType)
443  {
444  case Geometry::POINT: return 0;
445  case Geometry::SEGMENT: return 1;
446  case Geometry::TRIANGLE: return 0;
447  case Geometry::SQUARE: return 0;
448  default:
449  mfem_error ("RT0_2DFECollection: unknown geometry type.");
450  }
451  return 0; // Make some compilers happy
452 }
453
455 const
456 {
457  static int ind_pos[] = { 0 };
458  static int ind_neg[] = { -1 };
459
460  if (Or > 0)
461  {
462  return ind_pos;
463  }
464  return ind_neg;
465 }
466
467
468 const FiniteElement *
470 {
471  switch (GeomType)
472  {
473  case Geometry::SEGMENT: return &SegmentFE;
474  case Geometry::TRIANGLE: return &TriangleFE;
476  default:
477  mfem_error ("RT1_2DFECollection: unknown geometry type.");
478  }
479  return &SegmentFE; // Make some compilers happy
480 }
481
482 int RT1_2DFECollection::DofForGeometry(int GeomType) const
483 {
484  switch (GeomType)
485  {
486  case Geometry::POINT: return 0;
487  case Geometry::SEGMENT: return 2;
488  case Geometry::TRIANGLE: return 2;
489  case Geometry::SQUARE: return 4;
490  default:
491  mfem_error ("RT1_2DFECollection: unknown geometry type.");
492  }
493  return 0; // Make some compilers happy
494 }
495
497 const
498 {
499  static int ind_pos[] = { 0, 1 };
500  static int ind_neg[] = { -2, -1 };
501
502  if (Or > 0)
503  {
504  return ind_pos;
505  }
506  return ind_neg;
507 }
508
509 const FiniteElement *
511 {
512  switch (GeomType)
513  {
514  case Geometry::SEGMENT: return &SegmentFE;
515  case Geometry::TRIANGLE: return &TriangleFE;
517  default:
518  mfem_error ("RT2_2DFECollection: unknown geometry type.");
519  }
520  return &SegmentFE; // Make some compilers happy
521 }
522
523 int RT2_2DFECollection::DofForGeometry(int GeomType) const
524 {
525  switch (GeomType)
526  {
527  case Geometry::POINT: return 0;
528  case Geometry::SEGMENT: return 3;
529  case Geometry::TRIANGLE: return 6;
530  case Geometry::SQUARE: return 12;
531  default:
532  mfem_error ("RT2_2DFECollection: unknown geometry type.");
533  }
534  return 0; // Make some compilers happy
535 }
536
538 const
539 {
540  static int ind_pos[] = { 0, 1, 2 };
541  static int ind_neg[] = { -3, -2, -1 };
542
543  if (Or > 0)
544  {
545  return ind_pos;
546  }
547  return ind_neg;
548 }
549
550
551 const FiniteElement *
553 {
554  switch (GeomType)
555  {
556  case Geometry::TRIANGLE: return &TriangleFE;
558  default:
559  mfem_error ("Const2DFECollection: unknown geometry type.");
560  }
561  return &TriangleFE; // Make some compilers happy
562 }
563
564 int Const2DFECollection::DofForGeometry(int GeomType) const
565 {
566  switch (GeomType)
567  {
568  case Geometry::POINT: return 0;
569  case Geometry::SEGMENT: return 0;
570  case Geometry::TRIANGLE: return 1;
571  case Geometry::SQUARE: return 1;
572  default:
573  mfem_error ("Const2DFECollection: unknown geometry type.");
574  }
575  return 0; // Make some compilers happy
576 }
577
579 const
580 {
581  return NULL;
582 }
583
584
585 const FiniteElement *
587 {
588  switch (GeomType)
589  {
590  case Geometry::TRIANGLE: return &TriangleFE;
592  default:
593  mfem_error ("LinearDiscont2DFECollection: unknown geometry type.");
594  }
595  return &TriangleFE; // Make some compilers happy
596 }
597
599 {
600  switch (GeomType)
601  {
602  case Geometry::POINT: return 0;
603  case Geometry::SEGMENT: return 0;
604  case Geometry::TRIANGLE: return 3;
605  case Geometry::SQUARE: return 4;
606  default:
607  mfem_error ("LinearDiscont2DFECollection: unknown geometry type.");
608  }
609  return 0; // Make some compilers happy
610 }
611
613 const
614 {
615  return NULL;
616 }
617
618
619 const FiniteElement *
621 {
622  switch (GeomType)
623  {
624  case Geometry::TRIANGLE: return &TriangleFE;
626  default:
627  mfem_error ("GaussLinearDiscont2DFECollection:"
628  " unknown geometry type.");
629  }
630  return &TriangleFE; // Make some compilers happy
631 }
632
634 {
635  switch (GeomType)
636  {
637  case Geometry::POINT: return 0;
638  case Geometry::SEGMENT: return 0;
639  case Geometry::TRIANGLE: return 3;
640  case Geometry::SQUARE: return 4;
641  default:
642  mfem_error ("GaussLinearDiscont2DFECollection:"
643  " unknown geometry type.");
644  }
645  return 0; // Make some compilers happy
646 }
647
649  int GeomType, int Or) const
650 {
651  return NULL;
652 }
653
654
655 const FiniteElement *
657 {
658  if (GeomType != Geometry::SQUARE)
659  {
660  mfem_error ("P1OnQuadFECollection: unknown geometry type.");
661  }
663 }
664
666 {
667  switch (GeomType)
668  {
669  case Geometry::POINT: return 0;
670  case Geometry::SEGMENT: return 0;
671  case Geometry::SQUARE: return 3;
672  default:
673  mfem_error ("P1OnQuadFECollection: unknown geometry type.");
674  }
675  return 0; // Make some compilers happy
676 }
677
679  int GeomType, int Or) const
680 {
681  return NULL;
682 }
683
684
685 const FiniteElement *
687 {
688  switch (GeomType)
689  {
690  case Geometry::TRIANGLE: return &TriangleFE;
692  default:
693  mfem_error ("QuadraticDiscont2DFECollection: unknown geometry type.");
694  }
695  return &TriangleFE; // Make some compilers happy
696 }
697
699 {
700  switch (GeomType)
701  {
702  case Geometry::POINT: return 0;
703  case Geometry::SEGMENT: return 0;
704  case Geometry::TRIANGLE: return 6;
705  case Geometry::SQUARE: return 9;
706  default:
707  mfem_error ("QuadraticDiscont2DFECollection: unknown geometry type.");
708  }
709  return 0; // Make some compilers happy
710 }
711
713  int GeomType, int Or) const
714 {
715  return NULL;
716 }
717
718
719 const FiniteElement *
721 {
722  switch (GeomType)
723  {
725  default:
726  mfem_error ("QuadraticPosDiscont2DFECollection: unknown geometry type.");
727  }
728  return NULL; // Make some compilers happy
729 }
730
732 {
733  switch (GeomType)
734  {
735  case Geometry::POINT: return 0;
736  case Geometry::SEGMENT: return 0;
737  case Geometry::SQUARE: return 9;
738  default:
739  mfem_error ("QuadraticPosDiscont2DFECollection: unknown geometry type.");
740  }
741  return 0; // Make some compilers happy
742 }
743
744
745 const FiniteElement *
747 const
748 {
749  switch (GeomType)
750  {
751  case Geometry::TRIANGLE: return &TriangleFE;
753  default:
755  " unknown geometry type.");
756  }
757  return &QuadrilateralFE; // Make some compilers happy
758 }
759
761 {
762  switch (GeomType)
763  {
764  case Geometry::POINT: return 0;
765  case Geometry::SEGMENT: return 0;
766  case Geometry::TRIANGLE: return 6;
767  case Geometry::SQUARE: return 9;
768  default:
770  " unknown geometry type.");
771  }
772  return 0; // Make some compilers happy
773 }
774
776  int GeomType, int Or) const
777 {
778  return NULL;
779 }
780
781
782 const FiniteElement *
784 {
785  switch (GeomType)
786  {
787  case Geometry::TRIANGLE: return &TriangleFE;
789  default:
790  mfem_error ("CubicDiscont2DFECollection: unknown geometry type.");
791  }
792  return &TriangleFE; // Make some compilers happy
793 }
794
796 {
797  switch (GeomType)
798  {
799  case Geometry::POINT: return 0;
800  case Geometry::SEGMENT: return 0;
801  case Geometry::TRIANGLE: return 10;
802  case Geometry::SQUARE: return 16;
803  default:
804  mfem_error ("CubicDiscont2DFECollection: unknown geometry type.");
805  }
806  return 0; // Make some compilers happy
807 }
808
810 const
811 {
812  return NULL;
813 }
814
815
816 const FiniteElement *
818 {
819  switch (GeomType)
820  {
821  case Geometry::TRIANGLE: return &TriangleFE;
823  case Geometry::TETRAHEDRON: return &TetrahedronFE;
824  case Geometry::CUBE: return &ParallelepipedFE;
825  default:
826  mfem_error ("LinearNonConf3DFECollection: unknown geometry type.");
827  }
828  return &TriangleFE; // Make some compilers happy
829 }
830
832 {
833  switch (GeomType)
834  {
835  case Geometry::POINT: return 0;
836  case Geometry::SEGMENT: return 0;
837  case Geometry::TRIANGLE: return 1;
838  case Geometry::SQUARE: return 1;
839  case Geometry::TETRAHEDRON: return 0;
840  case Geometry::CUBE: return 0;
841  default:
842  mfem_error ("LinearNonConf3DFECollection: unknown geometry type.");
843  }
844  return 0; // Make some compilers happy
845 }
846
848 const
849 {
850  static int indexes[] = { 0 };
851
852  return indexes;
853 }
854
855
856 const FiniteElement *
858 {
859  switch (GeomType)
860  {
861  case Geometry::TETRAHEDRON: return &TetrahedronFE;
862  case Geometry::CUBE: return &ParallelepipedFE;
863  default:
864  mfem_error ("Const3DFECollection: unknown geometry type.");
865  }
866  return &TetrahedronFE; // Make some compilers happy
867 }
868
869 int Const3DFECollection::DofForGeometry(int GeomType) const
870 {
871  switch (GeomType)
872  {
873  case Geometry::POINT: return 0;
874  case Geometry::SEGMENT: return 0;
875  case Geometry::TRIANGLE: return 0;
876  case Geometry::TETRAHEDRON: return 1;
877  case Geometry::SQUARE: return 0;
878  case Geometry::CUBE: return 1;
879  default:
880  mfem_error ("Const3DFECollection: unknown geometry type.");
881  }
882  return 0; // Make some compilers happy
883 }
884
886 const
887 {
888  return NULL;
889 }
890
891
892 const FiniteElement *
894 {
895  switch (GeomType)
896  {
897  case Geometry::TETRAHEDRON: return &TetrahedronFE;
898  case Geometry::CUBE: return &ParallelepipedFE;
899  default:
900  mfem_error ("LinearDiscont3DFECollection: unknown geometry type.");
901  }
902  return &TetrahedronFE; // Make some compilers happy
903 }
904
906 {
907  switch (GeomType)
908  {
909  case Geometry::POINT: return 0;
910  case Geometry::SEGMENT: return 0;
911  case Geometry::TRIANGLE: return 0;
912  case Geometry::SQUARE: return 0;
913  case Geometry::TETRAHEDRON: return 4;
914  case Geometry::CUBE: return 8;
915  default:
916  mfem_error ("LinearDiscont3DFECollection: unknown geometry type.");
917  }
918  return 0; // Make some compilers happy
919 }
920
922 const
923 {
924  return NULL;
925 }
926
927
928 const FiniteElement *
930 {
931  switch (GeomType)
932  {
933  case Geometry::TETRAHEDRON: return &TetrahedronFE;
934  case Geometry::CUBE: return &ParallelepipedFE;
935  default:
936  mfem_error ("QuadraticDiscont3DFECollection: unknown geometry type.");
937  }
938  return &TetrahedronFE; // Make some compilers happy
939 }
940
942 {
943  switch (GeomType)
944  {
945  case Geometry::POINT: return 0;
946  case Geometry::SEGMENT: return 0;
947  case Geometry::TRIANGLE: return 0;
948  case Geometry::SQUARE: return 0;
949  case Geometry::TETRAHEDRON: return 10;
950  case Geometry::CUBE: return 27;
951  default:
952  mfem_error ("QuadraticDiscont3DFECollection: unknown geometry type.");
953  }
954  return 0; // Make some compilers happy
955 }
956
958  int GeomType, int Or) const
959 {
960  return NULL;
961 }
962
963 const FiniteElement *
965 {
966  switch (GeomType)
967  {
968  case Geometry::POINT: return &PointFE;
969  case Geometry::SEGMENT: return &SegmentFE;
970  case Geometry::TRIANGLE: return &TriangleFE;
972  case Geometry::TETRAHEDRON: return &TetrahedronFE;
973  case Geometry::CUBE: return &ParallelepipedFE;
974  default:
975  mfem_error ("RefinedLinearFECollection: unknown geometry type.");
976  }
977  return &SegmentFE; // Make some compilers happy
978 }
979
981 {
982  switch (GeomType)
983  {
984  case Geometry::POINT: return 1;
985  case Geometry::SEGMENT: return 1;
986  case Geometry::TRIANGLE: return 0;
987  case Geometry::SQUARE: return 1;
988  case Geometry::TETRAHEDRON: return 0;
989  case Geometry::CUBE: return 1;
990  default:
991  mfem_error ("RefinedLinearFECollection: unknown geometry type.");
992  }
993  return 0; // Make some compilers happy
994 }
995
997  int Or) const
998 {
999  static int indexes[] = { 0 };
1000
1001  return indexes;
1002 }
1003
1004
1005 const FiniteElement *
1007 {
1008  switch (GeomType)
1009  {
1010  case Geometry::CUBE: return &HexahedronFE;
1011  case Geometry::TETRAHEDRON: return &TetrahedronFE;
1012  default:
1013  mfem_error ("ND1_3DFECollection: unknown geometry type.");
1014  }
1015  return &HexahedronFE; // Make some compilers happy
1016 }
1017
1018 int ND1_3DFECollection::DofForGeometry(int GeomType) const
1019 {
1020  switch (GeomType)
1021  {
1022  case Geometry::POINT: return 0;
1023  case Geometry::SEGMENT: return 1;
1024  case Geometry::TRIANGLE: return 0;
1025  case Geometry::SQUARE: return 0;
1026  case Geometry::TETRAHEDRON: return 0;
1027  case Geometry::CUBE: return 0;
1028  default:
1029  mfem_error ("ND1_3DFECollection: unknown geometry type.");
1030  }
1031  return 0; // Make some compilers happy
1032 }
1033
1035 const
1036 {
1037  static int ind_pos[] = { 0 };
1038  static int ind_neg[] = { -1 };
1039
1040  if (Or > 0)
1041  {
1042  return ind_pos;
1043  }
1044  return ind_neg;
1045 }
1046
1047
1048 const FiniteElement *
1050 {
1051  switch (GeomType)
1052  {
1053  case Geometry::TRIANGLE: return &TriangleFE;
1055  case Geometry::CUBE: return &HexahedronFE;
1056  case Geometry::TETRAHEDRON: return &TetrahedronFE;
1057  default:
1058  mfem_error ("RT0_3DFECollection: unknown geometry type.");
1059  }
1060  return &HexahedronFE; // Make some compilers happy
1061 }
1062
1063 int RT0_3DFECollection::DofForGeometry(int GeomType) const
1064 {
1065  switch (GeomType)
1066  {
1067  case Geometry::POINT: return 0;
1068  case Geometry::SEGMENT: return 0;
1069  case Geometry::TRIANGLE: return 1;
1070  case Geometry::SQUARE: return 1;
1071  case Geometry::TETRAHEDRON: return 0;
1072  case Geometry::CUBE: return 0;
1073  default:
1074  mfem_error ("RT0_3DFECollection: unknown geometry type.");
1075  }
1076  return 0; // Make some compilers happy
1077 }
1078
1080 const
1081 {
1082  static int ind_pos[] = { 0 };
1083  static int ind_neg[] = { -1 };
1084
1085  if ((GeomType == Geometry::TRIANGLE) || (GeomType == Geometry::SQUARE))
1086  {
1087  if (Or % 2 == 0)
1088  {
1089  return ind_pos;
1090  }
1091  return ind_neg;
1092  }
1093  return NULL;
1094 }
1095
1096 const FiniteElement *
1098 {
1099  switch (GeomType)
1100  {
1101  case Geometry::TRIANGLE: return &TriangleFE;
1103  case Geometry::CUBE: return &HexahedronFE;
1104  default:
1105  mfem_error ("RT1_3DFECollection: unknown geometry type.");
1106  }
1107  return &HexahedronFE; // Make some compilers happy
1108 }
1109
1110 int RT1_3DFECollection::DofForGeometry(int GeomType) const
1111 {
1112  switch (GeomType)
1113  {
1114  case Geometry::POINT: return 0;
1115  case Geometry::SEGMENT: return 0;
1116  case Geometry::TRIANGLE: return 2;
1117  case Geometry::SQUARE: return 4;
1118  case Geometry::CUBE: return 12;
1119  default:
1120  mfem_error ("RT1_3DFECollection: unknown geometry type.");
1121  }
1122  return 0; // Make some compilers happy
1123 }
1124
1126 const
1127 {
1128  if (GeomType == Geometry::SQUARE)
1129  {
1130  static int sq_ind[8][4] =
1131  {
1132  {0, 1, 2, 3}, {-1, -3, -2, -4},
1133  {2, 0, 3, 1}, {-2, -1, -4, -3},
1134  {3, 2, 1, 0}, {-4, -2, -3, -1},
1135  {1, 3, 0, 2}, {-3, -4, -1, -2}
1136  };
1137
1138  return sq_ind[Or];
1139  }
1140  else
1141  {
1142  return NULL;
1143  }
1144 }
1145
1146
1147 H1_FECollection::H1_FECollection(const int p, const int dim, const int type)
1148 {
1149  const int pm1 = p - 1, pm2 = pm1 - 1, pm3 = pm2 - 1;
1150
1151  if (type == 0)
1152  {
1153  snprintf(h1_name, 32, "H1_%dD_P%d", dim, p);
1154  }
1155  else
1156  {
1157  snprintf(h1_name, 32, "H1Pos_%dD_P%d", dim, p);
1158  }
1159
1160  for (int g = 0; g < Geometry::NumGeom; g++)
1161  {
1162  H1_dof[g] = 0;
1163  H1_Elements[g] = NULL;
1164  }
1165  for (int i = 0; i < 2; i++)
1166  {
1167  SegDofOrd[i] = NULL;
1168  }
1169  for (int i = 0; i < 6; i++)
1170  {
1171  TriDofOrd[i] = NULL;
1172  }
1173  for (int i = 0; i < 8; i++)
1174  {
1176  }
1177
1178  H1_dof[Geometry::POINT] = 1;
1179  H1_Elements[Geometry::POINT] = new PointFiniteElement;
1180
1181  if (dim >= 1)
1182  {
1183  H1_dof[Geometry::SEGMENT] = pm1;
1184  if (type == 0)
1185  {
1186  H1_Elements[Geometry::SEGMENT] = new H1_SegmentElement(p);
1187  }
1188  else
1189  {
1190  H1_Elements[Geometry::SEGMENT] = new H1Pos_SegmentElement(p);
1191  }
1192
1193  SegDofOrd[0] = new int[2*pm1];
1194  SegDofOrd[1] = SegDofOrd[0] + pm1;
1195  for (int i = 0; i < pm1; i++)
1196  {
1197  SegDofOrd[0][i] = i;
1198  SegDofOrd[1][i] = pm2 - i;
1199  }
1200  }
1201
1202  if (dim >= 2)
1203  {
1204  H1_dof[Geometry::TRIANGLE] = (pm1*pm2)/2;
1205  H1_dof[Geometry::SQUARE] = pm1*pm1;
1206  if (type == 0)
1207  {
1208  H1_Elements[Geometry::TRIANGLE] = new H1_TriangleElement(p);
1210  }
1211  else
1212  {
1213  H1_Elements[Geometry::TRIANGLE] = NULL; // TODO
1215  }
1216
1217  const int &TriDof = H1_dof[Geometry::TRIANGLE];
1218  const int &QuadDof = H1_dof[Geometry::SQUARE];
1219  TriDofOrd[0] = new int[6*TriDof];
1220  for (int i = 1; i < 6; i++)
1221  {
1222  TriDofOrd[i] = TriDofOrd[i-1] + TriDof;
1223  }
1224  // see Mesh::GetTriOrientation in mesh/mesh.cpp
1225  for (int j = 0; j < pm2; j++)
1226  for (int i = 0; i + j < pm2; i++)
1227  {
1228  int o = TriDof - ((pm1 - j)*(pm2 - j))/2 + i;
1229  int k = pm3 - j - i;
1230  TriDofOrd[0][o] = o; // (0,1,2)
1231  TriDofOrd[1][o] = TriDof - ((pm1-j)*(pm2-j))/2 + k; // (1,0,2)
1232  TriDofOrd[2][o] = TriDof - ((pm1-i)*(pm2-i))/2 + k; // (2,0,1)
1233  TriDofOrd[3][o] = TriDof - ((pm1-k)*(pm2-k))/2 + i; // (2,1,0)
1234  TriDofOrd[4][o] = TriDof - ((pm1-k)*(pm2-k))/2 + j; // (1,2,0)
1235  TriDofOrd[5][o] = TriDof - ((pm1-i)*(pm2-i))/2 + j; // (0,2,1)
1236  }
1237
1239  for (int i = 1; i < 8; i++)
1240  {
1242  }
1243  // see Mesh::GetQuadOrientation in mesh/mesh.cpp
1244  for (int j = 0; j < pm1; j++)
1245  for (int i = 0; i < pm1; i++)
1246  {
1247  int o = i + j*pm1;
1248  QuadDofOrd[0][o] = i + j*pm1; // (0,1,2,3)
1249  QuadDofOrd[1][o] = j + i*pm1; // (0,3,2,1)
1250  QuadDofOrd[2][o] = j + (pm2 - i)*pm1; // (1,2,3,0)
1251  QuadDofOrd[3][o] = (pm2 - i) + j*pm1; // (1,0,3,2)
1252  QuadDofOrd[4][o] = (pm2 - i) + (pm2 - j)*pm1; // (2,3,0,1)
1253  QuadDofOrd[5][o] = (pm2 - j) + (pm2 - i)*pm1; // (2,1,0,3)
1254  QuadDofOrd[6][o] = (pm2 - j) + i*pm1; // (3,0,1,2)
1255  QuadDofOrd[7][o] = i + (pm2 - j)*pm1; // (3,2,1,0)
1256  }
1257
1258  if (dim >= 3)
1259  {
1260  H1_dof[Geometry::TETRAHEDRON] = (TriDof*pm3)/3;
1262  if (type == 0)
1263  {
1264  H1_Elements[Geometry::TETRAHEDRON] = new H1_TetrahedronElement(p);
1265  H1_Elements[Geometry::CUBE] = new H1_HexahedronElement(p);
1266  }
1267  else
1268  {
1269  H1_Elements[Geometry::TETRAHEDRON] = NULL; // TODO
1270  H1_Elements[Geometry::CUBE] = new H1Pos_HexahedronElement(p);
1271  }
1272  }
1273  }
1274 }
1275
1276 int *H1_FECollection::DofOrderForOrientation(int GeomType, int Or) const
1277 {
1278  if (GeomType == Geometry::SEGMENT)
1279  {
1280  if (Or > 0)
1281  {
1282  return SegDofOrd[0];
1283  }
1284  return SegDofOrd[1];
1285  }
1286  else if (GeomType == Geometry::TRIANGLE)
1287  {
1288  return TriDofOrd[Or%6];
1289  }
1290  else if (GeomType == Geometry::SQUARE)
1291  {
1293  }
1294  return NULL;
1295 }
1296
1298 {
1299  int p = H1_dof[Geometry::SEGMENT] + 1;
1300  if (!strncmp(h1_name, "H1_", 3))
1301  {
1302  return new H1_Trace_FECollection(p, atoi(h1_name + 3));
1303  }
1304  else if (!strncmp(h1_name, "H1Pos_", 6))
1305  {
1306  return new H1_Trace_FECollection(p, atoi(h1_name + 6), 1);
1307  }
1308  return NULL;
1309 }
1310
1312 {
1313  delete [] SegDofOrd[0];
1314  delete [] TriDofOrd[0];
1316  for (int g = 0; g < Geometry::NumGeom; g++)
1317  {
1318  delete H1_Elements[g];
1319  }
1320 }
1321
1322
1324  const int type)
1325  : H1_FECollection(p, dim-1, type)
1326 {
1327  if (type == 0)
1328  {
1329  snprintf(h1_name, 32, "H1_Trace_%dD_P%d", dim, p);
1330  }
1331  else
1332  {
1333  snprintf(h1_name, 32, "H1Pos_Trace_%dD_P%d", dim, p);
1334  }
1335 }
1336
1337
1338 L2_FECollection::L2_FECollection(const int p, const int dim, const int type)
1339 {
1340  if (type == 0)
1341  {
1342  snprintf(d_name, 32, "L2_%dD_P%d", dim, p);
1343  }
1344  else
1345  {
1346  snprintf(d_name, 32, "L2_T%d_%dD_P%d", type, dim, p);
1347  }
1348
1349  for (int g = 0; g < Geometry::NumGeom; g++)
1350  {
1351  L2_Elements[g] = NULL;
1352  Tr_Elements[g] = NULL;
1353  }
1354  for (int i = 0; i < 2; i++)
1355  {
1356  SegDofOrd[i] = NULL;
1357  }
1358  for (int i = 0; i < 6; i++)
1359  {
1360  TriDofOrd[i] = NULL;
1361  }
1362
1363  if (dim == 1)
1364  {
1365  if (type == 0 || type == 1)
1366  {
1367  L2_Elements[Geometry::SEGMENT] = new L2_SegmentElement(p, type);
1368  }
1369  else
1370  {
1371  L2_Elements[Geometry::SEGMENT] = new L2Pos_SegmentElement(p);
1372  }
1373
1374  Tr_Elements[Geometry::POINT] = new PointFiniteElement;
1375
1376  const int pp1 = p + 1;
1377  SegDofOrd[0] = new int[2*pp1];
1378  SegDofOrd[1] = SegDofOrd[0] + pp1;
1379  for (int i = 0; i <= p; i++)
1380  {
1381  SegDofOrd[0][i] = i;
1382  SegDofOrd[1][i] = p - i;
1383  }
1384  }
1385  else if (dim == 2)
1386  {
1387  if (type == 0 || type == 1)
1388  {
1389  L2_Elements[Geometry::TRIANGLE] = new L2_TriangleElement(p, type);
1390  L2_Elements[Geometry::SQUARE] = new L2_QuadrilateralElement(p, type);
1391  }
1392  else
1393  {
1394  L2_Elements[Geometry::TRIANGLE] = new L2Pos_TriangleElement(p);
1396  }
1397  Tr_Elements[Geometry::SEGMENT] = new L2_SegmentElement(p, 0);
1398
1399  const int TriDof = L2_Elements[Geometry::TRIANGLE]->GetDof();
1400  TriDofOrd[0] = new int[6*TriDof];
1401  for (int i = 1; i < 6; i++)
1402  {
1403  TriDofOrd[i] = TriDofOrd[i-1] + TriDof;
1404  }
1405  const int pp1 = p + 1, pp2 = pp1 + 1;
1406  for (int j = 0; j <= p; j++)
1407  for (int i = 0; i + j <= p; i++)
1408  {
1409  int o = TriDof - ((pp2 - j)*(pp1 - j))/2 + i;
1410  int k = p - j - i;
1411  TriDofOrd[0][o] = o; // (0,1,2)
1412  TriDofOrd[1][o] = TriDof - ((pp2-j)*(pp1-j))/2 + k; // (1,0,2)
1413  TriDofOrd[2][o] = TriDof - ((pp2-i)*(pp1-i))/2 + k; // (2,0,1)
1414  TriDofOrd[3][o] = TriDof - ((pp2-k)*(pp1-k))/2 + i; // (2,1,0)
1415  TriDofOrd[4][o] = TriDof - ((pp2-k)*(pp1-k))/2 + j; // (1,2,0)
1416  TriDofOrd[5][o] = TriDof - ((pp2-i)*(pp1-i))/2 + j; // (0,2,1)
1417  }
1418  }
1419  else if (dim == 3)
1420  {
1421  if (type == 0 || type == 1)
1422  {
1423  L2_Elements[Geometry::TETRAHEDRON] =
1424  new L2_TetrahedronElement(p, type);
1425  L2_Elements[Geometry::CUBE] = new L2_HexahedronElement(p, type);
1426  }
1427  else
1428  {
1429  L2_Elements[Geometry::TETRAHEDRON] = new L2Pos_TetrahedronElement(p);
1430  L2_Elements[Geometry::CUBE] = new L2Pos_HexahedronElement(p);
1431  }
1432  Tr_Elements[Geometry::TRIANGLE] = new L2_TriangleElement(p, 0);
1433  Tr_Elements[Geometry::SQUARE] = new L2_QuadrilateralElement(p, 0);
1434  }
1435  else
1436  {
1437  cerr << "L2_FECollection::L2_FECollection : dim = "
1438  << dim << endl;
1439  mfem_error();
1440  }
1441 }
1442
1443 int *L2_FECollection::DofOrderForOrientation(int GeomType, int Or) const
1444 {
1445  if (GeomType == Geometry::SEGMENT)
1446  {
1447  if (Or > 0)
1448  {
1449  return SegDofOrd[0];
1450  }
1451  return SegDofOrd[1];
1452  }
1453  else if (GeomType == Geometry::TRIANGLE)
1454  {
1455  return TriDofOrd[Or%6];
1456  }
1457  return NULL;
1458 }
1459
1461 {
1462  delete [] SegDofOrd[0];
1463  delete [] TriDofOrd[0];
1464  for (int i = 0; i < Geometry::NumGeom; i++)
1465  {
1466  delete L2_Elements[i];
1467  delete Tr_Elements[i];
1468  }
1469 }
1470
1471
1472 RT_FECollection::RT_FECollection(const int p, const int dim)
1473 {
1474  InitFaces(p, dim, FiniteElement::INTEGRAL, true);
1475
1476  snprintf(rt_name, 32, "RT_%dD_P%d", dim, p);
1477
1478  const int pp1 = p + 1;
1479  if (dim == 2)
1480  {
1482  RT_dof[Geometry::TRIANGLE] = p*pp1;
1483
1485  RT_dof[Geometry::SQUARE] = 2*p*pp1;
1486  }
1487  else if (dim == 3)
1488  {
1490  RT_dof[Geometry::TETRAHEDRON] = p*pp1*(p + 2)/2;
1491
1493  RT_dof[Geometry::CUBE] = 3*p*pp1*pp1;
1494  }
1495  else
1496  {
1497  cerr << "RT_FECollection::RT_FECollection : dim = " << dim << endl;
1498  mfem_error();
1499  }
1500 }
1501
1502 void RT_FECollection::InitFaces(const int p, const int dim, const int map_type,
1503  const bool signs)
1504 {
1505  const int pp1 = p + 1, pp2 = p + 2;
1506
1507  for (int g = 0; g < Geometry::NumGeom; g++)
1508  {
1509  RT_Elements[g] = NULL;
1510  RT_dof[g] = 0;
1511  }
1512  for (int i = 0; i < 2; i++)
1513  {
1514  SegDofOrd[i] = NULL;
1515  }
1516  for (int i = 0; i < 6; i++)
1517  {
1518  TriDofOrd[i] = NULL;
1519  }
1520  for (int i = 0; i < 8; i++)
1521  {
1523  }
1524
1525  if (dim == 2)
1526  {
1527  L2_SegmentElement *l2_seg = new L2_SegmentElement(p);
1528  l2_seg->SetMapType(map_type);
1529  RT_Elements[Geometry::SEGMENT] = l2_seg;
1530  RT_dof[Geometry::SEGMENT] = pp1;
1531
1532  SegDofOrd[0] = new int[2*pp1];
1533  SegDofOrd[1] = SegDofOrd[0] + pp1;
1534  for (int i = 0; i <= p; i++)
1535  {
1536  SegDofOrd[0][i] = i;
1537  SegDofOrd[1][i] = signs ? (-1 - (p - i)) : (p - i);
1538  }
1539  }
1540  else if (dim == 3)
1541  {
1542  L2_TriangleElement *l2_tri = new L2_TriangleElement(p);
1543  l2_tri->SetMapType(map_type);
1544  RT_Elements[Geometry::TRIANGLE] = l2_tri;
1545  RT_dof[Geometry::TRIANGLE] = pp1*pp2/2;
1546
1550  RT_dof[Geometry::SQUARE] = pp1*pp1;
1551
1552  int TriDof = RT_dof[Geometry::TRIANGLE];
1553  TriDofOrd[0] = new int[6*TriDof];
1554  for (int i = 1; i < 6; i++)
1555  {
1556  TriDofOrd[i] = TriDofOrd[i-1] + TriDof;
1557  }
1558  // see Mesh::GetTriOrientation in mesh/mesh.cpp,
1559  // the constructor of H1_FECollection
1560  for (int j = 0; j <= p; j++)
1561  for (int i = 0; i + j <= p; i++)
1562  {
1563  int o = TriDof - ((pp2 - j)*(pp1 - j))/2 + i;
1564  int k = p - j - i;
1565  TriDofOrd[0][o] = o; // (0,1,2)
1566  TriDofOrd[1][o] = -1-(TriDof-((pp2-j)*(pp1-j))/2+k); // (1,0,2)
1567  TriDofOrd[2][o] = TriDof-((pp2-i)*(pp1-i))/2+k; // (2,0,1)
1568  TriDofOrd[3][o] = -1-(TriDof-((pp2-k)*(pp1-k))/2+i); // (2,1,0)
1569  TriDofOrd[4][o] = TriDof-((pp2-k)*(pp1-k))/2+j; // (1,2,0)
1570  TriDofOrd[5][o] = -1-(TriDof-((pp2-i)*(pp1-i))/2+j); // (0,2,1)
1571  if (!signs)
1572  {
1573  for (int k = 1; k < 6; k += 2)
1574  {
1575  TriDofOrd[k][o] = -1 - TriDofOrd[k][o];
1576  }
1577  }
1578  }
1579
1582  for (int i = 1; i < 8; i++)
1583  {
1585  }
1586  // see Mesh::GetQuadOrientation in mesh/mesh.cpp
1587  for (int j = 0; j <= p; j++)
1588  for (int i = 0; i <= p; i++)
1589  {
1590  int o = i + j*pp1;
1591  QuadDofOrd[0][o] = i + j*pp1; // (0,1,2,3)
1592  QuadDofOrd[1][o] = -1 - (j + i*pp1); // (0,3,2,1)
1593  QuadDofOrd[2][o] = j + (p - i)*pp1; // (1,2,3,0)
1594  QuadDofOrd[3][o] = -1 - ((p - i) + j*pp1); // (1,0,3,2)
1595  QuadDofOrd[4][o] = (p - i) + (p - j)*pp1; // (2,3,0,1)
1596  QuadDofOrd[5][o] = -1 - ((p - j) + (p - i)*pp1); // (2,1,0,3)
1597  QuadDofOrd[6][o] = (p - j) + i*pp1; // (3,0,1,2)
1598  QuadDofOrd[7][o] = -1 - (i + (p - j)*pp1); // (3,2,1,0)
1599  if (!signs)
1600  {
1601  for (int k = 1; k < 8; k += 2)
1602  {
1604  }
1605  }
1606  }
1607  }
1608 }
1609
1610 int *RT_FECollection::DofOrderForOrientation(int GeomType, int Or) const
1611 {
1612  if (GeomType == Geometry::SEGMENT)
1613  {
1614  if (Or > 0)
1615  {
1616  return SegDofOrd[0];
1617  }
1618  return SegDofOrd[1];
1619  }
1620  else if (GeomType == Geometry::TRIANGLE)
1621  {
1622  return TriDofOrd[Or%6];
1623  }
1624  else if (GeomType == Geometry::SQUARE)
1625  {
1627  }
1628  return NULL;
1629 }
1630
1632 {
1633  return new RT_Trace_FECollection(atoi(rt_name + 7), atoi(rt_name + 3));
1634 }
1635
1637 {
1638  delete [] SegDofOrd[0];
1639  delete [] TriDofOrd[0];
1641  for (int g = 0; g < Geometry::NumGeom; g++)
1642  {
1643  delete RT_Elements[g];
1644  }
1645 }
1646
1648  const int map_type)
1649  : RT_FECollection(p, dim, map_type, true)
1650 {
1651  if (map_type == FiniteElement::INTEGRAL)
1652  {
1653  snprintf(rt_name, 32, "RT_Trace_%dD_P%d", dim, p);
1654  }
1655  else
1656  {
1657  snprintf(rt_name, 32, "RT_ValTrace_%dD_P%d", dim, p);
1658  }
1659
1660  MFEM_VERIFY(dim == 2 || dim == 3, "Wrong dimension, dim = " << dim);
1661 }
1662
1664  const int map_type)
1665  : RT_FECollection(p, dim, map_type, false)
1666 {
1667  if (map_type == FiniteElement::VALUE)
1668  {
1669  snprintf(rt_name, 32, "DG_Iface_%dD_P%d", dim, p);
1670  }
1671  else
1672  {
1673  snprintf(rt_name, 32, "DG_IntIface_%dD_P%d", dim, p);
1674  }
1675
1676  MFEM_VERIFY(dim == 2 || dim == 3, "Wrong dimension, dim = " << dim);
1677 }
1678
1679 ND_FECollection::ND_FECollection(const int p, const int dim)
1680 {
1681  const int pm1 = p - 1, pm2 = p - 2;
1682
1683  snprintf(nd_name, 32, "ND_%dD_P%d", dim, p);
1684
1685  for (int g = 0; g < Geometry::NumGeom; g++)
1686  {
1687  ND_Elements[g] = NULL;
1688  ND_dof[g] = 0;
1689  }
1690  for (int i = 0; i < 2; i++)
1691  {
1692  SegDofOrd[i] = NULL;
1693  }
1694  for (int i = 0; i < 6; i++)
1695  {
1696  TriDofOrd[i] = NULL;
1697  }
1698  for (int i = 0; i < 8; i++)
1699  {
1701  }
1702
1703  if (dim >= 1)
1704  {
1707
1708  SegDofOrd[0] = new int[2*p];
1709  SegDofOrd[1] = SegDofOrd[0] + p;
1710  for (int i = 0; i < p; i++)
1711  {
1712  SegDofOrd[0][i] = i;
1713  SegDofOrd[1][i] = -1 - (pm1 - i);
1714  }
1715  }
1716
1717  if (dim >= 2)
1718  {
1720  ND_dof[Geometry::SQUARE] = 2*p*pm1;
1721
1723  ND_dof[Geometry::TRIANGLE] = p*pm1;
1724
1727  for (int i = 1; i < 8; i++)
1728  {
1730  }
1731  // see Mesh::GetQuadOrientation in mesh/mesh.cpp
1732  for (int j = 0; j < pm1; j++)
1733  {
1734  for (int i = 0; i < p; i++)
1735  {
1736  int d1 = i + j*p; // x-component
1737  int d2 = p*pm1 + j + i*pm1; // y-component
1738  // (0,1,2,3)
1741  // (0,3,2,1)
1744  // (1,2,3,0)
1745  // QuadDofOrd[2][d1] = p*pm1 + (pm2 - j) + i*pm1;
1746  // QuadDofOrd[2][d2] = -1 - ((pm1 - i) + j*p);
1747  QuadDofOrd[2][d1] = -1 - (p*pm1 + j + (pm1 - i)*pm1);
1748  QuadDofOrd[2][d2] = i + (pm2 - j)*p;
1749  // (1,0,3,2)
1750  QuadDofOrd[3][d1] = -1 - ((pm1 - i) + j*p);
1751  QuadDofOrd[3][d2] = p*pm1 + (pm2 - j) + i*pm1;
1752  // (2,3,0,1)
1753  QuadDofOrd[4][d1] = -1 - ((pm1 - i) + (pm2 - j)*p);
1754  QuadDofOrd[4][d2] = -1 - (p*pm1 + (pm2 - j) + (pm1 - i)*pm1);
1755  // (2,1,0,3)
1756  QuadDofOrd[5][d1] = -1 - (p*pm1 + (pm2 - j) + (pm1 - i)*pm1);
1757  QuadDofOrd[5][d2] = -1 - ((pm1 - i) + (pm2 - j)*p);
1758  // (3,0,1,2)
1759  // QuadDofOrd[6][d1] = -1 - (p*pm1 + j + (pm1 - i)*pm1);
1760  // QuadDofOrd[6][d2] = i + (pm2 - j)*p;
1761  QuadDofOrd[6][d1] = p*pm1 + (pm2 - j) + i*pm1;
1762  QuadDofOrd[6][d2] = -1 - ((pm1 - i) + j*p);
1763  // (3,2,1,0)
1764  QuadDofOrd[7][d1] = i + (pm2 - j)*p;
1765  QuadDofOrd[7][d2] = -1 - (p*pm1 + j + (pm1 - i)*pm1);
1766  }
1767  }
1768
1769  int TriDof = ND_dof[Geometry::TRIANGLE];
1770  TriDofOrd[0] = new int[6*TriDof];
1771  for (int i = 1; i < 6; i++)
1772  {
1773  TriDofOrd[i] = TriDofOrd[i-1] + TriDof;
1774  }
1775  // see Mesh::GetTriOrientation in mesh/mesh.cpp,
1776  // the constructor of H1_FECollection
1777  for (int j = 0; j <= pm2; j++)
1778  {
1779  for (int i = 0; i + j <= pm2; i++)
1780  {
1781  int k1 = p*pm1 - (p - j)*(pm1 - j) + 2*i;
1782  int k2 = p*pm1 - (p - i)*(pm1 - i) + 2*j;
1783  // (0,1,2)
1784  TriDofOrd[0][k1 ] = k1;
1785  TriDofOrd[0][k1+1] = k1 + 1;
1786  // (0,2,1)
1787  TriDofOrd[5][k1 ] = k2 + 1;
1788  TriDofOrd[5][k1+1] = k2;
1789
1790  // The other orientations can not be supported with the current
1791  // interface. The method Mesh::ReorientTetMesh will ensure that
1792  // only orientations 0 and 5 are generated.
1793  }
1794  }
1795  }
1796
1797  if (dim >= 3)
1798  {
1800  ND_dof[Geometry::CUBE] = 3*p*pm1*pm1;
1801
1803  ND_dof[Geometry::TETRAHEDRON] = p*pm1*pm2/2;
1804  }
1805 }
1806
1807 int *ND_FECollection::DofOrderForOrientation(int GeomType, int Or) const
1808 {
1809  if (GeomType == Geometry::SEGMENT)
1810  {
1811  if (Or > 0)
1812  {
1813  return SegDofOrd[0];
1814  }
1815  return SegDofOrd[1];
1816  }
1817  else if (GeomType == Geometry::TRIANGLE)
1818  {
1819  if (Or != 0 && Or != 5)
1820  {
1821  MFEM_ABORT("ND_FECollection::DofOrderForOrientation: "
1822  "triangle face orientation " << Or << " is not supported! "
1823  "Use Mesh::ReorientTetMesh to fix it.");
1824  }
1825  return TriDofOrd[Or%6];
1826  }
1827  else if (GeomType == Geometry::SQUARE)
1828  {
1830  }
1831  return NULL;
1832 }
1833
1835 {
1837  atoi(nd_name + 3));
1838 }
1839
1841 {
1842  delete [] SegDofOrd[0];
1843  delete [] TriDofOrd[0];
1845  for (int g = 0; g < Geometry::NumGeom; g++)
1846  {
1847  delete ND_Elements[g];
1848  }
1849 }
1850
1851
1853  : ND_FECollection(p, dim-1)
1854 {
1855  snprintf(nd_name, 32, "ND_Trace_%dD_P%d", dim, p);
1856 }
1857
1858
1860 {
1861  snprintf(d_name, 32, "Local_%s", fe_name);
1862
1863  Local_Element = NULL;
1864
1865  if (!strcmp(fe_name, "BiCubic2DFiniteElement") ||
1867  {
1868  GeomType = Geometry::SQUARE;
1869  Local_Element = new BiCubic2DFiniteElement;
1870  }
1871  else if (!strcmp(fe_name, "Nedelec1HexFiniteElement") ||
1872  !strcmp(fe_name, "Hex_ND1"))
1873  {
1874  GeomType = Geometry::CUBE;
1875  Local_Element = new Nedelec1HexFiniteElement;
1876  }
1877  else if (!strncmp(fe_name, "H1_", 3))
1878  {
1879  GeomType = Geometry::SQUARE;
1880  Local_Element = new H1_QuadrilateralElement(atoi(fe_name + 7));
1881  }
1882  else if (!strncmp(fe_name, "L2_", 3))
1883  {
1884  GeomType = Geometry::SQUARE;
1885  Local_Element = new L2_QuadrilateralElement(atoi(fe_name + 7));
1886  }
1887  else
1888  {
1889  cerr << "Local_FECollection::Local_FECollection : fe_name = "
1890  << fe_name << endl;
1891  mfem_error();
1892  }
1893 }
1894
1895
1896 void NURBSFECollection::Allocate(int Order)
1897 {
1898  SegmentFE = new NURBS1DFiniteElement(Order);
1900  ParallelepipedFE = new NURBS3DFiniteElement(Order);
1901
1902  snprintf(name, 16, "NURBS%i", Order);
1903 }
1904
1905 void NURBSFECollection::Deallocate()
1906 {
1907  delete ParallelepipedFE;
1909  delete SegmentFE;
1910 }
1911
1912 const FiniteElement *
1914 {
1915  switch (GeomType)
1916  {
1917  case Geometry::SEGMENT: return SegmentFE;
1919  case Geometry::CUBE: return ParallelepipedFE;
1920  default:
1921  mfem_error ("NURBSFECollection: unknown geometry type.");
1922  }
1923  return SegmentFE; // Make some compilers happy
1924 }
1925
1926 int NURBSFECollection::DofForGeometry(int GeomType) const
1927 {
1928  mfem_error("NURBSFECollection::DofForGeometry");
1929  return 0; // Make some compilers happy
1930 }
1931
1932 int *NURBSFECollection::DofOrderForOrientation(int GeomType, int Or) const
1933 {
1934  mfem_error("NURBSFECollection::DofOrderForOrientation");
1935  return NULL;
1936 }
1937
1939 {
1940  MFEM_ABORT("NURBS finite elements can not be statically condensed!");
1941  return NULL;
1942 }
1943
1944 }
Abstract class for Finite Elements.
Definition: fe.hpp:44
Arbitrary order non-uniform rational B-splines (NURBS) finite elements.
Definition: fe_coll.hpp:217
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:929
Version of QuadraticDiscont2DFECollection with dofs in the Gaussian points.
Definition: fe_coll.hpp:570
int RT_dof[Geometry::NumGeom]
Definition: fe_coll.hpp:138
FiniteElement * ND_Elements[Geometry::NumGeom]
Definition: fe_coll.hpp:189
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:1034
virtual FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:38
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:301
FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:1938
static const int NumGeom
Definition: geom.hpp:34
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:1932
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:1097
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:857
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:775
virtual ~L2_FECollection()
Definition: fe_coll.cpp:1460
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:390
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:1018
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:620
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:678
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:656
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:586
Piecewise-(bi)linear continuous finite elements.
Definition: fe_coll.hpp:260
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:1125
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:957
ND_FECollection(const int p, const int dim)
Definition: fe_coll.cpp:1679
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:496
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:964
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:1063
FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:1834
L2_FECollection(const int p, const int dim, const int type=0)
Definition: fe_coll.cpp:1338
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:809
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:795
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:578
void InitFaces(const int p, const int dim, const int map_type, const bool signs)
Definition: fe_coll.cpp:1502
FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:1631
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:996
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:893
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:523
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:746
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:469
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:225
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:248
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:510
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:648
virtual ~RT_FECollection()
Definition: fe_coll.cpp:1636
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:564
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:1610
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:1079
Piecewise-(bi)cubic continuous finite elements.
Definition: fe_coll.hpp:327
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:612
PointFiniteElement PointFE
Definition: point.cpp:30
int dim
Definition: ex3.cpp:48
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:1110
FiniteElement * RT_Elements[Geometry::NumGeom]
Definition: fe_coll.hpp:137
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:314
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:598
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:241
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:665
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:869
ND_Trace_FECollection(const int p, const int dim)
Definition: fe_coll.cpp:1852
TriLinear3DFiniteElement HexahedronFE
Definition: hexahedron.cpp:68
Discontinuous collection defined locally by a given finite element.
Definition: fe_coll.hpp:763
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:289
Version of QuadraticDiscont2DFECollection with positive basis functions.
Definition: fe_coll.hpp:554
virtual ~ND_FECollection()
Definition: fe_coll.cpp:1840
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:885
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:427
RT_Trace_FECollection(const int p, const int dim, const int map_type=FiniteElement::INTEGRAL)
Definition: fe_coll.cpp:1647
H1_FECollection(const int p, const int dim=3, const int type=0)
Definition: fe_coll.cpp:1147
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:698
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:552
Version of QuadraticFECollection with positive basis functions.
Definition: fe_coll.hpp:307
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:817
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:980
virtual ~H1_FECollection()
Definition: fe_coll.cpp:1311
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:209
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:1049
FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:1297
Piecewise-linear nonconforming finite elements in 3D.
Definition: fe_coll.hpp:371
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:356
DG_Interface_FECollection(const int p, const int dim, const int map_type=FiniteElement::VALUE)
Definition: fe_coll.cpp:1663
Crouzeix-Raviart nonconforming elements in 2D.
Definition: fe_coll.hpp:351
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:133
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:417
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:403
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:482
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:905
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:712
int GetDof() const
Returns the degrees of freedom in the FE space.
Definition: fe.hpp:88
void mfem_error(const char *msg)
Definition: error.cpp:23
H1_Trace_FECollection(const int p, const int dim, const int type=0)
Definition: fe_coll.cpp:1323
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:720
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:280
Linear3DFiniteElement TetrahedronFE
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:633
Linear2DFiniteElement TriangleFE
Definition: triangle.cpp:134
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:941
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:783
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:760
Version of LinearDiscont2DFECollection with dofs in the Gaussian points.
Definition: fe_coll.hpp:498
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:1913
static FiniteElementCollection * New(const char *name)
Definition: fe_coll.cpp:44
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:921
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:1006
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:264
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:537
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:1807
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:1276
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:454
RT_FECollection(const int p, const int dim, const int map_type, const bool signs)
Definition: fe_coll.hpp:146
Definition: fe_coll.hpp:283
int HasFaceDofs(int GeomType) const
Definition: fe_coll.cpp:25
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:185
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:847
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:1443
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:54
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:1926
Local_FECollection(const char *fe_name)
Definition: fe_coll.cpp:1859
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:731
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:831
Linear (P1) finite elements on quadrilaterals.
Definition: fe_coll.hpp:519
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:340
Linear1DFiniteElement SegmentFE
Definition: segment.cpp:49
int ND_dof[Geometry::NumGeom]
Definition: fe_coll.hpp:190
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:440
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:324
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:95
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:686