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