MFEM  v4.6.0
Finite element discretization library
fe_coll.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
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 const FiniteElement *
27 {
28  ErrorMode save_error_mode = error_mode;
29  error_mode = RETURN_NULL;
30  const FiniteElement *fe = nullptr;
31  for (int g = Geometry::DimStart[dim]; g < Geometry::DimStart[dim+1]; g++)
32  {
33  fe = FiniteElementForGeometry((Geometry::Type)g);
34  if (fe != nullptr) { break; }
35  }
36  error_mode = save_error_mode;
37  return fe;
38 }
39 
41 {
42  const FiniteElement *fe = FiniteElementForDim(dim);
43  if (fe)
44  {
45  return fe->GetRangeType();
46  }
48 }
49 
51 {
52  const FiniteElement *fe = FiniteElementForDim(dim);
53  if (fe)
54  {
55  return fe->GetDerivRangeType();
56  }
58 }
59 
61 {
62  const FiniteElement *fe = FiniteElementForDim(dim);
63  if (fe)
64  {
65  return fe->GetMapType();
66  }
68 }
69 
71 {
72  const FiniteElement *fe = FiniteElementForDim(dim);
73  if (fe)
74  {
75  return fe->GetDerivType();
76  }
77  return FiniteElement::NONE;
78 }
79 
81 {
82  const FiniteElement *fe = FiniteElementForDim(dim);
83  if (fe)
84  {
85  return fe->GetDerivMapType();
86  }
88 }
89 
91 {
92  switch (geom)
93  {
95  return GetNumDof(Geometry::TRIANGLE, p);
96  case Geometry::CUBE:
97  return GetNumDof(Geometry::SQUARE, p);
98  case Geometry::PRISM:
99  return max(GetNumDof(Geometry::TRIANGLE, p),
100  GetNumDof(Geometry::SQUARE, p));
101  case Geometry::PYRAMID:
102  return max(GetNumDof(Geometry::TRIANGLE, p),
103  GetNumDof(Geometry::SQUARE, p));
104  default:
105  MFEM_ABORT("unknown geometry type");
106  }
107  return 0;
108 }
109 
111 {
112  MFEM_ABORT("this method is not implemented in this derived class!");
113  return NULL;
114 }
115 
117 {
118  FiniteElementCollection *fec = NULL;
119 
120  if (!strcmp(name, "Linear"))
121  {
122  fec = new LinearFECollection;
123  }
124  else if (!strcmp(name, "Quadratic"))
125  {
126  fec = new QuadraticFECollection;
127  }
128  else if (!strcmp(name, "QuadraticPos"))
129  {
130  fec = new QuadraticPosFECollection;
131  }
132  else if (!strcmp(name, "Cubic"))
133  {
134  fec = new CubicFECollection;
135  }
136  else if (!strcmp(name, "Const3D"))
137  {
138  fec = new Const3DFECollection;
139  }
140  else if (!strcmp(name, "Const2D"))
141  {
142  fec = new Const2DFECollection;
143  }
144  else if (!strcmp(name, "LinearDiscont2D"))
145  {
146  fec = new LinearDiscont2DFECollection;
147  }
148  else if (!strcmp(name, "GaussLinearDiscont2D"))
149  {
151  }
152  else if (!strcmp(name, "P1OnQuad"))
153  {
154  fec = new P1OnQuadFECollection;
155  }
156  else if (!strcmp(name, "QuadraticDiscont2D"))
157  {
159  }
160  else if (!strcmp(name, "QuadraticPosDiscont2D"))
161  {
163  }
164  else if (!strcmp(name, "GaussQuadraticDiscont2D"))
165  {
167  }
168  else if (!strcmp(name, "CubicDiscont2D"))
169  {
170  fec = new CubicDiscont2DFECollection;
171  }
172  else if (!strcmp(name, "LinearDiscont3D"))
173  {
174  fec = new LinearDiscont3DFECollection;
175  }
176  else if (!strcmp(name, "QuadraticDiscont3D"))
177  {
179  }
180  else if (!strcmp(name, "LinearNonConf3D"))
181  {
182  fec = new LinearNonConf3DFECollection;
183  }
184  else if (!strcmp(name, "CrouzeixRaviart"))
185  {
186  fec = new CrouzeixRaviartFECollection;
187  }
188  else if (!strcmp(name, "ND1_3D"))
189  {
190  fec = new ND1_3DFECollection;
191  }
192  else if (!strcmp(name, "RT0_2D"))
193  {
194  fec = new RT0_2DFECollection;
195  }
196  else if (!strcmp(name, "RT1_2D"))
197  {
198  fec = new RT1_2DFECollection;
199  }
200  else if (!strcmp(name, "RT2_2D"))
201  {
202  fec = new RT2_2DFECollection;
203  }
204  else if (!strcmp(name, "RT0_3D"))
205  {
206  fec = new RT0_3DFECollection;
207  }
208  else if (!strcmp(name, "RT1_3D"))
209  {
210  fec = new RT1_3DFECollection;
211  }
212  else if (!strncmp(name, "H1_Trace_", 9))
213  {
214  fec = new H1_Trace_FECollection(atoi(name + 13), atoi(name + 9));
215  }
216  else if (!strncmp(name, "H1_Trace@", 9))
217  {
218  fec = new H1_Trace_FECollection(atoi(name + 15), atoi(name + 11),
219  BasisType::GetType(name[9]));
220  }
221  else if (!strncmp(name, "H1_", 3))
222  {
223  fec = new H1_FECollection(atoi(name + 7), atoi(name + 3));
224  }
225  else if (!strncmp(name, "H1Pos_Trace_", 12))
226  {
227  fec = new H1_Trace_FECollection(atoi(name + 16), atoi(name + 12),
229  }
230  else if (!strncmp(name, "H1Pos_", 6))
231  {
232  fec = new H1Pos_FECollection(atoi(name + 10), atoi(name + 6));
233  }
234  else if (!strncmp(name, "H1Ser_", 6))
235  {
236  fec = new H1Ser_FECollection(atoi(name + 10), atoi(name + 6));
237  }
238  else if (!strncmp(name, "H1@", 3))
239  {
240  fec = new H1_FECollection(atoi(name + 9), atoi(name + 5),
241  BasisType::GetType(name[3]));
242  }
243  else if (!strncmp(name, "L2_T", 4))
244  fec = new L2_FECollection(atoi(name + 10), atoi(name + 6),
245  atoi(name + 4));
246  else if (!strncmp(name, "L2_", 3))
247  {
248  fec = new L2_FECollection(atoi(name + 7), atoi(name + 3));
249  }
250  else if (!strncmp(name, "L2Int_T", 7))
251  {
252  fec = new L2_FECollection(atoi(name + 13), atoi(name + 9),
253  atoi(name + 7), FiniteElement::INTEGRAL);
254  }
255  else if (!strncmp(name, "L2Int_", 6))
256  {
257  fec = new L2_FECollection(atoi(name + 10), atoi(name + 6),
260  }
261  else if (!strncmp(name, "RT_Trace_", 9))
262  {
263  fec = new RT_Trace_FECollection(atoi(name + 13), atoi(name + 9));
264  }
265  else if (!strncmp(name, "RT_ValTrace_", 12))
266  {
267  fec = new RT_Trace_FECollection(atoi(name + 16), atoi(name + 12),
269  }
270  else if (!strncmp(name, "RT_Trace@", 9))
271  {
272  fec = new RT_Trace_FECollection(atoi(name + 15), atoi(name + 11),
274  BasisType::GetType(name[9]));
275  }
276  else if (!strncmp(name, "RT_ValTrace@", 12))
277  {
278  fec = new RT_Trace_FECollection(atoi(name + 18), atoi(name + 14),
280  BasisType::GetType(name[12]));
281  }
282  else if (!strncmp(name, "DG_Iface_", 9))
283  {
284  fec = new DG_Interface_FECollection(atoi(name + 13), atoi(name + 9));
285  }
286  else if (!strncmp(name, "DG_Iface@", 9))
287  {
288  fec = new DG_Interface_FECollection(atoi(name + 15), atoi(name + 11),
290  BasisType::GetType(name[9]));
291  }
292  else if (!strncmp(name, "DG_IntIface_", 12))
293  {
294  fec = new DG_Interface_FECollection(atoi(name + 16), atoi(name + 12),
296  }
297  else if (!strncmp(name, "DG_IntIface@", 12))
298  {
299  fec = new DG_Interface_FECollection(atoi(name + 18), atoi(name + 14),
301  BasisType::GetType(name[12]));
302  }
303  else if (!strncmp(name, "RT_", 3))
304  {
305  fec = new RT_FECollection(atoi(name + 7), atoi(name + 3));
306  }
307  else if (!strncmp(name, "RT@", 3))
308  {
309  fec = new RT_FECollection(atoi(name + 10), atoi(name + 6),
310  BasisType::GetType(name[3]),
311  BasisType::GetType(name[4]));
312  }
313  else if (!strncmp(name, "ND_Trace_", 9))
314  {
315  fec = new ND_Trace_FECollection(atoi(name + 13), atoi(name + 9));
316  }
317  else if (!strncmp(name, "ND_Trace@", 9))
318  {
319  fec = new ND_Trace_FECollection(atoi(name + 16), atoi(name + 12),
320  BasisType::GetType(name[9]),
321  BasisType::GetType(name[10]));
322  }
323  else if (!strncmp(name, "ND_", 3))
324  {
325  fec = new ND_FECollection(atoi(name + 7), atoi(name + 3));
326  }
327  else if (!strncmp(name, "ND@", 3))
328  {
329  fec = new ND_FECollection(atoi(name + 10), atoi(name + 6),
330  BasisType::GetType(name[3]),
331  BasisType::GetType(name[4]));
332  }
333  else if (!strncmp(name, "Local_", 6))
334  {
335  fec = new Local_FECollection(name + 6);
336  }
337  else if (!strncmp(name, "NURBS", 5))
338  {
339  if (name[5] != '\0')
340  {
341  // "NURBS" + "number" --> fixed order nurbs collection
342  fec = new NURBSFECollection(atoi(name + 5));
343  }
344  else
345  {
346  // "NURBS" --> variable order nurbs collection
347  fec = new NURBSFECollection();
348  }
349  }
350  else
351  {
352  MFEM_ABORT("unknown FiniteElementCollection: " << name);
353  }
354  MFEM_VERIFY(!strcmp(fec->Name(), name), "input name: \"" << name
355  << "\" does not match the created collection name: \""
356  << fec->Name() << '"');
357 
358  return fec;
359 }
360 
362 {
363  // default implementation for collections that don't care about variable p
364  MFEM_ABORT("Collection " << Name() << " does not support variable orders.");
365  (void) p;
366  return NULL;
367 }
368 
370 {
371  if (p >= var_orders.Size())
372  {
373  var_orders.SetSize(p+1, NULL);
374  }
375  var_orders[p] = Clone(p);
376 }
377 
379 {
380  for (int i = 0; i < var_orders.Size(); i++)
381  {
382  delete var_orders[i];
383  }
384 }
385 
386 template <Geometry::Type geom>
387 inline void FiniteElementCollection::GetNVE(int &nv, int &ne)
388 {
389  typedef typename Geometry::Constants<geom> g_consts;
390 
391  nv = g_consts::NumVert;
392  ne = g_consts::NumEdges;
393 }
394 
395 template <Geometry::Type geom, typename v_t>
396 inline void FiniteElementCollection::
397 GetEdge(int &nv, v_t &v, int &ne, int &e, int &eo, const int edge_info)
398 {
399  typedef typename Geometry::Constants<Geometry::SEGMENT> e_consts;
400  typedef typename Geometry::Constants<geom> g_consts;
401 
402  nv = e_consts::NumVert;
403  ne = 1;
404  e = edge_info/64;
405  eo = edge_info%64;
406  MFEM_ASSERT(0 <= e && e < g_consts::NumEdges, "");
407  MFEM_ASSERT(0 <= eo && eo < e_consts::NumOrient, "");
408  v[0] = e_consts::Orient[eo][0];
409  v[1] = e_consts::Orient[eo][1];
410  v[0] = g_consts::Edges[e][v[0]];
411  v[1] = g_consts::Edges[e][v[1]];
412 }
413 
414 template <Geometry::Type geom, Geometry::Type f_geom,
415  typename v_t, typename e_t, typename eo_t>
416 inline void FiniteElementCollection::
417 GetFace(int &nv, v_t &v, int &ne, e_t &e, eo_t &eo,
418  int &nf, int &f, Geometry::Type &fg, int &fo, const int face_info)
419 {
420  typedef typename Geometry::Constants< geom> g_consts;
421  typedef typename Geometry::Constants<f_geom> f_consts;
422 
423  nv = f_consts::NumVert;
424  nf = 1;
425  f = face_info/64;
426  fg = f_geom;
427  fo = face_info%64;
428  MFEM_ASSERT(0 <= f && f < g_consts::NumFaces, "");
429  MFEM_ASSERT(0 <= fo && fo < f_consts::NumOrient, "");
430  for (int i = 0; i < f_consts::NumVert; i++)
431  {
432  v[i] = f_consts::Orient[fo][i];
433  v[i] = g_consts::FaceVert[f][v[i]];
434  }
435  ne = f_consts::NumEdges;
436  for (int i = 0; i < f_consts::NumEdges; i++)
437  {
438  int v0 = v[f_consts::Edges[i][0]];
439  int v1 = v[f_consts::Edges[i][1]];
440  int eor = 0;
441  if (v0 > v1) { swap(v0, v1); eor = 1; }
442  for (int j = g_consts::VertToVert::I[v0]; true; j++)
443  {
444  MFEM_ASSERT(j < g_consts::VertToVert::I[v0+1],
445  "internal error, edge not found");
446  if (v1 == g_consts::VertToVert::J[j][0])
447  {
448  int en = g_consts::VertToVert::J[j][1];
449  if (en < 0)
450  {
451  en = -1-en;
452  eor = 1-eor;
453  }
454  e[i] = en;
455  eo[i] = eor;
456  break;
457  }
458  }
459  }
460 }
461 
463  int Info,
464  Array<int> &dofs) const
465 {
466  // Info = 64 * SubIndex + SubOrientation
467  MFEM_ASSERT(0 <= Geom && Geom < Geometry::NumGeom,
468  "invalid Geom = " << Geom);
469  MFEM_ASSERT(0 <= SDim && SDim <= Geometry::Dimension[Geom],
470  "invalid SDim = " << SDim <<
471  " for Geom = " << Geometry::Name[Geom]);
472 
473  const int nvd = DofForGeometry(Geometry::POINT);
474  if (SDim == 0) // vertex
475  {
476  const int off = nvd*(Info/64);
477  dofs.SetSize(nvd);
478  for (int i = 0; i < nvd; i++)
479  {
480  dofs[i] = off + i;
481  }
482  }
483  else
484  {
485  int v[4], e[4], eo[4], f[1], fo[1];
486  int av = 0, nv = 0, ae = 0, ne = 0, nf = 0;
487  Geometry::Type fg[1];
488 
489  switch (Geom)
490  {
491  case Geometry::SEGMENT:
492  {
493  GetNVE<Geometry::SEGMENT>(av, ae);
494  GetEdge<Geometry::SEGMENT>(nv, v, ne, e[0], eo[0], Info);
495  break;
496  }
497 
498  case Geometry::TRIANGLE:
499  {
500  GetNVE<Geometry::TRIANGLE>(av, ae);
501  switch (SDim)
502  {
503  case 1:
504  GetEdge<Geometry::TRIANGLE>(nv, v, ne, e[0], eo[0], Info);
505  break;
506  case 2:
507  GetFace<Geometry::TRIANGLE,Geometry::TRIANGLE>(
508  nv, v, ne, e, eo, nf, f[0], fg[0], fo[0], Info);
509  break;
510  default:
511  goto not_supp;
512  }
513  break;
514  }
515 
516  case Geometry::SQUARE:
517  {
518  GetNVE<Geometry::SQUARE>(av, ae);
519  switch (SDim)
520  {
521  case 1:
522  GetEdge<Geometry::SQUARE>(nv, v, ne, e[0], eo[0], Info);
523  break;
524  case 2:
525  GetFace<Geometry::SQUARE,Geometry::SQUARE>(
526  nv, v, ne, e, eo, nf, f[0], fg[0], fo[0], Info);
527  break;
528  default:
529  goto not_supp;
530  }
531  break;
532  }
533 
535  {
536  GetNVE<Geometry::TETRAHEDRON>(av, ae);
537  switch (SDim)
538  {
539  case 1:
540  GetEdge<Geometry::TETRAHEDRON>(nv, v, ne, e[0], eo[0], Info);
541  break;
542  case 2:
543  GetFace<Geometry::TETRAHEDRON,Geometry::TRIANGLE>(
544  nv, v, ne, e, eo, nf, f[0], fg[0], fo[0], Info);
545  break;
546  default:
547  goto not_supp;
548  }
549  break;
550  }
551 
552  case Geometry::CUBE:
553  {
554  GetNVE<Geometry::CUBE>(av, ae);
555  switch (SDim)
556  {
557  case 1:
558  GetEdge<Geometry::CUBE>(nv, v, ne, e[0], eo[0], Info);
559  break;
560  case 2:
561  GetFace<Geometry::CUBE,Geometry::SQUARE>(
562  nv, v, ne, e, eo, nf, f[0], fg[0], fo[0], Info);
563  break;
564  default:
565  goto not_supp;
566  }
567  break;
568  }
569 
570  default:
571  MFEM_ABORT("invalid Geom = " << Geom);
572  }
573 
574  int ned = (ne > 0) ? DofForGeometry(Geometry::SEGMENT) : 0;
575 
576  // add vertex dofs
577  dofs.SetSize(nv*nvd+ne*ned);
578  for (int i = 0; i < nv; i++)
579  {
580  for (int j = 0; j < nvd; j++)
581  {
582  dofs[i*nvd+j] = v[i]*nvd+j;
583  }
584  }
585  int l_off = nv*nvd, g_off = av*nvd;
586 
587  // add edge dofs
588  if (ned > 0)
589  {
590  for (int i = 0; i < ne; i++)
591  {
592  const int *ed = DofOrderForOrientation(Geometry::SEGMENT,
593  eo[i] ? -1 : 1);
594  for (int j = 0; j < ned; j++)
595  {
596  dofs[l_off+i*ned+j] =
597  ed[j] >= 0 ?
598  g_off+e[i]*ned+ed[j] :
599  -1-(g_off+e[i]*ned+(-1-ed[j]));
600  }
601  }
602  l_off += ne*ned;
603  g_off += ae*ned;
604  }
605 
606  // add face dofs
607  if (nf > 0)
608  {
609  const int nfd = DofForGeometry(fg[0]); // assume same face geometry
610  dofs.SetSize(dofs.Size()+nf*nfd);
611  for (int i = 0; i < nf; i++)
612  {
613  const int *fd = DofOrderForOrientation(fg[i], fo[i]);
614  for (int j = 0; j < nfd; j++)
615  {
616  dofs[l_off+i*nfd+j] =
617  fd[j] >= 0 ?
618  g_off+f[i]*nfd+fd[j] :
619  -1-(g_off+f[i]*nfd+(-1-fd[j]));
620  }
621  }
622  }
623 
624  // add volume dofs ...
625  }
626  return;
627 
628 not_supp:
629  MFEM_ABORT("Geom = " << Geometry::Name[Geom] <<
630  ", SDim = " << SDim << " is not supported");
631 }
632 
633 const FiniteElement *
635 {
636  switch (GeomType)
637  {
638  case Geometry::POINT: return &PointFE;
639  case Geometry::SEGMENT: return &SegmentFE;
640  case Geometry::TRIANGLE: return &TriangleFE;
641  case Geometry::SQUARE: return &QuadrilateralFE;
642  case Geometry::TETRAHEDRON: return &TetrahedronFE;
643  case Geometry::CUBE: return &ParallelepipedFE;
644  case Geometry::PRISM: return &WedgeFE;
645  case Geometry::PYRAMID: return &PyramidFE;
646  default:
647  if (error_mode == RETURN_NULL) { return nullptr; }
648  mfem_error ("LinearFECollection: unknown geometry type.");
649  }
650  return &SegmentFE; // Make some compilers happy
651 }
652 
654 {
655  switch (GeomType)
656  {
657  case Geometry::POINT: return 1;
658  case Geometry::SEGMENT: return 0;
659  case Geometry::TRIANGLE: return 0;
660  case Geometry::SQUARE: return 0;
661  case Geometry::TETRAHEDRON: return 0;
662  case Geometry::CUBE: return 0;
663  case Geometry::PRISM: return 0;
664  case Geometry::PYRAMID: return 0;
665  default:
666  mfem_error ("LinearFECollection: unknown geometry type.");
667  }
668  return 0; // Make some compilers happy
669 }
670 
672  int Or) const
673 {
674  return NULL;
675 }
676 
677 
678 const FiniteElement *
680 {
681  switch (GeomType)
682  {
683  case Geometry::POINT: return &PointFE;
684  case Geometry::SEGMENT: return &SegmentFE;
685  case Geometry::TRIANGLE: return &TriangleFE;
686  case Geometry::SQUARE: return &QuadrilateralFE;
687  case Geometry::TETRAHEDRON: return &TetrahedronFE;
688  case Geometry::CUBE: return &ParallelepipedFE;
689  case Geometry::PRISM: return &WedgeFE;
690  default:
691  if (error_mode == RETURN_NULL) { return nullptr; }
692  mfem_error ("QuadraticFECollection: unknown geometry type.");
693  }
694  return &SegmentFE; // Make some compilers happy
695 }
696 
698 {
699  switch (GeomType)
700  {
701  case Geometry::POINT: return 1;
702  case Geometry::SEGMENT: return 1;
703  case Geometry::TRIANGLE: return 0;
704  case Geometry::SQUARE: return 1;
705  case Geometry::TETRAHEDRON: return 0;
706  case Geometry::CUBE: return 1;
707  case Geometry::PRISM: return 0;
708  default:
709  mfem_error ("QuadraticFECollection: unknown geometry type.");
710  }
711  return 0; // Make some compilers happy
712 }
713 
715  Geometry::Type GeomType, int Or) const
716 {
717  static int indexes[] = { 0 };
718 
719  return indexes;
720 }
721 
722 
723 const FiniteElement *
725  Geometry::Type GeomType) const
726 {
727  switch (GeomType)
728  {
729  case Geometry::SEGMENT: return &SegmentFE;
730  case Geometry::SQUARE: return &QuadrilateralFE;
731  default:
732  if (error_mode == RETURN_NULL) { return nullptr; }
733  mfem_error ("QuadraticPosFECollection: unknown geometry type.");
734  }
735  return NULL; // Make some compilers happy
736 }
737 
739 {
740  switch (GeomType)
741  {
742  case Geometry::POINT: return 1;
743  case Geometry::SEGMENT: return 1;
744  case Geometry::SQUARE: return 1;
745  default:
746  mfem_error ("QuadraticPosFECollection: unknown geometry type.");
747  }
748  return 0; // Make some compilers happy
749 }
750 
752  Geometry::Type GeomType, int Or) const
753 {
754  static int indexes[] = { 0 };
755 
756  return indexes;
757 }
758 
759 
760 const FiniteElement *
762 {
763  switch (GeomType)
764  {
765  case Geometry::POINT: return &PointFE;
766  case Geometry::SEGMENT: return &SegmentFE;
767  case Geometry::TRIANGLE: return &TriangleFE;
768  case Geometry::SQUARE: return &QuadrilateralFE;
769  case Geometry::TETRAHEDRON: return &TetrahedronFE;
770  case Geometry::CUBE: return &ParallelepipedFE;
771  case Geometry::PRISM: return &WedgeFE;
772  default:
773  if (error_mode == RETURN_NULL) { return nullptr; }
774  mfem_error ("CubicFECollection: unknown geometry type.");
775  }
776  return &SegmentFE; // Make some compilers happy
777 }
778 
780 {
781  switch (GeomType)
782  {
783  case Geometry::POINT: return 1;
784  case Geometry::SEGMENT: return 2;
785  case Geometry::TRIANGLE: return 1;
786  case Geometry::SQUARE: return 4;
787  case Geometry::TETRAHEDRON: return 0;
788  case Geometry::CUBE: return 8;
789  case Geometry::PRISM: return 2;
790  default:
791  mfem_error ("CubicFECollection: unknown geometry type.");
792  }
793  return 0; // Make some compilers happy
794 }
795 
797  int Or) const
798 {
799  if (GeomType == Geometry::SEGMENT)
800  {
801  static int ind_pos[] = { 0, 1 };
802  static int ind_neg[] = { 1, 0 };
803 
804  if (Or < 0)
805  {
806  return ind_neg;
807  }
808  return ind_pos;
809  }
810  else if (GeomType == Geometry::TRIANGLE)
811  {
812  static int indexes[] = { 0 };
813 
814  return indexes;
815  }
816  else if (GeomType == Geometry::SQUARE)
817  {
818  static int sq_ind[8][4] = {{0, 1, 2, 3}, {0, 2, 1, 3},
819  {2, 0, 3, 1}, {1, 0, 3, 2},
820  {3, 2, 1, 0}, {3, 1, 2, 0},
821  {1, 3, 0, 2}, {2, 3, 0, 1}
822  };
823  return sq_ind[Or];
824  }
825 
826  return NULL;
827 }
828 
829 
830 const FiniteElement *
832  Geometry::Type GeomType) const
833 {
834  switch (GeomType)
835  {
836  case Geometry::SEGMENT: return &SegmentFE;
837  case Geometry::TRIANGLE: return &TriangleFE;
838  case Geometry::SQUARE: return &QuadrilateralFE;
839  default:
840  if (error_mode == RETURN_NULL) { return nullptr; }
841  mfem_error ("CrouzeixRaviartFECollection: unknown geometry type.");
842  }
843  return &SegmentFE; // Make some compilers happy
844 }
845 
847 {
848  switch (GeomType)
849  {
850  case Geometry::POINT: return 0;
851  case Geometry::SEGMENT: return 1;
852  case Geometry::TRIANGLE: return 0;
853  case Geometry::SQUARE: return 0;
854  default:
855  mfem_error ("CrouzeixRaviartFECollection: unknown geometry type.");
856  }
857  return 0; // Make some compilers happy
858 }
859 
861  Geometry::Type GeomType, int Or) const
862 {
863  static int indexes[] = { 0 };
864 
865  return indexes;
866 }
867 
868 
869 const FiniteElement *
871 {
872  switch (GeomType)
873  {
874  case Geometry::SEGMENT: return &SegmentFE;
875  case Geometry::TRIANGLE: return &TriangleFE;
876  case Geometry::SQUARE: return &QuadrilateralFE;
877  default:
878  if (error_mode == RETURN_NULL) { return nullptr; }
879  mfem_error ("RT0_2DFECollection: unknown geometry type.");
880  }
881  return &SegmentFE; // Make some compilers happy
882 }
883 
885 {
886  switch (GeomType)
887  {
888  case Geometry::POINT: return 0;
889  case Geometry::SEGMENT: return 1;
890  case Geometry::TRIANGLE: return 0;
891  case Geometry::SQUARE: return 0;
892  default:
893  mfem_error ("RT0_2DFECollection: unknown geometry type.");
894  }
895  return 0; // Make some compilers happy
896 }
897 
899  int Or) const
900 {
901  static int ind_pos[] = { 0 };
902  static int ind_neg[] = { -1 };
903 
904  if (Or > 0)
905  {
906  return ind_pos;
907  }
908  return ind_neg;
909 }
910 
911 
912 const FiniteElement *
914 {
915  switch (GeomType)
916  {
917  case Geometry::SEGMENT: return &SegmentFE;
918  case Geometry::TRIANGLE: return &TriangleFE;
919  case Geometry::SQUARE: return &QuadrilateralFE;
920  default:
921  if (error_mode == RETURN_NULL) { return nullptr; }
922  mfem_error ("RT1_2DFECollection: unknown geometry type.");
923  }
924  return &SegmentFE; // Make some compilers happy
925 }
926 
928 {
929  switch (GeomType)
930  {
931  case Geometry::POINT: return 0;
932  case Geometry::SEGMENT: return 2;
933  case Geometry::TRIANGLE: return 2;
934  case Geometry::SQUARE: return 4;
935  default:
936  mfem_error ("RT1_2DFECollection: unknown geometry type.");
937  }
938  return 0; // Make some compilers happy
939 }
940 
942  int Or) const
943 {
944  static int ind_pos[] = { 0, 1 };
945  static int ind_neg[] = { -2, -1 };
946 
947  if (Or > 0)
948  {
949  return ind_pos;
950  }
951  return ind_neg;
952 }
953 
954 const FiniteElement *
956 {
957  switch (GeomType)
958  {
959  case Geometry::SEGMENT: return &SegmentFE;
960  case Geometry::TRIANGLE: return &TriangleFE;
961  case Geometry::SQUARE: return &QuadrilateralFE;
962  default:
963  if (error_mode == RETURN_NULL) { return nullptr; }
964  mfem_error ("RT2_2DFECollection: unknown geometry type.");
965  }
966  return &SegmentFE; // Make some compilers happy
967 }
968 
970 {
971  switch (GeomType)
972  {
973  case Geometry::POINT: return 0;
974  case Geometry::SEGMENT: return 3;
975  case Geometry::TRIANGLE: return 6;
976  case Geometry::SQUARE: return 12;
977  default:
978  mfem_error ("RT2_2DFECollection: unknown geometry type.");
979  }
980  return 0; // Make some compilers happy
981 }
982 
984  int Or) const
985 {
986  static int ind_pos[] = { 0, 1, 2 };
987  static int ind_neg[] = { -3, -2, -1 };
988 
989  if (Or > 0)
990  {
991  return ind_pos;
992  }
993  return ind_neg;
994 }
995 
996 
997 const FiniteElement *
999 {
1000  switch (GeomType)
1001  {
1002  case Geometry::TRIANGLE: return &TriangleFE;
1003  case Geometry::SQUARE: return &QuadrilateralFE;
1004  default:
1005  if (error_mode == RETURN_NULL) { return nullptr; }
1006  mfem_error ("Const2DFECollection: unknown geometry type.");
1007  }
1008  return &TriangleFE; // Make some compilers happy
1009 }
1010 
1012 {
1013  switch (GeomType)
1014  {
1015  case Geometry::POINT: return 0;
1016  case Geometry::SEGMENT: return 0;
1017  case Geometry::TRIANGLE: return 1;
1018  case Geometry::SQUARE: return 1;
1019  default:
1020  mfem_error ("Const2DFECollection: unknown geometry type.");
1021  }
1022  return 0; // Make some compilers happy
1023 }
1024 
1026  int Or) const
1027 {
1028  return NULL;
1029 }
1030 
1031 
1032 const FiniteElement *
1034  Geometry::Type GeomType) const
1035 {
1036  switch (GeomType)
1037  {
1038  case Geometry::TRIANGLE: return &TriangleFE;
1039  case Geometry::SQUARE: return &QuadrilateralFE;
1040  default:
1041  if (error_mode == RETURN_NULL) { return nullptr; }
1042  mfem_error ("LinearDiscont2DFECollection: unknown geometry type.");
1043  }
1044  return &TriangleFE; // Make some compilers happy
1045 }
1046 
1048 {
1049  switch (GeomType)
1050  {
1051  case Geometry::POINT: return 0;
1052  case Geometry::SEGMENT: return 0;
1053  case Geometry::TRIANGLE: return 3;
1054  case Geometry::SQUARE: return 4;
1055  default:
1056  mfem_error ("LinearDiscont2DFECollection: unknown geometry type.");
1057  }
1058  return 0; // Make some compilers happy
1059 }
1060 
1062  Geometry::Type GeomType, int Or) const
1063 {
1064  return NULL;
1065 }
1066 
1067 
1068 const FiniteElement *
1070  Geometry::Type GeomType) const
1071 {
1072  switch (GeomType)
1073  {
1074  case Geometry::TRIANGLE: return &TriangleFE;
1075  case Geometry::SQUARE: return &QuadrilateralFE;
1076  default:
1077  if (error_mode == RETURN_NULL) { return nullptr; }
1078  mfem_error ("GaussLinearDiscont2DFECollection:"
1079  " unknown geometry type.");
1080  }
1081  return &TriangleFE; // Make some compilers happy
1082 }
1083 
1085 const
1086 {
1087  switch (GeomType)
1088  {
1089  case Geometry::POINT: return 0;
1090  case Geometry::SEGMENT: return 0;
1091  case Geometry::TRIANGLE: return 3;
1092  case Geometry::SQUARE: return 4;
1093  default:
1094  mfem_error ("GaussLinearDiscont2DFECollection:"
1095  " unknown geometry type.");
1096  }
1097  return 0; // Make some compilers happy
1098 }
1099 
1101  Geometry::Type GeomType, int Or) const
1102 {
1103  return NULL;
1104 }
1105 
1106 
1107 const FiniteElement *
1109 {
1110  if (GeomType != Geometry::SQUARE)
1111  {
1112  if (error_mode == RETURN_NULL) { return nullptr; }
1113  mfem_error ("P1OnQuadFECollection: unknown geometry type.");
1114  }
1115  return &QuadrilateralFE;
1116 }
1117 
1119 {
1120  switch (GeomType)
1121  {
1122  case Geometry::POINT: return 0;
1123  case Geometry::SEGMENT: return 0;
1124  case Geometry::SQUARE: return 3;
1125  default:
1126  mfem_error ("P1OnQuadFECollection: unknown geometry type.");
1127  }
1128  return 0; // Make some compilers happy
1129 }
1130 
1132  Geometry::Type GeomType, int Or) const
1133 {
1134  return NULL;
1135 }
1136 
1137 
1138 const FiniteElement *
1140  Geometry::Type GeomType) const
1141 {
1142  switch (GeomType)
1143  {
1144  case Geometry::TRIANGLE: return &TriangleFE;
1145  case Geometry::SQUARE: return &QuadrilateralFE;
1146  default:
1147  if (error_mode == RETURN_NULL) { return nullptr; }
1148  mfem_error ("QuadraticDiscont2DFECollection: unknown geometry type.");
1149  }
1150  return &TriangleFE; // Make some compilers happy
1151 }
1152 
1154 const
1155 {
1156  switch (GeomType)
1157  {
1158  case Geometry::POINT: return 0;
1159  case Geometry::SEGMENT: return 0;
1160  case Geometry::TRIANGLE: return 6;
1161  case Geometry::SQUARE: return 9;
1162  default:
1163  mfem_error ("QuadraticDiscont2DFECollection: unknown geometry type.");
1164  }
1165  return 0; // Make some compilers happy
1166 }
1167 
1169  Geometry::Type GeomType, int Or) const
1170 {
1171  return NULL;
1172 }
1173 
1174 
1175 const FiniteElement *
1177  Geometry::Type GeomType) const
1178 {
1179  switch (GeomType)
1180  {
1181  case Geometry::SQUARE: return &QuadrilateralFE;
1182  default:
1183  if (error_mode == RETURN_NULL) { return nullptr; }
1184  mfem_error ("QuadraticPosDiscont2DFECollection: unknown geometry type.");
1185  }
1186  return NULL; // Make some compilers happy
1187 }
1188 
1190 const
1191 {
1192  switch (GeomType)
1193  {
1194  case Geometry::POINT: return 0;
1195  case Geometry::SEGMENT: return 0;
1196  case Geometry::SQUARE: return 9;
1197  default:
1198  mfem_error ("QuadraticPosDiscont2DFECollection: unknown geometry type.");
1199  }
1200  return 0; // Make some compilers happy
1201 }
1202 
1203 
1204 const FiniteElement *
1206  Geometry::Type GeomType)
1207 const
1208 {
1209  switch (GeomType)
1210  {
1211  case Geometry::TRIANGLE: return &TriangleFE;
1212  case Geometry::SQUARE: return &QuadrilateralFE;
1213  default:
1214  if (error_mode == RETURN_NULL) { return nullptr; }
1215  mfem_error ("GaussQuadraticDiscont2DFECollection:"
1216  " unknown geometry type.");
1217  }
1218  return &QuadrilateralFE; // Make some compilers happy
1219 }
1220 
1222  Geometry::Type GeomType) const
1223 {
1224  switch (GeomType)
1225  {
1226  case Geometry::POINT: return 0;
1227  case Geometry::SEGMENT: return 0;
1228  case Geometry::TRIANGLE: return 6;
1229  case Geometry::SQUARE: return 9;
1230  default:
1231  mfem_error ("GaussQuadraticDiscont2DFECollection:"
1232  " unknown geometry type.");
1233  }
1234  return 0; // Make some compilers happy
1235 }
1236 
1238  Geometry::Type GeomType, int Or) const
1239 {
1240  return NULL;
1241 }
1242 
1243 
1244 const FiniteElement *
1246  Geometry::Type GeomType) const
1247 {
1248  switch (GeomType)
1249  {
1250  case Geometry::TRIANGLE: return &TriangleFE;
1251  case Geometry::SQUARE: return &QuadrilateralFE;
1252  default:
1253  if (error_mode == RETURN_NULL) { return nullptr; }
1254  mfem_error ("CubicDiscont2DFECollection: unknown geometry type.");
1255  }
1256  return &TriangleFE; // Make some compilers happy
1257 }
1258 
1260 {
1261  switch (GeomType)
1262  {
1263  case Geometry::POINT: return 0;
1264  case Geometry::SEGMENT: return 0;
1265  case Geometry::TRIANGLE: return 10;
1266  case Geometry::SQUARE: return 16;
1267  default:
1268  mfem_error ("CubicDiscont2DFECollection: unknown geometry type.");
1269  }
1270  return 0; // Make some compilers happy
1271 }
1272 
1274  Geometry::Type GeomType, int Or) const
1275 {
1276  return NULL;
1277 }
1278 
1279 
1280 const FiniteElement *
1282  Geometry::Type GeomType) const
1283 {
1284  switch (GeomType)
1285  {
1286  case Geometry::TRIANGLE: return &TriangleFE;
1287  case Geometry::SQUARE: return &QuadrilateralFE;
1288  case Geometry::TETRAHEDRON: return &TetrahedronFE;
1289  case Geometry::CUBE: return &ParallelepipedFE;
1290  default:
1291  if (error_mode == RETURN_NULL) { return nullptr; }
1292  mfem_error ("LinearNonConf3DFECollection: unknown geometry type.");
1293  }
1294  return &TriangleFE; // Make some compilers happy
1295 }
1296 
1298 {
1299  switch (GeomType)
1300  {
1301  case Geometry::POINT: return 0;
1302  case Geometry::SEGMENT: return 0;
1303  case Geometry::TRIANGLE: return 1;
1304  case Geometry::SQUARE: return 1;
1305  case Geometry::TETRAHEDRON: return 0;
1306  case Geometry::CUBE: return 0;
1307  default:
1308  mfem_error ("LinearNonConf3DFECollection: unknown geometry type.");
1309  }
1310  return 0; // Make some compilers happy
1311 }
1312 
1314  Geometry::Type GeomType, int Or) const
1315 {
1316  static int indexes[] = { 0 };
1317 
1318  return indexes;
1319 }
1320 
1321 
1322 const FiniteElement *
1324 {
1325  switch (GeomType)
1326  {
1327  case Geometry::TETRAHEDRON: return &TetrahedronFE;
1328  case Geometry::CUBE: return &ParallelepipedFE;
1329  case Geometry::PRISM: return &WedgeFE;
1330  case Geometry::PYRAMID: return &PyramidFE;
1331  default:
1332  if (error_mode == RETURN_NULL) { return nullptr; }
1333  mfem_error ("Const3DFECollection: unknown geometry type.");
1334  }
1335  return &TetrahedronFE; // Make some compilers happy
1336 }
1337 
1339 {
1340  switch (GeomType)
1341  {
1342  case Geometry::POINT: return 0;
1343  case Geometry::SEGMENT: return 0;
1344  case Geometry::TRIANGLE: return 0;
1345  case Geometry::SQUARE: return 0;
1346  case Geometry::TETRAHEDRON: return 1;
1347  case Geometry::CUBE: return 1;
1348  case Geometry::PRISM: return 1;
1349  case Geometry::PYRAMID: return 1;
1350  default:
1351  mfem_error ("Const3DFECollection: unknown geometry type.");
1352  }
1353  return 0; // Make some compilers happy
1354 }
1355 
1357  int Or) const
1358 {
1359  return NULL;
1360 }
1361 
1362 
1363 const FiniteElement *
1365  Geometry::Type GeomType) const
1366 {
1367  switch (GeomType)
1368  {
1369  case Geometry::TETRAHEDRON: return &TetrahedronFE;
1370  case Geometry::PYRAMID: return &PyramidFE;
1371  case Geometry::PRISM: return &WedgeFE;
1372  case Geometry::CUBE: return &ParallelepipedFE;
1373  default:
1374  if (error_mode == RETURN_NULL) { return nullptr; }
1375  mfem_error ("LinearDiscont3DFECollection: unknown geometry type.");
1376  }
1377  return &TetrahedronFE; // Make some compilers happy
1378 }
1379 
1381 {
1382  switch (GeomType)
1383  {
1384  case Geometry::POINT: return 0;
1385  case Geometry::SEGMENT: return 0;
1386  case Geometry::TRIANGLE: return 0;
1387  case Geometry::SQUARE: return 0;
1388  case Geometry::TETRAHEDRON: return 4;
1389  case Geometry::PYRAMID: return 5;
1390  case Geometry::PRISM: return 6;
1391  case Geometry::CUBE: return 8;
1392  default:
1393  mfem_error ("LinearDiscont3DFECollection: unknown geometry type.");
1394  }
1395  return 0; // Make some compilers happy
1396 }
1397 
1399  Geometry::Type GeomType, int Or) const
1400 {
1401  return NULL;
1402 }
1403 
1404 
1405 const FiniteElement *
1407  Geometry::Type GeomType) const
1408 {
1409  switch (GeomType)
1410  {
1411  case Geometry::TETRAHEDRON: return &TetrahedronFE;
1412  case Geometry::CUBE: return &ParallelepipedFE;
1413  default:
1414  if (error_mode == RETURN_NULL) { return nullptr; }
1415  mfem_error ("QuadraticDiscont3DFECollection: unknown geometry type.");
1416  }
1417  return &TetrahedronFE; // Make some compilers happy
1418 }
1419 
1421 const
1422 {
1423  switch (GeomType)
1424  {
1425  case Geometry::POINT: return 0;
1426  case Geometry::SEGMENT: return 0;
1427  case Geometry::TRIANGLE: return 0;
1428  case Geometry::SQUARE: return 0;
1429  case Geometry::TETRAHEDRON: return 10;
1430  case Geometry::CUBE: return 27;
1431  default:
1432  mfem_error ("QuadraticDiscont3DFECollection: unknown geometry type.");
1433  }
1434  return 0; // Make some compilers happy
1435 }
1436 
1438  Geometry::Type GeomType, int Or) const
1439 {
1440  return NULL;
1441 }
1442 
1443 const FiniteElement *
1445  Geometry::Type GeomType) const
1446 {
1447  switch (GeomType)
1448  {
1449  case Geometry::POINT: return &PointFE;
1450  case Geometry::SEGMENT: return &SegmentFE;
1451  case Geometry::TRIANGLE: return &TriangleFE;
1452  case Geometry::SQUARE: return &QuadrilateralFE;
1453  case Geometry::TETRAHEDRON: return &TetrahedronFE;
1454  case Geometry::CUBE: return &ParallelepipedFE;
1455  default:
1456  if (error_mode == RETURN_NULL) { return nullptr; }
1457  mfem_error ("RefinedLinearFECollection: unknown geometry type.");
1458  }
1459  return &SegmentFE; // Make some compilers happy
1460 }
1461 
1463 {
1464  switch (GeomType)
1465  {
1466  case Geometry::POINT: return 1;
1467  case Geometry::SEGMENT: return 1;
1468  case Geometry::TRIANGLE: return 0;
1469  case Geometry::SQUARE: return 1;
1470  case Geometry::TETRAHEDRON: return 0;
1471  case Geometry::CUBE: return 1;
1472  default:
1473  mfem_error ("RefinedLinearFECollection: unknown geometry type.");
1474  }
1475  return 0; // Make some compilers happy
1476 }
1477 
1479  Geometry::Type GeomType, int Or) const
1480 {
1481  static int indexes[] = { 0 };
1482 
1483  return indexes;
1484 }
1485 
1486 
1487 const FiniteElement *
1489 {
1490  switch (GeomType)
1491  {
1492  case Geometry::CUBE: return &HexahedronFE;
1493  case Geometry::TETRAHEDRON: return &TetrahedronFE;
1494  case Geometry::PRISM: return &WedgeFE;
1495  case Geometry::PYRAMID: return &PyramidFE;
1496  default:
1497  if (error_mode == RETURN_NULL) { return nullptr; }
1498  mfem_error ("ND1_3DFECollection: unknown geometry type.");
1499  }
1500  return &HexahedronFE; // Make some compilers happy
1501 }
1502 
1504 {
1505  switch (GeomType)
1506  {
1507  case Geometry::POINT: return 0;
1508  case Geometry::SEGMENT: return 1;
1509  case Geometry::TRIANGLE: return 0;
1510  case Geometry::SQUARE: return 0;
1511  case Geometry::TETRAHEDRON: return 0;
1512  case Geometry::CUBE: return 0;
1513  case Geometry::PRISM: return 0;
1514  case Geometry::PYRAMID: return 0;
1515  default:
1516  mfem_error ("ND1_3DFECollection: unknown geometry type.");
1517  }
1518  return 0; // Make some compilers happy
1519 }
1520 
1522  int Or) const
1523 {
1524  static int ind_pos[] = { 0 };
1525  static int ind_neg[] = { -1 };
1526 
1527  if (Or > 0)
1528  {
1529  return ind_pos;
1530  }
1531  return ind_neg;
1532 }
1533 
1534 
1535 const FiniteElement *
1537 {
1538  switch (GeomType)
1539  {
1540  case Geometry::TRIANGLE: return &TriangleFE;
1541  case Geometry::SQUARE: return &QuadrilateralFE;
1542  case Geometry::CUBE: return &HexahedronFE;
1543  case Geometry::TETRAHEDRON: return &TetrahedronFE;
1544  case Geometry::PRISM: return &WedgeFE;
1545  case Geometry::PYRAMID: return &PyramidFE;
1546  default:
1547  if (error_mode == RETURN_NULL) { return nullptr; }
1548  mfem_error ("RT0_3DFECollection: unknown geometry type.");
1549  }
1550  return &HexahedronFE; // Make some compilers happy
1551 }
1552 
1554 {
1555  switch (GeomType)
1556  {
1557  case Geometry::POINT: return 0;
1558  case Geometry::SEGMENT: return 0;
1559  case Geometry::TRIANGLE: return 1;
1560  case Geometry::SQUARE: return 1;
1561  case Geometry::TETRAHEDRON: return 0;
1562  case Geometry::CUBE: return 0;
1563  case Geometry::PRISM: return 0;
1564  case Geometry::PYRAMID: return 0;
1565  default:
1566  mfem_error ("RT0_3DFECollection: unknown geometry type.");
1567  }
1568  return 0; // Make some compilers happy
1569 }
1570 
1572  int Or) const
1573 {
1574  static int ind_pos[] = { 0 };
1575  static int ind_neg[] = { -1 };
1576 
1577  if ((GeomType == Geometry::TRIANGLE) || (GeomType == Geometry::SQUARE))
1578  {
1579  if (Or % 2 == 0)
1580  {
1581  return ind_pos;
1582  }
1583  return ind_neg;
1584  }
1585  return NULL;
1586 }
1587 
1588 const FiniteElement *
1590 {
1591  switch (GeomType)
1592  {
1593  case Geometry::TRIANGLE: return &TriangleFE;
1594  case Geometry::SQUARE: return &QuadrilateralFE;
1595  case Geometry::CUBE: return &HexahedronFE;
1596  default:
1597  if (error_mode == RETURN_NULL) { return nullptr; }
1598  mfem_error ("RT1_3DFECollection: unknown geometry type.");
1599  }
1600  return &HexahedronFE; // Make some compilers happy
1601 }
1602 
1604 {
1605  switch (GeomType)
1606  {
1607  case Geometry::POINT: return 0;
1608  case Geometry::SEGMENT: return 0;
1609  case Geometry::TRIANGLE: return 2;
1610  case Geometry::SQUARE: return 4;
1611  case Geometry::CUBE: return 12;
1612  default:
1613  mfem_error ("RT1_3DFECollection: unknown geometry type.");
1614  }
1615  return 0; // Make some compilers happy
1616 }
1617 
1619  int Or) const
1620 {
1621  if (GeomType == Geometry::SQUARE)
1622  {
1623  static int sq_ind[8][4] =
1624  {
1625  {0, 1, 2, 3}, {-1, -3, -2, -4},
1626  {2, 0, 3, 1}, {-2, -1, -4, -3},
1627  {3, 2, 1, 0}, {-4, -2, -3, -1},
1628  {1, 3, 0, 2}, {-3, -4, -1, -2}
1629  };
1630 
1631  return sq_ind[Or];
1632  }
1633  else
1634  {
1635  return NULL;
1636  }
1637 }
1638 
1639 
1640 H1_FECollection::H1_FECollection(const int p, const int dim, const int btype)
1642  , dim(dim)
1643 {
1644  MFEM_VERIFY(p >= 1, "H1_FECollection requires order >= 1.");
1645  MFEM_VERIFY(dim >= 0 && dim <= 3, "H1_FECollection requires 0 <= dim <= 3.");
1646 
1647  const int pm1 = p - 1, pm2 = pm1 - 1, pm3 = pm2 - 1, pm4 = pm3 - 1;
1648 
1649  int pt_type = BasisType::GetQuadrature1D(btype);
1650  b_type = BasisType::Check(btype);
1651  switch (btype)
1652  {
1654  {
1655  snprintf(h1_name, 32, "H1_%dD_P%d", dim, p);
1656  break;
1657  }
1658  case BasisType::Positive:
1659  {
1660  snprintf(h1_name, 32, "H1Pos_%dD_P%d", dim, p);
1661  break;
1662  }
1664  {
1665  snprintf(h1_name, 32, "H1Ser_%dD_P%d", dim, p);
1666  break;
1667  }
1668  default:
1669  {
1670  MFEM_VERIFY(Quadrature1D::CheckClosed(pt_type) !=
1672  "unsupported BasisType: " << BasisType::Name(btype));
1673 
1674  snprintf(h1_name, 32, "H1@%c_%dD_P%d",
1675  (int)BasisType::GetChar(btype), dim, p);
1676  }
1677  }
1678 
1679  for (int g = 0; g < Geometry::NumGeom; g++)
1680  {
1681  H1_dof[g] = 0;
1682  H1_Elements[g] = NULL;
1683  }
1684  for (int i = 0; i < 2; i++)
1685  {
1686  SegDofOrd[i] = NULL;
1687  }
1688  for (int i = 0; i < 6; i++)
1689  {
1690  TriDofOrd[i] = NULL;
1691  }
1692  for (int i = 0; i < 8; i++)
1693  {
1694  QuadDofOrd[i] = NULL;
1695  }
1696  for (int i = 0; i < 24; i++)
1697  {
1698  TetDofOrd[i] = NULL;
1699  }
1700 
1701  H1_dof[Geometry::POINT] = 1;
1703 
1704  if (dim >= 1)
1705  {
1706  H1_dof[Geometry::SEGMENT] = pm1;
1707  if (b_type == BasisType::Positive)
1708  {
1710  }
1711  else
1712  {
1714  }
1715 
1716  SegDofOrd[0] = (pm1 > 0) ? new int[2*pm1] : nullptr;
1717  SegDofOrd[1] = SegDofOrd[0] + pm1;
1718  for (int i = 0; i < pm1; i++)
1719  {
1720  SegDofOrd[0][i] = i;
1721  SegDofOrd[1][i] = pm2 - i;
1722  }
1723  }
1724 
1725  if (dim >= 2)
1726  {
1727  H1_dof[Geometry::TRIANGLE] = (pm1*pm2)/2;
1728  H1_dof[Geometry::SQUARE] = pm1*pm1;
1729  if (b_type == BasisType::Positive)
1730  {
1733  }
1734  else if (b_type == BasisType::Serendipity)
1735  {
1736  // Note: in fe_coll.hpp the DofForGeometry(Geometry::Type) method
1737  // returns H1_dof[GeomType], so we need to fix the value of H1_dof here
1738  // for the serendipity case.
1739 
1740  // formula for number of interior serendipity DoFs (when p>1)
1741  H1_dof[Geometry::SQUARE] = (pm3*pm2)/2;
1743  // allows for mixed tri/quad meshes
1745  }
1746  else
1747  {
1750  }
1751 
1752  const int &TriDof = H1_dof[Geometry::TRIANGLE];
1753  const int &QuadDof = H1_dof[Geometry::SQUARE];
1754  TriDofOrd[0] = (TriDof > 0) ? new int[6*TriDof] : nullptr;
1755  for (int i = 1; i < 6; i++)
1756  {
1757  TriDofOrd[i] = TriDofOrd[i-1] + TriDof;
1758  }
1759  // see Mesh::GetTriOrientation in mesh/mesh.cpp
1760  for (int j = 0; j < pm2; j++)
1761  {
1762  for (int i = 0; i + j < pm2; i++)
1763  {
1764  int o = TriDof - ((pm1 - j)*(pm2 - j))/2 + i;
1765  int k = pm3 - j - i;
1766  TriDofOrd[0][o] = o; // (0,1,2)
1767  TriDofOrd[1][o] = TriDof - ((pm1-j)*(pm2-j))/2 + k; // (1,0,2)
1768  TriDofOrd[2][o] = TriDof - ((pm1-i)*(pm2-i))/2 + k; // (2,0,1)
1769  TriDofOrd[3][o] = TriDof - ((pm1-k)*(pm2-k))/2 + i; // (2,1,0)
1770  TriDofOrd[4][o] = TriDof - ((pm1-k)*(pm2-k))/2 + j; // (1,2,0)
1771  TriDofOrd[5][o] = TriDof - ((pm1-i)*(pm2-i))/2 + j; // (0,2,1)
1772  }
1773  }
1774 
1775  QuadDofOrd[0] = (QuadDof > 0) ? new int[8*QuadDof] : nullptr;
1776  for (int i = 1; i < 8; i++)
1777  {
1778  QuadDofOrd[i] = QuadDofOrd[i-1] + QuadDof;
1779  }
1780 
1781  // For serendipity order >=4, the QuadDofOrd array must be re-defined. We
1782  // do this by computing the corresponding tensor product QuadDofOrd array
1783  // or two orders less, which contains enough DoFs for their serendipity
1784  // basis. This could be optimized.
1786  {
1787  if (p < 4)
1788  {
1789  // no face dofs --> don't need to adjust QuadDofOrd
1790  }
1791  else // p >= 4 --> have face dofs
1792  {
1793  // Exactly the same as tensor product case, but with all orders
1794  // reduced by 2 e.g. in case p=5 it builds a 2x2 array, even though
1795  // there are only 3 serendipity dofs.
1796  // In the tensor product case, the i and j index tensor directions,
1797  // and o index from 0 to (pm1)^2,
1798 
1799  for (int j = 0; j < pm3; j++) // pm3 instead of pm1, etc
1800  {
1801  for (int i = 0; i < pm3; i++)
1802  {
1803  int o = i + j*pm3;
1804  QuadDofOrd[0][o] = i + j*pm3; // (0,1,2,3)
1805  QuadDofOrd[1][o] = j + i*pm3; // (0,3,2,1)
1806  QuadDofOrd[2][o] = j + (pm4 - i)*pm3; // (1,2,3,0)
1807  QuadDofOrd[3][o] = (pm4 - i) + j*pm3; // (1,0,3,2)
1808  QuadDofOrd[4][o] = (pm4 - i) + (pm4 - j)*pm3; // (2,3,0,1)
1809  QuadDofOrd[5][o] = (pm4 - j) + (pm4 - i)*pm3; // (2,1,0,3)
1810  QuadDofOrd[6][o] = (pm4 - j) + i*pm3; // (3,0,1,2)
1811  QuadDofOrd[7][o] = i + (pm4 - j)*pm3; // (3,2,1,0)
1812  }
1813  }
1814 
1815  }
1816  }
1817  else // not serendipity
1818  {
1819  for (int j = 0; j < pm1; j++)
1820  {
1821  for (int i = 0; i < pm1; i++)
1822  {
1823  int o = i + j*pm1;
1824  QuadDofOrd[0][o] = i + j*pm1; // (0,1,2,3)
1825  QuadDofOrd[1][o] = j + i*pm1; // (0,3,2,1)
1826  QuadDofOrd[2][o] = j + (pm2 - i)*pm1; // (1,2,3,0)
1827  QuadDofOrd[3][o] = (pm2 - i) + j*pm1; // (1,0,3,2)
1828  QuadDofOrd[4][o] = (pm2 - i) + (pm2 - j)*pm1; // (2,3,0,1)
1829  QuadDofOrd[5][o] = (pm2 - j) + (pm2 - i)*pm1; // (2,1,0,3)
1830  QuadDofOrd[6][o] = (pm2 - j) + i*pm1; // (3,0,1,2)
1831  QuadDofOrd[7][o] = i + (pm2 - j)*pm1; // (3,2,1,0)
1832  }
1833  }
1834  }
1835 
1836  if (dim >= 3)
1837  {
1838  H1_dof[Geometry::TETRAHEDRON] = (TriDof*pm3)/3;
1839  H1_dof[Geometry::CUBE] = QuadDof*pm1;
1840  H1_dof[Geometry::PRISM] = TriDof*pm1;
1842  if (b_type == BasisType::Positive)
1843  {
1847  }
1848  else
1849  {
1851  new H1_TetrahedronElement(p, btype);
1854  }
1856 
1857  const int &TetDof = H1_dof[Geometry::TETRAHEDRON];
1858  TetDofOrd[0] = (TetDof > 0) ? new int[24*TetDof] : nullptr;
1859  for (int i = 1; i < 24; i++)
1860  {
1861  TetDofOrd[i] = TetDofOrd[i-1] + TetDof;
1862  }
1863  // see Mesh::GetTetOrientation in mesh/mesh.cpp
1864  for (int k = 0; k < pm3; k++)
1865  {
1866  for (int j = 0; j + k < pm3; j++)
1867  {
1868  for (int i = 0; i + j + k < pm3; i++)
1869  {
1870  int l = pm4 - k - j - i;
1871  int o = TetDof - ((pm1 - k) * (pm2 - k) * (pm3 - k)) / 6
1872  + (j * (2 * p - 5 - j - 2 * k)) / 2 + i;
1873  int o1 = TetDof - ((pm1 - j) * (pm2 - j) * (pm3 - j)) / 6
1874  + (k * (2 * p - 5 - k - 2 * j)) / 2 + i;
1875  int o2 = TetDof - ((pm1 - i) * (pm2 - i) * (pm3 - i)) / 6
1876  + (k * (2 * p - 5 - k - 2 * i)) / 2 + j;
1877  int o3 = TetDof - ((pm1 - k) * (pm2 - k) * (pm3 - k)) / 6
1878  + (i * (2 * p - 5 - i - 2 * k)) / 2 + j;
1879  int o4 = TetDof - ((pm1 - j) * (pm2 - j) * (pm3 - j)) / 6
1880  + (i * (2 * p - 5 - i - 2 * j)) / 2 + k;
1881  int o5 = TetDof - ((pm1 - i) * (pm2 - i) * (pm3 - i)) / 6
1882  + (j * (2 * p - 5 - j - 2 * i)) / 2 + k;
1883  int o6 = TetDof - ((pm1 - k) * (pm2 - k) * (pm3 - k)) / 6
1884  + (l * (2 * p - 5 - l - 2 * k)) / 2 + j;
1885  int o7 = TetDof - ((pm1 - l) * (pm2 - l) * (pm3 - l)) / 6
1886  + (k * (2 * p - 5 - k - 2 * l)) / 2 + j;
1887  int o8 = TetDof - ((pm1 - l) * (pm2 - l) * (pm3 - l)) / 6
1888  + (j * (2 * p - 5 - j - 2 * l)) / 2 + k;
1889  int o9 = TetDof - ((pm1 - j) * (pm2 - j) * (pm3 - j)) / 6
1890  + (l * (2 * p - 5 - l - 2 * j)) / 2 + k;
1891  int o10 = TetDof - ((pm1 - j) * (pm2 - j) * (pm3 - j)) / 6
1892  + (k * (2 * p - 5 - k - 2 * j)) / 2 + l;
1893  int o11 = TetDof - ((pm1 - k) * (pm2 - k) * (pm3 - k)) / 6
1894  + (j * (2 * p - 5 - j - 2 * k)) / 2 + l;
1895  int o12 = TetDof - ((pm1 - i) * (pm2 - i) * (pm3 - i)) / 6
1896  + (l * (2 * p - 5 - l - 2 * i)) / 2 + k;
1897  int o13 = TetDof - ((pm1 - l) * (pm2 - l) * (pm3 - l)) / 6
1898  + (i * (2 * p - 5 - i - 2 * l)) / 2 + k;
1899  int o14 = TetDof - ((pm1 - k) * (pm2 - k) * (pm3 - k)) / 6
1900  + (i * (2 * p - 5 - i - 2 * k)) / 2 + l;
1901  int o15 = TetDof - ((pm1 - i) * (pm2 - i) * (pm3 - i)) / 6
1902  + (k * (2 * p - 5 - k - 2 * i)) / 2 + l;
1903  int o16 = TetDof - ((pm1 - l) * (pm2 - l) * (pm3 - l)) / 6
1904  + (k * (2 * p - 5 - k - 2 * l)) / 2 + i;
1905  int o17 = TetDof - ((pm1 - k) * (pm2 - k) * (pm3 - k)) / 6
1906  + (l * (2 * p - 5 - l - 2 * k)) / 2 + i;
1907  int o18 = TetDof - ((pm1 - i) * (pm2 - i) * (pm3 - i)) / 6
1908  + (j * (2 * p - 5 - j - 2 * i)) / 2 + l;
1909  int o19 = TetDof - ((pm1 - j) * (pm2 - j) * (pm3 - j)) / 6
1910  + (i * (2 * p - 5 - i - 2 * j)) / 2 + l;
1911  int o20 = TetDof - ((pm1 - j) * (pm2 - j) * (pm3 - j)) / 6
1912  + (l * (2 * p - 5 - l - 2 * j)) / 2 + i;
1913  int o21 = TetDof - ((pm1 - l) * (pm2 - l) * (pm3 - l)) / 6
1914  + (j * (2 * p - 5 - j - 2 * l)) / 2 + i;
1915  int o22 = TetDof - ((pm1 - l) * (pm2 - l) * (pm3 - l)) / 6
1916  + (i * (2 * p - 5 - i - 2 * l)) / 2 + j;
1917  int o23 = TetDof - ((pm1 - i) * (pm2 - i) * (pm3 - i)) / 6
1918  + (l * (2 * p - 5 - l - 2 * i)) / 2 + j;
1919  TetDofOrd[ 0][o] = o; // (0,1,2,3)
1920  TetDofOrd[ 1][o] = o1; // (0,1,3,2)
1921  TetDofOrd[ 2][o] = o2; // (0,2,3,1)
1922  TetDofOrd[ 3][o] = o3; // (0,2,1,3)
1923  TetDofOrd[ 4][o] = o4; // (0,3,1,2)
1924  TetDofOrd[ 5][o] = o5; // (0,3,2,1)
1925  TetDofOrd[ 6][o] = o6; // (1,2,0,3)
1926  TetDofOrd[ 7][o] = o7; // (1,2,3,0)
1927  TetDofOrd[ 8][o] = o8; // (1,3,2,0)
1928  TetDofOrd[ 9][o] = o9; // (1,3,0,2)
1929  TetDofOrd[10][o] = o10; // (1,0,3,2)
1930  TetDofOrd[11][o] = o11; // (1,0,2,3)
1931  TetDofOrd[12][o] = o12; // (2,3,0,1)
1932  TetDofOrd[13][o] = o13; // (2,3,1,0)
1933  TetDofOrd[14][o] = o14; // (2,0,1,3)
1934  TetDofOrd[15][o] = o15; // (2,0,3,1)
1935  TetDofOrd[16][o] = o16; // (2,1,3,0)
1936  TetDofOrd[17][o] = o17; // (2,1,0,3)
1937  TetDofOrd[18][o] = o18; // (3,0,2,1)
1938  TetDofOrd[19][o] = o19; // (3,0,1,2)
1939  TetDofOrd[20][o] = o20; // (3,1,0,2)
1940  TetDofOrd[21][o] = o21; // (3,1,2,0)
1941  TetDofOrd[22][o] = o22; // (3,2,1,0)
1942  TetDofOrd[23][o] = o23; // (3,2,0,1)
1943  }
1944  }
1945  }
1946  }
1947  }
1948 }
1949 
1950 const FiniteElement *
1952 {
1953  if (GeomType != Geometry::PYRAMID || this->GetOrder() == 1)
1954  {
1955  return H1_Elements[GeomType];
1956  }
1957  else
1958  {
1959  if (error_mode == RETURN_NULL) { return nullptr; }
1960  MFEM_ABORT("H1 Pyramid basis functions are not yet supported "
1961  "for order > 1.");
1962  return NULL;
1963  }
1964 }
1965 
1967  int Or) const
1968 {
1969  if (GeomType == Geometry::SEGMENT)
1970  {
1971  return (Or > 0) ? SegDofOrd[0] : SegDofOrd[1];
1972  }
1973  else if (GeomType == Geometry::TRIANGLE)
1974  {
1975  return TriDofOrd[Or%6];
1976  }
1977  else if (GeomType == Geometry::SQUARE)
1978  {
1979  return QuadDofOrd[Or%8];
1980  }
1981  else if (GeomType == Geometry::TETRAHEDRON)
1982  {
1983  return TetDofOrd[Or%24];
1984  }
1985  return NULL;
1986 }
1987 
1989 {
1990  int tr_p = H1_dof[Geometry::SEGMENT] + 1;
1991  int tr_dim = -1;
1992  if (!strncmp(h1_name, "H1_", 3))
1993  {
1994  tr_dim = atoi(h1_name + 3);
1995  }
1996  else if (!strncmp(h1_name, "H1Pos_", 6))
1997  {
1998  tr_dim = atoi(h1_name + 6);
1999  }
2000  else if (!strncmp(h1_name, "H1@", 3))
2001  {
2002  tr_dim = atoi(h1_name + 5);
2003  }
2004  return (dim < 0) ? NULL : new H1_Trace_FECollection(tr_p, tr_dim, b_type);
2005 }
2006 
2007 const int *H1_FECollection::GetDofMap(Geometry::Type GeomType) const
2008 {
2009  const int *dof_map = NULL;
2010  const FiniteElement *fe = H1_Elements[GeomType];
2011  const NodalFiniteElement *nodal_fe =
2012  dynamic_cast<const NodalFiniteElement*>(fe);
2013  if (nodal_fe)
2014  {
2015  dof_map = nodal_fe->GetLexicographicOrdering().GetData();
2016  }
2017  else
2018  {
2019  MFEM_ABORT("Geometry type " << Geometry::Name[GeomType] << " is not "
2020  "implemented");
2021  }
2022  return dof_map;
2023 }
2024 
2025 const int *H1_FECollection::GetDofMap(Geometry::Type GeomType, int p) const
2026 {
2027  if (p == base_p) { return GetDofMap(GeomType); }
2028  if (p >= var_orders.Size() || !var_orders[p]) { InitVarOrder(p); }
2029  return ((H1_FECollection*) var_orders[p])->GetDofMap(GeomType);
2030 }
2031 
2033 {
2034  delete [] SegDofOrd[0];
2035  delete [] TriDofOrd[0];
2036  delete [] QuadDofOrd[0];
2037  delete [] TetDofOrd[0];
2038  for (int g = 0; g < Geometry::NumGeom; g++)
2039  {
2040  delete H1_Elements[g];
2041  }
2042 }
2043 
2044 
2046  const int btype)
2047  : H1_FECollection(p, dim-1, btype)
2048 {
2049  if (btype == BasisType::GaussLobatto)
2050  {
2051  snprintf(h1_name, 32, "H1_Trace_%dD_P%d", dim, p);
2052  }
2053  else if (btype == BasisType::Positive)
2054  {
2055  snprintf(h1_name, 32, "H1Pos_Trace_%dD_P%d", dim, p);
2056  }
2057  else // base class checks that type is closed
2058  {
2059  snprintf(h1_name, 32, "H1_Trace@%c_%dD_P%d",
2060  (int)BasisType::GetChar(btype), dim, p);
2061  }
2062 }
2063 
2064 
2065 L2_FECollection::L2_FECollection(const int p, const int dim, const int btype,
2066  const int map_type)
2068  , dim(dim)
2069  , m_type(map_type)
2070 {
2071  MFEM_VERIFY(p >= 0, "L2_FECollection requires order >= 0.");
2072 
2073  b_type = BasisType::Check(btype);
2074  const char *prefix = NULL;
2075  switch (map_type)
2076  {
2077  case FiniteElement::VALUE: prefix = "L2"; break;
2078  case FiniteElement::INTEGRAL: prefix = "L2Int"; break;
2079  default:
2080  MFEM_ABORT("invalid map_type: " << map_type);
2081  }
2082  switch (btype)
2083  {
2085  snprintf(d_name, 32, "%s_%dD_P%d", prefix, dim, p);
2086  break;
2087  default:
2088  snprintf(d_name, 32, "%s_T%d_%dD_P%d", prefix, btype, dim, p);
2089  }
2090 
2091  for (int g = 0; g < Geometry::NumGeom; g++)
2092  {
2093  L2_Elements[g] = NULL;
2094  Tr_Elements[g] = NULL;
2095  }
2096  for (int i = 0; i < 2; i++)
2097  {
2098  SegDofOrd[i] = NULL;
2099  }
2100  for (int i = 0; i < 6; i++)
2101  {
2102  TriDofOrd[i] = NULL;
2103  }
2104  for (int i = 0; i < 24; i++)
2105  {
2106  TetDofOrd[i] = NULL;
2107  }
2108  OtherDofOrd = NULL;
2109 
2110  if (dim == 0)
2111  {
2112  L2_Elements[Geometry::POINT] = new PointFiniteElement;
2113  }
2114  else if (dim == 1)
2115  {
2116  if (b_type == BasisType::Positive)
2117  {
2118  L2_Elements[Geometry::SEGMENT] = new L2Pos_SegmentElement(p);
2119  }
2120  else
2121  {
2122  L2_Elements[Geometry::SEGMENT] = new L2_SegmentElement(p, btype);
2123  }
2124  L2_Elements[Geometry::SEGMENT]->SetMapType(map_type);
2125 
2126  Tr_Elements[Geometry::POINT] = new PointFiniteElement;
2127  // No need to set the map_type for Tr_Elements.
2128 
2129  const int pp1 = p + 1;
2130  SegDofOrd[0] = (pp1 > 0) ? new int[2*pp1] : nullptr;
2131  SegDofOrd[1] = SegDofOrd[0] + pp1;
2132  for (int i = 0; i <= p; i++)
2133  {
2134  SegDofOrd[0][i] = i;
2135  SegDofOrd[1][i] = p - i;
2136  }
2137  }
2138  else if (dim == 2)
2139  {
2140  if (b_type == BasisType::Positive)
2141  {
2142  L2_Elements[Geometry::TRIANGLE] = new L2Pos_TriangleElement(p);
2143  L2_Elements[Geometry::SQUARE] = new L2Pos_QuadrilateralElement(p);
2144  }
2145  else
2146  {
2147  L2_Elements[Geometry::TRIANGLE] = new L2_TriangleElement(p, btype);
2148  L2_Elements[Geometry::SQUARE] = new L2_QuadrilateralElement(p, btype);
2149  }
2150  L2_Elements[Geometry::TRIANGLE]->SetMapType(map_type);
2151  L2_Elements[Geometry::SQUARE]->SetMapType(map_type);
2152  // Trace element use the default Gauss-Legendre nodal points for positive basis
2153  if (b_type == BasisType::Positive)
2154  {
2155  Tr_Elements[Geometry::SEGMENT] = new L2Pos_SegmentElement(p);
2156  }
2157  else
2158  {
2159  Tr_Elements[Geometry::SEGMENT] = new L2_SegmentElement(p, btype);
2160  }
2161 
2162  const int TriDof = L2_Elements[Geometry::TRIANGLE]->GetDof();
2163  TriDofOrd[0] = (TriDof > 0) ? new int[6*TriDof] : nullptr;
2164  for (int i = 1; i < 6; i++)
2165  {
2166  TriDofOrd[i] = TriDofOrd[i-1] + TriDof;
2167  }
2168  const int pp1 = p + 1, pp2 = pp1 + 1;
2169  for (int j = 0; j <= p; j++)
2170  {
2171  for (int i = 0; i + j <= p; i++)
2172  {
2173  int o = TriDof - ((pp2 - j)*(pp1 - j))/2 + i;
2174  int k = p - j - i;
2175  TriDofOrd[0][o] = o; // (0,1,2)
2176  TriDofOrd[1][o] = TriDof - ((pp2-j)*(pp1-j))/2 + k; // (1,0,2)
2177  TriDofOrd[2][o] = TriDof - ((pp2-i)*(pp1-i))/2 + k; // (2,0,1)
2178  TriDofOrd[3][o] = TriDof - ((pp2-k)*(pp1-k))/2 + i; // (2,1,0)
2179  TriDofOrd[4][o] = TriDof - ((pp2-k)*(pp1-k))/2 + j; // (1,2,0)
2180  TriDofOrd[5][o] = TriDof - ((pp2-i)*(pp1-i))/2 + j; // (0,2,1)
2181  }
2182  }
2183  const int QuadDof = L2_Elements[Geometry::SQUARE]->GetDof();
2184  OtherDofOrd = (QuadDof > 0) ? new int[QuadDof] : nullptr;
2185  for (int j = 0; j < QuadDof; j++)
2186  {
2187  OtherDofOrd[j] = j; // for Or == 0
2188  }
2189  }
2190  else if (dim == 3)
2191  {
2192  if (b_type == BasisType::Positive)
2193  {
2194  L2_Elements[Geometry::TETRAHEDRON] = new L2Pos_TetrahedronElement(p);
2195  L2_Elements[Geometry::CUBE] = new L2Pos_HexahedronElement(p);
2196  L2_Elements[Geometry::PRISM] = new L2Pos_WedgeElement(p);
2197  }
2198  else
2199  {
2200  L2_Elements[Geometry::TETRAHEDRON] =
2201  new L2_TetrahedronElement(p, btype);
2202  L2_Elements[Geometry::CUBE] = new L2_HexahedronElement(p, btype);
2203  L2_Elements[Geometry::PRISM] = new L2_WedgeElement(p, btype);
2204  }
2205  L2_Elements[Geometry::PYRAMID] = new P0PyrFiniteElement;
2206 
2207  L2_Elements[Geometry::TETRAHEDRON]->SetMapType(map_type);
2208  L2_Elements[Geometry::CUBE]->SetMapType(map_type);
2209  L2_Elements[Geometry::PRISM]->SetMapType(map_type);
2210  L2_Elements[Geometry::PYRAMID]->SetMapType(map_type);
2211  // Trace element use the default Gauss-Legendre nodal points for positive basis
2212  if (b_type == BasisType::Positive)
2213  {
2214  Tr_Elements[Geometry::TRIANGLE] = new L2Pos_TriangleElement(p);
2215  Tr_Elements[Geometry::SQUARE] = new L2Pos_QuadrilateralElement(p);
2216  }
2217  else
2218  {
2219  Tr_Elements[Geometry::TRIANGLE] = new L2_TriangleElement(p, btype);
2220  Tr_Elements[Geometry::SQUARE] = new L2_QuadrilateralElement(p, btype);
2221  }
2222 
2223  const int TetDof = L2_Elements[Geometry::TETRAHEDRON]->GetDof();
2224  const int HexDof = L2_Elements[Geometry::CUBE]->GetDof();
2225  const int PriDof = L2_Elements[Geometry::PRISM]->GetDof();
2226  const int MaxDof = std::max(TetDof, std::max(PriDof, HexDof));
2227 
2228  TetDofOrd[0] = (TetDof > 0) ? new int[24*TetDof] : nullptr;
2229  for (int i = 1; i < 24; i++)
2230  {
2231  TetDofOrd[i] = TetDofOrd[i-1] + TetDof;
2232  }
2233  // see Mesh::GetTetOrientation in mesh/mesh.cpp
2234  const int pp1 = p + 1, pp2 = pp1 + 1, pp3 = pp2 + 1;
2235  for (int k = 0; k <= p; k++)
2236  {
2237  for (int j = 0; j + k <= p; j++)
2238  {
2239  for (int i = 0; i + j + k <= p; i++)
2240  {
2241  int l = p - k - j - i;
2242  int o = TetDof - ((pp1 - k) * (pp2 - k) * (pp3 - k)) / 6
2243  + (j * (2 * p + 3 - j - 2 * k)) / 2 + i;
2244  int o1 = TetDof - ((pp1 - j) * (pp2 - j) * (pp3 - j)) / 6
2245  + (k * (2 * p + 3 - k - 2 * j)) / 2 + i;
2246  int o2 = TetDof - ((pp1 - i) * (pp2 - i) * (pp3 - i)) / 6
2247  + (k * (2 * p + 3 - k - 2 * i)) / 2 + j;
2248  int o3 = TetDof - ((pp1 - k) * (pp2 - k) * (pp3 - k)) / 6
2249  + (i * (2 * p + 3 - i - 2 * k)) / 2 + j;
2250  int o4 = TetDof - ((pp1 - j) * (pp2 - j) * (pp3 - j)) / 6
2251  + (i * (2 * p + 3 - i - 2 * j)) / 2 + k;
2252  int o5 = TetDof - ((pp1 - i) * (pp2 - i) * (pp3 - i)) / 6
2253  + (j * (2 * p + 3 - j - 2 * i)) / 2 + k;
2254  int o6 = TetDof - ((pp1 - k) * (pp2 - k) * (pp3 - k)) / 6
2255  + (l * (2 * p + 3 - l - 2 * k)) / 2 + j;
2256  int o7 = TetDof - ((pp1 - l) * (pp2 - l) * (pp3 - l)) / 6
2257  + (k * (2 * p + 3 - k - 2 * l)) / 2 + j;
2258  int o8 = TetDof - ((pp1 - l) * (pp2 - l) * (pp3 - l)) / 6
2259  + (j * (2 * p + 3 - j - 2 * l)) / 2 + k;
2260  int o9 = TetDof - ((pp1 - j) * (pp2 - j) * (pp3 - j)) / 6
2261  + (l * (2 * p + 3 - l - 2 * j)) / 2 + k;
2262  int o10 = TetDof - ((pp1 - j) * (pp2 - j) * (pp3 - j)) / 6
2263  + (k * (2 * p + 3 - k - 2 * j)) / 2 + l;
2264  int o11 = TetDof - ((pp1 - k) * (pp2 - k) * (pp3 - k)) / 6
2265  + (j * (2 * p + 3 - j - 2 * k)) / 2 + l;
2266  int o12 = TetDof - ((pp1 - i) * (pp2 - i) * (pp3 - i)) / 6
2267  + (l * (2 * p + 3 - l - 2 * i)) / 2 + k;
2268  int o13 = TetDof - ((pp1 - l) * (pp2 - l) * (pp3 - l)) / 6
2269  + (i * (2 * p + 3 - i - 2 * l)) / 2 + k;
2270  int o14 = TetDof - ((pp1 - k) * (pp2 - k) * (pp3 - k)) / 6
2271  + (i * (2 * p + 3 - i - 2 * k)) / 2 + l;
2272  int o15 = TetDof - ((pp1 - i) * (pp2 - i) * (pp3 - i)) / 6
2273  + (k * (2 * p + 3 - k - 2 * i)) / 2 + l;
2274  int o16 = TetDof - ((pp1 - l) * (pp2 - l) * (pp3 - l)) / 6
2275  + (k * (2 * p + 3 - k - 2 * l)) / 2 + i;
2276  int o17 = TetDof - ((pp1 - k) * (pp2 - k) * (pp3 - k)) / 6
2277  + (l * (2 * p + 3 - l - 2 * k)) / 2 + i;
2278  int o18 = TetDof - ((pp1 - i) * (pp2 - i) * (pp3 - i)) / 6
2279  + (j * (2 * p + 3 - j - 2 * i)) / 2 + l;
2280  int o19 = TetDof - ((pp1 - j) * (pp2 - j) * (pp3 - j)) / 6
2281  + (i * (2 * p + 3 - i - 2 * j)) / 2 + l;
2282  int o20 = TetDof - ((pp1 - j) * (pp2 - j) * (pp3 - j)) / 6
2283  + (l * (2 * p + 3 - l - 2 * j)) / 2 + i;
2284  int o21 = TetDof - ((pp1 - l) * (pp2 - l) * (pp3 - l)) / 6
2285  + (j * (2 * p + 3 - j - 2 * l)) / 2 + i;
2286  int o22 = TetDof - ((pp1 - l) * (pp2 - l) * (pp3 - l)) / 6
2287  + (i * (2 * p + 3 - i - 2 * l)) / 2 + j;
2288  int o23 = TetDof - ((pp1 - i) * (pp2 - i) * (pp3 - i)) / 6
2289  + (l * (2 * p + 3 - l - 2 * i)) / 2 + j;
2290  TetDofOrd[ 0][o] = o; // (0,1,2,3)
2291  TetDofOrd[ 1][o] = o1; // (0,1,3,2)
2292  TetDofOrd[ 2][o] = o2; // (0,2,3,1)
2293  TetDofOrd[ 3][o] = o3; // (0,2,1,3)
2294  TetDofOrd[ 4][o] = o4; // (0,3,1,2)
2295  TetDofOrd[ 5][o] = o5; // (0,3,2,1)
2296  TetDofOrd[ 6][o] = o6; // (1,2,0,3)
2297  TetDofOrd[ 7][o] = o7; // (1,2,3,0)
2298  TetDofOrd[ 8][o] = o8; // (1,3,2,0)
2299  TetDofOrd[ 9][o] = o9; // (1,3,0,2)
2300  TetDofOrd[10][o] = o10; // (1,0,3,2)
2301  TetDofOrd[11][o] = o11; // (1,0,2,3)
2302  TetDofOrd[12][o] = o12; // (2,3,0,1)
2303  TetDofOrd[13][o] = o13; // (2,3,1,0)
2304  TetDofOrd[14][o] = o14; // (2,0,1,3)
2305  TetDofOrd[15][o] = o15; // (2,0,3,1)
2306  TetDofOrd[16][o] = o16; // (2,1,3,0)
2307  TetDofOrd[17][o] = o17; // (2,1,0,3)
2308  TetDofOrd[18][o] = o18; // (3,0,2,1)
2309  TetDofOrd[19][o] = o19; // (3,0,1,2)
2310  TetDofOrd[20][o] = o20; // (3,1,0,2)
2311  TetDofOrd[21][o] = o21; // (3,1,2,0)
2312  TetDofOrd[22][o] = o22; // (3,2,1,0)
2313  TetDofOrd[23][o] = o23; // (3,2,0,1)
2314  }
2315  }
2316  }
2317  OtherDofOrd = (MaxDof > 0) ? new int[MaxDof] : nullptr;
2318  for (int j = 0; j < MaxDof; j++)
2319  {
2320  OtherDofOrd[j] = j; // for Or == 0
2321  }
2322  }
2323  else
2324  {
2325  mfem::err << "L2_FECollection::L2_FECollection : dim = "
2326  << dim << endl;
2327  mfem_error();
2328  }
2329 }
2330 
2331 const FiniteElement *
2333 {
2334  if (GeomType != Geometry::PYRAMID || this->GetOrder() == 0)
2335  {
2336  return L2_Elements[GeomType];
2337  }
2338  else
2339  {
2340  if (error_mode == RETURN_NULL) { return nullptr; }
2341  MFEM_ABORT("L2 Pyramid basis functions are not yet supported "
2342  "for order > 0.");
2343  return NULL;
2344  }
2345 }
2346 
2348  int Or) const
2349 {
2350  switch (GeomType)
2351  {
2352  case Geometry::SEGMENT:
2353  return (Or > 0) ? SegDofOrd[0] : SegDofOrd[1];
2354 
2355  case Geometry::TRIANGLE:
2356  return TriDofOrd[Or%6];
2357 
2358  case Geometry::TETRAHEDRON:
2359  return TetDofOrd[Or%24];
2360 
2361  default:
2362  return (Or == 0) ? OtherDofOrd : NULL;
2363  }
2364 }
2365 
2367 {
2368  delete [] OtherDofOrd;
2369  delete [] SegDofOrd[0];
2370  delete [] TriDofOrd[0];
2371  delete [] TetDofOrd[0];
2372  for (int i = 0; i < Geometry::NumGeom; i++)
2373  {
2374  delete L2_Elements[i];
2375  delete Tr_Elements[i];
2376  }
2377 }
2378 
2379 
2380 RT_FECollection::RT_FECollection(const int order, const int dim,
2381  const int cb_type, const int ob_type)
2382  : FiniteElementCollection(order + 1)
2383  , dim(dim)
2384  , cb_type(cb_type)
2385  , ob_type(ob_type)
2386 {
2387  int p = order;
2388  MFEM_VERIFY(p >= 0, "RT_FECollection requires order >= 0.");
2389 
2390  int cp_type = BasisType::GetQuadrature1D(cb_type);
2391  int op_type = BasisType::GetQuadrature1D(ob_type);
2392 
2394  {
2395  const char *cb_name = BasisType::Name(cb_type); // this may abort
2396  MFEM_ABORT("unknown closed BasisType: " << cb_name);
2397  }
2400  {
2401  const char *ob_name = BasisType::Name(ob_type); // this may abort
2402  MFEM_ABORT("unknown open BasisType: " << ob_name);
2403  }
2404 
2406 
2409  {
2410  snprintf(rt_name, 32, "RT_%dD_P%d", dim, p);
2411  }
2412  else
2413  {
2414  snprintf(rt_name, 32, "RT@%c%c_%dD_P%d", (int)BasisType::GetChar(cb_type),
2415  (int)BasisType::GetChar(ob_type), dim, p);
2416  }
2417 
2418  const int pp1 = p + 1;
2419  if (dim == 2)
2420  {
2421  // TODO: cb_type, ob_type for triangles
2423  RT_dof[Geometry::TRIANGLE] = p*pp1;
2424 
2426  ob_type);
2427  // two vector components * n_unk_face *
2428  RT_dof[Geometry::SQUARE] = 2*p*pp1;
2429  }
2430  else if (dim == 3)
2431  {
2432  // TODO: cb_type, ob_type for tets
2434  RT_dof[Geometry::TETRAHEDRON] = p*pp1*(p + 2)/2;
2435 
2437  RT_dof[Geometry::CUBE] = 3*p*pp1*pp1;
2438 
2440  RT_dof[Geometry::PRISM] = p*pp1*(3*p + 4)/2;
2441 
2444  }
2445  else
2446  {
2447  MFEM_ABORT("invalid dim = " << dim);
2448  }
2449 }
2450 
2451 // This is a special protected constructor only used by RT_Trace_FECollection
2452 // and DG_Interface_FECollection
2454  const int map_type, const bool signs,
2455  const int ob_type)
2456  : FiniteElementCollection(p + 1)
2457  , ob_type(ob_type)
2458 {
2461  {
2462  const char *ob_name = BasisType::Name(ob_type); // this may abort
2463  MFEM_ABORT("Invalid open basis type: " << ob_name);
2464  }
2465  InitFaces(p, dim, map_type, signs);
2466 }
2467 
2468 void RT_FECollection::InitFaces(const int p, const int dim_,
2469  const int map_type,
2470  const bool signs)
2471 {
2472  int op_type = BasisType::GetQuadrature1D(ob_type);
2473 
2474  MFEM_VERIFY(Quadrature1D::CheckOpen(op_type) != Quadrature1D::Invalid,
2475  "invalid open point type");
2476 
2477  const int pp1 = p + 1, pp2 = p + 2;
2478 
2479  for (int g = 0; g < Geometry::NumGeom; g++)
2480  {
2481  RT_Elements[g] = NULL;
2482  RT_dof[g] = 0;
2483  }
2484  // Degree of Freedom orderings
2485  for (int i = 0; i < 2; i++)
2486  {
2487  SegDofOrd[i] = NULL;
2488  }
2489  for (int i = 0; i < 6; i++)
2490  {
2491  TriDofOrd[i] = NULL;
2492  }
2493  for (int i = 0; i < 8; i++)
2494  {
2495  QuadDofOrd[i] = NULL;
2496  }
2497 
2498  if (dim_ == 2)
2499  {
2500  L2_SegmentElement *l2_seg = new L2_SegmentElement(p, ob_type);
2501  l2_seg->SetMapType(map_type);
2502  RT_Elements[Geometry::SEGMENT] = l2_seg;
2503  RT_dof[Geometry::SEGMENT] = pp1;
2504 
2505  SegDofOrd[0] = (pp1 > 0) ? new int[2*pp1] : nullptr;
2506  SegDofOrd[1] = SegDofOrd[0] + pp1;
2507  for (int i = 0; i <= p; i++)
2508  {
2509  SegDofOrd[0][i] = i;
2510  SegDofOrd[1][i] = signs ? (-1 - (p - i)) : (p - i);
2511  }
2512  }
2513  else if (dim_ == 3)
2514  {
2516  l2_tri->SetMapType(map_type);
2517  RT_Elements[Geometry::TRIANGLE] = l2_tri;
2518  RT_dof[Geometry::TRIANGLE] = pp1*pp2/2;
2519 
2521  l2_quad->SetMapType(map_type);
2522  RT_Elements[Geometry::SQUARE] = l2_quad;
2523  RT_dof[Geometry::SQUARE] = pp1*pp1;
2524 
2525  int TriDof = RT_dof[Geometry::TRIANGLE];
2526  TriDofOrd[0] = (TriDof > 0) ? new int[6*TriDof] : nullptr;
2527  for (int i = 1; i < 6; i++)
2528  {
2529  TriDofOrd[i] = TriDofOrd[i-1] + TriDof;
2530  }
2531  // see Mesh::GetTriOrientation in mesh/mesh.cpp,
2532  // the constructor of H1_FECollection
2533  for (int j = 0; j <= p; j++)
2534  {
2535  for (int i = 0; i + j <= p; i++)
2536  {
2537  int o = TriDof - ((pp2 - j)*(pp1 - j))/2 + i;
2538  int k = p - j - i;
2539  TriDofOrd[0][o] = o; // (0,1,2)
2540  TriDofOrd[1][o] = -1-(TriDof-((pp2-j)*(pp1-j))/2+k); // (1,0,2)
2541  TriDofOrd[2][o] = TriDof-((pp2-i)*(pp1-i))/2+k; // (2,0,1)
2542  TriDofOrd[3][o] = -1-(TriDof-((pp2-k)*(pp1-k))/2+i); // (2,1,0)
2543  TriDofOrd[4][o] = TriDof-((pp2-k)*(pp1-k))/2+j; // (1,2,0)
2544  TriDofOrd[5][o] = -1-(TriDof-((pp2-i)*(pp1-i))/2+j); // (0,2,1)
2545  if (!signs)
2546  {
2547  for (int kk = 1; kk < 6; kk += 2)
2548  {
2549  TriDofOrd[kk][o] = -1 - TriDofOrd[kk][o];
2550  }
2551  }
2552  }
2553  }
2554 
2555  int QuadDof = RT_dof[Geometry::SQUARE];
2556  QuadDofOrd[0] = (QuadDof > 0) ? new int[8*QuadDof] : nullptr;
2557  for (int i = 1; i < 8; i++)
2558  {
2559  QuadDofOrd[i] = QuadDofOrd[i-1] + QuadDof;
2560  }
2561  // see Mesh::GetQuadOrientation in mesh/mesh.cpp
2562  for (int j = 0; j <= p; j++)
2563  {
2564  for (int i = 0; i <= p; i++)
2565  {
2566  int o = i + j*pp1;
2567  QuadDofOrd[0][o] = i + j*pp1; // (0,1,2,3)
2568  QuadDofOrd[1][o] = -1 - (j + i*pp1); // (0,3,2,1)
2569  QuadDofOrd[2][o] = j + (p - i)*pp1; // (1,2,3,0)
2570  QuadDofOrd[3][o] = -1 - ((p - i) + j*pp1); // (1,0,3,2)
2571  QuadDofOrd[4][o] = (p - i) + (p - j)*pp1; // (2,3,0,1)
2572  QuadDofOrd[5][o] = -1 - ((p - j) + (p - i)*pp1); // (2,1,0,3)
2573  QuadDofOrd[6][o] = (p - j) + i*pp1; // (3,0,1,2)
2574  QuadDofOrd[7][o] = -1 - (i + (p - j)*pp1); // (3,2,1,0)
2575  if (!signs)
2576  {
2577  for (int k = 1; k < 8; k += 2)
2578  {
2579  QuadDofOrd[k][o] = -1 - QuadDofOrd[k][o];
2580  }
2581  }
2582  }
2583  }
2584  }
2585 }
2586 
2587 const FiniteElement *
2589 {
2590  if (GeomType != Geometry::PYRAMID || this->GetOrder() == 1)
2591  {
2592  return RT_Elements[GeomType];
2593  }
2594  else
2595  {
2596  if (error_mode == RETURN_NULL) { return nullptr; }
2597  MFEM_ABORT("RT Pyramid basis functions are not yet supported "
2598  "for order > 0.");
2599  return NULL;
2600  }
2601 }
2602 
2604  int Or) const
2605 {
2606  if (GeomType == Geometry::SEGMENT)
2607  {
2608  return (Or > 0) ? SegDofOrd[0] : SegDofOrd[1];
2609  }
2610  else if (GeomType == Geometry::TRIANGLE)
2611  {
2612  return TriDofOrd[Or%6];
2613  }
2614  else if (GeomType == Geometry::SQUARE)
2615  {
2616  return QuadDofOrd[Or%8];
2617  }
2618  return NULL;
2619 }
2620 
2622 {
2623  int tr_dim, tr_p;
2624  if (!strncmp(rt_name, "RT_", 3))
2625  {
2626  tr_dim = atoi(rt_name + 3);
2627  tr_p = atoi(rt_name + 7);
2628  }
2629  else // rt_name = RT@.._.D_P*
2630  {
2631  tr_dim = atoi(rt_name + 6);
2632  tr_p = atoi(rt_name + 10);
2633  }
2634  return new RT_Trace_FECollection(tr_p, tr_dim, FiniteElement::INTEGRAL,
2635  ob_type);
2636 }
2637 
2639 {
2640  delete [] SegDofOrd[0];
2641  delete [] TriDofOrd[0];
2642  delete [] QuadDofOrd[0];
2643  for (int g = 0; g < Geometry::NumGeom; g++)
2644  {
2645  delete RT_Elements[g];
2646  }
2647 }
2648 
2649 
2651  const int map_type,
2652  const int ob_type)
2653  : RT_FECollection(p, dim, map_type, true, ob_type)
2654 {
2655  const char *prefix =
2656  (map_type == FiniteElement::INTEGRAL) ? "RT_Trace" : "RT_ValTrace";
2657  char ob_str[3] = { '\0', '\0', '\0' };
2658 
2660  {
2661  ob_str[0] = '@';
2662  ob_str[1] = BasisType::GetChar(ob_type);
2663  }
2664  snprintf(rt_name, 32, "%s%s_%dD_P%d", prefix, ob_str, dim, p);
2665 
2666  MFEM_VERIFY(dim == 2 || dim == 3, "Wrong dimension, dim = " << dim);
2667 }
2668 
2669 
2671  const int map_type,
2672  const int ob_type)
2673  : RT_FECollection(p, dim, map_type, false, ob_type)
2674 {
2675  MFEM_VERIFY(dim == 2 || dim == 3, "Wrong dimension, dim = " << dim);
2676 
2677  const char *prefix =
2678  (map_type == FiniteElement::VALUE) ? "DG_Iface" : "DG_IntIface";
2680  {
2681  snprintf(rt_name, 32, "%s_%dD_P%d", prefix, dim, p);
2682  }
2683  else
2684  {
2685  snprintf(rt_name, 32, "%s@%c_%dD_P%d", prefix,
2686  (int)BasisType::GetChar(ob_type), dim, p);
2687  }
2688 }
2689 
2691  const int cb_type, const int ob_type)
2692  : FiniteElementCollection(dim > 1 ? p : p - 1)
2693  , dim(dim)
2694  , cb_type(cb_type)
2695  , ob_type(ob_type)
2696 {
2697  MFEM_VERIFY(p >= 1, "ND_FECollection requires order >= 1.");
2698  MFEM_VERIFY(dim >= 1 && dim <= 3, "ND_FECollection requires 1 <= dim <= 3.");
2699 
2700  const int pm1 = p - 1, pm2 = p - 2;
2701 
2704  {
2705  snprintf(nd_name, 32, "ND_%dD_P%d", dim, p);
2706  }
2707  else
2708  {
2709  snprintf(nd_name, 32, "ND@%c%c_%dD_P%d", (int)BasisType::GetChar(cb_type),
2710  (int)BasisType::GetChar(ob_type), dim, p);
2711  }
2712 
2713  for (int g = 0; g < Geometry::NumGeom; g++)
2714  {
2715  ND_Elements[g] = NULL;
2716  ND_dof[g] = 0;
2717  }
2718  for (int i = 0; i < 2; i++)
2719  {
2720  SegDofOrd[i] = NULL;
2721  }
2722  for (int i = 0; i < 6; i++)
2723  {
2724  TriDofOrd[i] = NULL;
2725  }
2726  for (int i = 0; i < 8; i++)
2727  {
2728  QuadDofOrd[i] = NULL;
2729  }
2730 
2731  int op_type = BasisType::GetQuadrature1D(ob_type);
2732  int cp_type = BasisType::GetQuadrature1D(cb_type);
2733 
2734  // Error checking
2737  {
2738  const char *ob_name = BasisType::Name(ob_type);
2739  MFEM_ABORT("Invalid open basis point type: " << ob_name);
2740  }
2742  {
2743  const char *cb_name = BasisType::Name(cb_type);
2744  MFEM_ABORT("Invalid closed basis point type: " << cb_name);
2745  }
2746 
2747  if (dim >= 1)
2748  {
2751 
2752  SegDofOrd[0] = (p > 0) ? new int[2*p] : nullptr;
2753  SegDofOrd[1] = SegDofOrd[0] + p;
2754  for (int i = 0; i < p; i++)
2755  {
2756  SegDofOrd[0][i] = i;
2757  SegDofOrd[1][i] = -1 - (pm1 - i);
2758  }
2759  }
2760 
2761  if (dim >= 2)
2762  {
2764  ob_type);
2765  ND_dof[Geometry::SQUARE] = 2*p*pm1;
2766 
2767  // TODO: cb_type and ob_type for triangles
2769  ND_dof[Geometry::TRIANGLE] = p*pm1;
2770 
2771  int QuadDof = ND_dof[Geometry::SQUARE];
2772  QuadDofOrd[0] = (QuadDof > 0) ? new int[8*QuadDof] : nullptr;
2773  for (int i = 1; i < 8; i++)
2774  {
2775  QuadDofOrd[i] = QuadDofOrd[i-1] + QuadDof;
2776  }
2777  // see Mesh::GetQuadOrientation in mesh/mesh.cpp
2778  for (int j = 0; j < pm1; j++)
2779  {
2780  for (int i = 0; i < p; i++)
2781  {
2782  int d1 = i + j*p; // x-component
2783  int d2 = p*pm1 + j + i*pm1; // y-component
2784  // (0,1,2,3)
2785  QuadDofOrd[0][d1] = d1;
2786  QuadDofOrd[0][d2] = d2;
2787  // (0,3,2,1)
2788  QuadDofOrd[1][d1] = d2;
2789  QuadDofOrd[1][d2] = d1;
2790  // (1,2,3,0)
2791  // QuadDofOrd[2][d1] = p*pm1 + (pm2 - j) + i*pm1;
2792  // QuadDofOrd[2][d2] = -1 - ((pm1 - i) + j*p);
2793  QuadDofOrd[2][d1] = -1 - (p*pm1 + j + (pm1 - i)*pm1);
2794  QuadDofOrd[2][d2] = i + (pm2 - j)*p;
2795  // (1,0,3,2)
2796  QuadDofOrd[3][d1] = -1 - ((pm1 - i) + j*p);
2797  QuadDofOrd[3][d2] = p*pm1 + (pm2 - j) + i*pm1;
2798  // (2,3,0,1)
2799  QuadDofOrd[4][d1] = -1 - ((pm1 - i) + (pm2 - j)*p);
2800  QuadDofOrd[4][d2] = -1 - (p*pm1 + (pm2 - j) + (pm1 - i)*pm1);
2801  // (2,1,0,3)
2802  QuadDofOrd[5][d1] = -1 - (p*pm1 + (pm2 - j) + (pm1 - i)*pm1);
2803  QuadDofOrd[5][d2] = -1 - ((pm1 - i) + (pm2 - j)*p);
2804  // (3,0,1,2)
2805  // QuadDofOrd[6][d1] = -1 - (p*pm1 + j + (pm1 - i)*pm1);
2806  // QuadDofOrd[6][d2] = i + (pm2 - j)*p;
2807  QuadDofOrd[6][d1] = p*pm1 + (pm2 - j) + i*pm1;
2808  QuadDofOrd[6][d2] = -1 - ((pm1 - i) + j*p);
2809  // (3,2,1,0)
2810  QuadDofOrd[7][d1] = i + (pm2 - j)*p;
2811  QuadDofOrd[7][d2] = -1 - (p*pm1 + j + (pm1 - i)*pm1);
2812  }
2813  }
2814 
2815  int TriDof = ND_dof[Geometry::TRIANGLE];
2816  TriDofOrd[0] = (TriDof > 0) ? new int[6*TriDof] : nullptr;
2817  for (int i = 1; i < 6; i++)
2818  {
2819  TriDofOrd[i] = TriDofOrd[i-1] + TriDof;
2820  }
2821  // see Mesh::GetTriOrientation in mesh/mesh.cpp,
2822  // the constructor of H1_FECollection
2823  for (int j = 0; j <= pm2; j++)
2824  {
2825  for (int i = 0; i + j <= pm2; i++)
2826  {
2827  int k0 = p*pm1 - (p - j)*(pm1 - j) + 2*i;
2828  int k1 = 2*pm2 - 2*i + ((2*p-3)-j)*j;
2829  int k2 = 2*pm2 - 2*j + ((2*p-3)-i)*i;
2830  int k3 = p*pm1 - 2 - 3*j - i - (i+j)*(i+j);
2831  int k4 = p*pm1 - 2 - 3*i - j - (i+j)*(i+j);
2832  int k5 = p*pm1 - (p - i)*(pm1 - i) + 2*j;
2833 
2834  // (0,1,2)
2835  TriDofOrd[0][k0 ] = k0;
2836  TriDofOrd[0][k0+1] = k0 + 1;
2837  // (1,0,2)
2838  TriDofOrd[1][k0 ] = k1;
2839  TriDofOrd[1][k0+1] = k1 + 1;
2840  // (2,0,1)
2841  TriDofOrd[2][k0 ] = k2;
2842  TriDofOrd[2][k0+1] = k2 + 1;
2843  // (2,1,0)
2844  TriDofOrd[3][k0 ] = k3;
2845  TriDofOrd[3][k0+1] = k3 + 1;
2846  // (1,2,0)
2847  TriDofOrd[4][k0 ] = k4;
2848  TriDofOrd[4][k0+1] = k4 + 1;
2849  // (0,2,1)
2850  TriDofOrd[5][k0 ] = k5;
2851  TriDofOrd[5][k0+1] = k5 + 1;
2852  }
2853  }
2854  }
2855 
2856  if (dim >= 3)
2857  {
2859  ND_dof[Geometry::CUBE] = 3*p*pm1*pm1;
2860 
2861  // TODO: cb_type and ob_type for tets
2863  ND_dof[Geometry::TETRAHEDRON] = p*pm1*pm2/2;
2864 
2866  ND_dof[Geometry::PRISM] = p*pm1*(3*p-4)/2;
2867 
2870  }
2871 }
2872 
2873 const FiniteElement *
2875 {
2876  if (GeomType != Geometry::PYRAMID || this->GetOrder() == 1)
2877  {
2878  return ND_Elements[GeomType];
2879  }
2880  else
2881  {
2882  if (error_mode == RETURN_NULL) { return nullptr; }
2883  MFEM_ABORT("ND Pyramid basis functions are not yet supported "
2884  "for order > 1.");
2885  return NULL;
2886  }
2887 }
2888 
2891 {
2892  if (!Geometry::IsTensorProduct(GeomType) && this->GetOrder() > 1)
2893  {
2894  return FiniteElementForGeometry(GeomType)->GetDofTransformation();
2895  }
2896  else
2897  {
2898  return NULL;
2899  }
2900 }
2901 
2903  int Or) const
2904 {
2905  if (GeomType == Geometry::SEGMENT)
2906  {
2907  return (Or > 0) ? SegDofOrd[0] : SegDofOrd[1];
2908  }
2909  else if (GeomType == Geometry::TRIANGLE)
2910  {
2911  return TriDofOrd[Or%6];
2912  }
2913  else if (GeomType == Geometry::SQUARE)
2914  {
2915  return QuadDofOrd[Or%8];
2916  }
2917  return NULL;
2918 }
2919 
2921 {
2922  int tr_p, tr_dim, tr_cb_type, tr_ob_type;
2923 
2924  tr_p = ND_dof[Geometry::SEGMENT];
2925  if (nd_name[2] == '_') // ND_
2926  {
2927  tr_dim = atoi(nd_name + 3);
2928  tr_cb_type = BasisType::GaussLobatto;
2929  tr_ob_type = BasisType::GaussLegendre;
2930  }
2931  else // ND@
2932  {
2933  tr_dim = atoi(nd_name + 6);
2934  tr_cb_type = BasisType::GetType(nd_name[3]);
2935  tr_ob_type = BasisType::GetType(nd_name[4]);
2936  }
2937  return new ND_Trace_FECollection(tr_p, tr_dim, tr_cb_type, tr_ob_type);
2938 }
2939 
2941 {
2942  delete [] SegDofOrd[0];
2943  delete [] TriDofOrd[0];
2944  delete [] QuadDofOrd[0];
2945  for (int g = 0; g < Geometry::NumGeom; g++)
2946  {
2947  delete ND_Elements[g];
2948  }
2949 }
2950 
2951 
2953  const int cb_type,
2954  const int ob_type)
2955  : ND_FECollection(p, dim-1, cb_type, ob_type)
2956 {
2959  {
2960  snprintf(nd_name, 32, "ND_Trace_%dD_P%d", dim, p);
2961  }
2962  else
2963  {
2964  snprintf(nd_name, 32, "ND_Trace@%c%c_%dD_P%d",
2966  (int)BasisType::GetChar(ob_type), dim, p);
2967  }
2968 }
2969 
2970 
2972  const int cb_type, const int ob_type)
2974 {
2975  MFEM_VERIFY(p >= 1, "ND_R1D_FECollection requires order >= 1.");
2976  MFEM_VERIFY(dim == 1, "ND_R1D_FECollection requires dim == 1.");
2977 
2978  if (cb_type == BasisType::GaussLobatto &&
2979  ob_type == BasisType::GaussLegendre)
2980  {
2981  snprintf(nd_name, 32, "ND_R1D_%dD_P%d", dim, p);
2982  }
2983  else
2984  {
2985  snprintf(nd_name, 32, "ND_R1D@%c%c_%dD_P%d",
2986  (int)BasisType::GetChar(cb_type),
2987  (int)BasisType::GetChar(ob_type), dim, p);
2988  }
2989 
2990  for (int g = 0; g < Geometry::NumGeom; g++)
2991  {
2992  ND_Elements[g] = NULL;
2993  ND_dof[g] = 0;
2994  }
2995 
2996  int op_type = BasisType::GetQuadrature1D(ob_type);
2997  int cp_type = BasisType::GetQuadrature1D(cb_type);
2998 
2999  // Error checking
3001  {
3002  const char *ob_name = BasisType::Name(ob_type);
3003  MFEM_ABORT("Invalid open basis point type: " << ob_name);
3004  }
3006  {
3007  const char *cb_name = BasisType::Name(cb_type);
3008  MFEM_ABORT("Invalid closed basis point type: " << cb_name);
3009  }
3010 
3012  ND_dof[Geometry::POINT] = 2;
3013 
3015  cb_type,
3016  ob_type);
3017  ND_dof[Geometry::SEGMENT] = 3 * p - 2;
3018 }
3019 
3021  int Or) const
3022 {
3023  return NULL;
3024 }
3025 
3027 {
3028  return NULL;
3029 }
3030 
3032 {
3033  for (int g = 0; g < Geometry::NumGeom; g++)
3034  {
3035  delete ND_Elements[g];
3036  }
3037 }
3038 
3039 
3041  const int cb_type, const int ob_type)
3042  : FiniteElementCollection(p + 1)
3043 {
3044  MFEM_VERIFY(p >= 0, "RT_R1D_FECollection requires order >= 0.");
3045  MFEM_VERIFY(dim == 1, "RT_R1D_FECollection requires dim == 1.");
3046 
3047  if (cb_type == BasisType::GaussLobatto &&
3048  ob_type == BasisType::GaussLegendre)
3049  {
3050  snprintf(rt_name, 32, "RT_R1D_%dD_P%d", dim, p);
3051  }
3052  else
3053  {
3054  snprintf(rt_name, 32, "RT_R1D@%c%c_%dD_P%d",
3055  (int)BasisType::GetChar(cb_type),
3056  (int)BasisType::GetChar(ob_type), dim, p);
3057  }
3058 
3059  for (int g = 0; g < Geometry::NumGeom; g++)
3060  {
3061  RT_Elements[g] = NULL;
3062  RT_dof[g] = 0;
3063  }
3064 
3065  int op_type = BasisType::GetQuadrature1D(ob_type);
3066  int cp_type = BasisType::GetQuadrature1D(cb_type);
3067 
3068  // Error checking
3070  {
3071  const char *ob_name = BasisType::Name(ob_type);
3072  MFEM_ABORT("Invalid open basis point type: " << ob_name);
3073  }
3075  {
3076  const char *cb_name = BasisType::Name(cb_type);
3077  MFEM_ABORT("Invalid closed basis point type: " << cb_name);
3078  }
3079 
3081  RT_dof[Geometry::POINT] = 1;
3082 
3084  cb_type,
3085  ob_type);
3086  RT_dof[Geometry::SEGMENT] = 3 * p + 2;
3087 }
3088 
3090  int Or) const
3091 {
3092  return NULL;
3093 }
3094 
3096 {
3097  MFEM_ABORT("this method is not implemented in RT_R1D_FECollection!");
3098  return NULL;
3099 }
3100 
3102 {
3103  for (int g = 0; g < Geometry::NumGeom; g++)
3104  {
3105  delete RT_Elements[g];
3106  }
3107 }
3108 
3109 
3111  const int cb_type, const int ob_type)
3113 {
3114  MFEM_VERIFY(p >= 1, "ND_R2D_FECollection requires order >= 1.");
3115  MFEM_VERIFY(dim >= 1 && dim <= 2,
3116  "ND_R2D_FECollection requires 1 <= dim <= 2.");
3117 
3118  const int pm1 = p - 1, pm2 = p - 2;
3119 
3120  if (cb_type == BasisType::GaussLobatto &&
3121  ob_type == BasisType::GaussLegendre)
3122  {
3123  snprintf(nd_name, 32, "ND_R2D_%dD_P%d", dim, p);
3124  }
3125  else
3126  {
3127  snprintf(nd_name, 32, "ND_R2D@%c%c_%dD_P%d",
3128  (int)BasisType::GetChar(cb_type),
3129  (int)BasisType::GetChar(ob_type), dim, p);
3130  }
3131 
3132  for (int g = 0; g < Geometry::NumGeom; g++)
3133  {
3134  ND_Elements[g] = NULL;
3135  ND_dof[g] = 0;
3136  }
3137  for (int i = 0; i < 2; i++)
3138  {
3139  SegDofOrd[i] = NULL;
3140  }
3141 
3142  int op_type = BasisType::GetQuadrature1D(ob_type);
3143  int cp_type = BasisType::GetQuadrature1D(cb_type);
3144 
3145  // Error checking
3147  {
3148  const char *ob_name = BasisType::Name(ob_type);
3149  MFEM_ABORT("Invalid open basis point type: " << ob_name);
3150  }
3152  {
3153  const char *cb_name = BasisType::Name(cb_type);
3154  MFEM_ABORT("Invalid closed basis point type: " << cb_name);
3155  }
3156 
3157  ND_dof[Geometry::POINT] = 1;
3158 
3159  if (dim >= 1)
3160  {
3162  cb_type,
3163  ob_type);
3164  ND_dof[Geometry::SEGMENT] = 2 * p - 1;
3165 
3166  SegDofOrd[0] = (4*p > 2) ? new int[4 * p - 2] : nullptr;
3167  SegDofOrd[1] = SegDofOrd[0] + 2 * p - 1;
3168  for (int i = 0; i < p; i++)
3169  {
3170  SegDofOrd[0][i] = i;
3171  SegDofOrd[1][i] = -1 - (pm1 - i);
3172  }
3173  for (int i = 0; i < pm1; i++)
3174  {
3175  SegDofOrd[0][p+i] = p + i;
3176  SegDofOrd[1][p+i] = 2 * pm1 - i;
3177  }
3178  }
3179 
3180  if (dim >= 2)
3181  {
3183  cb_type,
3184  ob_type);
3185  ND_dof[Geometry::SQUARE] = 2*p*pm1 + pm1*pm1;
3186 
3187  // TODO: cb_type and ob_type for triangles
3189  ND_dof[Geometry::TRIANGLE] = p*pm1 + (pm1*pm2)/2;
3190  }
3191 }
3192 
3194  int Or) const
3195 {
3196  if (GeomType == Geometry::SEGMENT)
3197  {
3198  return (Or > 0) ? SegDofOrd[0] : SegDofOrd[1];
3199  }
3200  return NULL;
3201 }
3202 
3204 {
3205  int p, dim, cb_type, ob_type;
3206 
3208  if (nd_name[5] == '_') // ND_R2D_
3209  {
3210  dim = atoi(nd_name + 6);
3211  cb_type = BasisType::GaussLobatto;
3212  ob_type = BasisType::GaussLegendre;
3213  }
3214  else // ND_R2D@
3215  {
3216  dim = atoi(nd_name + 9);
3217  cb_type = BasisType::GetType(nd_name[6]);
3218  ob_type = BasisType::GetType(nd_name[7]);
3219  }
3220  return new ND_R2D_Trace_FECollection(p, dim, cb_type, ob_type);
3221 }
3222 
3224 {
3225  delete [] SegDofOrd[0];
3226  for (int g = 0; g < Geometry::NumGeom; g++)
3227  {
3228  delete ND_Elements[g];
3229  }
3230 }
3231 
3232 
3234  const int cb_type,
3235  const int ob_type)
3236  : ND_R2D_FECollection(p, dim-1, cb_type, ob_type)
3237 {
3238  if (cb_type == BasisType::GaussLobatto &&
3239  ob_type == BasisType::GaussLegendre)
3240  {
3241  snprintf(nd_name, 32, "ND_R2D_Trace_%dD_P%d", dim, p);
3242  }
3243  else
3244  {
3245  snprintf(nd_name, 32, "ND_R2D_Trace@%c%c_%dD_P%d",
3246  (int)BasisType::GetChar(cb_type),
3247  (int)BasisType::GetChar(ob_type), dim, p);
3248  }
3249 }
3250 
3251 
3253  const int cb_type, const int ob_type)
3254  : FiniteElementCollection(p + 1),
3255  ob_type(ob_type)
3256 {
3257  MFEM_VERIFY(p >= 0, "RT_R2D_FECollection requires order >= 0.");
3258  MFEM_VERIFY(dim >= 1 && dim <= 2,
3259  "RT_R2D_FECollection requires 1 <= dim <= 2.");
3260 
3261  int cp_type = BasisType::GetQuadrature1D(cb_type);
3262  int op_type = BasisType::GetQuadrature1D(ob_type);
3263 
3265  {
3266  const char *cb_name = BasisType::Name(cb_type); // this may abort
3267  MFEM_ABORT("unknown closed BasisType: " << cb_name);
3268  }
3270  {
3271  const char *ob_name = BasisType::Name(ob_type); // this may abort
3272  MFEM_ABORT("unknown open BasisType: " << ob_name);
3273  }
3274 
3276 
3277  if (cb_type == BasisType::GaussLobatto &&
3279  {
3280  snprintf(rt_name, 32, "RT_R2D_%dD_P%d", dim, p);
3281  }
3282  else
3283  {
3284  snprintf(rt_name, 32, "RT_R2D@%c%c_%dD_P%d",
3285  (int)BasisType::GetChar(cb_type),
3286  (int)BasisType::GetChar(ob_type), dim, p);
3287  }
3288 
3289  const int pp1 = p + 1;
3290  const int pp2 = p + 2;
3291  if (dim == 2)
3292  {
3293  // TODO: cb_type, ob_type for triangles
3295  RT_dof[Geometry::TRIANGLE] = p*pp1 + (pp1 * pp2) / 2;
3296 
3298  cb_type,
3299  ob_type);
3300  // two vector components * n_unk_face *
3301  RT_dof[Geometry::SQUARE] = 2*p*pp1 + pp1*pp1;
3302  }
3303 }
3304 
3305 // This is a special protected constructor only used by RT_Trace_FECollection
3306 // and DG_Interface_FECollection
3308  const int map_type,
3309  const bool signs, const int ob_type)
3310  : ob_type(ob_type)
3311 {
3314  {
3315  const char *ob_name = BasisType::Name(ob_type); // this may abort
3316  MFEM_ABORT("Invalid open basis type: " << ob_name);
3317  }
3318  InitFaces(p, dim, map_type, signs);
3319 }
3320 
3321 void RT_R2D_FECollection::InitFaces(const int p, const int dim,
3322  const int map_type,
3323  const bool signs)
3324 {
3325  int op_type = BasisType::GetQuadrature1D(ob_type);
3326 
3327  MFEM_VERIFY(Quadrature1D::CheckOpen(op_type) != Quadrature1D::Invalid,
3328  "invalid open point type");
3329 
3330  const int pp1 = p + 1;
3331 
3332  for (int g = 0; g < Geometry::NumGeom; g++)
3333  {
3334  RT_Elements[g] = NULL;
3335  RT_dof[g] = 0;
3336  }
3337  // Degree of Freedom orderings
3338  for (int i = 0; i < 2; i++)
3339  {
3340  SegDofOrd[i] = NULL;
3341  }
3342 
3343  if (dim == 2)
3344  {
3345  L2_SegmentElement *l2_seg = new L2_SegmentElement(p, ob_type);
3346  l2_seg->SetMapType(map_type);
3347  RT_Elements[Geometry::SEGMENT] = l2_seg;
3348  RT_dof[Geometry::SEGMENT] = pp1;
3349 
3350  SegDofOrd[0] = (pp1 > 0) ? new int[2*pp1] : nullptr;
3351  SegDofOrd[1] = SegDofOrd[0] + pp1;
3352  for (int i = 0; i <= p; i++)
3353  {
3354  SegDofOrd[0][i] = i;
3355  SegDofOrd[1][i] = signs ? (-1 - (p - i)) : (p - i);
3356  }
3357  }
3358 }
3359 
3361  int Or) const
3362 {
3363  if (GeomType == Geometry::SEGMENT)
3364  {
3365  return (Or > 0) ? SegDofOrd[0] : SegDofOrd[1];
3366  }
3367  return NULL;
3368 }
3369 
3371 {
3372  int dim, p;
3373  if (!strncmp(rt_name, "RT_R2D_", 7))
3374  {
3375  dim = atoi(rt_name + 7);
3376  p = atoi(rt_name + 11);
3377  }
3378  else // rt_name = RT_R2D@.._.D_P*
3379  {
3380  dim = atoi(rt_name + 10);
3381  p = atoi(rt_name + 14);
3382  }
3384 }
3385 
3387 {
3388  delete [] SegDofOrd[0];
3389  for (int g = 0; g < Geometry::NumGeom; g++)
3390  {
3391  delete RT_Elements[g];
3392  }
3393 }
3394 
3395 
3397  const int map_type,
3398  const int ob_type)
3399  : RT_R2D_FECollection(p, dim-1, map_type, true, ob_type)
3400 {
3401  const char *prefix =
3402  (map_type == FiniteElement::INTEGRAL) ? "RT_R2D_Trace" : "RT_R2D_ValTrace";
3403  char ob_str[3] = { '\0', '\0', '\0' };
3404 
3406  {
3407  ob_str[0] = '@';
3408  ob_str[1] = BasisType::GetChar(ob_type);
3409  }
3410  snprintf(rt_name, 32, "%s%s_%dD_P%d", prefix, ob_str, dim, p);
3411 
3412  MFEM_VERIFY(dim == 2, "Wrong dimension, dim = " << dim);
3413 }
3414 
3415 
3417 {
3418  snprintf(d_name, 32, "Local_%s", fe_name);
3419 
3420  Local_Element = NULL;
3421 
3422  if (!strcmp(fe_name, "BiCubic2DFiniteElement") ||
3423  !strcmp(fe_name, "Quad_Q3"))
3424  {
3425  GeomType = Geometry::SQUARE;
3426  Local_Element = new BiCubic2DFiniteElement;
3427  }
3428  else if (!strcmp(fe_name, "Nedelec1HexFiniteElement") ||
3429  !strcmp(fe_name, "Hex_ND1"))
3430  {
3431  GeomType = Geometry::CUBE;
3432  Local_Element = new Nedelec1HexFiniteElement;
3433  }
3434  else if (!strncmp(fe_name, "H1_", 3))
3435  {
3436  GeomType = Geometry::SQUARE;
3437  Local_Element = new H1_QuadrilateralElement(atoi(fe_name + 7));
3438  }
3439  else if (!strncmp(fe_name, "H1Pos_", 6))
3440  {
3441  GeomType = Geometry::SQUARE;
3442  Local_Element = new H1Pos_QuadrilateralElement(atoi(fe_name + 10));
3443  }
3444  else if (!strncmp(fe_name, "L2_", 3))
3445  {
3446  GeomType = Geometry::SQUARE;
3447  Local_Element = new L2_QuadrilateralElement(atoi(fe_name + 7));
3448  }
3449  else
3450  {
3451  mfem::err << "Local_FECollection::Local_FECollection : fe_name = "
3452  << fe_name << endl;
3453  mfem_error();
3454  }
3455 }
3456 
3457 
3459  : FiniteElementCollection((Order == VariableOrder) ? 1 : Order)
3460 {
3461  const int order = (Order == VariableOrder) ? 1 : Order;
3462  PointFE = new PointFiniteElement();
3463  SegmentFE = new NURBS1DFiniteElement(order);
3464  QuadrilateralFE = new NURBS2DFiniteElement(order);
3465  ParallelepipedFE = new NURBS3DFiniteElement(order);
3466 
3467  SetOrder(Order);
3468 }
3469 
3470 void NURBSFECollection::SetOrder(int Order) const
3471 {
3472  mOrder = Order;
3473  if (Order != VariableOrder)
3474  {
3475  snprintf(name, 16, "NURBS%i", Order);
3476  }
3477  else
3478  {
3479  snprintf(name, 16, "NURBS");
3480  }
3481 }
3482 
3484 {
3485  delete PointFE;
3486  delete SegmentFE;
3487  delete QuadrilateralFE;
3488  delete ParallelepipedFE;
3489 }
3490 
3491 const FiniteElement *
3493 {
3494  switch (GeomType)
3495  {
3496  case Geometry::POINT: return PointFE;
3497  case Geometry::SEGMENT: return SegmentFE;
3498  case Geometry::SQUARE: return QuadrilateralFE;
3499  case Geometry::CUBE: return ParallelepipedFE;
3500  default:
3501  if (error_mode == RETURN_NULL) { return nullptr; }
3502  mfem_error ("NURBSFECollection: unknown geometry type.");
3503  }
3504  return SegmentFE; // Make some compilers happy
3505 }
3506 
3508 {
3509  mfem_error("NURBSFECollection::DofForGeometry");
3510  return 0; // Make some compilers happy
3511 }
3512 
3514  int Or) const
3515 {
3516  mfem_error("NURBSFECollection::DofOrderForOrientation");
3517  return NULL;
3518 }
3519 
3521 {
3522  MFEM_ABORT("NURBS finite elements can not be statically condensed!");
3523  return NULL;
3524 }
3525 
3526 }
Abstract class for all finite elements.
Definition: fe_base.hpp:233
Arbitrary order non-uniform rational B-splines (NURBS) finite elements.
Definition: fe_coll.hpp:649
A 0D Nedelec finite element for the boundary of a 1D domain.
Definition: fe_nd.hpp:416
H1_Trace_FECollection(const int p, const int dim, const int btype=BasisType::GaussLobatto)
Definition: fe_coll.cpp:2045
Arbitrary order L2 elements in 3D on a wedge.
Definition: fe_l2.hpp:165
Arbitrary order Nedelec elements in 1D on a segment.
Definition: fe_nd.hpp:291
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1176
Piecewise-cubic discontinuous finite elements in 2D. This class is kept only for backward compatibili...
Definition: fe_coll.hpp:1095
Version of QuadraticDiscont2DFECollection with dofs in the Gaussian points.
Definition: fe_coll.hpp:1070
BiLinear2DFiniteElement QuadrilateralFE
int RT_dof[Geometry::NumGeom]
Definition: fe_coll.hpp:388
FiniteElement * ND_Elements[Geometry::NumGeom]
Definition: fe_coll.hpp:461
void SubDofOrder(Geometry::Type Geom, int SDim, int Info, Array< int > &dofs) const
Get the local dofs for a given sub-manifold.
Definition: fe_coll.cpp:462
FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:3095
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1462
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1420
static int Check(int b_type)
If the input does not represents a valid BasisType, abort with an error; otherwise return the input...
Definition: fe_base.hpp:45
int RT_dof[Geometry::NumGeom]
Definition: fe_coll.hpp:539
Arbitrary order L2 elements in 2D utilizing the Bernstein basis on a triangle.
Definition: fe_pos.hpp:314
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1084
virtual ~NURBSFECollection()
Definition: fe_coll.cpp:3483
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:870
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:1618
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:71
RT_FECollection(const int p, const int dim, const int map_type, const bool signs, const int ob_type=BasisType::GaussLegendre)
Definition: fe_coll.cpp:2453
Arbitrary order L2 elements in 2D on a triangle.
Definition: fe_l2.hpp:108
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:983
ErrorMode
How to treat errors in FiniteElementForGeometry() calls.
Definition: fe_coll.hpp:245
Linear1DFiniteElement SegmentFE
Definition: segment.cpp:49
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1603
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:3193
Arbitrary order Nedelec 3D elements in 2D on a square.
Definition: fe_nd.hpp:646
class LinearPyramidFiniteElement PyramidFE
Definition: fe.cpp:44
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:884
static const int NumGeom
Definition: geom.hpp:42
A 3D 1st order Nedelec element on a pyramid.
Arbitrary order L2 elements in 3D on a cube.
Definition: fe_l2.hpp:78
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:2874
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:927
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:1437
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1488
virtual ~L2_FECollection()
Definition: fe_coll.cpp:2366
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:1478
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:1356
Arbitrary order L2 elements in 1D utilizing the Bernstein basis on a segment.
Definition: fe_pos.hpp:260
int GetOrder() const
Return the order (polynomial degree) of the FE collection, corresponding to the order/degree returned...
Definition: fe_coll.hpp:225
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:679
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1951
T * GetData()
Returns the data.
Definition: array.hpp:115
Arbitrary order H1 elements in 2D on a square.
Definition: fe_h1.hpp:41
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:3360
Piecewise-(bi/tri)linear continuous finite elements.
Definition: fe_coll.hpp:707
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:2603
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:3020
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:634
A 0D point finite element.
Arbitrary order 3D H(curl)-trace finite elements in 2D defined on the interface between mesh elements...
Definition: fe_coll.hpp:591
static int GetQuadrature1D(int b_type)
Get the corresponding Quadrature1D constant, when that makes sense; otherwise return Quadrature1D::In...
Definition: fe_base.hpp:61
Arbitrary order Nedelec elements in 3D on a tetrahedron.
Definition: fe_nd.hpp:170
Arbitrary order 3D H(curl)-conforming Nedelec finite elements in 2D.
Definition: fe_coll.hpp:561
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:1237
Arbitrary order, three component, Raviart-Thomas elements in 1D on a segment.
Definition: fe_rt.hpp:338
virtual ~RT_R2D_FECollection()
Definition: fe_coll.cpp:3386
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1069
TriLinear3DFiniteElement HexahedronFE
Definition: hexahedron.cpp:52
const int base_p
Order as returned by GetOrder().
Definition: fe_coll.hpp:235
Second order Raviart-Thomas finite elements in 2D. This class is kept only for backward compatibility...
Definition: fe_coll.hpp:892
void InitFaces(const int p, const int dim, const int map_type, const bool signs)
Definition: fe_coll.cpp:2468
STL namespace.
Arbitrary order, three component, Nedelec elements in 1D on a segment.
Definition: fe_nd.hpp:437
DG_Interface_FECollection(const int p, const int dim, const int map_type=FiniteElement::VALUE, const int ob_type=BasisType::GaussLegendre)
Definition: fe_coll.cpp:2670
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1553
RT_R2D_FECollection(const int p, const int dim, const int map_type, const bool signs, const int ob_type=BasisType::GaussLegendre)
Definition: fe_coll.cpp:3307
static const char * Name(int b_type)
Check and convert a BasisType constant to a string identifier.
Definition: fe_base.hpp:92
A 2D bi-cubic element on a square with uniformly spaces nodes.
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:796
Arbitrary order Nedelec elements in 2D on a triangle.
Definition: fe_nd.hpp:232
static int GetType(char b_ident)
Convert char basis identifier to a BasisType constant.
Definition: fe_base.hpp:111
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:1100
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1118
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:1313
NURBSFECollection(int Order=VariableOrder)
The parameter Order must be either a positive number, for fixed order, or VariableOrder (default)...
Definition: fe_coll.cpp:3458
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1406
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:779
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:738
std::function< double(const Vector &)> f(double mass_coeff)
Definition: lor_mms.hpp:30
Arbitrary order Raviart-Thomas 3D elements in 2D on a triangle.
Definition: fe_rt.hpp:472
virtual int GetRangeType(int dim) const
Definition: fe_coll.cpp:40
virtual ~RT_FECollection()
Definition: fe_coll.cpp:2638
Arbitrary order Raviart-Thomas 3D elements in 2D on a square.
Definition: fe_rt.hpp:500
L2_FECollection(const int p, const int dim, const int btype=BasisType::GaussLegendre, const int map_type=FiniteElement::VALUE)
Definition: fe_coll.cpp:2065
virtual const FiniteElement * FiniteElementForDim(int dim) const
Returns the first non-NULL FiniteElement for the given dimension.
Definition: fe_coll.cpp:26
virtual FiniteElementCollection * Clone(int p) const
Instantiate a new collection of the same type with a different order.
Definition: fe_coll.cpp:361
Class for standard nodal finite elements.
Definition: fe_base.hpp:708
Piecewise-(bi)cubic continuous finite elements.
Definition: fe_coll.hpp:787
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:3089
Arbitrary order L2 elements in 2D utilizing the Bernstein basis on a square.
Definition: fe_pos.hpp:278
FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:1988
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:2332
virtual void SetMapType(int M)
Set the FiniteElement::MapType of the element to either VALUE or INTEGRAL. Also sets the FiniteElemen...
Definition: fe_base.hpp:675
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1323
int H1_dof[Geometry::NumGeom]
Definition: fe_coll.hpp:266
Arbitrary order L2 elements in 1D on a segment.
Definition: fe_l2.hpp:21
static void GetFace(int &nv, v_t &v, int &ne, e_t &e, eo_t &eo, int &nf, int &f, Geometry::Type &fg, int &fo, const int face_info)
Definition: fe_coll.cpp:417
FiniteElement * RT_Elements[Geometry::NumGeom]
Definition: fe_coll.hpp:387
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1297
static const int Dimension[NumGeom]
Definition: geom.hpp:47
static void GetNVE(int &nv, int &ne)
Definition: fe_coll.cpp:387
FiniteElement * ND_Elements[Geometry::NumGeom]
Definition: fe_coll.hpp:511
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:969
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1205
virtual int GetDerivType(int dim) const
Definition: fe_coll.cpp:70
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:860
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:913
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:653
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1221
Piecewise-linear discontinuous finite elements in 2D. This class is kept only for backward compatibil...
Definition: fe_coll.hpp:965
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:2902
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1139
FiniteElement * ND_Elements[Geometry::NumGeom]
Definition: fe_coll.hpp:565
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
H1_FECollection(const int p, const int dim=3, const int btype=BasisType::GaussLobatto)
Definition: fe_coll.cpp:1640
ND_Trace_FECollection(const int p, const int dim, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Definition: fe_coll.cpp:2952
Arbitrary order "H^{1/2}-conforming" trace finite elements defined on the interface between mesh elem...
Definition: fe_coll.hpp:319
ND_R2D_Trace_FECollection(const int p, const int dim, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Definition: fe_coll.cpp:3233
ND_FECollection(const int p, const int dim, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Definition: fe_coll.cpp:2690
ErrorMode error_mode
How to treat errors in FiniteElementForGeometry() calls.
Definition: fe_coll.hpp:255
Discontinuous collection defined locally by a given finite element.
Definition: fe_coll.hpp:1295
Arbitrary order H(curl)-trace finite elements defined on the interface between mesh elements (faces...
Definition: fe_coll.hpp:498
Version of QuadraticDiscont2DFECollection with positive basis functions.
Definition: fe_coll.hpp:1052
virtual ~ND_FECollection()
Definition: fe_coll.cpp:2940
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:941
Bernstein polynomials.
Definition: fe_base.hpp:33
RT_R1D_FECollection(const int p, const int dim, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Definition: fe_coll.cpp:3040
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1536
virtual FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:110
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1108
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:697
Array< FiniteElementCollection * > var_orders
Definition: fe_coll.hpp:242
Third order Raviart-Thomas finite elements in 2D. This class is kept only for backward compatibility...
Definition: fe_coll.hpp:916
PointFiniteElement PointFE
Definition: point.cpp:30
First order Raviart-Thomas finite elements in 3D. This class is kept only for backward compatibility...
Definition: fe_coll.hpp:1247
virtual const char * Name() const
Definition: fe_coll.hpp:80
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:3492
Version of QuadraticFECollection with positive basis functions.
Definition: fe_coll.hpp:764
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1338
virtual ~H1_FECollection()
Definition: fe_coll.cpp:2032
Piecewise-constant discontinuous finite elements in 3D. This class is kept only for backward compatib...
Definition: fe_coll.hpp:1119
Arbitrary order H1 elements in 3D on a tetrahedron.
Definition: fe_h1.hpp:105
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:1025
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:898
int ND_dof[Geometry::NumGeom]
Definition: fe_coll.hpp:512
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:1521
Piecewise-constant discontinuous finite elements in 2D. This class is kept only for backward compatib...
Definition: fe_coll.hpp:941
FiniteElement * H1_Elements[Geometry::NumGeom]
Definition: fe_coll.hpp:265
static const char * Name[NumGeom]
Definition: geom.hpp:45
Piecewise-linear nonconforming finite elements in 3D.
Definition: fe_coll.hpp:841
Arbitrary order H1 elements in 3D utilizing the Bernstein basis on a wedge.
Definition: fe_pos.hpp:237
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1189
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:1273
Crouzeix-Raviart nonconforming elements in 2D.
Definition: fe_coll.hpp:818
FiniteElement * RT_Elements[Geometry::NumGeom]
Definition: fe_coll.hpp:605
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1503
const Array< int > & GetLexicographicOrdering() const
Get an Array<int> that maps lexicographically ordered indices to the indices of the respective nodes/...
Definition: fe_base.hpp:783
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:380
Arbitrary order L2 elements in 3D on a tetrahedron.
Definition: fe_l2.hpp:138
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:3513
static void GetEdge(int &nv, v_t &v, int &ne, int &e, int &eo, const int edge_info)
Definition: fe_coll.cpp:397
static int CheckClosed(int type)
If the Quadrature1D type is not closed return Invalid; otherwise return type.
Definition: intrules.cpp:911
Arbitrary order Raviart-Thomas elements in 3D on a tetrahedron.
Definition: fe_rt.hpp:223
int GetDerivType() const
Returns the FiniteElement::DerivType of the element describing the spatial derivative method implemen...
Definition: fe_base.hpp:354
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1380
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:2588
Arbitrary order H1 elements in 3D utilizing the Bernstein basis on a cube.
Definition: fe_pos.hpp:161
ND_R1D_FECollection(const int p, const int dim, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Definition: fe_coll.cpp:2971
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:955
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:3520
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
Piecewise-quadratic discontinuous finite elements in 3D. This class is kept only for backward compati...
Definition: fe_coll.hpp:1171
Arbitrary order L2 elements in 3D utilizing the Bernstein basis on a wedge.
Definition: fe_pos.hpp:351
Arbitrary order H1 elements in 3D on a cube.
Definition: fe_h1.hpp:62
FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:3370
Serendipity basis (squares / cubes)
Definition: fe_base.hpp:37
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
virtual int GetDerivMapType(int dim) const
Definition: fe_coll.cpp:80
MFEM_EXPORT class Linear3DFiniteElement TetrahedronFE
Definition: fe.cpp:36
Arbitrary order Raviart-Thomas elements in 2D on a triangle.
Definition: fe_rt.hpp:163
MFEM_EXPORT Linear2DFiniteElement TriangleFE
Definition: fe.cpp:32
static int CheckOpen(int type)
If the Quadrature1D type is not open return Invalid; otherwise return type.
Definition: intrules.cpp:923
void SetMapType(const int map_type_) override
Set the FiniteElement::MapType of the element to either VALUE or INTEGRAL. Also sets the FiniteElemen...
Definition: fe_base.cpp:2511
Integrated GLL indicator functions.
Definition: fe_base.hpp:39
Version of LinearDiscont2DFECollection with dofs in the Gaussian points.
Definition: fe_coll.hpp:988
int GetDerivMapType() const
Returns the FiniteElement::DerivType of the element describing how reference function derivatives are...
Definition: fe_base.hpp:359
ND_R2D_FECollection(const int p, const int dim, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Definition: fe_coll.cpp:3110
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name...
Definition: fe_coll.cpp:116
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1033
static char GetChar(int b_type)
Check and convert a BasisType constant to a char basis identifier.
Definition: fe_base.hpp:104
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:831
Arbitrary order H1-conforming (continuous) finite elements with positive basis functions.
Definition: fe_coll.hpp:299
Lowest order Nedelec finite elements in 3D. This class is kept only for backward compatibility, consider using the new ND_FECollection instead.
Definition: fe_coll.hpp:1222
FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:3203
Arbitrary order Nedelec elements in 2D on a square.
Definition: fe_nd.hpp:104
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1589
static bool IsTensorProduct(Type geom)
Definition: geom.hpp:108
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1153
FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:2920
int HasFaceDofs(Geometry::Type geom, int p) const
Definition: fe_coll.cpp:90
virtual StatelessDofTransformation * GetDofTransformation() const
Return a DoF transformation object for this particular type of basis.
Definition: fe_base.hpp:598
Arbitrary order H1 serendipity elements in 2D on a quad.
Definition: fe_ser.hpp:21
Arbitrary order Raviart-Thomas elements in 2D on a square.
Definition: fe_rt.hpp:23
static const int DimStart[MaxDim+2]
Definition: geom.hpp:48
virtual int GetMapType(int dim) const
Definition: fe_coll.cpp:60
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:1131
MFEM_EXPORT class LinearWedgeFiniteElement WedgeFE
Definition: fe.cpp:40
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:751
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1364
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:1061
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:714
Arbitrary order H1 elements in 2D on a triangle.
Definition: fe_h1.hpp:83
int dim
Definition: ex24.cpp:53
void InitFaces(const int p, const int dim, const int map_type, const bool signs)
Definition: fe_coll.cpp:3321
Piecewise-linear discontinuous finite elements in 3D. This class is kept only for backward compatibil...
Definition: fe_coll.hpp:1145
RT_R2D_Trace_FECollection(const int p, const int dim, const int map_type=FiniteElement::INTEGRAL, const int ob_type=BasisType::GaussLegendre)
Definition: fe_coll.cpp:3396
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:724
Arbitrary order L2 elements in 2D on a square.
Definition: fe_l2.hpp:44
int ND_dof[Geometry::NumGeom]
Definition: fe_coll.hpp:566
Arbitrary order H1 elements in 1D.
Definition: fe_h1.hpp:21
int GetDerivRangeType() const
Returns the FiniteElement::RangeType of the element derivative, either SCALAR or VECTOR.
Definition: fe_base.hpp:344
int RT_dof[Geometry::NumGeom]
Definition: fe_coll.hpp:606
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:1398
RT_Trace_FECollection(const int p, const int dim, const int map_type=FiniteElement::INTEGRAL, const int ob_type=BasisType::GaussLegendre)
Definition: fe_coll.cpp:2650
Arbitrary order H1 elements in 1D utilizing the Bernstein basis.
Definition: fe_pos.hpp:119
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition: fe_base.hpp:349
Arbitrary order Nedelec 3D elements in 2D on a triangle.
Definition: fe_nd.hpp:615
Piecewise-(bi)quadratic continuous finite elements.
Definition: fe_coll.hpp:735
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1444
First order Raviart-Thomas finite elements in 2D. This class is kept only for backward compatibility...
Definition: fe_coll.hpp:868
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:1966
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
Arbitrary order Nedelec elements in 3D on a cube.
Definition: fe_nd.hpp:22
Arbitrary order H1 elements in 2D utilizing the Bernstein basis on a square.
Definition: fe_pos.hpp:142
A 3D constant element on a pyramid.
Arbitrary order L2 elements in 3D utilizing the Bernstein basis on a cube.
Definition: fe_pos.hpp:296
virtual ~RT_R1D_FECollection()
Definition: fe_coll.cpp:3101
A linear element defined on a square pyramid.
const int * GetDofMap(Geometry::Type GeomType) const
Get the Cartesian to local H1 dof map.
Definition: fe_coll.cpp:2007
virtual int GetDerivRangeType(int dim) const
Definition: fe_coll.cpp:50
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:846
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:454
Arbitrary order "H^{-1/2}-conforming" face finite elements defined on the interface between mesh elem...
Definition: fe_coll.hpp:434
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:3507
Arbitrary order 3D H(div)-conforming Raviart-Thomas finite elements in 2D.
Definition: fe_coll.hpp:600
virtual ~ND_R1D_FECollection()
Definition: fe_coll.cpp:3031
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
Definition: fe_base.hpp:340
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:2347
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:671
Second order Raviart-Thomas finite elements in 3D. This class is kept only for backward compatibility...
Definition: fe_coll.hpp:1273
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1011
A 3D 0th order Raviert-Thomas element on a pyramid.
FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:3026
Arbitrary order H1 elements in 3D on a wedge.
Definition: fe_h1.hpp:130
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1047
Piecewise-quadratic discontinuous finite elements in 2D. This class is kept only for backward compati...
Definition: fe_coll.hpp:1030
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1245
Arbitrary order 3D "H^{-1/2}-conforming" face finite elements defined on the interface between mesh e...
Definition: fe_coll.hpp:640
Local_FECollection(const char *fe_name)
Definition: fe_coll.cpp:3416
FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:2621
Arbitrary order Raviart-Thomas elements in 3D on a cube.
Definition: fe_rt.hpp:94
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1281
void SetOrder(int Order) const
Set the order and the name, based on the given Order: either a positive number for fixed order...
Definition: fe_coll.cpp:3470
Linear (P1) finite elements on quadrilaterals.
Definition: fe_coll.hpp:1012
void InitVarOrder(int p) const
Definition: fe_coll.cpp:369
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:761
An arbitrary order 3D NURBS element on a cube.
Definition: fe_nurbs.hpp:119
An arbitrary order 1D NURBS element on a segment.
Definition: fe_nurbs.hpp:67
int ND_dof[Geometry::NumGeom]
Definition: fe_coll.hpp:462
FiniteElement * RT_Elements[Geometry::NumGeom]
Definition: fe_coll.hpp:538
An arbitrary order 2D NURBS element on a square.
Definition: fe_nurbs.hpp:87
Arbitrary order H1 elements in 2D utilizing the Bernstein basis on a triangle.
Definition: fe_pos.hpp:180
virtual StatelessDofTransformation * DofTransformationForGeometry(Geometry::Type GeomType) const
Returns a DoF transformation object compatible with this basis and geometry type. ...
Definition: fe_coll.cpp:2890
A 3D 1st order Nedelec element on a cube.
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:1571
virtual ~ND_R2D_FECollection()
Definition: fe_coll.cpp:3223
No derivatives implemented.
Definition: fe_base.hpp:294
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
Definition: fe_coll.cpp:1168
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:327
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:998
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1259