MFEM  v3.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
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.googlecode.com.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
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  FiniteElementCollection *fec = NULL;
41 
42  if (!strcmp(name, "Linear"))
43  fec = new LinearFECollection;
44  else if (!strcmp(name, "Quadratic"))
45  fec = new QuadraticFECollection;
46  else if (!strcmp(name, "QuadraticPos"))
47  fec = new QuadraticPosFECollection;
48  else if (!strcmp(name, "Cubic"))
49  fec = new CubicFECollection;
50  else if (!strcmp(name, "Const3D"))
51  fec = new Const3DFECollection;
52  else if (!strcmp(name, "Const2D"))
53  fec = new Const2DFECollection;
54  else if (!strcmp(name, "LinearDiscont2D"))
56  else if (!strcmp(name, "GaussLinearDiscont2D"))
58  else if (!strcmp(name, "P1OnQuad"))
59  fec = new P1OnQuadFECollection;
60  else if (!strcmp(name, "QuadraticDiscont2D"))
62  else if (!strcmp(name, "QuadraticPosDiscont2D"))
64  else if (!strcmp(name, "GaussQuadraticDiscont2D"))
66  else if (!strcmp(name, "CubicDiscont2D"))
68  else if (!strcmp(name, "LinearDiscont3D"))
70  else if (!strcmp(name, "QuadraticDiscont3D"))
72  else if (!strcmp(name, "LinearNonConf3D"))
74  else if (!strcmp(name, "CrouzeixRaviart"))
76  else if (!strcmp(name, "ND1_3D"))
77  fec = new ND1_3DFECollection;
78  else if (!strcmp(name, "RT0_2D"))
79  fec = new RT0_2DFECollection;
80  else if (!strcmp(name, "RT1_2D"))
81  fec = new RT1_2DFECollection;
82  else if (!strcmp(name, "RT2_2D"))
83  fec = new RT2_2DFECollection;
84  else if (!strcmp(name, "RT0_3D"))
85  fec = new RT0_3DFECollection;
86  else if (!strcmp(name, "RT1_3D"))
87  fec = new RT1_3DFECollection;
88  else if (!strncmp(name, "H1_", 3))
89  fec = new H1_FECollection(atoi(name + 7), atoi(name + 3));
90  else if (!strncmp(name, "H1Pos_", 6))
91  fec = new H1Pos_FECollection(atoi(name + 10), atoi(name + 6));
92  else if (!strncmp(name, "L2_T", 4))
93  fec = new L2_FECollection(atoi(name + 10), atoi(name + 6),
94  atoi(name + 4));
95  else if (!strncmp(name, "L2_", 3))
96  fec = new L2_FECollection(atoi(name + 7), atoi(name + 3));
97  else if (!strncmp(name, "RT_", 3))
98  fec = new RT_FECollection(atoi(name + 7), atoi(name + 3));
99  else if (!strncmp(name, "RT_Trace_", 9))
100  fec = new RT_Trace_FECollection(atoi(name + 13), atoi(name + 9));
101  else if (!strncmp(name, "RT_ValTrace_", 12))
102  fec = new RT_Trace_FECollection(atoi(name + 16), atoi(name + 12),
104  else if (!strncmp(name, "ND_", 3))
105  fec = new ND_FECollection(atoi(name + 7), atoi(name + 3));
106  else if (!strncmp(name, "Local_", 6))
107  fec = new Local_FECollection(name + 6);
108  else if (!strncmp(name, "NURBS", 5))
109  fec = new NURBSFECollection(atoi(name + 5));
110  else
111  mfem_error ("FiniteElementCollection::New : "
112  "Unknown FiniteElementCollection!");
113 
114  return fec;
115 }
116 
117 const FiniteElement *
119 {
120  switch (GeomType)
121  {
122  case Geometry::POINT: return &PointFE;
123  case Geometry::SEGMENT: return &SegmentFE;
124  case Geometry::TRIANGLE: return &TriangleFE;
125  case Geometry::SQUARE: return &QuadrilateralFE;
126  case Geometry::TETRAHEDRON: return &TetrahedronFE;
127  case Geometry::CUBE: return &ParallelepipedFE;
128  default:
129  mfem_error ("LinearFECollection: unknown geometry type.");
130  }
131  return &SegmentFE; // Make some compilers happy
132 }
133 
134 int LinearFECollection::DofForGeometry(int GeomType) const
135 {
136  switch (GeomType)
137  {
138  case Geometry::POINT: return 1;
139  case Geometry::SEGMENT: return 0;
140  case Geometry::TRIANGLE: return 0;
141  case Geometry::SQUARE: return 0;
142  case Geometry::TETRAHEDRON: return 0;
143  case Geometry::CUBE: return 0;
144  default:
145  mfem_error ("LinearFECollection: unknown geometry type.");
146  }
147  return 0; // Make some compilers happy
148 }
149 
150 int * LinearFECollection::DofOrderForOrientation(int GeomType, int Or) const
151 {
152  return NULL;
153 }
154 
155 
156 const FiniteElement *
158 {
159  switch (GeomType)
160  {
161  case Geometry::POINT: return &PointFE;
162  case Geometry::SEGMENT: return &SegmentFE;
163  case Geometry::TRIANGLE: return &TriangleFE;
164  case Geometry::SQUARE: return &QuadrilateralFE;
165  case Geometry::TETRAHEDRON: return &TetrahedronFE;
166  case Geometry::CUBE: return &ParallelepipedFE;
167  default:
168  mfem_error ("QuadraticFECollection: unknown geometry type.");
169  }
170  return &SegmentFE; // Make some compilers happy
171 }
172 
174 {
175  switch (GeomType)
176  {
177  case Geometry::POINT: return 1;
178  case Geometry::SEGMENT: return 1;
179  case Geometry::TRIANGLE: return 0;
180  case Geometry::SQUARE: return 1;
181  case Geometry::TETRAHEDRON: return 0;
182  case Geometry::CUBE: return 1;
183  default:
184  mfem_error ("QuadraticFECollection: unknown geometry type.");
185  }
186  return 0; // Make some compilers happy
187 }
188 
189 int * QuadraticFECollection::DofOrderForOrientation(int GeomType, int Or) const
190 {
191  static int indexes[] = { 0 };
192 
193  return indexes;
194 }
195 
196 
197 const FiniteElement *
199 {
200  switch (GeomType)
201  {
202  case Geometry::SEGMENT: return &SegmentFE;
203  case Geometry::SQUARE: return &QuadrilateralFE;
204  default:
205  mfem_error ("QuadraticPosFECollection: unknown geometry type.");
206  }
207  return NULL; // Make some compilers happy
208 }
209 
211 {
212  switch (GeomType)
213  {
214  case Geometry::POINT: return 1;
215  case Geometry::SEGMENT: return 1;
216  case Geometry::SQUARE: return 1;
217  default:
218  mfem_error ("QuadraticPosFECollection: unknown geometry type.");
219  }
220  return 0; // Make some compilers happy
221 }
222 
224  const
225 {
226  static int indexes[] = { 0 };
227 
228  return indexes;
229 }
230 
231 
232 const FiniteElement *
234 {
235  switch (GeomType)
236  {
237  case Geometry::POINT: return &PointFE;
238  case Geometry::SEGMENT: return &SegmentFE;
239  case Geometry::TRIANGLE: return &TriangleFE;
240  case Geometry::SQUARE: return &QuadrilateralFE;
241  case Geometry::TETRAHEDRON: return &TetrahedronFE;
242  case Geometry::CUBE: return &ParallelepipedFE;
243  default:
244  mfem_error ("CubicFECollection: unknown geometry type.");
245  }
246  return &SegmentFE; // Make some compilers happy
247 }
248 
249 int CubicFECollection::DofForGeometry(int GeomType) const
250 {
251  switch (GeomType)
252  {
253  case Geometry::POINT: return 1;
254  case Geometry::SEGMENT: return 2;
255  case Geometry::TRIANGLE: return 1;
256  case Geometry::SQUARE: return 4;
257  case Geometry::TETRAHEDRON: return 0;
258  case Geometry::CUBE: return 8;
259  default:
260  mfem_error ("CubicFECollection: unknown geometry type.");
261  }
262  return 0; // Make some compilers happy
263 }
264 
265 int * CubicFECollection::DofOrderForOrientation(int GeomType, int Or) const
266 {
267  if (GeomType == Geometry::SEGMENT)
268  {
269  static int ind_pos[] = { 0, 1 };
270  static int ind_neg[] = { 1, 0 };
271 
272  if (Or < 0)
273  return ind_neg;
274  return ind_pos;
275  }
276  else if (GeomType == Geometry::TRIANGLE)
277  {
278  static int indexes[] = { 0 };
279 
280  return indexes;
281  }
282  else if (GeomType == Geometry::SQUARE)
283  {
284  static int sq_ind[8][4] = {{0, 1, 2, 3}, {0, 2, 1, 3},
285  {2, 0, 3, 1}, {1, 0, 3, 2},
286  {3, 2, 1, 0}, {3, 1, 2, 0},
287  {1, 3, 0, 2}, {2, 3, 0, 1}};
288  return sq_ind[Or];
289  }
290 
291  return NULL;
292 }
293 
294 
295 const FiniteElement *
297 {
298  switch (GeomType)
299  {
300  case Geometry::SEGMENT: return &SegmentFE;
301  case Geometry::TRIANGLE: return &TriangleFE;
302  case Geometry::SQUARE: return &QuadrilateralFE;
303  default:
304  mfem_error ("CrouzeixRaviartFECollection: unknown geometry type.");
305  }
306  return &SegmentFE; // Make some compilers happy
307 }
308 
310 {
311  switch (GeomType)
312  {
313  case Geometry::POINT: return 0;
314  case Geometry::SEGMENT: return 1;
315  case Geometry::TRIANGLE: return 0;
316  case Geometry::SQUARE: return 0;
317  default:
318  mfem_error ("CrouzeixRaviartFECollection: unknown geometry type.");
319  }
320  return 0; // Make some compilers happy
321 }
322 
324  const
325 {
326  static int indexes[] = { 0 };
327 
328  return indexes;
329 }
330 
331 
332 const FiniteElement *
334 {
335  switch (GeomType)
336  {
337  case Geometry::SEGMENT: return &SegmentFE;
338  case Geometry::TRIANGLE: return &TriangleFE;
339  case Geometry::SQUARE: return &QuadrilateralFE;
340  default:
341  mfem_error ("RT0_2DFECollection: unknown geometry type.");
342  }
343  return &SegmentFE; // Make some compilers happy
344 }
345 
346 int RT0_2DFECollection::DofForGeometry(int GeomType) const
347 {
348  switch (GeomType)
349  {
350  case Geometry::POINT: return 0;
351  case Geometry::SEGMENT: return 1;
352  case Geometry::TRIANGLE: return 0;
353  case Geometry::SQUARE: return 0;
354  default:
355  mfem_error ("RT0_2DFECollection: unknown geometry type.");
356  }
357  return 0; // Make some compilers happy
358 }
359 
361  const
362 {
363  static int ind_pos[] = { 0 };
364  static int ind_neg[] = { -1 };
365 
366  if (Or > 0)
367  return ind_pos;
368  return ind_neg;
369 }
370 
371 
372 const FiniteElement *
374 {
375  switch (GeomType)
376  {
377  case Geometry::SEGMENT: return &SegmentFE;
378  case Geometry::TRIANGLE: return &TriangleFE;
379  case Geometry::SQUARE: return &QuadrilateralFE;
380  default:
381  mfem_error ("RT1_2DFECollection: unknown geometry type.");
382  }
383  return &SegmentFE; // Make some compilers happy
384 }
385 
386 int RT1_2DFECollection::DofForGeometry(int GeomType) const
387 {
388  switch (GeomType)
389  {
390  case Geometry::POINT: return 0;
391  case Geometry::SEGMENT: return 2;
392  case Geometry::TRIANGLE: return 2;
393  case Geometry::SQUARE: return 4;
394  default:
395  mfem_error ("RT1_2DFECollection: unknown geometry type.");
396  }
397  return 0; // Make some compilers happy
398 }
399 
401  const
402 {
403  static int ind_pos[] = { 0, 1 };
404  static int ind_neg[] = { -2, -1 };
405 
406  if (Or > 0)
407  return ind_pos;
408  return ind_neg;
409 }
410 
411 const FiniteElement *
413 {
414  switch (GeomType)
415  {
416  case Geometry::SEGMENT: return &SegmentFE;
417  case Geometry::TRIANGLE: return &TriangleFE;
418  case Geometry::SQUARE: return &QuadrilateralFE;
419  default:
420  mfem_error ("RT2_2DFECollection: unknown geometry type.");
421  }
422  return &SegmentFE; // Make some compilers happy
423 }
424 
425 int RT2_2DFECollection::DofForGeometry(int GeomType) const
426 {
427  switch (GeomType)
428  {
429  case Geometry::POINT: return 0;
430  case Geometry::SEGMENT: return 3;
431  case Geometry::TRIANGLE: return 6;
432  case Geometry::SQUARE: return 12;
433  default:
434  mfem_error ("RT2_2DFECollection: unknown geometry type.");
435  }
436  return 0; // Make some compilers happy
437 }
438 
440  const
441 {
442  static int ind_pos[] = { 0, 1, 2 };
443  static int ind_neg[] = { -3, -2, -1 };
444 
445  if (Or > 0)
446  return ind_pos;
447  return ind_neg;
448 }
449 
450 
451 const FiniteElement *
453 {
454  switch (GeomType)
455  {
456  case Geometry::TRIANGLE: return &TriangleFE;
457  case Geometry::SQUARE: return &QuadrilateralFE;
458  default:
459  mfem_error ("Const2DFECollection: unknown geometry type.");
460  }
461  return &TriangleFE; // Make some compilers happy
462 }
463 
464 int Const2DFECollection::DofForGeometry(int GeomType) const
465 {
466  switch (GeomType)
467  {
468  case Geometry::POINT: return 0;
469  case Geometry::SEGMENT: return 0;
470  case Geometry::TRIANGLE: return 1;
471  case Geometry::SQUARE: return 1;
472  default:
473  mfem_error ("Const2DFECollection: unknown geometry type.");
474  }
475  return 0; // Make some compilers happy
476 }
477 
479  const
480 {
481  return NULL;
482 }
483 
484 
485 const FiniteElement *
487 {
488  switch (GeomType)
489  {
490  case Geometry::TRIANGLE: return &TriangleFE;
491  case Geometry::SQUARE: return &QuadrilateralFE;
492  default:
493  mfem_error ("LinearDiscont2DFECollection: unknown geometry type.");
494  }
495  return &TriangleFE; // Make some compilers happy
496 }
497 
499 {
500  switch (GeomType)
501  {
502  case Geometry::POINT: return 0;
503  case Geometry::SEGMENT: return 0;
504  case Geometry::TRIANGLE: return 3;
505  case Geometry::SQUARE: return 4;
506  default:
507  mfem_error ("LinearDiscont2DFECollection: unknown geometry type.");
508  }
509  return 0; // Make some compilers happy
510 }
511 
513  const
514 {
515  return NULL;
516 }
517 
518 
519 const FiniteElement *
521 {
522  switch (GeomType)
523  {
524  case Geometry::TRIANGLE: return &TriangleFE;
525  case Geometry::SQUARE: return &QuadrilateralFE;
526  default:
527  mfem_error ("GaussLinearDiscont2DFECollection:"
528  " unknown geometry type.");
529  }
530  return &TriangleFE; // Make some compilers happy
531 }
532 
534 {
535  switch (GeomType)
536  {
537  case Geometry::POINT: return 0;
538  case Geometry::SEGMENT: return 0;
539  case Geometry::TRIANGLE: return 3;
540  case Geometry::SQUARE: return 4;
541  default:
542  mfem_error ("GaussLinearDiscont2DFECollection:"
543  " unknown geometry type.");
544  }
545  return 0; // Make some compilers happy
546 }
547 
549  int GeomType, int Or) const
550 {
551  return NULL;
552 }
553 
554 
555 const FiniteElement *
557 {
558  if (GeomType != Geometry::SQUARE)
559  {
560  mfem_error ("P1OnQuadFECollection: unknown geometry type.");
561  }
562  return &QuadrilateralFE;
563 }
564 
566 {
567  switch (GeomType)
568  {
569  case Geometry::POINT: return 0;
570  case Geometry::SEGMENT: return 0;
571  case Geometry::SQUARE: return 3;
572  default:
573  mfem_error ("P1OnQuadFECollection: unknown geometry type.");
574  }
575  return 0; // Make some compilers happy
576 }
577 
579  int GeomType, int Or) const
580 {
581  return NULL;
582 }
583 
584 
585 const FiniteElement *
587 {
588  switch (GeomType)
589  {
590  case Geometry::TRIANGLE: return &TriangleFE;
591  case Geometry::SQUARE: return &QuadrilateralFE;
592  default:
593  mfem_error ("QuadraticDiscont2DFECollection: 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 6;
605  case Geometry::SQUARE: return 9;
606  default:
607  mfem_error ("QuadraticDiscont2DFECollection: unknown geometry type.");
608  }
609  return 0; // Make some compilers happy
610 }
611 
613  int GeomType, int Or) const
614 {
615  return NULL;
616 }
617 
618 
619 const FiniteElement *
621 {
622  switch (GeomType)
623  {
624  case Geometry::SQUARE: return &QuadrilateralFE;
625  default:
626  mfem_error ("QuadraticPosDiscont2DFECollection: unknown geometry type.");
627  }
628  return NULL; // Make some compilers happy
629 }
630 
632 {
633  switch (GeomType)
634  {
635  case Geometry::POINT: return 0;
636  case Geometry::SEGMENT: return 0;
637  case Geometry::SQUARE: return 9;
638  default:
639  mfem_error ("QuadraticPosDiscont2DFECollection: unknown geometry type.");
640  }
641  return 0; // Make some compilers happy
642 }
643 
644 
645 const FiniteElement *
647  const
648 {
649  switch (GeomType)
650  {
651  case Geometry::TRIANGLE: return &TriangleFE;
652  case Geometry::SQUARE: return &QuadrilateralFE;
653  default:
654  mfem_error ("GaussQuadraticDiscont2DFECollection:"
655  " unknown geometry type.");
656  }
657  return &QuadrilateralFE; // Make some compilers happy
658 }
659 
661 {
662  switch (GeomType)
663  {
664  case Geometry::POINT: return 0;
665  case Geometry::SEGMENT: return 0;
666  case Geometry::TRIANGLE: return 6;
667  case Geometry::SQUARE: return 9;
668  default:
669  mfem_error ("GaussQuadraticDiscont2DFECollection:"
670  " unknown geometry type.");
671  }
672  return 0; // Make some compilers happy
673 }
674 
676  int GeomType, int Or) const
677 {
678  return NULL;
679 }
680 
681 
682 const FiniteElement *
684 {
685  switch (GeomType)
686  {
687  case Geometry::TRIANGLE: return &TriangleFE;
688  case Geometry::SQUARE: return &QuadrilateralFE;
689  default:
690  mfem_error ("CubicDiscont2DFECollection: unknown geometry type.");
691  }
692  return &TriangleFE; // Make some compilers happy
693 }
694 
696 {
697  switch (GeomType)
698  {
699  case Geometry::POINT: return 0;
700  case Geometry::SEGMENT: return 0;
701  case Geometry::TRIANGLE: return 10;
702  case Geometry::SQUARE: return 16;
703  default:
704  mfem_error ("CubicDiscont2DFECollection: unknown geometry type.");
705  }
706  return 0; // Make some compilers happy
707 }
708 
710  const
711 {
712  return NULL;
713 }
714 
715 
716 const FiniteElement *
718 {
719  switch (GeomType)
720  {
721  case Geometry::TRIANGLE: return &TriangleFE;
722  case Geometry::SQUARE: return &QuadrilateralFE;
723  case Geometry::TETRAHEDRON: return &TetrahedronFE;
724  case Geometry::CUBE: return &ParallelepipedFE;
725  default:
726  mfem_error ("LinearNonConf3DFECollection: unknown geometry type.");
727  }
728  return &TriangleFE; // 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::TRIANGLE: return 1;
738  case Geometry::SQUARE: return 1;
739  case Geometry::TETRAHEDRON: return 0;
740  case Geometry::CUBE: return 0;
741  default:
742  mfem_error ("LinearNonConf3DFECollection: unknown geometry type.");
743  }
744  return 0; // Make some compilers happy
745 }
746 
748  const
749 {
750  static int indexes[] = { 0 };
751 
752  return indexes;
753 }
754 
755 
756 const FiniteElement *
758 {
759  switch (GeomType)
760  {
761  case Geometry::TETRAHEDRON: return &TetrahedronFE;
762  case Geometry::CUBE: return &ParallelepipedFE;
763  default:
764  mfem_error ("Const3DFECollection: unknown geometry type.");
765  }
766  return &TetrahedronFE; // Make some compilers happy
767 }
768 
769 int Const3DFECollection::DofForGeometry(int GeomType) const
770 {
771  switch (GeomType)
772  {
773  case Geometry::POINT: return 0;
774  case Geometry::SEGMENT: return 0;
775  case Geometry::TRIANGLE: return 0;
776  case Geometry::TETRAHEDRON: return 1;
777  case Geometry::SQUARE: return 0;
778  case Geometry::CUBE: return 1;
779  default:
780  mfem_error ("Const3DFECollection: unknown geometry type.");
781  }
782  return 0; // Make some compilers happy
783 }
784 
786  const
787 {
788  return NULL;
789 }
790 
791 
792 const FiniteElement *
794 {
795  switch (GeomType)
796  {
797  case Geometry::TETRAHEDRON: return &TetrahedronFE;
798  case Geometry::CUBE: return &ParallelepipedFE;
799  default:
800  mfem_error ("LinearDiscont3DFECollection: unknown geometry type.");
801  }
802  return &TetrahedronFE; // Make some compilers happy
803 }
804 
806 {
807  switch (GeomType)
808  {
809  case Geometry::POINT: return 0;
810  case Geometry::SEGMENT: return 0;
811  case Geometry::TRIANGLE: return 0;
812  case Geometry::SQUARE: return 0;
813  case Geometry::TETRAHEDRON: return 4;
814  case Geometry::CUBE: return 8;
815  default:
816  mfem_error ("LinearDiscont3DFECollection: unknown geometry type.");
817  }
818  return 0; // Make some compilers happy
819 }
820 
822  const
823 {
824  return NULL;
825 }
826 
827 
828 const FiniteElement *
830 {
831  switch (GeomType)
832  {
833  case Geometry::TETRAHEDRON: return &TetrahedronFE;
834  case Geometry::CUBE: return &ParallelepipedFE;
835  default:
836  mfem_error ("QuadraticDiscont3DFECollection: unknown geometry type.");
837  }
838  return &TetrahedronFE; // Make some compilers happy
839 }
840 
842 {
843  switch (GeomType)
844  {
845  case Geometry::POINT: return 0;
846  case Geometry::SEGMENT: return 0;
847  case Geometry::TRIANGLE: return 0;
848  case Geometry::SQUARE: return 0;
849  case Geometry::TETRAHEDRON: return 10;
850  case Geometry::CUBE: return 27;
851  default:
852  mfem_error ("QuadraticDiscont3DFECollection: unknown geometry type.");
853  }
854  return 0; // Make some compilers happy
855 }
856 
858  int GeomType, int Or) const
859 {
860  return NULL;
861 }
862 
863 const FiniteElement *
865 {
866  switch (GeomType)
867  {
868  case Geometry::POINT: return &PointFE;
869  case Geometry::SEGMENT: return &SegmentFE;
870  case Geometry::TRIANGLE: return &TriangleFE;
871  case Geometry::SQUARE: return &QuadrilateralFE;
872  case Geometry::TETRAHEDRON: return &TetrahedronFE;
873  case Geometry::CUBE: return &ParallelepipedFE;
874  default:
875  mfem_error ("RefinedLinearFECollection: unknown geometry type.");
876  }
877  return &SegmentFE; // Make some compilers happy
878 }
879 
881 {
882  switch (GeomType)
883  {
884  case Geometry::POINT: return 1;
885  case Geometry::SEGMENT: return 1;
886  case Geometry::TRIANGLE: return 0;
887  case Geometry::SQUARE: return 1;
888  case Geometry::TETRAHEDRON: return 0;
889  case Geometry::CUBE: return 1;
890  default:
891  mfem_error ("RefinedLinearFECollection: unknown geometry type.");
892  }
893  return 0; // Make some compilers happy
894 }
895 
896 int * RefinedLinearFECollection::DofOrderForOrientation(int GeomType, int Or) const
897 {
898  static int indexes[] = { 0 };
899 
900  return indexes;
901 }
902 
903 
904 const FiniteElement *
906 {
907  switch (GeomType)
908  {
909  case Geometry::CUBE: return &HexahedronFE;
910  case Geometry::TETRAHEDRON: return &TetrahedronFE;
911  default:
912  mfem_error ("ND1_3DFECollection: unknown geometry type.");
913  }
914  return &HexahedronFE; // Make some compilers happy
915 }
916 
917 int ND1_3DFECollection::DofForGeometry(int GeomType) const
918 {
919  switch (GeomType)
920  {
921  case Geometry::POINT: return 0;
922  case Geometry::SEGMENT: return 1;
923  case Geometry::TRIANGLE: return 0;
924  case Geometry::SQUARE: return 0;
925  case Geometry::TETRAHEDRON: return 0;
926  case Geometry::CUBE: return 0;
927  default:
928  mfem_error ("ND1_3DFECollection: unknown geometry type.");
929  }
930  return 0; // Make some compilers happy
931 }
932 
934  const
935 {
936  static int ind_pos[] = { 0 };
937  static int ind_neg[] = { -1 };
938 
939  if (Or > 0)
940  return ind_pos;
941  return ind_neg;
942 }
943 
944 
945 const FiniteElement *
947 {
948  switch (GeomType)
949  {
950  case Geometry::TRIANGLE: return &TriangleFE;
951  case Geometry::SQUARE: return &QuadrilateralFE;
952  case Geometry::CUBE: return &HexahedronFE;
953  case Geometry::TETRAHEDRON: return &TetrahedronFE;
954  default:
955  mfem_error ("RT0_3DFECollection: unknown geometry type.");
956  }
957  return &HexahedronFE; // Make some compilers happy
958 }
959 
960 int RT0_3DFECollection::DofForGeometry(int GeomType) const
961 {
962  switch (GeomType)
963  {
964  case Geometry::POINT: return 0;
965  case Geometry::SEGMENT: return 0;
966  case Geometry::TRIANGLE: return 1;
967  case Geometry::SQUARE: return 1;
968  case Geometry::TETRAHEDRON: return 0;
969  case Geometry::CUBE: return 0;
970  default:
971  mfem_error ("RT0_3DFECollection: unknown geometry type.");
972  }
973  return 0; // Make some compilers happy
974 }
975 
977  const
978 {
979  static int ind_pos[] = { 0 };
980  static int ind_neg[] = { -1 };
981 
982  if ((GeomType == Geometry::TRIANGLE) || (GeomType == Geometry::SQUARE))
983  {
984  if (Or % 2 == 0)
985  return ind_pos;
986  return ind_neg;
987  }
988  return NULL;
989 }
990 
991 const FiniteElement *
993 {
994  switch (GeomType)
995  {
996  case Geometry::TRIANGLE: return &TriangleFE;
997  case Geometry::SQUARE: return &QuadrilateralFE;
998  case Geometry::CUBE: return &HexahedronFE;
999  default:
1000  mfem_error ("RT1_3DFECollection: unknown geometry type.");
1001  }
1002  return &HexahedronFE; // Make some compilers happy
1003 }
1004 
1005 int RT1_3DFECollection::DofForGeometry(int GeomType) const
1006 {
1007  switch (GeomType)
1008  {
1009  case Geometry::POINT: return 0;
1010  case Geometry::SEGMENT: return 0;
1011  case Geometry::TRIANGLE: return 2;
1012  case Geometry::SQUARE: return 4;
1013  case Geometry::CUBE: return 12;
1014  default:
1015  mfem_error ("RT1_3DFECollection: unknown geometry type.");
1016  }
1017  return 0; // Make some compilers happy
1018 }
1019 
1021  const
1022 {
1023  if (GeomType == Geometry::SQUARE)
1024  {
1025  static int sq_ind[8][4] = {
1026  {0, 1, 2, 3}, {-1, -3, -2, -4},
1027  {2, 0, 3, 1}, {-2, -1, -4, -3},
1028  {3, 2, 1, 0}, {-4, -2, -3, -1},
1029  {1, 3, 0, 2}, {-3, -4, -1, -2}
1030  };
1031 
1032  return sq_ind[Or];
1033  }
1034  else
1035  return NULL;
1036 }
1037 
1038 
1039 H1_FECollection::H1_FECollection(const int p, const int dim, const int type)
1040 {
1041  const int pm1 = p - 1, pm2 = pm1 - 1, pm3 = pm2 - 1;
1042 
1043  if (type == 0)
1044  snprintf(h1_name, 32, "H1_%dD_P%d", dim, p);
1045  else
1046  snprintf(h1_name, 32, "H1Pos_%dD_P%d", dim, p);
1047 
1048  for (int g = 0; g < Geometry::NumGeom; g++)
1049  {
1050  H1_dof[g] = 0;
1051  H1_Elements[g] = NULL;
1052  }
1053  for (int i = 0; i < 2; i++)
1054  SegDofOrd[i] = NULL;
1055  for (int i = 0; i < 6; i++)
1056  TriDofOrd[i] = NULL;
1057  for (int i = 0; i < 8; i++)
1058  QuadDofOrd[i] = NULL;
1059 
1060  H1_dof[Geometry::POINT] = 1;
1061  H1_Elements[Geometry::POINT] = new PointFiniteElement;
1062  H1_dof[Geometry::SEGMENT] = pm1;
1063  if (type == 0)
1064  H1_Elements[Geometry::SEGMENT] = new H1_SegmentElement(p);
1065  else
1066  H1_Elements[Geometry::SEGMENT] = new H1Pos_SegmentElement(p);
1067 
1068  SegDofOrd[0] = new int[2*pm1];
1069  SegDofOrd[1] = SegDofOrd[0] + pm1;
1070  for (int i = 0; i < pm1; i++)
1071  {
1072  SegDofOrd[0][i] = i;
1073  SegDofOrd[1][i] = pm2 - i;
1074  }
1075 
1076  if (dim >= 2)
1077  {
1078  H1_dof[Geometry::TRIANGLE] = (pm1*pm2)/2;
1079  H1_dof[Geometry::SQUARE] = pm1*pm1;
1080  if (type == 0)
1081  {
1082  H1_Elements[Geometry::TRIANGLE] = new H1_TriangleElement(p);
1083  H1_Elements[Geometry::SQUARE] = new H1_QuadrilateralElement(p);
1084  }
1085  else
1086  {
1087  H1_Elements[Geometry::TRIANGLE] = NULL; // TODO
1088  H1_Elements[Geometry::SQUARE] = new H1Pos_QuadrilateralElement(p);
1089  }
1090 
1091  const int &TriDof = H1_dof[Geometry::TRIANGLE];
1092  const int &QuadDof = H1_dof[Geometry::SQUARE];
1093  TriDofOrd[0] = new int[6*TriDof];
1094  for (int i = 1; i < 6; i++)
1095  TriDofOrd[i] = TriDofOrd[i-1] + TriDof;
1096  // see Mesh::GetTriOrientation in mesh/mesh.cpp
1097  for (int j = 0; j < pm2; j++)
1098  for (int i = 0; i + j < pm2; i++)
1099  {
1100  int o = TriDof - ((pm1 - j)*(pm2 - j))/2 + i;
1101  int k = pm3 - j - i;
1102  TriDofOrd[0][o] = o; // (0,1,2)
1103  TriDofOrd[1][o] = TriDof - ((pm1-j)*(pm2-j))/2 + k; // (1,0,2)
1104  TriDofOrd[2][o] = TriDof - ((pm1-i)*(pm2-i))/2 + k; // (2,0,1)
1105  TriDofOrd[3][o] = TriDof - ((pm1-k)*(pm2-k))/2 + i; // (2,1,0)
1106  TriDofOrd[4][o] = TriDof - ((pm1-k)*(pm2-k))/2 + j; // (1,2,0)
1107  TriDofOrd[5][o] = TriDof - ((pm1-i)*(pm2-i))/2 + j; // (0,2,1)
1108  }
1109 
1110  QuadDofOrd[0] = new int[8*QuadDof];
1111  for (int i = 1; i < 8; i++)
1112  QuadDofOrd[i] = QuadDofOrd[i-1] + QuadDof;
1113  // see Mesh::GetQuadOrientation in mesh/mesh.cpp
1114  for (int j = 0; j < pm1; j++)
1115  for (int i = 0; i < pm1; i++)
1116  {
1117  int o = i + j*pm1;
1118  QuadDofOrd[0][o] = i + j*pm1; // (0,1,2,3)
1119  QuadDofOrd[1][o] = j + i*pm1; // (0,3,2,1)
1120  QuadDofOrd[2][o] = j + (pm2 - i)*pm1; // (1,2,3,0)
1121  QuadDofOrd[3][o] = (pm2 - i) + j*pm1; // (1,0,3,2)
1122  QuadDofOrd[4][o] = (pm2 - i) + (pm2 - j)*pm1; // (2,3,0,1)
1123  QuadDofOrd[5][o] = (pm2 - j) + (pm2 - i)*pm1; // (2,1,0,3)
1124  QuadDofOrd[6][o] = (pm2 - j) + i*pm1; // (3,0,1,2)
1125  QuadDofOrd[7][o] = i + (pm2 - j)*pm1; // (3,2,1,0)
1126  }
1127 
1128  if (dim >= 3)
1129  {
1130  H1_dof[Geometry::TETRAHEDRON] = (TriDof*pm3)/3;
1131  H1_dof[Geometry::CUBE] = QuadDof*pm1;
1132  if (type == 0)
1133  {
1134  H1_Elements[Geometry::TETRAHEDRON] = new H1_TetrahedronElement(p);
1135  H1_Elements[Geometry::CUBE] = new H1_HexahedronElement(p);
1136  }
1137  else
1138  {
1139  H1_Elements[Geometry::TETRAHEDRON] = NULL; // TODO
1140  H1_Elements[Geometry::CUBE] = new H1Pos_HexahedronElement(p);
1141  }
1142  }
1143  }
1144 }
1145 
1146 int *H1_FECollection::DofOrderForOrientation(int GeomType, int Or) const
1147 {
1148  if (GeomType == Geometry::SEGMENT)
1149  {
1150  if (Or > 0)
1151  return SegDofOrd[0];
1152  return SegDofOrd[1];
1153  }
1154  else if (GeomType == Geometry::TRIANGLE)
1155  {
1156  return TriDofOrd[Or%6];
1157  }
1158  else if (GeomType == Geometry::SQUARE)
1159  {
1160  return QuadDofOrd[Or%8];
1161  }
1162  return NULL;
1163 }
1164 
1166 {
1167  delete [] SegDofOrd[0];
1168  delete [] TriDofOrd[0];
1169  delete [] QuadDofOrd[0];
1170  for (int g = 0; g < Geometry::NumGeom; g++)
1171  delete H1_Elements[g];
1172 }
1173 
1174 
1175 L2_FECollection::L2_FECollection(const int p, const int dim, const int type)
1176 {
1177  if (type == 0)
1178  snprintf(d_name, 32, "L2_%dD_P%d", dim, p);
1179  else
1180  snprintf(d_name, 32, "L2_T%d_%dD_P%d", type, dim, p);
1181 
1182  for (int g = 0; g < Geometry::NumGeom; g++)
1183  {
1184  L2_Elements[g] = NULL;
1185  Tr_Elements[g] = NULL;
1186  }
1187  for (int i = 0; i < 2; i++)
1188  SegDofOrd[i] = NULL;
1189  for (int i = 0; i < 6; i++)
1190  TriDofOrd[i] = NULL;
1191 
1192  if (dim == 1)
1193  {
1194  if (type == 0 || type == 1)
1195  L2_Elements[Geometry::SEGMENT] = new L2_SegmentElement(p, type);
1196  else
1197  L2_Elements[Geometry::SEGMENT] = new L2Pos_SegmentElement(p);
1198 
1199  Tr_Elements[Geometry::POINT] = new PointFiniteElement;
1200 
1201  const int pp1 = p + 1;
1202  SegDofOrd[0] = new int[2*pp1];
1203  SegDofOrd[1] = SegDofOrd[0] + pp1;
1204  for (int i = 0; i <= p; i++)
1205  {
1206  SegDofOrd[0][i] = i;
1207  SegDofOrd[1][i] = p - i;
1208  }
1209  }
1210  else if (dim == 2)
1211  {
1212  if (type == 0 || type == 1)
1213  {
1214  L2_Elements[Geometry::TRIANGLE] = new L2_TriangleElement(p, type);
1215  L2_Elements[Geometry::SQUARE] = new L2_QuadrilateralElement(p, type);
1216  }
1217  else
1218  {
1219  L2_Elements[Geometry::TRIANGLE] = new L2Pos_TriangleElement(p);
1220  L2_Elements[Geometry::SQUARE] = new L2Pos_QuadrilateralElement(p);
1221  }
1222  Tr_Elements[Geometry::SEGMENT] = new L2_SegmentElement(p, 0);
1223 
1224  const int TriDof = L2_Elements[Geometry::TRIANGLE]->GetDof();
1225  TriDofOrd[0] = new int[6*TriDof];
1226  for (int i = 1; i < 6; i++)
1227  TriDofOrd[i] = TriDofOrd[i-1] + TriDof;
1228  const int pp1 = p + 1, pp2 = pp1 + 1;
1229  for (int j = 0; j <= p; j++)
1230  for (int i = 0; i + j <= p; i++)
1231  {
1232  int o = TriDof - ((pp2 - j)*(pp1 - j))/2 + i;
1233  int k = p - j - i;
1234  TriDofOrd[0][o] = o; // (0,1,2)
1235  TriDofOrd[1][o] = TriDof - ((pp2-j)*(pp1-j))/2 + k; // (1,0,2)
1236  TriDofOrd[2][o] = TriDof - ((pp2-i)*(pp1-i))/2 + k; // (2,0,1)
1237  TriDofOrd[3][o] = TriDof - ((pp2-k)*(pp1-k))/2 + i; // (2,1,0)
1238  TriDofOrd[4][o] = TriDof - ((pp2-k)*(pp1-k))/2 + j; // (1,2,0)
1239  TriDofOrd[5][o] = TriDof - ((pp2-i)*(pp1-i))/2 + j; // (0,2,1)
1240  }
1241  }
1242  else if (dim == 3)
1243  {
1244  if (type == 0 || type == 1)
1245  {
1246  L2_Elements[Geometry::TETRAHEDRON] =
1247  new L2_TetrahedronElement(p, type);
1248  L2_Elements[Geometry::CUBE] = new L2_HexahedronElement(p, type);
1249  }
1250  else
1251  {
1252  L2_Elements[Geometry::TETRAHEDRON] = new L2Pos_TetrahedronElement(p);
1253  L2_Elements[Geometry::CUBE] = new L2Pos_HexahedronElement(p);
1254  }
1255  Tr_Elements[Geometry::TRIANGLE] = new L2_TriangleElement(p, 0);
1256  Tr_Elements[Geometry::SQUARE] = new L2_QuadrilateralElement(p, 0);
1257  }
1258  else
1259  {
1260  cerr << "L2_FECollection::L2_FECollection : dim = "
1261  << dim << endl;
1262  mfem_error();
1263  }
1264 }
1265 
1266 int *L2_FECollection::DofOrderForOrientation(int GeomType, int Or) const
1267 {
1268  if (GeomType == Geometry::SEGMENT)
1269  {
1270  if (Or > 0)
1271  return SegDofOrd[0];
1272  return SegDofOrd[1];
1273  }
1274  else if (GeomType == Geometry::TRIANGLE)
1275  {
1276  return TriDofOrd[Or%6];
1277  }
1278  return NULL;
1279 }
1280 
1282 {
1283  delete [] SegDofOrd[0];
1284  delete [] TriDofOrd[0];
1285  for (int i = 0; i < Geometry::NumGeom; i++)
1286  {
1287  delete L2_Elements[i];
1288  delete Tr_Elements[i];
1289  }
1290 }
1291 
1292 
1293 RT_FECollection::RT_FECollection(const int p, const int dim)
1294 {
1295  InitFaces(p, dim, FiniteElement::INTEGRAL);
1296 
1297  snprintf(rt_name, 32, "RT_%dD_P%d", dim, p);
1298 
1299  const int pp1 = p + 1;
1300  if (dim == 2)
1301  {
1302  RT_Elements[Geometry::TRIANGLE] = new RT_TriangleElement(p);
1303  RT_dof[Geometry::TRIANGLE] = p*pp1;
1304 
1305  RT_Elements[Geometry::SQUARE] = new RT_QuadrilateralElement(p);
1306  RT_dof[Geometry::SQUARE] = 2*p*pp1;
1307  }
1308  else if (dim == 3)
1309  {
1310  RT_Elements[Geometry::TETRAHEDRON] = new RT_TetrahedronElement(p);
1311  RT_dof[Geometry::TETRAHEDRON] = p*pp1*(p + 2)/2;
1312 
1313  RT_Elements[Geometry::CUBE] = new RT_HexahedronElement(p);
1314  RT_dof[Geometry::CUBE] = 3*p*pp1*pp1;
1315  }
1316  else
1317  {
1318  cerr << "RT_FECollection::RT_FECollection : dim = " << dim << endl;
1319  mfem_error();
1320  }
1321 }
1322 
1323 void RT_FECollection::InitFaces(const int p, const int dim, const int map_type)
1324 {
1325  const int pp1 = p + 1, pp2 = p + 2;
1326 
1327  for (int g = 0; g < Geometry::NumGeom; g++)
1328  {
1329  RT_Elements[g] = NULL;
1330  RT_dof[g] = 0;
1331  }
1332  for (int i = 0; i < 2; i++)
1333  SegDofOrd[i] = NULL;
1334  for (int i = 0; i < 6; i++)
1335  TriDofOrd[i] = NULL;
1336  for (int i = 0; i < 8; i++)
1337  QuadDofOrd[i] = NULL;
1338 
1339  if (dim == 2)
1340  {
1341  L2_SegmentElement *l2_seg = new L2_SegmentElement(p);
1342  l2_seg->SetMapType(map_type);
1343  RT_Elements[Geometry::SEGMENT] = l2_seg;
1344  RT_dof[Geometry::SEGMENT] = pp1;
1345 
1346  SegDofOrd[0] = new int[2*pp1];
1347  SegDofOrd[1] = SegDofOrd[0] + pp1;
1348  for (int i = 0; i <= p; i++)
1349  {
1350  SegDofOrd[0][i] = i;
1351  SegDofOrd[1][i] = -1 - (p - i);
1352  }
1353  }
1354  else if (dim == 3)
1355  {
1356  L2_TriangleElement *l2_tri = new L2_TriangleElement(p);
1357  l2_tri->SetMapType(map_type);
1358  RT_Elements[Geometry::TRIANGLE] = l2_tri;
1359  RT_dof[Geometry::TRIANGLE] = pp1*pp2/2;
1360 
1362  l2_quad->SetMapType(map_type);
1363  RT_Elements[Geometry::SQUARE] = l2_quad;
1364  RT_dof[Geometry::SQUARE] = pp1*pp1;
1365 
1366  int TriDof = RT_dof[Geometry::TRIANGLE];
1367  TriDofOrd[0] = new int[6*TriDof];
1368  for (int i = 1; i < 6; i++)
1369  TriDofOrd[i] = TriDofOrd[i-1] + TriDof;
1370  // see Mesh::GetTriOrientation in mesh/mesh.cpp,
1371  // the constructor of H1_FECollection
1372  for (int j = 0; j <= p; j++)
1373  for (int i = 0; i + j <= p; i++)
1374  {
1375  int o = TriDof - ((pp2 - j)*(pp1 - j))/2 + i;
1376  int k = p - j - i;
1377  TriDofOrd[0][o] = o; // (0,1,2)
1378  TriDofOrd[1][o] = -1-(TriDof-((pp2-j)*(pp1-j))/2+k); // (1,0,2)
1379  TriDofOrd[2][o] = TriDof-((pp2-i)*(pp1-i))/2+k; // (2,0,1)
1380  TriDofOrd[3][o] = -1-(TriDof-((pp2-k)*(pp1-k))/2+i); // (2,1,0)
1381  TriDofOrd[4][o] = TriDof-((pp2-k)*(pp1-k))/2+j; // (1,2,0)
1382  TriDofOrd[5][o] = -1-(TriDof-((pp2-i)*(pp1-i))/2+j); // (0,2,1)
1383  }
1384 
1385  int QuadDof = RT_dof[Geometry::SQUARE];
1386  QuadDofOrd[0] = new int[8*QuadDof];
1387  for (int i = 1; i < 8; i++)
1388  QuadDofOrd[i] = QuadDofOrd[i-1] + QuadDof;
1389  // see Mesh::GetQuadOrientation in mesh/mesh.cpp
1390  for (int j = 0; j <= p; j++)
1391  for (int i = 0; i <= p; i++)
1392  {
1393  int o = i + j*pp1;
1394  QuadDofOrd[0][o] = i + j*pp1; // (0,1,2,3)
1395  QuadDofOrd[1][o] = -1 - (j + i*pp1); // (0,3,2,1)
1396  QuadDofOrd[2][o] = j + (p - i)*pp1; // (1,2,3,0)
1397  QuadDofOrd[3][o] = -1 - ((p - i) + j*pp1); // (1,0,3,2)
1398  QuadDofOrd[4][o] = (p - i) + (p - j)*pp1; // (2,3,0,1)
1399  QuadDofOrd[5][o] = -1 - ((p - j) + (p - i)*pp1); // (2,1,0,3)
1400  QuadDofOrd[6][o] = (p - j) + i*pp1; // (3,0,1,2)
1401  QuadDofOrd[7][o] = -1 - (i + (p - j)*pp1); // (3,2,1,0)
1402  }
1403  }
1404 }
1405 
1406 int *RT_FECollection::DofOrderForOrientation(int GeomType, int Or) const
1407 {
1408  if (GeomType == Geometry::SEGMENT)
1409  {
1410  if (Or > 0)
1411  return SegDofOrd[0];
1412  return SegDofOrd[1];
1413  }
1414  else if (GeomType == Geometry::TRIANGLE)
1415  {
1416  return TriDofOrd[Or%6];
1417  }
1418  else if (GeomType == Geometry::SQUARE)
1419  {
1420  return QuadDofOrd[Or%8];
1421  }
1422  return NULL;
1423 }
1424 
1426 {
1427  delete [] SegDofOrd[0];
1428  delete [] TriDofOrd[0];
1429  delete [] QuadDofOrd[0];
1430  for (int g = 0; g < Geometry::NumGeom; g++)
1431  delete RT_Elements[g];
1432 }
1433 
1435  const int map_type)
1436  : RT_FECollection(p, dim, map_type)
1437 {
1438  if (map_type == FiniteElement::INTEGRAL)
1439  snprintf(rt_name, 32, "RT_Trace_%dD_P%d", dim, p);
1440  else
1441  snprintf(rt_name, 32, "RT_ValTrace_%dD_P%d", dim, p);
1442 
1443  MFEM_VERIFY(dim == 2 || dim == 3, "Wrong dimension, dim = " << dim);
1444 }
1445 
1446 ND_FECollection::ND_FECollection(const int p, const int dim)
1447 {
1448  const int pm1 = p - 1, pm2 = p - 2;
1449 
1450  snprintf(nd_name, 32, "ND_%dD_P%d", dim, p);
1451 
1452  for (int g = 0; g < Geometry::NumGeom; g++)
1453  {
1454  ND_Elements[g] = NULL;
1455  ND_dof[g] = 0;
1456  }
1457  for (int i = 0; i < 2; i++)
1458  SegDofOrd[i] = NULL;
1459  for (int i = 0; i < 6; i++)
1460  TriDofOrd[i] = NULL;
1461  for (int i = 0; i < 8; i++)
1462  QuadDofOrd[i] = NULL;
1463 
1464  if (dim == 2 || dim == 3)
1465  {
1466  ND_Elements[Geometry::SQUARE] = new ND_QuadrilateralElement(p);
1467  ND_dof[Geometry::SQUARE] = 2*p*pm1;
1468 
1469  ND_Elements[Geometry::TRIANGLE] = new ND_TriangleElement(p);
1470  ND_dof[Geometry::TRIANGLE] = p*pm1;
1471 
1472  L2_SegmentElement *l2_seg = new L2_SegmentElement(p-1);
1473  l2_seg->SetMapType(FiniteElement::INTEGRAL);
1474  ND_Elements[Geometry::SEGMENT] = l2_seg;
1475  ND_dof[Geometry::SEGMENT] = p;
1476 
1477  SegDofOrd[0] = new int[2*p];
1478  SegDofOrd[1] = SegDofOrd[0] + p;
1479  for (int i = 0; i < p; i++)
1480  {
1481  SegDofOrd[0][i] = i;
1482  SegDofOrd[1][i] = -1 - (pm1 - i);
1483  }
1484  }
1485  else
1486  {
1487  mfem_error("ND_FECollection::ND_FECollection : dim != 2 or 3");
1488  }
1489 
1490  if (dim == 3)
1491  {
1492  ND_Elements[Geometry::CUBE] = new ND_HexahedronElement(p);
1493  ND_dof[Geometry::CUBE] = 3*p*pm1*pm1;
1494 
1495  ND_Elements[Geometry::TETRAHEDRON] = new ND_TetrahedronElement(p);
1496  ND_dof[Geometry::TETRAHEDRON] = p*pm1*pm2/2;
1497 
1498  int QuadDof = ND_dof[Geometry::SQUARE];
1499  QuadDofOrd[0] = new int[8*QuadDof];
1500  for (int i = 1; i < 8; i++)
1501  QuadDofOrd[i] = QuadDofOrd[i-1] + QuadDof;
1502  // see Mesh::GetQuadOrientation in mesh/mesh.cpp
1503  for (int j = 0; j < pm1; j++)
1504  for (int i = 0; i < p; i++)
1505  {
1506  int d1 = i + j*p; // x-component
1507  int d2 = p*pm1 + j + i*pm1; // y-component
1508  // (0,1,2,3)
1509  QuadDofOrd[0][d1] = d1;
1510  QuadDofOrd[0][d2] = d2;
1511  // (0,3,2,1)
1512  QuadDofOrd[1][d1] = d2;
1513  QuadDofOrd[1][d2] = d1;
1514  // (1,2,3,0)
1515  // QuadDofOrd[2][d1] = p*pm1 + (pm2 - j) + i*pm1;
1516  // QuadDofOrd[2][d2] = -1 - ((pm1 - i) + j*p);
1517  QuadDofOrd[2][d1] = -1 - (p*pm1 + j + (pm1 - i)*pm1);
1518  QuadDofOrd[2][d2] = i + (pm2 - j)*p;
1519  // (1,0,3,2)
1520  QuadDofOrd[3][d1] = -1 - ((pm1 - i) + j*p);
1521  QuadDofOrd[3][d2] = p*pm1 + (pm2 - j) + i*pm1;
1522  // (2,3,0,1)
1523  QuadDofOrd[4][d1] = -1 - ((pm1 - i) + (pm2 - j)*p);
1524  QuadDofOrd[4][d2] = -1 - (p*pm1 + (pm2 - j) + (pm1 - i)*pm1);
1525  // (2,1,0,3)
1526  QuadDofOrd[5][d1] = -1 - (p*pm1 + (pm2 - j) + (pm1 - i)*pm1);
1527  QuadDofOrd[5][d2] = -1 - ((pm1 - i) + (pm2 - j)*p);
1528  // (3,0,1,2)
1529  // QuadDofOrd[6][d1] = -1 - (p*pm1 + j + (pm1 - i)*pm1);
1530  // QuadDofOrd[6][d2] = i + (pm2 - j)*p;
1531  QuadDofOrd[6][d1] = p*pm1 + (pm2 - j) + i*pm1;
1532  QuadDofOrd[6][d2] = -1 - ((pm1 - i) + j*p);
1533  // (3,2,1,0)
1534  QuadDofOrd[7][d1] = i + (pm2 - j)*p;
1535  QuadDofOrd[7][d2] = -1 - (p*pm1 + j + (pm1 - i)*pm1);
1536  }
1537 
1538  int TriDof = ND_dof[Geometry::TRIANGLE];
1539  TriDofOrd[0] = new int[6*TriDof];
1540  for (int i = 1; i < 6; i++)
1541  TriDofOrd[i] = TriDofOrd[i-1] + TriDof;
1542  // see Mesh::GetTriOrientation in mesh/mesh.cpp,
1543  // the constructor of H1_FECollection
1544  for (int j = 0; j <= pm2; j++)
1545  for (int i = 0; i + j <= pm2; i++)
1546  {
1547  int k1 = p*pm1 - (p - j)*(pm1 - j) + 2*i;
1548  int k2 = p*pm1 - (p - i)*(pm1 - i) + 2*j;
1549  // (0,1,2)
1550  TriDofOrd[0][k1 ] = k1;
1551  TriDofOrd[0][k1+1] = k1 + 1;
1552  // (0,2,1)
1553  TriDofOrd[5][k1 ] = k2 + 1;
1554  TriDofOrd[5][k1+1] = k2;
1555 
1556  // The other orientations can not be supported with the current
1557  // interface. The method Mesh::ReorientTetMesh will ensure that
1558  // only orientations 0 and 5 are generated.
1559  }
1560  }
1561 }
1562 
1563 int *ND_FECollection::DofOrderForOrientation(int GeomType, int Or) const
1564 {
1565  if (GeomType == Geometry::SEGMENT)
1566  {
1567  if (Or > 0)
1568  return SegDofOrd[0];
1569  return SegDofOrd[1];
1570  }
1571  else if (GeomType == Geometry::TRIANGLE)
1572  {
1573  if (Or != 0 && Or != 5)
1574  {
1575  cerr <<
1576  "ND_FECollection::DofOrderForOrientation :\n"
1577  " triangle face orientation " << Or << " is not supported!\n"
1578  " Use Mesh::ReorientTetMesh to fix it." << endl;
1579  mfem_error();
1580  }
1581  return TriDofOrd[Or%6];
1582  }
1583  else if (GeomType == Geometry::SQUARE)
1584  {
1585  return QuadDofOrd[Or%8];
1586  }
1587  return NULL;
1588 }
1589 
1591 {
1592  delete [] SegDofOrd[0];
1593  delete [] TriDofOrd[0];
1594  delete [] QuadDofOrd[0];
1595  for (int g = 0; g < Geometry::NumGeom; g++)
1596  delete ND_Elements[g];
1597 }
1598 
1599 
1601 {
1602  snprintf(d_name, 32, "Local_%s", fe_name);
1603 
1604  Local_Element = NULL;
1605 
1606  if (!strcmp(fe_name, "BiCubic2DFiniteElement") ||
1607  !strcmp(fe_name, "Quad_Q3"))
1608  {
1609  GeomType = Geometry::SQUARE;
1610  Local_Element = new BiCubic2DFiniteElement;
1611  }
1612  else if (!strcmp(fe_name, "Nedelec1HexFiniteElement") ||
1613  !strcmp(fe_name, "Hex_ND1"))
1614  {
1615  GeomType = Geometry::CUBE;
1616  Local_Element = new Nedelec1HexFiniteElement;
1617  }
1618  else if (!strncmp(fe_name, "H1_", 3))
1619  {
1620  GeomType = Geometry::SQUARE;
1621  Local_Element = new H1_QuadrilateralElement(atoi(fe_name + 7));
1622  }
1623  else if (!strncmp(fe_name, "L2_", 3))
1624  {
1625  GeomType = Geometry::SQUARE;
1626  Local_Element = new L2_QuadrilateralElement(atoi(fe_name + 7));
1627  }
1628  else
1629  {
1630  cerr << "Local_FECollection::Local_FECollection : fe_name = "
1631  << fe_name << endl;
1632  mfem_error();
1633  }
1634 }
1635 
1636 
1637 void NURBSFECollection::Allocate(int Order)
1638 {
1639  SegmentFE = new NURBS1DFiniteElement(Order);
1640  QuadrilateralFE = new NURBS2DFiniteElement(Order);
1641  ParallelepipedFE = new NURBS3DFiniteElement(Order);
1642 
1643  snprintf(name, 16, "NURBS%i", Order);
1644 }
1645 
1646 void NURBSFECollection::Deallocate()
1647 {
1648  delete ParallelepipedFE;
1649  delete QuadrilateralFE;
1650  delete SegmentFE;
1651 }
1652 
1653 const FiniteElement *
1655 {
1656  switch (GeomType)
1657  {
1658  case Geometry::SEGMENT: return SegmentFE;
1659  case Geometry::SQUARE: return QuadrilateralFE;
1660  case Geometry::CUBE: return ParallelepipedFE;
1661  default:
1662  mfem_error ("NURBSFECollection: unknown geometry type.");
1663  }
1664  return SegmentFE; // Make some compilers happy
1665 }
1666 
1667 int NURBSFECollection::DofForGeometry(int GeomType) const
1668 {
1669  mfem_error("NURBSFECollection::DofForGeometry");
1670  return 0; // Make some compilers happy
1671 }
1672 
1673 int *NURBSFECollection::DofOrderForOrientation(int GeomType, int Or) const
1674 {
1675  mfem_error("NURBSFECollection::DofOrderForOrientation");
1676  return NULL;
1677 }
1678 
1679 }
Abstract class for Finite Elements.
Definition: fe.hpp:42
Arbitrary order non-uniform rational B-splines (NURBS) finite elements.
Definition: fe_coll.hpp:180
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:829
Version of QuadraticDiscont2DFECollection with dofs in the Gaussian points.
Definition: fe_coll.hpp:531
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:933
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:210
static const int NumGeom
Definition: geom.hpp:34
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:1673
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:992
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:757
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:675
virtual ~L2_FECollection()
Definition: fe_coll.cpp:1281
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:296
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:917
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:520
void InitFaces(const int p, const int dim, const int map_type)
Definition: fe_coll.cpp:1323
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:578
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:556
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:486
Piecewise-(bi)linear continuous finite elements.
Definition: fe_coll.hpp:221
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:1020
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:857
ND_FECollection(const int p, const int dim)
Definition: fe_coll.cpp:1446
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:400
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:864
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:960
L2_FECollection(const int p, const int dim, const int type=0)
Definition: fe_coll.cpp:1175
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:709
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:695
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:478
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:896
RT_FECollection(const int p, const int dim, const int map_type)
Definition: fe_coll.hpp:131
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:793
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:425
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:646
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:373
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:134
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:157
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:412
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:548
virtual ~RT_FECollection()
Definition: fe_coll.cpp:1425
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:464
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:1406
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:976
Piecewise-(bi)cubic continuous finite elements.
Definition: fe_coll.hpp:288
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:512
PointFiniteElement PointFE
Definition: point.cpp:28
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:1005
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:223
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:498
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:150
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:565
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:769
TriLinear3DFiniteElement HexahedronFE
Definition: hexahedron.cpp:60
Discontinuous collection defined locally by a given finite element.
Definition: fe_coll.hpp:724
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:198
Version of QuadraticDiscont2DFECollection with positive basis functions.
Definition: fe_coll.hpp:515
virtual ~ND_FECollection()
Definition: fe_coll.cpp:1590
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:785
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:333
RT_Trace_FECollection(const int p, const int dim, const int map_type=FiniteElement::INTEGRAL)
Definition: fe_coll.cpp:1434
H1_FECollection(const int p, const int dim=3, const int type=0)
Definition: fe_coll.cpp:1039
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:598
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:452
Version of QuadraticFECollection with positive basis functions.
Definition: fe_coll.hpp:268
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:717
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:880
virtual ~H1_FECollection()
Definition: fe_coll.cpp:1165
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:118
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:946
Piecewise-linear nonconforming finite elements in 3D.
Definition: fe_coll.hpp:332
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:265
Crouzeix-Raviart nonconforming elements in 2D.
Definition: fe_coll.hpp:312
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:119
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:323
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:309
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:386
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:805
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:612
void mfem_error(const char *msg)
Definition: error.cpp:23
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:620
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:189
Linear3DFiniteElement TetrahedronFE
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:533
Linear2DFiniteElement TriangleFE
Definition: triangle.cpp:127
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:841
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:683
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:660
Version of LinearDiscont2DFECollection with dofs in the Gaussian points.
Definition: fe_coll.hpp:459
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:1654
static FiniteElementCollection * New(const char *name)
Definition: fe_coll.cpp:38
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:821
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:905
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:173
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:439
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:1563
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:1146
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:360
Piecewise-(bi)quadratic continuous finite elements.
Definition: fe_coll.hpp:244
int HasFaceDofs(int GeomType) const
Definition: fe_coll.cpp:25
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:158
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:747
virtual int * DofOrderForOrientation(int GeomType, int Or) const
Definition: fe_coll.cpp:1266
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:52
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:1667
Local_FECollection(const char *fe_name)
Definition: fe_coll.cpp:1600
BiLinear2DFiniteElement QuadrilateralFE
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:631
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:731
Linear (P1) finite elements on quadrilaterals.
Definition: fe_coll.hpp:480
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:249
Linear1DFiniteElement SegmentFE
Definition: segment.cpp:41
virtual int DofForGeometry(int GeomType) const
Definition: fe_coll.cpp:346
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:233
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:83
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.cpp:586