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