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