MFEM  v4.1.0 Finite element discretization library
fe_coll.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
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
26 {
27  switch (GeomType)
28  {
29  case Geometry::TETRAHEDRON: return DofForGeometry (Geometry::TRIANGLE);
30  case Geometry::CUBE: return DofForGeometry (Geometry::SQUARE);
31  case Geometry::PRISM:
32  return max(DofForGeometry (Geometry::TRIANGLE),
33  DofForGeometry (Geometry::SQUARE));
34  default:
35  mfem_error ("FiniteElementCollection::HasFaceDofs:"
36  " unknown geometry type.");
37  }
38  return 0;
39 }
40
42 {
43  MFEM_ABORT("this method is not implemented in this derived class!");
44  return NULL;
45 }
46
48 {
49  FiniteElementCollection *fec = NULL;
50
51  if (!strcmp(name, "Linear"))
52  {
53  fec = new LinearFECollection;
54  }
56  {
58  }
60  {
62  }
63  else if (!strcmp(name, "Cubic"))
64  {
65  fec = new CubicFECollection;
66  }
67  else if (!strcmp(name, "Const3D"))
68  {
69  fec = new Const3DFECollection;
70  }
71  else if (!strcmp(name, "Const2D"))
72  {
73  fec = new Const2DFECollection;
74  }
75  else if (!strcmp(name, "LinearDiscont2D"))
76  {
78  }
79  else if (!strcmp(name, "GaussLinearDiscont2D"))
80  {
82  }
84  {
86  }
88  {
90  }
92  {
94  }
96  {
98  }
99  else if (!strcmp(name, "CubicDiscont2D"))
100  {
101  fec = new CubicDiscont2DFECollection;
102  }
103  else if (!strcmp(name, "LinearDiscont3D"))
104  {
105  fec = new LinearDiscont3DFECollection;
106  }
108  {
110  }
111  else if (!strcmp(name, "LinearNonConf3D"))
112  {
113  fec = new LinearNonConf3DFECollection;
114  }
115  else if (!strcmp(name, "CrouzeixRaviart"))
116  {
117  fec = new CrouzeixRaviartFECollection;
118  }
119  else if (!strcmp(name, "ND1_3D"))
120  {
121  fec = new ND1_3DFECollection;
122  }
123  else if (!strcmp(name, "RT0_2D"))
124  {
125  fec = new RT0_2DFECollection;
126  }
127  else if (!strcmp(name, "RT1_2D"))
128  {
129  fec = new RT1_2DFECollection;
130  }
131  else if (!strcmp(name, "RT2_2D"))
132  {
133  fec = new RT2_2DFECollection;
134  }
135  else if (!strcmp(name, "RT0_3D"))
136  {
137  fec = new RT0_3DFECollection;
138  }
139  else if (!strcmp(name, "RT1_3D"))
140  {
141  fec = new RT1_3DFECollection;
142  }
143  else if (!strncmp(name, "H1_Trace_", 9))
144  {
145  fec = new H1_Trace_FECollection(atoi(name + 13), atoi(name + 9));
146  }
147  else if (!strncmp(name, "H1_Trace@", 9))
148  {
149  fec = new H1_Trace_FECollection(atoi(name + 15), atoi(name + 11),
150  BasisType::GetType(name[9]));
151  }
152  else if (!strncmp(name, "H1_", 3))
153  {
154  fec = new H1_FECollection(atoi(name + 7), atoi(name + 3));
155  }
156  else if (!strncmp(name, "H1Pos_Trace_", 12))
157  {
158  fec = new H1_Trace_FECollection(atoi(name + 16), atoi(name + 12),
160  }
161  else if (!strncmp(name, "H1Pos_", 6))
162  {
163  fec = new H1Pos_FECollection(atoi(name + 10), atoi(name + 6));
164  }
165  else if (!strncmp(name, "H1Ser_", 6))
166  {
167  fec = new H1Ser_FECollection(atoi(name + 10), atoi(name + 6));
168  }
169  else if (!strncmp(name, "H1@", 3))
170  {
171  fec = new H1_FECollection(atoi(name + 9), atoi(name + 5),
172  BasisType::GetType(name[3]));
173  }
174  else if (!strncmp(name, "L2_T", 4))
175  fec = new L2_FECollection(atoi(name + 10), atoi(name + 6),
176  atoi(name + 4));
177  else if (!strncmp(name, "L2_", 3))
178  {
179  fec = new L2_FECollection(atoi(name + 7), atoi(name + 3));
180  }
181  else if (!strncmp(name, "L2Int_T", 7))
182  {
183  fec = new L2_FECollection(atoi(name + 13), atoi(name + 9),
184  atoi(name + 7), FiniteElement::INTEGRAL);
185  }
186  else if (!strncmp(name, "L2Int_", 6))
187  {
188  fec = new L2_FECollection(atoi(name + 10), atoi(name + 6),
191  }
192  else if (!strncmp(name, "RT_Trace_", 9))
193  {
194  fec = new RT_Trace_FECollection(atoi(name + 13), atoi(name + 9));
195  }
196  else if (!strncmp(name, "RT_ValTrace_", 12))
197  {
198  fec = new RT_Trace_FECollection(atoi(name + 16), atoi(name + 12),
200  }
201  else if (!strncmp(name, "RT_Trace@", 9))
202  {
203  fec = new RT_Trace_FECollection(atoi(name + 15), atoi(name + 11),
205  BasisType::GetType(name[9]));
206  }
207  else if (!strncmp(name, "RT_ValTrace@", 12))
208  {
209  fec = new RT_Trace_FECollection(atoi(name + 18), atoi(name + 14),
211  BasisType::GetType(name[12]));
212  }
213  else if (!strncmp(name, "DG_Iface_", 9))
214  {
215  fec = new DG_Interface_FECollection(atoi(name + 13), atoi(name + 9));
216  }
217  else if (!strncmp(name, "DG_Iface@", 9))
218  {
219  fec = new DG_Interface_FECollection(atoi(name + 15), atoi(name + 11),
221  BasisType::GetType(name[9]));
222  }
223  else if (!strncmp(name, "DG_IntIface_", 12))
224  {
225  fec = new DG_Interface_FECollection(atoi(name + 16), atoi(name + 12),
227  }
228  else if (!strncmp(name, "DG_IntIface@", 12))
229  {
230  fec = new DG_Interface_FECollection(atoi(name + 18), atoi(name + 14),
232  BasisType::GetType(name[12]));
233  }
234  else if (!strncmp(name, "RT_", 3))
235  {
236  fec = new RT_FECollection(atoi(name + 7), atoi(name + 3));
237  }
238  else if (!strncmp(name, "RT@", 3))
239  {
240  fec = new RT_FECollection(atoi(name + 10), atoi(name + 6),
241  BasisType::GetType(name[3]),
242  BasisType::GetType(name[4]));
243  }
244  else if (!strncmp(name, "ND_Trace_", 9))
245  {
246  fec = new ND_Trace_FECollection(atoi(name + 13), atoi(name + 9));
247  }
248  else if (!strncmp(name, "ND_Trace@", 9))
249  {
250  fec = new ND_Trace_FECollection(atoi(name + 16), atoi(name + 12),
251  BasisType::GetType(name[9]),
252  BasisType::GetType(name[10]));
253  }
254  else if (!strncmp(name, "ND_", 3))
255  {
256  fec = new ND_FECollection(atoi(name + 7), atoi(name + 3));
257  }
258  else if (!strncmp(name, "ND@", 3))
259  {
260  fec = new ND_FECollection(atoi(name + 10), atoi(name + 6),
261  BasisType::GetType(name[3]),
262  BasisType::GetType(name[4]));
263  }
264  else if (!strncmp(name, "Local_", 6))
265  {
266  fec = new Local_FECollection(name + 6);
267  }
268  else if (!strncmp(name, "NURBS", 5))
269  {
270  if (name[5] != '\0')
271  {
272  // "NURBS" + "number" --> fixed order nurbs collection
273  fec = new NURBSFECollection(atoi(name + 5));
274  }
275  else
276  {
277  // "NURBS" --> variable order nurbs collection
278  fec = new NURBSFECollection();
279  }
280  }
281  else
282  {
283  MFEM_ABORT("unknown FiniteElementCollection: " << name);
284  }
285  MFEM_VERIFY(!strcmp(fec->Name(), name), "input name: \"" << name
286  << "\" does not match the created collection name: \""
287  << fec->Name() << '"');
288
289  return fec;
290 }
291
292 template <Geometry::Type geom>
293 inline void FiniteElementCollection::GetNVE(int &nv, int &ne)
294 {
295  typedef typename Geometry::Constants<geom> g_consts;
296
297  nv = g_consts::NumVert;
298  ne = g_consts::NumEdges;
299 }
300
301 template <Geometry::Type geom, typename v_t>
302 inline void FiniteElementCollection::
303 GetEdge(int &nv, v_t &v, int &ne, int &e, int &eo, const int edge_info)
304 {
305  typedef typename Geometry::Constants<Geometry::SEGMENT> e_consts;
306  typedef typename Geometry::Constants<geom> g_consts;
307
308  nv = e_consts::NumVert;
309  ne = 1;
310  e = edge_info/64;
311  eo = edge_info%64;
312  MFEM_ASSERT(0 <= e && e < g_consts::NumEdges, "");
313  MFEM_ASSERT(0 <= eo && eo < e_consts::NumOrient, "");
314  v[0] = g_consts::Edges[e][0];
315  v[1] = g_consts::Edges[e][1];
316  v[0] = e_consts::Orient[eo][v[0]];
317  v[1] = e_consts::Orient[eo][v[1]];
318 }
319
320 template <Geometry::Type geom, Geometry::Type f_geom,
321  typename v_t, typename e_t, typename eo_t>
322 inline void FiniteElementCollection::
323 GetFace(int &nv, v_t &v, int &ne, e_t &e, eo_t &eo,
324  int &nf, int &f, Geometry::Type &fg, int &fo, const int face_info)
325 {
326  typedef typename Geometry::Constants< geom> g_consts;
327  typedef typename Geometry::Constants<f_geom> f_consts;
328
329  nv = f_consts::NumVert;
330  nf = 1;
331  f = face_info/64;
332  fg = f_geom;
333  fo = face_info%64;
334  MFEM_ASSERT(0 <= f && f < g_consts::NumFaces, "");
335  MFEM_ASSERT(0 <= fo && fo < f_consts::NumOrient, "");
336  for (int i = 0; i < f_consts::NumVert; i++)
337  {
338  v[i] = f_consts::Orient[fo][i];
339  v[i] = g_consts::FaceVert[f][v[i]];
340  }
341  ne = f_consts::NumEdges;
342  for (int i = 0; i < f_consts::NumEdges; i++)
343  {
344  int v0 = v[f_consts::Edges[i][0]];
345  int v1 = v[f_consts::Edges[i][1]];
346  int eor = 0;
347  if (v0 > v1) { swap(v0, v1); eor = 1; }
348  for (int j = g_consts::VertToVert::I[v0]; true; j++)
349  {
350  MFEM_ASSERT(j < g_consts::VertToVert::I[v0+1],
352  if (v1 == g_consts::VertToVert::J[j][0])
353  {
354  int en = g_consts::VertToVert::J[j][1];
355  if (en < 0)
356  {
357  en = -1-en;
358  eor = 1-eor;
359  }
360  e[i] = en;
361  eo[i] = eor;
362  break;
363  }
364  }
365  }
366 }
367
369  int Info,
370  Array<int> &dofs) const
371 {
372  // Info = 64 * SubIndex + SubOrientation
373  MFEM_ASSERT(0 <= Geom && Geom < Geometry::NumGeom,
374  "invalid Geom = " << Geom);
375  MFEM_ASSERT(0 <= SDim && SDim <= Geometry::Dimension[Geom],
376  "invalid SDim = " << SDim <<
377  " for Geom = " << Geometry::Name[Geom]);
378
379  const int nvd = DofForGeometry(Geometry::POINT);
380  if (SDim == 0) // vertex
381  {
382  const int off = nvd*(Info/64);
383  dofs.SetSize(nvd);
384  for (int i = 0; i < nvd; i++)
385  {
386  dofs[i] = off + i;
387  }
388  }
389  else
390  {
391  int v[4], e[4], eo[4], f[1], fo[1];
392  int av = 0, nv = 0, ae = 0, ne = 0, nf = 0;
393  Geometry::Type fg[1];
394
395  switch (Geom)
396  {
397  case Geometry::SEGMENT:
398  {
399  GetNVE<Geometry::SEGMENT>(av, ae);
400  GetEdge<Geometry::SEGMENT>(nv, v, ne, e[0], eo[0], Info);
401  break;
402  }
403
404  case Geometry::TRIANGLE:
405  {
406  GetNVE<Geometry::TRIANGLE>(av, ae);
407  switch (SDim)
408  {
409  case 1:
410  GetEdge<Geometry::TRIANGLE>(nv, v, ne, e[0], eo[0], Info);
411  break;
412  case 2:
413  GetFace<Geometry::TRIANGLE,Geometry::TRIANGLE>(
414  nv, v, ne, e, eo, nf, f[0], fg[0], fo[0], Info);
415  break;
416  default:
417  goto not_supp;
418  }
419  break;
420  }
421
422  case Geometry::SQUARE:
423  {
424  GetNVE<Geometry::SQUARE>(av, ae);
425  switch (SDim)
426  {
427  case 1:
428  GetEdge<Geometry::SQUARE>(nv, v, ne, e[0], eo[0], Info);
429  break;
430  case 2:
431  GetFace<Geometry::SQUARE,Geometry::SQUARE>(
432  nv, v, ne, e, eo, nf, f[0], fg[0], fo[0], Info);
433  break;
434  default:
435  goto not_supp;
436  }
437  break;
438  }
439
441  {
442  GetNVE<Geometry::TETRAHEDRON>(av, ae);
443  switch (SDim)
444  {
445  case 1:
446  GetEdge<Geometry::TETRAHEDRON>(nv, v, ne, e[0], eo[0], Info);
447  break;
448  case 2:
449  GetFace<Geometry::TETRAHEDRON,Geometry::TRIANGLE>(
450  nv, v, ne, e, eo, nf, f[0], fg[0], fo[0], Info);
451  break;
452  default:
453  goto not_supp;
454  }
455  break;
456  }
457
458  case Geometry::CUBE:
459  {
460  GetNVE<Geometry::CUBE>(av, ae);
461  switch (SDim)
462  {
463  case 1:
464  GetEdge<Geometry::CUBE>(nv, v, ne, e[0], eo[0], Info);
465  break;
466  case 2:
467  GetFace<Geometry::CUBE,Geometry::SQUARE>(
468  nv, v, ne, e, eo, nf, f[0], fg[0], fo[0], Info);
469  break;
470  default:
471  goto not_supp;
472  }
473  break;
474  }
475
476  default:
477  MFEM_ABORT("invalid Geom = " << Geom);
478  }
479
480  int ned = (ne > 0) ? DofForGeometry(Geometry::SEGMENT) : 0;
481
483  dofs.SetSize(nv*nvd+ne*ned);
484  for (int i = 0; i < nv; i++)
485  {
486  for (int j = 0; j < nvd; j++)
487  {
488  dofs[i*nvd+j] = v[i]*nvd+j;
489  }
490  }
491  int l_off = nv*nvd, g_off = av*nvd;
492
494  if (ned > 0)
495  {
496  for (int i = 0; i < ne; i++)
497  {
498  const int *ed = DofOrderForOrientation(Geometry::SEGMENT,
499  eo[i] ? -1 : 1);
500  for (int j = 0; j < ned; j++)
501  {
502  dofs[l_off+i*ned+j] =
503  ed[j] >= 0 ?
504  g_off+e[i]*ned+ed[j] :
505  -1-(g_off+e[i]*ned+(-1-ed[j]));
506  }
507  }
508  l_off += ne*ned;
509  g_off += ae*ned;
510  }
511
513  if (nf > 0)
514  {
515  const int nfd = DofForGeometry(fg[0]); // assume same face geometry
516  dofs.SetSize(dofs.Size()+nf*nfd);
517  for (int i = 0; i < nf; i++)
518  {
519  const int *fd = DofOrderForOrientation(fg[i], fo[i]);
520  for (int j = 0; j < nfd; j++)
521  {
522  dofs[l_off+i*nfd+j] =
523  fd[j] >= 0 ?
524  g_off+f[i]*nfd+fd[j] :
525  -1-(g_off+f[i]*nfd+(-1-fd[j]));
526  }
527  }
528  }
529
530  // add volume dofs ...
531  }
532  return;
533
534 not_supp:
535  MFEM_ABORT("Geom = " << Geometry::Name[Geom] <<
536  ", SDim = " << SDim << " is not supported");
537 }
538
539 const FiniteElement *
541 {
542  switch (GeomType)
543  {
544  case Geometry::POINT: return &PointFE;
545  case Geometry::SEGMENT: return &SegmentFE;
546  case Geometry::TRIANGLE: return &TriangleFE;
548  case Geometry::TETRAHEDRON: return &TetrahedronFE;
549  case Geometry::CUBE: return &ParallelepipedFE;
550  case Geometry::PRISM: return &WedgeFE;
551  default:
552  mfem_error ("LinearFECollection: unknown geometry type.");
553  }
554  return &SegmentFE; // Make some compilers happy
555 }
556
558 {
559  switch (GeomType)
560  {
561  case Geometry::POINT: return 1;
562  case Geometry::SEGMENT: return 0;
563  case Geometry::TRIANGLE: return 0;
564  case Geometry::SQUARE: return 0;
565  case Geometry::TETRAHEDRON: return 0;
566  case Geometry::CUBE: return 0;
567  case Geometry::PRISM: return 0;
568  default:
569  mfem_error ("LinearFECollection: unknown geometry type.");
570  }
571  return 0; // Make some compilers happy
572 }
573
575  int Or) const
576 {
577  return NULL;
578 }
579
580
581 const FiniteElement *
583 {
584  switch (GeomType)
585  {
586  case Geometry::POINT: return &PointFE;
587  case Geometry::SEGMENT: return &SegmentFE;
588  case Geometry::TRIANGLE: return &TriangleFE;
590  case Geometry::TETRAHEDRON: return &TetrahedronFE;
591  case Geometry::CUBE: return &ParallelepipedFE;
592  case Geometry::PRISM: return &WedgeFE;
593  default:
594  mfem_error ("QuadraticFECollection: unknown geometry type.");
595  }
596  return &SegmentFE; // Make some compilers happy
597 }
598
600 {
601  switch (GeomType)
602  {
603  case Geometry::POINT: return 1;
604  case Geometry::SEGMENT: return 1;
605  case Geometry::TRIANGLE: return 0;
606  case Geometry::SQUARE: return 1;
607  case Geometry::TETRAHEDRON: return 0;
608  case Geometry::CUBE: return 1;
609  case Geometry::PRISM: return 0;
610  default:
611  mfem_error ("QuadraticFECollection: unknown geometry type.");
612  }
613  return 0; // Make some compilers happy
614 }
615
617  Geometry::Type GeomType, int Or) const
618 {
619  static int indexes[] = { 0 };
620
621  return indexes;
622 }
623
624
625 const FiniteElement *
627  Geometry::Type GeomType) const
628 {
629  switch (GeomType)
630  {
631  case Geometry::SEGMENT: return &SegmentFE;
633  default:
634  mfem_error ("QuadraticPosFECollection: unknown geometry type.");
635  }
636  return NULL; // Make some compilers happy
637 }
638
640 {
641  switch (GeomType)
642  {
643  case Geometry::POINT: return 1;
644  case Geometry::SEGMENT: return 1;
645  case Geometry::SQUARE: return 1;
646  default:
647  mfem_error ("QuadraticPosFECollection: unknown geometry type.");
648  }
649  return 0; // Make some compilers happy
650 }
651
653  Geometry::Type GeomType, int Or) const
654 {
655  static int indexes[] = { 0 };
656
657  return indexes;
658 }
659
660
661 const FiniteElement *
663 {
664  switch (GeomType)
665  {
666  case Geometry::POINT: return &PointFE;
667  case Geometry::SEGMENT: return &SegmentFE;
668  case Geometry::TRIANGLE: return &TriangleFE;
670  case Geometry::TETRAHEDRON: return &TetrahedronFE;
671  case Geometry::CUBE: return &ParallelepipedFE;
672  case Geometry::PRISM: return &WedgeFE;
673  default:
674  mfem_error ("CubicFECollection: unknown geometry type.");
675  }
676  return &SegmentFE; // Make some compilers happy
677 }
678
680 {
681  switch (GeomType)
682  {
683  case Geometry::POINT: return 1;
684  case Geometry::SEGMENT: return 2;
685  case Geometry::TRIANGLE: return 1;
686  case Geometry::SQUARE: return 4;
687  case Geometry::TETRAHEDRON: return 0;
688  case Geometry::CUBE: return 8;
689  case Geometry::PRISM: return 2;
690  default:
691  mfem_error ("CubicFECollection: unknown geometry type.");
692  }
693  return 0; // Make some compilers happy
694 }
695
697  int Or) const
698 {
699  if (GeomType == Geometry::SEGMENT)
700  {
701  static int ind_pos[] = { 0, 1 };
702  static int ind_neg[] = { 1, 0 };
703
704  if (Or < 0)
705  {
706  return ind_neg;
707  }
708  return ind_pos;
709  }
710  else if (GeomType == Geometry::TRIANGLE)
711  {
712  static int indexes[] = { 0 };
713
714  return indexes;
715  }
716  else if (GeomType == Geometry::SQUARE)
717  {
718  static int sq_ind[8][4] = {{0, 1, 2, 3}, {0, 2, 1, 3},
719  {2, 0, 3, 1}, {1, 0, 3, 2},
720  {3, 2, 1, 0}, {3, 1, 2, 0},
721  {1, 3, 0, 2}, {2, 3, 0, 1}
722  };
723  return sq_ind[Or];
724  }
725
726  return NULL;
727 }
728
729
730 const FiniteElement *
732  Geometry::Type GeomType) const
733 {
734  switch (GeomType)
735  {
736  case Geometry::SEGMENT: return &SegmentFE;
737  case Geometry::TRIANGLE: return &TriangleFE;
739  default:
740  mfem_error ("CrouzeixRaviartFECollection: unknown geometry type.");
741  }
742  return &SegmentFE; // Make some compilers happy
743 }
744
746 {
747  switch (GeomType)
748  {
749  case Geometry::POINT: return 0;
750  case Geometry::SEGMENT: return 1;
751  case Geometry::TRIANGLE: return 0;
752  case Geometry::SQUARE: return 0;
753  default:
754  mfem_error ("CrouzeixRaviartFECollection: unknown geometry type.");
755  }
756  return 0; // Make some compilers happy
757 }
758
760  Geometry::Type GeomType, int Or) const
761 {
762  static int indexes[] = { 0 };
763
764  return indexes;
765 }
766
767
768 const FiniteElement *
770 {
771  switch (GeomType)
772  {
773  case Geometry::SEGMENT: return &SegmentFE;
774  case Geometry::TRIANGLE: return &TriangleFE;
776  default:
777  mfem_error ("RT0_2DFECollection: unknown geometry type.");
778  }
779  return &SegmentFE; // Make some compilers happy
780 }
781
783 {
784  switch (GeomType)
785  {
786  case Geometry::POINT: return 0;
787  case Geometry::SEGMENT: return 1;
788  case Geometry::TRIANGLE: return 0;
789  case Geometry::SQUARE: return 0;
790  default:
791  mfem_error ("RT0_2DFECollection: unknown geometry type.");
792  }
793  return 0; // Make some compilers happy
794 }
795
797  int Or) const
798 {
799  static int ind_pos[] = { 0 };
800  static int ind_neg[] = { -1 };
801
802  if (Or > 0)
803  {
804  return ind_pos;
805  }
806  return ind_neg;
807 }
808
809
810 const FiniteElement *
812 {
813  switch (GeomType)
814  {
815  case Geometry::SEGMENT: return &SegmentFE;
816  case Geometry::TRIANGLE: return &TriangleFE;
818  default:
819  mfem_error ("RT1_2DFECollection: unknown geometry type.");
820  }
821  return &SegmentFE; // Make some compilers happy
822 }
823
825 {
826  switch (GeomType)
827  {
828  case Geometry::POINT: return 0;
829  case Geometry::SEGMENT: return 2;
830  case Geometry::TRIANGLE: return 2;
831  case Geometry::SQUARE: return 4;
832  default:
833  mfem_error ("RT1_2DFECollection: unknown geometry type.");
834  }
835  return 0; // Make some compilers happy
836 }
837
839  int Or) const
840 {
841  static int ind_pos[] = { 0, 1 };
842  static int ind_neg[] = { -2, -1 };
843
844  if (Or > 0)
845  {
846  return ind_pos;
847  }
848  return ind_neg;
849 }
850
851 const FiniteElement *
853 {
854  switch (GeomType)
855  {
856  case Geometry::SEGMENT: return &SegmentFE;
857  case Geometry::TRIANGLE: return &TriangleFE;
859  default:
860  mfem_error ("RT2_2DFECollection: unknown geometry type.");
861  }
862  return &SegmentFE; // Make some compilers happy
863 }
864
866 {
867  switch (GeomType)
868  {
869  case Geometry::POINT: return 0;
870  case Geometry::SEGMENT: return 3;
871  case Geometry::TRIANGLE: return 6;
872  case Geometry::SQUARE: return 12;
873  default:
874  mfem_error ("RT2_2DFECollection: unknown geometry type.");
875  }
876  return 0; // Make some compilers happy
877 }
878
880  int Or) const
881 {
882  static int ind_pos[] = { 0, 1, 2 };
883  static int ind_neg[] = { -3, -2, -1 };
884
885  if (Or > 0)
886  {
887  return ind_pos;
888  }
889  return ind_neg;
890 }
891
892
893 const FiniteElement *
895 {
896  switch (GeomType)
897  {
898  case Geometry::TRIANGLE: return &TriangleFE;
900  default:
901  mfem_error ("Const2DFECollection: unknown geometry type.");
902  }
903  return &TriangleFE; // Make some compilers happy
904 }
905
907 {
908  switch (GeomType)
909  {
910  case Geometry::POINT: return 0;
911  case Geometry::SEGMENT: return 0;
912  case Geometry::TRIANGLE: return 1;
913  case Geometry::SQUARE: return 1;
914  default:
915  mfem_error ("Const2DFECollection: unknown geometry type.");
916  }
917  return 0; // Make some compilers happy
918 }
919
921  int Or) const
922 {
923  return NULL;
924 }
925
926
927 const FiniteElement *
929  Geometry::Type GeomType) const
930 {
931  switch (GeomType)
932  {
933  case Geometry::TRIANGLE: return &TriangleFE;
935  default:
936  mfem_error ("LinearDiscont2DFECollection: unknown geometry type.");
937  }
938  return &TriangleFE; // Make some compilers happy
939 }
940
942 {
943  switch (GeomType)
944  {
945  case Geometry::POINT: return 0;
946  case Geometry::SEGMENT: return 0;
947  case Geometry::TRIANGLE: return 3;
948  case Geometry::SQUARE: return 4;
949  default:
950  mfem_error ("LinearDiscont2DFECollection: unknown geometry type.");
951  }
952  return 0; // Make some compilers happy
953 }
954
956  Geometry::Type GeomType, int Or) const
957 {
958  return NULL;
959 }
960
961
962 const FiniteElement *
964  Geometry::Type GeomType) const
965 {
966  switch (GeomType)
967  {
968  case Geometry::TRIANGLE: return &TriangleFE;
970  default:
971  mfem_error ("GaussLinearDiscont2DFECollection:"
972  " unknown geometry type.");
973  }
974  return &TriangleFE; // Make some compilers happy
975 }
976
978 const
979 {
980  switch (GeomType)
981  {
982  case Geometry::POINT: return 0;
983  case Geometry::SEGMENT: return 0;
984  case Geometry::TRIANGLE: return 3;
985  case Geometry::SQUARE: return 4;
986  default:
987  mfem_error ("GaussLinearDiscont2DFECollection:"
988  " unknown geometry type.");
989  }
990  return 0; // Make some compilers happy
991 }
992
994  Geometry::Type GeomType, int Or) const
995 {
996  return NULL;
997 }
998
999
1000 const FiniteElement *
1002 {
1003  if (GeomType != Geometry::SQUARE)
1004  {
1005  mfem_error ("P1OnQuadFECollection: unknown geometry type.");
1006  }
1008 }
1009
1011 {
1012  switch (GeomType)
1013  {
1014  case Geometry::POINT: return 0;
1015  case Geometry::SEGMENT: return 0;
1016  case Geometry::SQUARE: return 3;
1017  default:
1018  mfem_error ("P1OnQuadFECollection: unknown geometry type.");
1019  }
1020  return 0; // Make some compilers happy
1021 }
1022
1024  Geometry::Type GeomType, int Or) const
1025 {
1026  return NULL;
1027 }
1028
1029
1030 const FiniteElement *
1032  Geometry::Type GeomType) const
1033 {
1034  switch (GeomType)
1035  {
1036  case Geometry::TRIANGLE: return &TriangleFE;
1038  default:
1039  mfem_error ("QuadraticDiscont2DFECollection: unknown geometry type.");
1040  }
1041  return &TriangleFE; // Make some compilers happy
1042 }
1043
1045 const
1046 {
1047  switch (GeomType)
1048  {
1049  case Geometry::POINT: return 0;
1050  case Geometry::SEGMENT: return 0;
1051  case Geometry::TRIANGLE: return 6;
1052  case Geometry::SQUARE: return 9;
1053  default:
1054  mfem_error ("QuadraticDiscont2DFECollection: unknown geometry type.");
1055  }
1056  return 0; // Make some compilers happy
1057 }
1058
1060  Geometry::Type GeomType, int Or) const
1061 {
1062  return NULL;
1063 }
1064
1065
1066 const FiniteElement *
1068  Geometry::Type GeomType) const
1069 {
1070  switch (GeomType)
1071  {
1073  default:
1074  mfem_error ("QuadraticPosDiscont2DFECollection: unknown geometry type.");
1075  }
1076  return NULL; // Make some compilers happy
1077 }
1078
1080 const
1081 {
1082  switch (GeomType)
1083  {
1084  case Geometry::POINT: return 0;
1085  case Geometry::SEGMENT: return 0;
1086  case Geometry::SQUARE: return 9;
1087  default:
1088  mfem_error ("QuadraticPosDiscont2DFECollection: unknown geometry type.");
1089  }
1090  return 0; // Make some compilers happy
1091 }
1092
1093
1094 const FiniteElement *
1096  Geometry::Type GeomType)
1097 const
1098 {
1099  switch (GeomType)
1100  {
1101  case Geometry::TRIANGLE: return &TriangleFE;
1103  default:
1105  " unknown geometry type.");
1106  }
1107  return &QuadrilateralFE; // Make some compilers happy
1108 }
1109
1111  Geometry::Type GeomType) const
1112 {
1113  switch (GeomType)
1114  {
1115  case Geometry::POINT: return 0;
1116  case Geometry::SEGMENT: return 0;
1117  case Geometry::TRIANGLE: return 6;
1118  case Geometry::SQUARE: return 9;
1119  default:
1121  " unknown geometry type.");
1122  }
1123  return 0; // Make some compilers happy
1124 }
1125
1127  Geometry::Type GeomType, int Or) const
1128 {
1129  return NULL;
1130 }
1131
1132
1133 const FiniteElement *
1135  Geometry::Type GeomType) const
1136 {
1137  switch (GeomType)
1138  {
1139  case Geometry::TRIANGLE: return &TriangleFE;
1141  default:
1142  mfem_error ("CubicDiscont2DFECollection: unknown geometry type.");
1143  }
1144  return &TriangleFE; // Make some compilers happy
1145 }
1146
1148 {
1149  switch (GeomType)
1150  {
1151  case Geometry::POINT: return 0;
1152  case Geometry::SEGMENT: return 0;
1153  case Geometry::TRIANGLE: return 10;
1154  case Geometry::SQUARE: return 16;
1155  default:
1156  mfem_error ("CubicDiscont2DFECollection: unknown geometry type.");
1157  }
1158  return 0; // Make some compilers happy
1159 }
1160
1162  Geometry::Type GeomType, int Or) const
1163 {
1164  return NULL;
1165 }
1166
1167
1168 const FiniteElement *
1170  Geometry::Type GeomType) const
1171 {
1172  switch (GeomType)
1173  {
1174  case Geometry::TRIANGLE: return &TriangleFE;
1176  case Geometry::TETRAHEDRON: return &TetrahedronFE;
1177  case Geometry::CUBE: return &ParallelepipedFE;
1178  default:
1179  mfem_error ("LinearNonConf3DFECollection: unknown geometry type.");
1180  }
1181  return &TriangleFE; // Make some compilers happy
1182 }
1183
1185 {
1186  switch (GeomType)
1187  {
1188  case Geometry::POINT: return 0;
1189  case Geometry::SEGMENT: return 0;
1190  case Geometry::TRIANGLE: return 1;
1191  case Geometry::SQUARE: return 1;
1192  case Geometry::TETRAHEDRON: return 0;
1193  case Geometry::CUBE: return 0;
1194  default:
1195  mfem_error ("LinearNonConf3DFECollection: unknown geometry type.");
1196  }
1197  return 0; // Make some compilers happy
1198 }
1199
1201  Geometry::Type GeomType, int Or) const
1202 {
1203  static int indexes[] = { 0 };
1204
1205  return indexes;
1206 }
1207
1208
1209 const FiniteElement *
1211 {
1212  switch (GeomType)
1213  {
1214  case Geometry::TETRAHEDRON: return &TetrahedronFE;
1215  case Geometry::CUBE: return &ParallelepipedFE;
1216  case Geometry::PRISM: return &WedgeFE;
1217  default:
1218  mfem_error ("Const3DFECollection: unknown geometry type.");
1219  }
1220  return &TetrahedronFE; // Make some compilers happy
1221 }
1222
1224 {
1225  switch (GeomType)
1226  {
1227  case Geometry::POINT: return 0;
1228  case Geometry::SEGMENT: return 0;
1229  case Geometry::TRIANGLE: return 0;
1230  case Geometry::SQUARE: return 0;
1231  case Geometry::TETRAHEDRON: return 1;
1232  case Geometry::CUBE: return 1;
1233  case Geometry::PRISM: return 1;
1234  default:
1235  mfem_error ("Const3DFECollection: unknown geometry type.");
1236  }
1237  return 0; // Make some compilers happy
1238 }
1239
1241  int Or) const
1242 {
1243  return NULL;
1244 }
1245
1246
1247 const FiniteElement *
1249  Geometry::Type GeomType) const
1250 {
1251  switch (GeomType)
1252  {
1253  case Geometry::TETRAHEDRON: return &TetrahedronFE;
1254  case Geometry::CUBE: return &ParallelepipedFE;
1255  default:
1256  mfem_error ("LinearDiscont3DFECollection: unknown geometry type.");
1257  }
1258  return &TetrahedronFE; // Make some compilers happy
1259 }
1260
1262 {
1263  switch (GeomType)
1264  {
1265  case Geometry::POINT: return 0;
1266  case Geometry::SEGMENT: return 0;
1267  case Geometry::TRIANGLE: return 0;
1268  case Geometry::SQUARE: return 0;
1269  case Geometry::TETRAHEDRON: return 4;
1270  case Geometry::CUBE: return 8;
1271  default:
1272  mfem_error ("LinearDiscont3DFECollection: unknown geometry type.");
1273  }
1274  return 0; // Make some compilers happy
1275 }
1276
1278  Geometry::Type GeomType, int Or) const
1279 {
1280  return NULL;
1281 }
1282
1283
1284 const FiniteElement *
1286  Geometry::Type GeomType) const
1287 {
1288  switch (GeomType)
1289  {
1290  case Geometry::TETRAHEDRON: return &TetrahedronFE;
1291  case Geometry::CUBE: return &ParallelepipedFE;
1292  default:
1293  mfem_error ("QuadraticDiscont3DFECollection: unknown geometry type.");
1294  }
1295  return &TetrahedronFE; // Make some compilers happy
1296 }
1297
1299 const
1300 {
1301  switch (GeomType)
1302  {
1303  case Geometry::POINT: return 0;
1304  case Geometry::SEGMENT: return 0;
1305  case Geometry::TRIANGLE: return 0;
1306  case Geometry::SQUARE: return 0;
1307  case Geometry::TETRAHEDRON: return 10;
1308  case Geometry::CUBE: return 27;
1309  default:
1310  mfem_error ("QuadraticDiscont3DFECollection: unknown geometry type.");
1311  }
1312  return 0; // Make some compilers happy
1313 }
1314
1316  Geometry::Type GeomType, int Or) const
1317 {
1318  return NULL;
1319 }
1320
1321 const FiniteElement *
1323  Geometry::Type GeomType) const
1324 {
1325  switch (GeomType)
1326  {
1327  case Geometry::POINT: return &PointFE;
1328  case Geometry::SEGMENT: return &SegmentFE;
1329  case Geometry::TRIANGLE: return &TriangleFE;
1331  case Geometry::TETRAHEDRON: return &TetrahedronFE;
1332  case Geometry::CUBE: return &ParallelepipedFE;
1333  default:
1334  mfem_error ("RefinedLinearFECollection: unknown geometry type.");
1335  }
1336  return &SegmentFE; // Make some compilers happy
1337 }
1338
1340 {
1341  switch (GeomType)
1342  {
1343  case Geometry::POINT: return 1;
1344  case Geometry::SEGMENT: return 1;
1345  case Geometry::TRIANGLE: return 0;
1346  case Geometry::SQUARE: return 1;
1347  case Geometry::TETRAHEDRON: return 0;
1348  case Geometry::CUBE: return 1;
1349  default:
1350  mfem_error ("RefinedLinearFECollection: unknown geometry type.");
1351  }
1352  return 0; // Make some compilers happy
1353 }
1354
1356  Geometry::Type GeomType, int Or) const
1357 {
1358  static int indexes[] = { 0 };
1359
1360  return indexes;
1361 }
1362
1363
1364 const FiniteElement *
1366 {
1367  switch (GeomType)
1368  {
1369  case Geometry::CUBE: return &HexahedronFE;
1370  case Geometry::TETRAHEDRON: return &TetrahedronFE;
1371  default:
1372  mfem_error ("ND1_3DFECollection: unknown geometry type.");
1373  }
1374  return &HexahedronFE; // Make some compilers happy
1375 }
1376
1378 {
1379  switch (GeomType)
1380  {
1381  case Geometry::POINT: return 0;
1382  case Geometry::SEGMENT: return 1;
1383  case Geometry::TRIANGLE: return 0;
1384  case Geometry::SQUARE: return 0;
1385  case Geometry::TETRAHEDRON: return 0;
1386  case Geometry::CUBE: return 0;
1387  default:
1388  mfem_error ("ND1_3DFECollection: unknown geometry type.");
1389  }
1390  return 0; // Make some compilers happy
1391 }
1392
1394  int Or) const
1395 {
1396  static int ind_pos[] = { 0 };
1397  static int ind_neg[] = { -1 };
1398
1399  if (Or > 0)
1400  {
1401  return ind_pos;
1402  }
1403  return ind_neg;
1404 }
1405
1406
1407 const FiniteElement *
1409 {
1410  switch (GeomType)
1411  {
1412  case Geometry::TRIANGLE: return &TriangleFE;
1414  case Geometry::CUBE: return &HexahedronFE;
1415  case Geometry::TETRAHEDRON: return &TetrahedronFE;
1416  default:
1417  mfem_error ("RT0_3DFECollection: unknown geometry type.");
1418  }
1419  return &HexahedronFE; // Make some compilers happy
1420 }
1421
1423 {
1424  switch (GeomType)
1425  {
1426  case Geometry::POINT: return 0;
1427  case Geometry::SEGMENT: return 0;
1428  case Geometry::TRIANGLE: return 1;
1429  case Geometry::SQUARE: return 1;
1430  case Geometry::TETRAHEDRON: return 0;
1431  case Geometry::CUBE: return 0;
1432  default:
1433  mfem_error ("RT0_3DFECollection: unknown geometry type.");
1434  }
1435  return 0; // Make some compilers happy
1436 }
1437
1439  int Or) const
1440 {
1441  static int ind_pos[] = { 0 };
1442  static int ind_neg[] = { -1 };
1443
1444  if ((GeomType == Geometry::TRIANGLE) || (GeomType == Geometry::SQUARE))
1445  {
1446  if (Or % 2 == 0)
1447  {
1448  return ind_pos;
1449  }
1450  return ind_neg;
1451  }
1452  return NULL;
1453 }
1454
1455 const FiniteElement *
1457 {
1458  switch (GeomType)
1459  {
1460  case Geometry::TRIANGLE: return &TriangleFE;
1462  case Geometry::CUBE: return &HexahedronFE;
1463  default:
1464  mfem_error ("RT1_3DFECollection: unknown geometry type.");
1465  }
1466  return &HexahedronFE; // Make some compilers happy
1467 }
1468
1470 {
1471  switch (GeomType)
1472  {
1473  case Geometry::POINT: return 0;
1474  case Geometry::SEGMENT: return 0;
1475  case Geometry::TRIANGLE: return 2;
1476  case Geometry::SQUARE: return 4;
1477  case Geometry::CUBE: return 12;
1478  default:
1479  mfem_error ("RT1_3DFECollection: unknown geometry type.");
1480  }
1481  return 0; // Make some compilers happy
1482 }
1483
1485  int Or) const
1486 {
1487  if (GeomType == Geometry::SQUARE)
1488  {
1489  static int sq_ind[8][4] =
1490  {
1491  {0, 1, 2, 3}, {-1, -3, -2, -4},
1492  {2, 0, 3, 1}, {-2, -1, -4, -3},
1493  {3, 2, 1, 0}, {-4, -2, -3, -1},
1494  {1, 3, 0, 2}, {-3, -4, -1, -2}
1495  };
1496
1497  return sq_ind[Or];
1498  }
1499  else
1500  {
1501  return NULL;
1502  }
1503 }
1504
1505
1506 H1_FECollection::H1_FECollection(const int p, const int dim, const int btype)
1507 {
1508  MFEM_VERIFY(p >= 1, "H1_FECollection requires order >= 1.");
1509  MFEM_VERIFY(dim >= 0 && dim <= 3, "H1_FECollection requires 0 <= dim <= 3.");
1510
1511  const int pm1 = p - 1, pm2 = pm1 - 1, pm3 = pm2 - 1;
1512
1514  b_type = BasisType::Check(btype);
1515  switch (btype)
1516  {
1518  {
1519  snprintf(h1_name, 32, "H1_%dD_P%d", dim, p);
1520  break;
1521  }
1522  case BasisType::Positive:
1523  {
1524  snprintf(h1_name, 32, "H1Pos_%dD_P%d", dim, p);
1525  break;
1526  }
1528  {
1529  snprintf(h1_name, 32, "H1Ser_%dD_P%d", dim, p);
1530  break;
1531  }
1532  default:
1533  {
1536  "unsupported BasisType: " << BasisType::Name(btype));
1537
1538  snprintf(h1_name, 32, "H1@%c_%dD_P%d",
1539  (int)BasisType::GetChar(btype), dim, p);
1540  }
1541  }
1542
1543  for (int g = 0; g < Geometry::NumGeom; g++)
1544  {
1545  H1_dof[g] = 0;
1546  H1_Elements[g] = NULL;
1547  }
1548  for (int i = 0; i < 2; i++)
1549  {
1550  SegDofOrd[i] = NULL;
1551  }
1552  for (int i = 0; i < 6; i++)
1553  {
1554  TriDofOrd[i] = NULL;
1555  }
1556  for (int i = 0; i < 8; i++)
1557  {
1559  }
1560
1561  H1_dof[Geometry::POINT] = 1;
1562  H1_Elements[Geometry::POINT] = new PointFiniteElement;
1563
1564  if (dim >= 1)
1565  {
1566  H1_dof[Geometry::SEGMENT] = pm1;
1567  if (b_type == BasisType::Positive)
1568  {
1569  H1_Elements[Geometry::SEGMENT] = new H1Pos_SegmentElement(p);
1570  }
1571  else
1572  {
1573  H1_Elements[Geometry::SEGMENT] = new H1_SegmentElement(p, btype);
1574  }
1575
1576  SegDofOrd[0] = new int[2*pm1];
1577  SegDofOrd[1] = SegDofOrd[0] + pm1;
1578  for (int i = 0; i < pm1; i++)
1579  {
1580  SegDofOrd[0][i] = i;
1581  SegDofOrd[1][i] = pm2 - i;
1582  }
1583  }
1584
1585  if (dim >= 2)
1586  {
1587  H1_dof[Geometry::TRIANGLE] = (pm1*pm2)/2;
1588  H1_dof[Geometry::SQUARE] = pm1*pm1;
1589  if (b_type == BasisType::Positive)
1590  {
1591  H1_Elements[Geometry::TRIANGLE] = new H1Pos_TriangleElement(p);
1593  }
1594  else if (b_type == BasisType::Serendipity)
1595  {
1596  // Note: in fe_coll.hpp the DofForGeometry(Geometry::Type) method
1597  // returns H1_dof[GeomType], so we need to fix the value of H1_dof here
1598  // for the serendipity case.
1599
1600  // formula for number of interior serendipity DoFs (when p>1)
1601  H1_dof[Geometry::SQUARE] = (pm3*pm2)/2;
1603  // allows for mixed tri/quad meshes
1604  H1_Elements[Geometry::TRIANGLE] = new H1Pos_TriangleElement(p);
1605  }
1606  else
1607  {
1608  H1_Elements[Geometry::TRIANGLE] = new H1_TriangleElement(p, btype);
1609  H1_Elements[Geometry::SQUARE] = new H1_QuadrilateralElement(p, btype);
1610  }
1611
1612  const int &TriDof = H1_dof[Geometry::TRIANGLE];
1613  const int &QuadDof = H1_dof[Geometry::SQUARE];
1614  TriDofOrd[0] = new int[6*TriDof];
1615  for (int i = 1; i < 6; i++)
1616  {
1617  TriDofOrd[i] = TriDofOrd[i-1] + TriDof;
1618  }
1619  // see Mesh::GetTriOrientation in mesh/mesh.cpp
1620  for (int j = 0; j < pm2; j++)
1621  {
1622  for (int i = 0; i + j < pm2; i++)
1623  {
1624  int o = TriDof - ((pm1 - j)*(pm2 - j))/2 + i;
1625  int k = pm3 - j - i;
1626  TriDofOrd[0][o] = o; // (0,1,2)
1627  TriDofOrd[1][o] = TriDof - ((pm1-j)*(pm2-j))/2 + k; // (1,0,2)
1628  TriDofOrd[2][o] = TriDof - ((pm1-i)*(pm2-i))/2 + k; // (2,0,1)
1629  TriDofOrd[3][o] = TriDof - ((pm1-k)*(pm2-k))/2 + i; // (2,1,0)
1630  TriDofOrd[4][o] = TriDof - ((pm1-k)*(pm2-k))/2 + j; // (1,2,0)
1631  TriDofOrd[5][o] = TriDof - ((pm1-i)*(pm2-i))/2 + j; // (0,2,1)
1632  }
1633  }
1634
1636  for (int i = 1; i < 8; i++)
1637  {
1639  }
1640
1641  // For serendipity order >=4, the QuadDofOrd array must be re-defined. We
1642  // do this by computing the corresponding tensor product QuadDofOrd array
1643  // or two orders less, which contains enough DoFs for their serendipity
1644  // basis. This could be optimized.
1645  if (b_type == BasisType::Serendipity)
1646  {
1647  if (p < 4)
1648  {
1650  }
1651  else // p >= 4 --> have face dofs
1652  {
1653  // Exactly the same as tensor product case, but with all orders
1654  // reduced by 2 e.g. in case p=5 it builds a 2x2 array, even though
1655  // there are only 3 serendipity dofs.
1656  // In the tensor product case, the i and j index tensor directions,
1657  // and o index from 0 to (pm1)^2,
1658  const int pm4 = pm3 -1;
1659
1660  for (int j = 0; j < pm3; j++) // pm3 instead of pm1, etc
1661  {
1662  for (int i = 0; i < pm3; i++)
1663  {
1664  int o = i + j*pm3;
1665  QuadDofOrd[0][o] = i + j*pm3; // (0,1,2,3)
1666  QuadDofOrd[1][o] = j + i*pm3; // (0,3,2,1)
1667  QuadDofOrd[2][o] = j + (pm4 - i)*pm3; // (1,2,3,0)
1668  QuadDofOrd[3][o] = (pm4 - i) + j*pm3; // (1,0,3,2)
1669  QuadDofOrd[4][o] = (pm4 - i) + (pm4 - j)*pm3; // (2,3,0,1)
1670  QuadDofOrd[5][o] = (pm4 - j) + (pm4 - i)*pm3; // (2,1,0,3)
1671  QuadDofOrd[6][o] = (pm4 - j) + i*pm3; // (3,0,1,2)
1672  QuadDofOrd[7][o] = i + (pm4 - j)*pm3; // (3,2,1,0)
1673  }
1674  }
1675
1676  }
1677  }
1678  else // not serendipity
1679  {
1680  for (int j = 0; j < pm1; j++)
1681  {
1682  for (int i = 0; i < pm1; i++)
1683  {
1684  int o = i + j*pm1;
1685  QuadDofOrd[0][o] = i + j*pm1; // (0,1,2,3)
1686  QuadDofOrd[1][o] = j + i*pm1; // (0,3,2,1)
1687  QuadDofOrd[2][o] = j + (pm2 - i)*pm1; // (1,2,3,0)
1688  QuadDofOrd[3][o] = (pm2 - i) + j*pm1; // (1,0,3,2)
1689  QuadDofOrd[4][o] = (pm2 - i) + (pm2 - j)*pm1; // (2,3,0,1)
1690  QuadDofOrd[5][o] = (pm2 - j) + (pm2 - i)*pm1; // (2,1,0,3)
1691  QuadDofOrd[6][o] = (pm2 - j) + i*pm1; // (3,0,1,2)
1692  QuadDofOrd[7][o] = i + (pm2 - j)*pm1; // (3,2,1,0)
1693  }
1694  }
1695  }
1696
1697  if (dim >= 3)
1698  {
1699  H1_dof[Geometry::TETRAHEDRON] = (TriDof*pm3)/3;
1701  H1_dof[Geometry::PRISM] = TriDof*pm1;
1702  if (b_type == BasisType::Positive)
1703  {
1704  H1_Elements[Geometry::TETRAHEDRON] = new H1Pos_TetrahedronElement(p);
1705  H1_Elements[Geometry::CUBE] = new H1Pos_HexahedronElement(p);
1706  H1_Elements[Geometry::PRISM] = new H1Pos_WedgeElement(p);
1707  }
1708  else
1709  {
1710  H1_Elements[Geometry::TETRAHEDRON] =
1711  new H1_TetrahedronElement(p, btype);
1712  H1_Elements[Geometry::CUBE] = new H1_HexahedronElement(p, btype);
1713  H1_Elements[Geometry::PRISM] = new H1_WedgeElement(p, btype);
1714  }
1715  }
1716  }
1717 }
1718
1720  int Or) const
1721 {
1722  if (GeomType == Geometry::SEGMENT)
1723  {
1724  return (Or > 0) ? SegDofOrd[0] : SegDofOrd[1];
1725  }
1726  else if (GeomType == Geometry::TRIANGLE)
1727  {
1728  return TriDofOrd[Or%6];
1729  }
1730  else if (GeomType == Geometry::SQUARE)
1731  {
1733  }
1734  return NULL;
1735 }
1736
1738 {
1739  int p = H1_dof[Geometry::SEGMENT] + 1;
1740  int dim = -1;
1741  if (!strncmp(h1_name, "H1_", 3))
1742  {
1743  dim = atoi(h1_name + 3);
1744  }
1745  else if (!strncmp(h1_name, "H1Pos_", 6))
1746  {
1747  dim = atoi(h1_name + 6);
1748  }
1749  else if (!strncmp(h1_name, "H1@", 3))
1750  {
1751  dim = atoi(h1_name + 5);
1752  }
1753  return (dim < 0) ? NULL : new H1_Trace_FECollection(p, dim, b_type);
1754 }
1755
1756 const int *H1_FECollection::GetDofMap(Geometry::Type GeomType) const
1757 {
1758  const int *dof_map = NULL;
1759  const FiniteElement *fe = H1_Elements[GeomType];
1760  switch (GeomType)
1761  {
1762  case Geometry::SEGMENT:
1763  case Geometry::SQUARE:
1764  case Geometry::CUBE:
1765  dof_map = dynamic_cast<const TensorBasisElement *>(fe)
1766  ->GetDofMap().GetData();
1767  break;
1768  default:
1769  MFEM_ABORT("Geometry type " << Geometry::Name[GeomType] << " is not "
1770  "implemented");
1771  // The "Cartesian" ordering for other geometries is defined by the
1772  // class GeometryRefiner.
1773  }
1774  return dof_map;
1775 }
1776
1778 {
1779  delete [] SegDofOrd[0];
1780  delete [] TriDofOrd[0];
1782  for (int g = 0; g < Geometry::NumGeom; g++)
1783  {
1784  delete H1_Elements[g];
1785  }
1786 }
1787
1788
1790  const int btype)
1791  : H1_FECollection(p, dim-1, btype)
1792 {
1793  if (btype == BasisType::GaussLobatto)
1794  {
1795  snprintf(h1_name, 32, "H1_Trace_%dD_P%d", dim, p);
1796  }
1797  else if (btype == BasisType::Positive)
1798  {
1799  snprintf(h1_name, 32, "H1Pos_Trace_%dD_P%d", dim, p);
1800  }
1801  else // base class checks that type is closed
1802  {
1803  snprintf(h1_name, 32, "H1_Trace@%c_%dD_P%d",
1804  (int)BasisType::GetChar(btype), dim, p);
1805  }
1806 }
1807
1808
1809 L2_FECollection::L2_FECollection(const int p, const int dim, const int btype,
1810  const int map_type)
1811 {
1812  MFEM_VERIFY(p >= 0, "L2_FECollection requires order >= 0.");
1813
1814  b_type = BasisType::Check(btype);
1815  const char *prefix = NULL;
1816  switch (map_type)
1817  {
1818  case FiniteElement::VALUE: prefix = "L2"; break;
1819  case FiniteElement::INTEGRAL: prefix = "L2Int"; break;
1820  default:
1821  MFEM_ABORT("invalid map_type: " << map_type);
1822  }
1823  switch (btype)
1824  {
1826  snprintf(d_name, 32, "%s_%dD_P%d", prefix, dim, p);
1827  break;
1828  default:
1829  snprintf(d_name, 32, "%s_T%d_%dD_P%d", prefix, btype, dim, p);
1830  }
1831
1832  for (int g = 0; g < Geometry::NumGeom; g++)
1833  {
1834  L2_Elements[g] = NULL;
1835  Tr_Elements[g] = NULL;
1836  }
1837  for (int i = 0; i < 2; i++)
1838  {
1839  SegDofOrd[i] = NULL;
1840  }
1841  for (int i = 0; i < 6; i++)
1842  {
1843  TriDofOrd[i] = NULL;
1844  }
1845  OtherDofOrd = NULL;
1846
1847  if (dim == 0)
1848  {
1849  L2_Elements[Geometry::POINT] = new PointFiniteElement;
1850  }
1851  else if (dim == 1)
1852  {
1853  if (b_type == BasisType::Positive)
1854  {
1855  L2_Elements[Geometry::SEGMENT] = new L2Pos_SegmentElement(p);
1856  }
1857  else
1858  {
1859  L2_Elements[Geometry::SEGMENT] = new L2_SegmentElement(p, btype);
1860  }
1861  L2_Elements[Geometry::SEGMENT]->SetMapType(map_type);
1862
1863  Tr_Elements[Geometry::POINT] = new PointFiniteElement;
1864  // No need to set the map_type for Tr_Elements.
1865
1866  const int pp1 = p + 1;
1867  SegDofOrd[0] = new int[2*pp1];
1868  SegDofOrd[1] = SegDofOrd[0] + pp1;
1869  for (int i = 0; i <= p; i++)
1870  {
1871  SegDofOrd[0][i] = i;
1872  SegDofOrd[1][i] = p - i;
1873  }
1874  }
1875  else if (dim == 2)
1876  {
1877  if (b_type == BasisType::Positive)
1878  {
1879  L2_Elements[Geometry::TRIANGLE] = new L2Pos_TriangleElement(p);
1881  }
1882  else
1883  {
1884  L2_Elements[Geometry::TRIANGLE] = new L2_TriangleElement(p, btype);
1885  L2_Elements[Geometry::SQUARE] = new L2_QuadrilateralElement(p, btype);
1886  }
1887  L2_Elements[Geometry::TRIANGLE]->SetMapType(map_type);
1888  L2_Elements[Geometry::SQUARE]->SetMapType(map_type);
1889  // Trace element use the default Gauss-Legendre nodal points for positive basis
1890  if (b_type == BasisType::Positive)
1891  {
1892  Tr_Elements[Geometry::SEGMENT] = new L2_SegmentElement(p);
1893  }
1894  else
1895  {
1896  Tr_Elements[Geometry::SEGMENT] = new L2_SegmentElement(p, btype);
1897  }
1898
1899  const int TriDof = L2_Elements[Geometry::TRIANGLE]->GetDof();
1900  TriDofOrd[0] = new int[6*TriDof];
1901  for (int i = 1; i < 6; i++)
1902  {
1903  TriDofOrd[i] = TriDofOrd[i-1] + TriDof;
1904  }
1905  const int pp1 = p + 1, pp2 = pp1 + 1;
1906  for (int j = 0; j <= p; j++)
1907  {
1908  for (int i = 0; i + j <= p; i++)
1909  {
1910  int o = TriDof - ((pp2 - j)*(pp1 - j))/2 + i;
1911  int k = p - j - i;
1912  TriDofOrd[0][o] = o; // (0,1,2)
1913  TriDofOrd[1][o] = TriDof - ((pp2-j)*(pp1-j))/2 + k; // (1,0,2)
1914  TriDofOrd[2][o] = TriDof - ((pp2-i)*(pp1-i))/2 + k; // (2,0,1)
1915  TriDofOrd[3][o] = TriDof - ((pp2-k)*(pp1-k))/2 + i; // (2,1,0)
1916  TriDofOrd[4][o] = TriDof - ((pp2-k)*(pp1-k))/2 + j; // (1,2,0)
1917  TriDofOrd[5][o] = TriDof - ((pp2-i)*(pp1-i))/2 + j; // (0,2,1)
1918  }
1919  }
1920  const int QuadDof = L2_Elements[Geometry::SQUARE]->GetDof();
1922  for (int j = 0; j < QuadDof; j++)
1923  {
1924  OtherDofOrd[j] = j; // for Or == 0
1925  }
1926  }
1927  else if (dim == 3)
1928  {
1929  if (b_type == BasisType::Positive)
1930  {
1931  L2_Elements[Geometry::TETRAHEDRON] = new L2Pos_TetrahedronElement(p);
1932  L2_Elements[Geometry::CUBE] = new L2Pos_HexahedronElement(p);
1933  L2_Elements[Geometry::PRISM] = new L2Pos_WedgeElement(p);
1934  }
1935  else
1936  {
1937  L2_Elements[Geometry::TETRAHEDRON] =
1938  new L2_TetrahedronElement(p, btype);
1939  L2_Elements[Geometry::CUBE] = new L2_HexahedronElement(p, btype);
1940  L2_Elements[Geometry::PRISM] = new L2_WedgeElement(p, btype);
1941  }
1942  L2_Elements[Geometry::TETRAHEDRON]->SetMapType(map_type);
1943  L2_Elements[Geometry::CUBE]->SetMapType(map_type);
1944  L2_Elements[Geometry::PRISM]->SetMapType(map_type);
1945  // Trace element use the default Gauss-Legendre nodal points for positive basis
1946  if (b_type == BasisType::Positive)
1947  {
1948  Tr_Elements[Geometry::TRIANGLE] = new L2_TriangleElement(p);
1950  }
1951  else
1952  {
1953  Tr_Elements[Geometry::TRIANGLE] = new L2_TriangleElement(p, btype);
1954  Tr_Elements[Geometry::SQUARE] = new L2_QuadrilateralElement(p, btype);
1955  }
1956
1957  const int TetDof = L2_Elements[Geometry::TETRAHEDRON]->GetDof();
1958  const int HexDof = L2_Elements[Geometry::CUBE]->GetDof();
1959  const int PriDof = L2_Elements[Geometry::PRISM]->GetDof();
1960  const int MaxDof = std::max(TetDof, std::max(PriDof, HexDof));
1961  OtherDofOrd = new int[MaxDof];
1962  for (int j = 0; j < MaxDof; j++)
1963  {
1964  OtherDofOrd[j] = j; // for Or == 0
1965  }
1966  }
1967  else
1968  {
1969  mfem::err << "L2_FECollection::L2_FECollection : dim = "
1970  << dim << endl;
1971  mfem_error();
1972  }
1973 }
1974
1976  int Or) const
1977 {
1978  switch (GeomType)
1979  {
1980  case Geometry::SEGMENT:
1981  return (Or > 0) ? SegDofOrd[0] : SegDofOrd[1];
1982
1983  case Geometry::TRIANGLE:
1984  return TriDofOrd[Or%6];
1985
1986  default:
1987  return (Or == 0) ? OtherDofOrd : NULL;
1988  }
1989 }
1990
1992 {
1993  delete [] OtherDofOrd;
1994  delete [] SegDofOrd[0];
1995  delete [] TriDofOrd[0];
1996  for (int i = 0; i < Geometry::NumGeom; i++)
1997  {
1998  delete L2_Elements[i];
1999  delete Tr_Elements[i];
2000  }
2001 }
2002
2003
2004 RT_FECollection::RT_FECollection(const int p, const int dim,
2005  const int cb_type, const int ob_type)
2006  : ob_type(ob_type)
2007 {
2008  MFEM_VERIFY(p >= 0, "RT_FECollection requires order >= 0.");
2009
2012
2014  {
2015  const char *cb_name = BasisType::Name(cb_type); // this may abort
2016  MFEM_ABORT("unknown closed BasisType: " << cb_name);
2017  }
2019  {
2020  const char *ob_name = BasisType::Name(ob_type); // this may abort
2021  MFEM_ABORT("unknown open BasisType: " << ob_name);
2022  }
2023
2024  InitFaces(p, dim, FiniteElement::INTEGRAL, true);
2025
2026  if (cb_type == BasisType::GaussLobatto &&
2027  ob_type == BasisType::GaussLegendre)
2028  {
2029  snprintf(rt_name, 32, "RT_%dD_P%d", dim, p);
2030  }
2031  else
2032  {
2033  snprintf(rt_name, 32, "RT@%c%c_%dD_P%d", (int)BasisType::GetChar(cb_type),
2034  (int)BasisType::GetChar(ob_type), dim, p);
2035  }
2036
2037  const int pp1 = p + 1;
2038  if (dim == 2)
2039  {
2040  // TODO: cb_type, ob_type for triangles
2042  RT_dof[Geometry::TRIANGLE] = p*pp1;
2043
2045  ob_type);
2046  // two vector components * n_unk_face *
2047  RT_dof[Geometry::SQUARE] = 2*p*pp1;
2048  }
2049  else if (dim == 3)
2050  {
2051  // TODO: cb_type, ob_type for tets
2053  RT_dof[Geometry::TETRAHEDRON] = p*pp1*(p + 2)/2;
2054
2055  RT_Elements[Geometry::CUBE] = new RT_HexahedronElement(p, cb_type, ob_type);
2056  RT_dof[Geometry::CUBE] = 3*p*pp1*pp1;
2057  }
2058  else
2059  {
2060  MFEM_ABORT("invalid dim = " << dim);
2061  }
2062 }
2063
2064 // This is a special protected constructor only used by RT_Trace_FECollection
2065 // and DG_Interface_FECollection
2066 RT_FECollection::RT_FECollection(const int p, const int dim, const int map_type,
2067  const bool signs, const int ob_type)
2068  : ob_type(ob_type)
2069 {
2072  {
2073  const char *ob_name = BasisType::Name(ob_type); // this may abort
2074  MFEM_ABORT("Invalid open basis type: " << ob_name);
2075  }
2076  InitFaces(p, dim, map_type, signs);
2077 }
2078
2079 void RT_FECollection::InitFaces(const int p, const int dim, const int map_type,
2080  const bool signs)
2081 {
2083
2085  "invalid open point type");
2086
2087  const int pp1 = p + 1, pp2 = p + 2;
2088
2089  for (int g = 0; g < Geometry::NumGeom; g++)
2090  {
2091  RT_Elements[g] = NULL;
2092  RT_dof[g] = 0;
2093  }
2094  // Degree of Freedom orderings
2095  for (int i = 0; i < 2; i++)
2096  {
2097  SegDofOrd[i] = NULL;
2098  }
2099  for (int i = 0; i < 6; i++)
2100  {
2101  TriDofOrd[i] = NULL;
2102  }
2103  for (int i = 0; i < 8; i++)
2104  {
2106  }
2107
2108  if (dim == 2)
2109  {
2110  L2_SegmentElement *l2_seg = new L2_SegmentElement(p, ob_type);
2111  l2_seg->SetMapType(map_type);
2112  RT_Elements[Geometry::SEGMENT] = l2_seg;
2113  RT_dof[Geometry::SEGMENT] = pp1;
2114
2115  SegDofOrd[0] = new int[2*pp1];
2116  SegDofOrd[1] = SegDofOrd[0] + pp1;
2117  for (int i = 0; i <= p; i++)
2118  {
2119  SegDofOrd[0][i] = i;
2120  SegDofOrd[1][i] = signs ? (-1 - (p - i)) : (p - i);
2121  }
2122  }
2123  else if (dim == 3)
2124  {
2125  L2_TriangleElement *l2_tri = new L2_TriangleElement(p, ob_type);
2126  l2_tri->SetMapType(map_type);
2127  RT_Elements[Geometry::TRIANGLE] = l2_tri;
2128  RT_dof[Geometry::TRIANGLE] = pp1*pp2/2;
2129
2133  RT_dof[Geometry::SQUARE] = pp1*pp1;
2134
2135  int TriDof = RT_dof[Geometry::TRIANGLE];
2136  TriDofOrd[0] = new int[6*TriDof];
2137  for (int i = 1; i < 6; i++)
2138  {
2139  TriDofOrd[i] = TriDofOrd[i-1] + TriDof;
2140  }
2141  // see Mesh::GetTriOrientation in mesh/mesh.cpp,
2142  // the constructor of H1_FECollection
2143  for (int j = 0; j <= p; j++)
2144  {
2145  for (int i = 0; i + j <= p; i++)
2146  {
2147  int o = TriDof - ((pp2 - j)*(pp1 - j))/2 + i;
2148  int k = p - j - i;
2149  TriDofOrd[0][o] = o; // (0,1,2)
2150  TriDofOrd[1][o] = -1-(TriDof-((pp2-j)*(pp1-j))/2+k); // (1,0,2)
2151  TriDofOrd[2][o] = TriDof-((pp2-i)*(pp1-i))/2+k; // (2,0,1)
2152  TriDofOrd[3][o] = -1-(TriDof-((pp2-k)*(pp1-k))/2+i); // (2,1,0)
2153  TriDofOrd[4][o] = TriDof-((pp2-k)*(pp1-k))/2+j; // (1,2,0)
2154  TriDofOrd[5][o] = -1-(TriDof-((pp2-i)*(pp1-i))/2+j); // (0,2,1)
2155  if (!signs)
2156  {
2157  for (int k = 1; k < 6; k += 2)
2158  {
2159  TriDofOrd[k][o] = -1 - TriDofOrd[k][o];
2160  }
2161  }
2162  }
2163  }
2164
2167  for (int i = 1; i < 8; i++)
2168  {
2170  }
2171  // see Mesh::GetQuadOrientation in mesh/mesh.cpp
2172  for (int j = 0; j <= p; j++)
2173  {
2174  for (int i = 0; i <= p; i++)
2175  {
2176  int o = i + j*pp1;
2177  QuadDofOrd[0][o] = i + j*pp1; // (0,1,2,3)
2178  QuadDofOrd[1][o] = -1 - (j + i*pp1); // (0,3,2,1)
2179  QuadDofOrd[2][o] = j + (p - i)*pp1; // (1,2,3,0)
2180  QuadDofOrd[3][o] = -1 - ((p - i) + j*pp1); // (1,0,3,2)
2181  QuadDofOrd[4][o] = (p - i) + (p - j)*pp1; // (2,3,0,1)
2182  QuadDofOrd[5][o] = -1 - ((p - j) + (p - i)*pp1); // (2,1,0,3)
2183  QuadDofOrd[6][o] = (p - j) + i*pp1; // (3,0,1,2)
2184  QuadDofOrd[7][o] = -1 - (i + (p - j)*pp1); // (3,2,1,0)
2185  if (!signs)
2186  {
2187  for (int k = 1; k < 8; k += 2)
2188  {
2190  }
2191  }
2192  }
2193  }
2194  }
2195 }
2196
2198  int Or) const
2199 {
2200  if (GeomType == Geometry::SEGMENT)
2201  {
2202  return (Or > 0) ? SegDofOrd[0] : SegDofOrd[1];
2203  }
2204  else if (GeomType == Geometry::TRIANGLE)
2205  {
2206  return TriDofOrd[Or%6];
2207  }
2208  else if (GeomType == Geometry::SQUARE)
2209  {
2211  }
2212  return NULL;
2213 }
2214
2216 {
2217  int dim, p;
2218  if (!strncmp(rt_name, "RT_", 3))
2219  {
2220  dim = atoi(rt_name + 3);
2221  p = atoi(rt_name + 7);
2222  }
2223  else // rt_name = RT@.._.D_P*
2224  {
2225  dim = atoi(rt_name + 6);
2226  p = atoi(rt_name + 10);
2227  }
2229 }
2230
2232 {
2233  delete [] SegDofOrd[0];
2234  delete [] TriDofOrd[0];
2236  for (int g = 0; g < Geometry::NumGeom; g++)
2237  {
2238  delete RT_Elements[g];
2239  }
2240 }
2241
2242
2244  const int map_type,
2245  const int ob_type)
2246  : RT_FECollection(p, dim, map_type, true, ob_type)
2247 {
2248  const char *prefix =
2249  (map_type == FiniteElement::INTEGRAL) ? "RT_Trace" : "RT_ValTrace";
2250  char ob_str[3] = { '\0', '\0', '\0' };
2251
2252  if (ob_type != BasisType::GaussLegendre)
2253  {
2254  ob_str[0] = '@';
2255  ob_str[1] = BasisType::GetChar(ob_type);
2256  }
2257  snprintf(rt_name, 32, "%s%s_%dD_P%d", prefix, ob_str, dim, p);
2258
2259  MFEM_VERIFY(dim == 2 || dim == 3, "Wrong dimension, dim = " << dim);
2260 }
2261
2262
2264  const int map_type,
2265  const int ob_type)
2266  : RT_FECollection(p, dim, map_type, false, ob_type)
2267 {
2268  MFEM_VERIFY(dim == 2 || dim == 3, "Wrong dimension, dim = " << dim);
2269
2270  const char *prefix =
2271  (map_type == FiniteElement::VALUE) ? "DG_Iface" : "DG_IntIface";
2272  if (ob_type == BasisType::GaussLegendre)
2273  {
2274  snprintf(rt_name, 32, "%s_%dD_P%d", prefix, dim, p);
2275  }
2276  else
2277  {
2278  snprintf(rt_name, 32, "%s@%c_%dD_P%d", prefix,
2279  (int)BasisType::GetChar(ob_type), dim, p);
2280  }
2281 }
2282
2283 ND_FECollection::ND_FECollection(const int p, const int dim,
2284  const int cb_type, const int ob_type)
2285 {
2286  MFEM_VERIFY(p >= 1, "ND_FECollection requires order >= 1.");
2287  MFEM_VERIFY(dim >= 1 && dim <= 3, "ND_FECollection requires 1 <= dim <= 3.");
2288
2289  const int pm1 = p - 1, pm2 = p - 2;
2290
2291  if (cb_type == BasisType::GaussLobatto &&
2292  ob_type == BasisType::GaussLegendre)
2293  {
2294  snprintf(nd_name, 32, "ND_%dD_P%d", dim, p);
2295  }
2296  else
2297  {
2298  snprintf(nd_name, 32, "ND@%c%c_%dD_P%d", (int)BasisType::GetChar(cb_type),
2299  (int)BasisType::GetChar(ob_type), dim, p);
2300  }
2301
2302  for (int g = 0; g < Geometry::NumGeom; g++)
2303  {
2304  ND_Elements[g] = NULL;
2305  ND_dof[g] = 0;
2306  }
2307  for (int i = 0; i < 2; i++)
2308  {
2309  SegDofOrd[i] = NULL;
2310  }
2311  for (int i = 0; i < 6; i++)
2312  {
2313  TriDofOrd[i] = NULL;
2314  }
2315  for (int i = 0; i < 8; i++)
2316  {
2318  }
2319
2322
2323  // Error checking
2325  {
2326  const char *ob_name = BasisType::Name(ob_type);
2327  MFEM_ABORT("Invalid open basis point type: " << ob_name);
2328  }
2330  {
2331  const char *cb_name = BasisType::Name(cb_type);
2332  MFEM_ABORT("Invalid closed basis point type: " << cb_name);
2333  }
2334
2335  if (dim >= 1)
2336  {
2337  ND_Elements[Geometry::SEGMENT] = new ND_SegmentElement(p, ob_type);
2339
2340  SegDofOrd[0] = new int[2*p];
2341  SegDofOrd[1] = SegDofOrd[0] + p;
2342  for (int i = 0; i < p; i++)
2343  {
2344  SegDofOrd[0][i] = i;
2345  SegDofOrd[1][i] = -1 - (pm1 - i);
2346  }
2347  }
2348
2349  if (dim >= 2)
2350  {
2352  ob_type);
2353  ND_dof[Geometry::SQUARE] = 2*p*pm1;
2354
2355  // TODO: cb_type and ob_type for triangles
2357  ND_dof[Geometry::TRIANGLE] = p*pm1;
2358
2361  for (int i = 1; i < 8; i++)
2362  {
2364  }
2365  // see Mesh::GetQuadOrientation in mesh/mesh.cpp
2366  for (int j = 0; j < pm1; j++)
2367  {
2368  for (int i = 0; i < p; i++)
2369  {
2370  int d1 = i + j*p; // x-component
2371  int d2 = p*pm1 + j + i*pm1; // y-component
2372  // (0,1,2,3)
2375  // (0,3,2,1)
2378  // (1,2,3,0)
2379  // QuadDofOrd[2][d1] = p*pm1 + (pm2 - j) + i*pm1;
2380  // QuadDofOrd[2][d2] = -1 - ((pm1 - i) + j*p);
2381  QuadDofOrd[2][d1] = -1 - (p*pm1 + j + (pm1 - i)*pm1);
2382  QuadDofOrd[2][d2] = i + (pm2 - j)*p;
2383  // (1,0,3,2)
2384  QuadDofOrd[3][d1] = -1 - ((pm1 - i) + j*p);
2385  QuadDofOrd[3][d2] = p*pm1 + (pm2 - j) + i*pm1;
2386  // (2,3,0,1)
2387  QuadDofOrd[4][d1] = -1 - ((pm1 - i) + (pm2 - j)*p);
2388  QuadDofOrd[4][d2] = -1 - (p*pm1 + (pm2 - j) + (pm1 - i)*pm1);
2389  // (2,1,0,3)
2390  QuadDofOrd[5][d1] = -1 - (p*pm1 + (pm2 - j) + (pm1 - i)*pm1);
2391  QuadDofOrd[5][d2] = -1 - ((pm1 - i) + (pm2 - j)*p);
2392  // (3,0,1,2)
2393  // QuadDofOrd[6][d1] = -1 - (p*pm1 + j + (pm1 - i)*pm1);
2394  // QuadDofOrd[6][d2] = i + (pm2 - j)*p;
2395  QuadDofOrd[6][d1] = p*pm1 + (pm2 - j) + i*pm1;
2396  QuadDofOrd[6][d2] = -1 - ((pm1 - i) + j*p);
2397  // (3,2,1,0)
2398  QuadDofOrd[7][d1] = i + (pm2 - j)*p;
2399  QuadDofOrd[7][d2] = -1 - (p*pm1 + j + (pm1 - i)*pm1);
2400  }
2401  }
2402
2403  int TriDof = ND_dof[Geometry::TRIANGLE];
2404  TriDofOrd[0] = new int[6*TriDof];
2405  for (int i = 1; i < 6; i++)
2406  {
2407  TriDofOrd[i] = TriDofOrd[i-1] + TriDof;
2408  }
2409  // see Mesh::GetTriOrientation in mesh/mesh.cpp,
2410  // the constructor of H1_FECollection
2411  for (int j = 0; j <= pm2; j++)
2412  {
2413  for (int i = 0; i + j <= pm2; i++)
2414  {
2415  int k1 = p*pm1 - (p - j)*(pm1 - j) + 2*i;
2416  int k2 = p*pm1 - (p - i)*(pm1 - i) + 2*j;
2417  // (0,1,2)
2418  TriDofOrd[0][k1 ] = k1;
2419  TriDofOrd[0][k1+1] = k1 + 1;
2420  // (0,2,1)
2421  TriDofOrd[5][k1 ] = k2 + 1;
2422  TriDofOrd[5][k1+1] = k2;
2423
2424  // The other orientations can not be supported with the current
2425  // interface. The method Mesh::ReorientTetMesh will ensure that
2426  // only orientations 0 and 5 are generated.
2427  }
2428  }
2429  }
2430
2431  if (dim >= 3)
2432  {
2433  ND_Elements[Geometry::CUBE] = new ND_HexahedronElement(p, cb_type, ob_type);
2434  ND_dof[Geometry::CUBE] = 3*p*pm1*pm1;
2435
2436  // TODO: cb_type and ob_type for tets
2438  ND_dof[Geometry::TETRAHEDRON] = p*pm1*pm2/2;
2439  }
2440 }
2441
2443  int Or) const
2444 {
2445  if (GeomType == Geometry::SEGMENT)
2446  {
2447  return (Or > 0) ? SegDofOrd[0] : SegDofOrd[1];
2448  }
2449  else if (GeomType == Geometry::TRIANGLE)
2450  {
2451  if (Or != 0 && Or != 5)
2452  {
2453  MFEM_ABORT("triangle face orientation " << Or << " is not supported! "
2454  "Use Mesh::ReorientTetMesh to fix it.");
2455  }
2456  return TriDofOrd[Or%6];
2457  }
2458  else if (GeomType == Geometry::SQUARE)
2459  {
2461  }
2462  return NULL;
2463 }
2464
2466 {
2467  int p, dim, cb_type, ob_type;
2468
2470  if (nd_name[2] == '_') // ND_
2471  {
2472  dim = atoi(nd_name + 3);
2473  cb_type = BasisType::GaussLobatto;
2474  ob_type = BasisType::GaussLegendre;
2475  }
2476  else // ND@
2477  {
2478  dim = atoi(nd_name + 6);
2479  cb_type = BasisType::GetType(nd_name[3]);
2480  ob_type = BasisType::GetType(nd_name[4]);
2481  }
2482  return new ND_Trace_FECollection(p, dim, cb_type, ob_type);
2483 }
2484
2486 {
2487  delete [] SegDofOrd[0];
2488  delete [] TriDofOrd[0];
2490  for (int g = 0; g < Geometry::NumGeom; g++)
2491  {
2492  delete ND_Elements[g];
2493  }
2494 }
2495
2496
2498  const int cb_type,
2499  const int ob_type)
2500  : ND_FECollection(p, dim-1, cb_type, ob_type)
2501 {
2502  if (cb_type == BasisType::GaussLobatto &&
2503  ob_type == BasisType::GaussLegendre)
2504  {
2505  snprintf(nd_name, 32, "ND_Trace_%dD_P%d", dim, p);
2506  }
2507  else
2508  {
2509  snprintf(nd_name, 32, "ND_Trace@%c%c_%dD_P%d",
2510  (int)BasisType::GetChar(cb_type),
2511  (int)BasisType::GetChar(ob_type), dim, p);
2512  }
2513 }
2514
2515
2517 {
2518  snprintf(d_name, 32, "Local_%s", fe_name);
2519
2520  Local_Element = NULL;
2521
2522  if (!strcmp(fe_name, "BiCubic2DFiniteElement") ||
2524  {
2525  GeomType = Geometry::SQUARE;
2526  Local_Element = new BiCubic2DFiniteElement;
2527  }
2528  else if (!strcmp(fe_name, "Nedelec1HexFiniteElement") ||
2529  !strcmp(fe_name, "Hex_ND1"))
2530  {
2531  GeomType = Geometry::CUBE;
2532  Local_Element = new Nedelec1HexFiniteElement;
2533  }
2534  else if (!strncmp(fe_name, "H1_", 3))
2535  {
2536  GeomType = Geometry::SQUARE;
2537  Local_Element = new H1_QuadrilateralElement(atoi(fe_name + 7));
2538  }
2539  else if (!strncmp(fe_name, "H1Pos_", 6))
2540  {
2541  GeomType = Geometry::SQUARE;
2542  Local_Element = new H1Pos_QuadrilateralElement(atoi(fe_name + 10));
2543  }
2544  else if (!strncmp(fe_name, "L2_", 3))
2545  {
2546  GeomType = Geometry::SQUARE;
2547  Local_Element = new L2_QuadrilateralElement(atoi(fe_name + 7));
2548  }
2549  else
2550  {
2551  mfem::err << "Local_FECollection::Local_FECollection : fe_name = "
2552  << fe_name << endl;
2553  mfem_error();
2554  }
2555 }
2556
2557
2559 {
2560  const int order = (Order == VariableOrder) ? 1 : Order;
2561  SegmentFE = new NURBS1DFiniteElement(order);
2563  ParallelepipedFE = new NURBS3DFiniteElement(order);
2564
2565  SetOrder(Order);
2566 }
2567
2568 void NURBSFECollection::SetOrder(int Order) const
2569 {
2570  mOrder = Order;
2571  if (Order != VariableOrder)
2572  {
2573  snprintf(name, 16, "NURBS%i", Order);
2574  }
2575  else
2576  {
2577  snprintf(name, 16, "NURBS");
2578  }
2579 }
2580
2582 {
2583  delete ParallelepipedFE;
2585  delete SegmentFE;
2586 }
2587
2588 const FiniteElement *
2590 {
2591  switch (GeomType)
2592  {
2593  case Geometry::SEGMENT: return SegmentFE;
2595  case Geometry::CUBE: return ParallelepipedFE;
2596  default:
2597  mfem_error ("NURBSFECollection: unknown geometry type.");
2598  }
2599  return SegmentFE; // Make some compilers happy
2600 }
2601
2603 {
2604  mfem_error("NURBSFECollection::DofForGeometry");
2605  return 0; // Make some compilers happy
2606 }
2607
2609  int Or) const
2610 {
2611  mfem_error("NURBSFECollection::DofOrderForOrientation");
2612  return NULL;
2613 }
2614
2616 {
2617  MFEM_ABORT("NURBS finite elements can not be statically condensed!");
2618  return NULL;
2619 }
2620
2621 }
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:838
Abstract class for Finite Elements.
Definition: fe.hpp:232
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:955
Arbitrary order non-uniform rational B-splines (NURBS) finite elements.
Definition: fe_coll.hpp:289
H1_Trace_FECollection(const int p, const int dim, const int btype=BasisType::GaussLobatto)
Definition: fe_coll.cpp:1789
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1169
int Size() const
Logical size of the array.
Definition: array.hpp:124
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:1161
For scalar fields; preserves point values.
Definition: fe.hpp:276
Version of QuadraticDiscont2DFECollection with dofs in the Gaussian points.
Definition: fe_coll.hpp:672
int RT_dof[Geometry::NumGeom]
Definition: fe_coll.hpp:197
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1031
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:977
FiniteElement * ND_Elements[Geometry::NumGeom]
Definition: fe_coll.hpp:255
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1210
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:1126
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:1240
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.hpp:45
virtual FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:41
virtual ~NURBSFECollection()
Definition: fe_coll.cpp:2581
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:1975
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:2066
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:574
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:2197
FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:2615
static const int NumGeom
Definition: geom.hpp:41
const Geometry::Type geom
Definition: ex1.cpp:40
virtual ~L2_FECollection()
Definition: fe_coll.cpp:1991
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:2589
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:759
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:616
Piecewise-(bi)linear continuous finite elements.
Definition: fe_coll.hpp:343
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:1719
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:928
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1377
Get the corresponding Quadrature1D constant, when that makes sense; otherwise return Quadrature1D::In...
Definition: fe.hpp:61
FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:2465
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:993
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:662
void InitFaces(const int p, const int dim, const int map_type, const bool signs)
Definition: fe_coll.cpp:2079
FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:2215
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:2263
int HasFaceDofs(Geometry::Type GeomType) const
Definition: fe_coll.cpp:25
static const char * Name(int b_type)
Check and convert a BasisType constant to a string identifier.
Definition: fe.hpp:89
static int GetType(char b_ident)
Convert char basis identifier to a BasisType constant.
Definition: fe.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:920
NURBSFECollection(int Order=VariableOrder)
The parameter Order must be either a positive number, for fixed order, or VariableOrder (default)...
Definition: fe_coll.cpp:2558
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1184
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1248
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:626
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:652
void SetMapType(int M)
Definition: fe.hpp:603
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:696
virtual ~RT_FECollection()
Definition: fe_coll.cpp:2231
L2_FECollection(const int p, const int dim, const int btype=BasisType::GaussLegendre, const int map_type=FiniteElement::VALUE)
Definition: fe_coll.cpp:1809
Piecewise-(bi)cubic continuous finite elements.
Definition: fe_coll.hpp:415
PointFiniteElement PointFE
Definition: point.cpp:30
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:824
const int * GetDofMap(Geometry::Type GeomType) const
Get the Cartesian to local H1 dof map.
Definition: fe_coll.cpp:1756
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:879
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:540
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:323
FiniteElement * RT_Elements[Geometry::NumGeom]
Definition: fe_coll.hpp:196
static const int Dimension[NumGeom]
Definition: geom.hpp:46
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:906
static void GetNVE(int &nv, int &ne)
Definition: fe_coll.cpp:293
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1365
For scalar fields; preserves volume integrals.
Definition: fe.hpp:277
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:1355
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:153
H1_FECollection(const int p, const int dim=3, const int btype=BasisType::GaussLobatto)
Definition: fe_coll.cpp:1506
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:2497
TriLinear3DFiniteElement HexahedronFE
Definition: hexahedron.cpp:52
ND_FECollection(const int p, const int dim, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Definition: fe_coll.cpp:2283
Discontinuous collection defined locally by a given finite element.
Definition: fe_coll.hpp:875
Version of QuadraticDiscont2DFECollection with positive basis functions.
Definition: fe_coll.hpp:655
virtual ~ND_FECollection()
Definition: fe_coll.cpp:2485
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:811
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:963
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:782
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:368
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:1059
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1261
Version of QuadraticFECollection with positive basis functions.
Definition: fe_coll.hpp:394
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:582
virtual ~H1_FECollection()
Definition: fe_coll.cpp:1777
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:1315
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:1277
FiniteElementCollection * GetTraceCollection() const
Definition: fe_coll.cpp:1737
static const char * Name[NumGeom]
Definition: geom.hpp:44
Piecewise-linear nonconforming finite elements in 3D.
Definition: fe_coll.hpp:463
Crouzeix-Raviart nonconforming elements in 2D.
Definition: fe_coll.hpp:442
class H1_WedgeElement WedgeFE
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:769
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1298
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1110
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:191
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1001
static void GetEdge(int &nv, v_t &v, int &ne, int &e, int &eo, const int edge_info)
Definition: fe_coll.cpp:303
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:1438
static int CheckClosed(int type)
If the Quadrature1D type is not closed return Invalid; otherwise return type.
Definition: intrules.cpp:824
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:2568
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:865
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:314
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
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:2442
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:745
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:635
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:599
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1285
class Linear3DFiniteElement TetrahedronFE
Definition: fe.cpp:12625
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1408
Linear2DFiniteElement TriangleFE
Definition: fe.cpp:12621
static int CheckOpen(int type)
If the Quadrature1D type is not open return Invalid; otherwise return type.
Definition: intrules.cpp:836
Serendipity basis (squares / cubes)
Definition: fe.hpp:39
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1339
Version of LinearDiscont2DFECollection with dofs in the Gaussian points.
Definition: fe_coll.hpp:596
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name...
Definition: fe_coll.cpp:47
static char GetChar(int b_type)
Check and convert a BasisType constant to a char basis identifier.
Definition: fe.hpp:99
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1322
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:1023
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:731
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1010
virtual const char * Name() const
Definition: fe_coll.hpp:53
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1469
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:1484
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:2608
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:639
int dim
Definition: ex24.cpp:43
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1067
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1079
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1223
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
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:679
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1422
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:2243
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:557
Definition: fe_coll.hpp:368
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:1200
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:852
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:251
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1044
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:83
Bernstein polynomials.
Definition: fe.hpp:35
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1095
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:2602
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:941
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:894
Local_FECollection(const char *fe_name)
Definition: fe_coll.cpp:2516
virtual int DofForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1147
Linear (P1) finite elements on quadrilaterals.
Definition: fe_coll.hpp:618
Linear1DFiniteElement SegmentFE
Definition: segment.cpp:49
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:1393
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1134
int ND_dof[Geometry::NumGeom]
Definition: fe_coll.hpp:256
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const
Definition: fe_coll.cpp:1456
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:143