MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
fe_base.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12// Finite Element Base classes
13
14#include "fe_base.hpp"
15#include "face_map_utils.hpp"
16#include "../coefficient.hpp"
17
18namespace mfem
19{
20
21using namespace std;
22
24 int Do, int O, int F)
25 : Nodes(Do)
26{
27 dim = D ; geom_type = G ; dof = Do ; order = O ; func_space = F;
28 vdim = 0 ; cdim = 0;
34 for (int i = 0; i < Geometry::MaxDim; i++) { orders[i] = -1; }
37#endif
38}
39
41 const IntegrationPoint &ip, DenseMatrix &shape) const
42{
43 MFEM_ABORT("method is not implemented for this class");
44}
45
47 ElementTransformation &Trans, DenseMatrix &shape) const
48{
49 MFEM_ABORT("method is not implemented for this class");
50}
51
53 const IntegrationPoint &ip, Vector &divshape) const
54{
55 MFEM_ABORT("method is not implemented for this class");
56}
57
59 ElementTransformation &Trans, Vector &div_shape) const
60{
61 CalcDivShape(Trans.GetIntPoint(), div_shape);
62 div_shape *= (1.0 / Trans.Weight());
63}
64
66 DenseMatrix &curl_shape) const
67{
68 MFEM_ABORT("method is not implemented for this class");
69}
70
72 DenseMatrix &curl_shape) const
73{
74 switch (dim)
75 {
76 case 3:
77 {
80#endif
82 MultABt(vshape, Trans.Jacobian(), curl_shape);
83 curl_shape *= (1.0 / Trans.Weight());
84 break;
85 }
86 case 2:
87 // This is valid for both 2x2 and 3x2 Jacobians
88 CalcCurlShape(Trans.GetIntPoint(), curl_shape);
89 curl_shape *= (1.0 / Trans.Weight());
90 break;
91 default:
92 MFEM_ABORT("Invalid dimension, Dim = " << dim);
93 }
94}
95
96void FiniteElement::GetFaceDofs(int face, int **dofs, int *ndofs) const
97{
99}
100
102 DenseMatrix &h) const
103{
105}
106
108 DenseMatrix &I) const
109{
111}
112
114 DenseMatrix &) const
115{
117}
118
121 DenseMatrix &I) const
122{
124}
125
127 Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
128{
130}
131
133 VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
134{
136}
137
139 Vector &dofs) const
140{
141 mfem_error("FiniteElement::ProjectFromNodes() (vector) is not overloaded!");
142}
143
146{
148}
149
150void FiniteElement::ProjectDelta(int vertex, Vector &dofs) const
151{
152 MFEM_ABORT("method is not implemented for this element");
153}
154
156 const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &I) const
157{
158 MFEM_ABORT("method is not implemented for this element");
159}
160
162 const FiniteElement &fe, ElementTransformation &Trans,
164{
165 MFEM_ABORT("method is not implemented for this element");
166}
167
169 const FiniteElement &fe, ElementTransformation &Trans,
170 DenseMatrix &curl) const
171{
172 MFEM_ABORT("method is not implemented for this element");
173}
174
176 const FiniteElement &fe, ElementTransformation &Trans,
177 DenseMatrix &div) const
178{
179 MFEM_ABORT("method is not implemented for this element");
180}
181
183 Vector &shape) const
184{
185 CalcShape(Trans.GetIntPoint(), shape);
186 if (map_type == INTEGRAL)
187 {
188 shape /= Trans.Weight();
189 }
190}
191
193 DenseMatrix &dshape) const
194{
195 MFEM_ASSERT(map_type == VALUE, "");
198#endif
199 CalcDShape(Trans.GetIntPoint(), vshape);
200 Mult(vshape, Trans.InverseJacobian(), dshape);
201}
202
204 Vector &Laplacian) const
205{
206 MFEM_ASSERT(map_type == VALUE, "");
207
208 // Simpler routine if mapping is affine
209 if (Trans.Hessian().FNorm2() < 1e-20)
210 {
211 CalcPhysLinLaplacian(Trans, Laplacian);
212 return;
213 }
214
215 // Compute full Hessian first if non-affine
216 int size = (dim*(dim+1))/2;
217 DenseMatrix hess(dof, size);
218 CalcPhysHessian(Trans,hess);
219
220 if (dim == 3)
221 {
222 for (int nd = 0; nd < dof; nd++)
223 {
224 Laplacian[nd] = hess(nd,0) + hess(nd,4) + hess(nd,5);
225 }
226 }
227 else if (dim == 2)
228 {
229 for (int nd = 0; nd < dof; nd++)
230 {
231 Laplacian[nd] = hess(nd,0) + hess(nd,2);
232 }
233 }
234 else
235 {
236 for (int nd = 0; nd < dof; nd++)
237 {
238 Laplacian[nd] = hess(nd,0);
239 }
240 }
241}
242
243// Assume a linear mapping
245 Vector &Laplacian) const
246{
247 MFEM_ASSERT(map_type == VALUE, "");
248 int size = (dim*(dim+1))/2;
249 DenseMatrix hess(dof, size);
250 DenseMatrix Gij(dim,dim);
251 Vector scale(size);
252
253 CalcHessian(Trans.GetIntPoint(), hess);
254 MultAAt(Trans.InverseJacobian(), Gij);
255
256 if (dim == 3)
257 {
258 scale[0] = Gij(0,0);
259 scale[1] = 2*Gij(0,1);
260 scale[2] = 2*Gij(0,2);
261
262 scale[3] = 2*Gij(1,2);
263 scale[4] = Gij(2,2);
264
265 scale[5] = Gij(1,1);
266 }
267 else if (dim == 2)
268 {
269 scale[0] = Gij(0,0);
270 scale[1] = 2*Gij(0,1);
271 scale[2] = Gij(1,1);
272 }
273 else
274 {
275 scale[0] = Gij(0,0);
276 }
277
278 for (int nd = 0; nd < dof; nd++)
279 {
280 Laplacian[nd] = 0.0;
281 for (int ii = 0; ii < size; ii++)
282 {
283 Laplacian[nd] += hess(nd,ii)*scale[ii];
284 }
285 }
286}
287
289 DenseMatrix& Hessian) const
290{
291 MFEM_ASSERT(map_type == VALUE, "");
292
293 // Roll 2-Tensors in vectors and 4-Tensor in Matrix, exploiting symmetry
294 Array<int> map(dim*dim);
295 if (dim == 3)
296 {
297 map[0] = 0;
298 map[1] = 1;
299 map[2] = 2;
300
301 map[3] = 1;
302 map[4] = 5;
303 map[5] = 3;
304
305 map[6] = 2;
306 map[7] = 3;
307 map[8] = 4;
308 }
309 else if (dim == 2)
310 {
311 map[0] = 0;
312 map[1] = 1;
313
314 map[2] = 1;
315 map[3] = 2;
316 }
317 else
318 {
319 map[0] = 0;
320 }
321
322 // Hessian in ref coords
323 int size = (dim*(dim+1))/2;
324 DenseMatrix hess(dof, size);
325 CalcHessian(Trans.GetIntPoint(), hess);
326
327 // Gradient in physical coords
328 if (Trans.Hessian().FNorm2() > 1e-10)
329 {
332 DenseMatrix gmap(dof, size);
334 hess -= gmap;
335 }
336
337 // LHM
338 DenseMatrix lhm(size,size);
339 DenseMatrix invJ = Trans.Jacobian();
340 lhm = 0.0;
341 for (int i = 0; i < dim; i++)
342 {
343 for (int j = 0; j < dim; j++)
344 {
345 for (int k = 0; k < dim; k++)
346 {
347 for (int l = 0; l < dim; l++)
348 {
349 lhm(map[i*dim+j],map[k*dim+l]) += invJ(i,k)*invJ(j,l);
350 }
351 }
352 }
353 }
354 // Correct multiplicity
355 Vector mult(size);
356 mult = 0.0;
357 for (int i = 0; i < dim*dim; i++) { mult[map[i]]++; }
358 lhm.InvRightScaling(mult);
359
360 // Hessian in physical coords
361 lhm.Invert();
362 Mult(hess, lhm, Hessian);
363}
364
367{
369 MFEM_VERIFY(mode == DofToQuad::FULL, "invalid mode requested");
370
373#endif
374 {
375 for (int i = 0; i < dof2quad_array.Size(); i++)
376 {
378 if (d2q->IntRule != &ir || d2q->mode != mode) { d2q = nullptr; }
379 }
380 if (!d2q)
381 {
384#endif
386 const int nqpt = ir.GetNPoints();
387 d2q->FE = this;
388 d2q->IntRule = &ir;
389 d2q->mode = mode;
390 d2q->ndof = dof;
391 d2q->nqpt = nqpt;
392 switch (range_type)
393 {
394 case SCALAR:
395 {
396 d2q->B.SetSize(nqpt*dof);
397 d2q->Bt.SetSize(dof*nqpt);
398
399 Vector shape;
400 vshape.GetColumnReference(0, shape);
401 for (int i = 0; i < nqpt; i++)
402 {
403 const IntegrationPoint &ip = ir.IntPoint(i);
404 CalcShape(ip, shape);
405 for (int j = 0; j < dof; j++)
406 {
407 d2q->B[i+nqpt*j] = d2q->Bt[j+dof*i] = shape(j);
408 }
409 }
410 break;
411 }
412 case VECTOR:
413 {
414 d2q->B.SetSize(nqpt*dim*dof);
415 d2q->Bt.SetSize(dof*nqpt*dim);
416
417 for (int i = 0; i < nqpt; i++)
418 {
419 const IntegrationPoint &ip = ir.IntPoint(i);
420 CalcVShape(ip, vshape);
421 for (int d = 0; d < dim; d++)
422 {
423 for (int j = 0; j < dof; j++)
424 {
425 d2q->B[i+nqpt*(d+dim*j)] =
426 d2q->Bt[j+dof*(i+nqpt*d)] = vshape(j, d);
427 }
428 }
429 }
430 break;
431 }
433 // Skip B and Bt for unknown range type
434 break;
435 }
436 switch (deriv_type)
437 {
439 {
440 d2q->G.SetSize(nqpt*dim*dof);
441 d2q->Gt.SetSize(dof*nqpt*dim);
442
443 for (int i = 0; i < nqpt; i++)
444 {
445 const IntegrationPoint &ip = ir.IntPoint(i);
446 CalcDShape(ip, vshape);
447 for (int d = 0; d < dim; d++)
448 {
449 for (int j = 0; j < dof; j++)
450 {
451 d2q->G[i+nqpt*(d+dim*j)] =
452 d2q->Gt[j+dof*(i+nqpt*d)] = vshape(j, d);
453 }
454 }
455 }
456 break;
457 }
458 case DIV:
459 {
460 d2q->G.SetSize(nqpt*dof);
461 d2q->Gt.SetSize(dof*nqpt);
462
463 Vector divshape;
464 vshape.GetColumnReference(0, divshape);
465 for (int i = 0; i < nqpt; i++)
466 {
467 const IntegrationPoint &ip = ir.IntPoint(i);
468 CalcDivShape(ip, divshape);
469 for (int j = 0; j < dof; j++)
470 {
471 d2q->G[i+nqpt*j] = d2q->Gt[j+dof*i] = divshape(j);
472 }
473 }
474 break;
475 }
476 case CURL:
477 {
478 d2q->G.SetSize(nqpt*cdim*dof);
479 d2q->Gt.SetSize(dof*nqpt*cdim);
480
481 DenseMatrix curlshape(vshape.GetData(), dof, cdim); // cdim <= dim
482 for (int i = 0; i < nqpt; i++)
483 {
484 const IntegrationPoint &ip = ir.IntPoint(i);
485 CalcCurlShape(ip, curlshape);
486 for (int d = 0; d < cdim; d++)
487 {
488 for (int j = 0; j < dof; j++)
489 {
490 d2q->G[i+nqpt*(d+cdim*j)] =
491 d2q->Gt[j+dof*(i+nqpt*d)] = curlshape(j, d);
492 }
493 }
494 }
495 break;
496 }
497 case NONE:
498 // Skip G and Gt for unknown derivative type
499 break;
500 }
502 }
503 }
504 return *d2q;
505}
506
507void FiniteElement::GetFaceMap(const int face_id,
508 Array<int> &face_map) const
509{
510 MFEM_ABORT("method is not implemented for this element");
511}
512
514{
515 for (int i = 0; i < dof2quad_array.Size(); i++)
516 {
518 }
519}
520
521
524 const ScalarFiniteElement &fine_fe) const
525{
527 Vector vv(v, dim);
528 IntegrationPoint f_ip;
529
531 Vector shape(dof);
532#else
533 Vector shape;
534 vshape.GetColumnReference(0, shape);
535#endif
536
537 MFEM_ASSERT(map_type == fine_fe.GetMapType(), "");
538
539 I.SetSize(fine_fe.dof, dof);
540 for (int i = 0; i < fine_fe.dof; i++)
541 {
542 Trans.Transform(fine_fe.Nodes.IntPoint(i), vv);
543 f_ip.Set(v, dim);
544 CalcShape(f_ip, shape);
545 for (int j = 0; j < dof; j++)
546 {
547 if (fabs(I(i,j) = shape(j)) < 1.0e-12)
548 {
549 I(i,j) = 0.0;
550 }
551 }
552 }
553 if (map_type == INTEGRAL)
554 {
555 // assuming Trans is linear; this should be ok for all refinement types
557 I *= Trans.Weight();
558 }
559}
560
563 const ScalarFiniteElement &fine_fe) const
564{
565 // General "interpolation", defined by L2 projection
566
568 Vector vv(v, dim);
569 IntegrationPoint f_ip;
570
571 const int fs = fine_fe.GetDof(), cs = this->GetDof();
572 I.SetSize(fs, cs);
573 Vector fine_shape(fs), coarse_shape(cs);
574 DenseMatrix fine_mass(fs), fine_coarse_mass(fs, cs); // initialized with 0
575 const int ir_order =
576 std::max(GetOrder(), fine_fe.GetOrder()) + fine_fe.GetOrder();
577 const IntegrationRule &ir = IntRules.Get(fine_fe.GetGeomType(), ir_order);
578
579 for (int i = 0; i < ir.GetNPoints(); i++)
580 {
581 const IntegrationPoint &ip = ir.IntPoint(i);
582 fine_fe.CalcShape(ip, fine_shape);
583 Trans.Transform(ip, vv);
584 f_ip.Set(v, dim);
585 this->CalcShape(f_ip, coarse_shape);
586
589 }
590
591 DenseMatrixInverse fine_mass_inv(fine_mass);
592 fine_mass_inv.Mult(fine_coarse_mass, I);
593
594 if (map_type == INTEGRAL)
595 {
596 // assuming Trans is linear; this should be ok for all refinement types
598 I *= Trans.Weight();
599 }
600}
601
604 const ScalarFiniteElement &coarse_fe) const
605{
606 // General "restriction", defined by L2 projection
608 Vector vv(v, dim);
609
610 const int cs = coarse_fe.GetDof(), fs = this->GetDof();
611 R.SetSize(cs, fs);
612 Vector fine_shape(fs), coarse_shape(cs);
613 DenseMatrix coarse_mass(cs), coarse_fine_mass(cs, fs); // initialized with 0
614 const int ir_order = GetOrder() + coarse_fe.GetOrder();
615 const IntegrationRule &ir = IntRules.Get(coarse_fe.GetGeomType(), ir_order);
616
617 // integrate coarse_mass in the coarse space
618 for (int i = 0; i < ir.GetNPoints(); i++)
619 {
620 const IntegrationPoint &c_ip = ir.IntPoint(i);
621 coarse_fe.CalcShape(c_ip, coarse_shape);
623 }
624
625 // integrate coarse_fine_mass in the fine space
627 for (int i = 0; i < ir.GetNPoints(); i++)
628 {
629 const IntegrationPoint &f_ip = ir.IntPoint(i);
630 this->CalcShape(f_ip, fine_shape);
631 Trans.Transform(f_ip, vv);
632
633 IntegrationPoint c_ip;
634 c_ip.Set(v, dim);
635 coarse_fe.CalcShape(c_ip, coarse_shape);
637 coarse_fine_mass);
638 }
639
640 DenseMatrixInverse coarse_mass_inv(coarse_mass);
641 coarse_mass_inv.Mult(coarse_fine_mass, R);
642
643 if (map_type == INTEGRAL)
644 {
645 // assuming Trans is linear; this should be ok for all refinement types
647 R *= 1.0 / Trans.Weight();
648 }
649}
650
651void NodalFiniteElement::CreateLexicographicFullMap(const IntegrationRule &ir)
652const
653{
654 // Get the FULL version of the map.
656 //Undo the native ordering which is what FiniteElement::GetDofToQuad returns.
657 auto *d2q_new = new DofToQuad(d2q);
659 const int nqpt = ir.GetNPoints();
660
661 const int b_dim = (range_type == VECTOR) ? dim : 1;
662
663 for (int i = 0; i < nqpt; i++)
664 {
665 for (int d = 0; d < b_dim; d++)
666 {
667 for (int j = 0; j < dof; j++)
668 {
669 const double val = d2q.B[i + nqpt*(d+b_dim*lex_ordering[j])];
670 d2q_new->B[i+nqpt*(d+b_dim*j)] = val;
671 d2q_new->Bt[j+dof*(i+nqpt*d)] = val;
672 }
673 }
674 }
675
676 const int g_dim = [this]()
677 {
678 switch (deriv_type)
679 {
681 case DIV: return 1;
682 case CURL: return cdim;
683 default: return 0;
684 }
685 }();
686
687 for (int i = 0; i < nqpt; i++)
688 {
689 for (int d = 0; d < g_dim; d++)
690 {
691 for (int j = 0; j < dof; j++)
692 {
693 const double val = d2q.G[i + nqpt*(d+g_dim*lex_ordering[j])];
694 d2q_new->G[i+nqpt*(d+g_dim*j)] = val;
695 d2q_new->Gt[j+dof*(i+nqpt*d)] = val;
696 }
697 }
698 }
699
701}
702
705{
706 //Should make this loop a function of FiniteElement
707 for (int i = 0; i < dof2quad_array.Size(); i++)
708 {
710 if (d2q.IntRule == &ir && d2q.mode == mode) { return d2q; }
711 }
712
714 {
716 }
717 else
718 {
719 CreateLexicographicFullMap(ir);
721 }
722}
723
725 const FiniteElement &fe, ElementTransformation &Trans,
726 DenseMatrix &curl) const
727{
728 DenseMatrix curl_shape(fe.GetDof(), 1);
729
730 curl.SetSize(dof, fe.GetDof());
731 for (int i = 0; i < dof; i++)
732 {
733 fe.CalcCurlShape(Nodes.IntPoint(i), curl_shape);
734
735 real_t w = 1.0;
737 {
738 Trans.SetIntPoint(&Nodes.IntPoint(i));
739 w /= Trans.Weight();
740 }
741 for (int j = 0; j < fe.GetDof(); j++)
742 {
743 curl(i,j) = w * curl_shape(j,0);
744 }
745 }
746}
747
749 const IntegrationPoint &pt, Vector &x)
750{
751 // invert a linear transform with one Newton step
753 p0.Set3(0, 0, 0);
754 trans.Transform(p0, x);
755
756 real_t store[3];
757 Vector v(store, x.Size());
758 pt.Get(store, x.Size());
759 v -= x;
760
761 trans.InverseJacobian().Mult(v, x);
762}
763
765 DenseMatrix &R) const
766{
768 Vector pt(&ipt.x, dim);
769
771 Vector shape(dof);
772#else
773 Vector shape;
774 vshape.GetColumnReference(0, shape);
775#endif
776
777 Trans.SetIntPoint(&Nodes[0]);
778
779 for (int j = 0; j < dof; j++)
780 {
781 InvertLinearTrans(Trans, Nodes[j], pt);
782 if (Geometries.CheckPoint(geom_type, ipt)) // do we need an epsilon here?
783 {
784 CalcShape(ipt, shape);
785 R.SetRow(j, shape);
786 }
787 else
788 {
789 // Set the whole row to avoid valgrind warnings in R.Threshold().
790 R.SetRow(j, infinity());
791 }
792 }
793 R.Threshold(1e-12);
794}
795
797 Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
798{
799 for (int i = 0; i < dof; i++)
800 {
801 const IntegrationPoint &ip = Nodes.IntPoint(i);
802 // some coefficients expect that Trans.IntPoint is the same
803 // as the second argument of Eval
804 Trans.SetIntPoint(&ip);
805 dofs(i) = coeff.Eval(Trans, ip);
806 if (map_type == INTEGRAL)
807 {
808 dofs(i) *= Trans.Weight();
809 }
810 }
811}
812
814 VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
815{
816 MFEM_ASSERT(dofs.Size() == vc.GetVDim()*dof, "");
817 Vector x(vc.GetVDim());
818
819 for (int i = 0; i < dof; i++)
820 {
821 const IntegrationPoint &ip = Nodes.IntPoint(i);
822 Trans.SetIntPoint(&ip);
823 vc.Eval (x, Trans, ip);
824 if (map_type == INTEGRAL)
825 {
826 x *= Trans.Weight();
827 }
828 for (int j = 0; j < x.Size(); j++)
829 {
830 dofs(dof*j+i) = x(j);
831 }
832 }
833}
834
837{
838 // (mc.height x mc.width) @ DOFs -> (dof x mc.width x mc.height) in dofs
839 MFEM_ASSERT(dofs.Size() == mc.GetHeight()*mc.GetWidth()*dof, "");
840 DenseMatrix MQ(mc.GetHeight(), mc.GetWidth());
841
842 for (int k = 0; k < dof; k++)
843 {
845 mc.Eval(MQ, T, Nodes.IntPoint(k));
846 if (map_type == INTEGRAL) { MQ *= T.Weight(); }
847 for (int r = 0; r < MQ.Height(); r++)
848 {
849 for (int d = 0; d < MQ.Width(); d++)
850 {
851 dofs(k+dof*(d+MQ.Width()*r)) = MQ(r,d);
852 }
853 }
854 }
855}
856
858 const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &I) const
859{
860 if (fe.GetRangeType() == SCALAR)
861 {
862 Vector shape(fe.GetDof());
863
864 I.SetSize(dof, fe.GetDof());
865 if (map_type == fe.GetMapType())
866 {
867 for (int k = 0; k < dof; k++)
868 {
869 fe.CalcShape(Nodes.IntPoint(k), shape);
870 for (int j = 0; j < shape.Size(); j++)
871 {
872 I(k,j) = (fabs(shape(j)) < 1e-12) ? 0.0 : shape(j);
873 }
874 }
875 }
876 else
877 {
878 for (int k = 0; k < dof; k++)
879 {
880 Trans.SetIntPoint(&Nodes.IntPoint(k));
881 fe.CalcPhysShape(Trans, shape);
882 if (map_type == INTEGRAL)
883 {
884 shape *= Trans.Weight();
885 }
886 for (int j = 0; j < shape.Size(); j++)
887 {
888 I(k,j) = (fabs(shape(j)) < 1e-12) ? 0.0 : shape(j);
889 }
890 }
891 }
892 }
893 else
894 {
895 DenseMatrix vshape(fe.GetDof(), std::max(Trans.GetSpaceDim(),
896 fe.GetRangeDim()));
897
898 I.SetSize(vshape.Width()*dof, fe.GetDof());
899 for (int k = 0; k < dof; k++)
900 {
901 Trans.SetIntPoint(&Nodes.IntPoint(k));
902 fe.CalcVShape(Trans, vshape);
903 if (map_type == INTEGRAL)
904 {
905 vshape *= Trans.Weight();
906 }
907 for (int j = 0; j < vshape.Height(); j++)
908 for (int d = 0; d < vshape.Width(); d++)
909 {
910 I(k+d*dof,j) = vshape(j,d);
911 }
912 }
913 }
914}
915
917 const FiniteElement &fe, ElementTransformation &Trans,
919{
920 MFEM_ASSERT(fe.GetMapType() == VALUE, "");
921 MFEM_ASSERT(Trans.GetSpaceDim() == dim, "")
922
923 DenseMatrix dshape(fe.GetDof(), dim), grad_k(fe.GetDof(), dim), Jinv(dim);
924
926 for (int k = 0; k < dof; k++)
927 {
928 const IntegrationPoint &ip = Nodes.IntPoint(k);
929 fe.CalcDShape(ip, dshape);
930 Trans.SetIntPoint(&ip);
931 CalcInverse(Trans.Jacobian(), Jinv);
933 if (map_type == INTEGRAL)
934 {
936 }
937 for (int j = 0; j < grad_k.Height(); j++)
938 for (int d = 0; d < dim; d++)
939 {
941 }
942 }
943}
944
946 const FiniteElement &fe, ElementTransformation &Trans,
947 DenseMatrix &div) const
948{
949 real_t detJ;
950 Vector div_shape(fe.GetDof());
951
952 div.SetSize(dof, fe.GetDof());
953 for (int k = 0; k < dof; k++)
954 {
955 const IntegrationPoint &ip = Nodes.IntPoint(k);
956 fe.CalcDivShape(ip, div_shape);
957 if (map_type == VALUE)
958 {
959 Trans.SetIntPoint(&ip);
960 detJ = Trans.Weight();
961 for (int j = 0; j < div_shape.Size(); j++)
962 {
963 div(k,j) = (fabs(div_shape(j)) < 1e-12) ? 0.0 : div_shape(j)/detJ;
964 }
965 }
966 else
967 {
968 for (int j = 0; j < div_shape.Size(); j++)
969 {
970 div(k,j) = (fabs(div_shape(j)) < 1e-12) ? 0.0 : div_shape(j);
971 }
972 }
973 }
974}
975
976
978 int Do, int O, int M, int F)
979 : FiniteElement(D, G, Do, O, F)
980{
982 map_type = M;
984 is_nodal = true;
985 vdim = dim;
986 if (map_type == H_CURL)
987 {
988 cdim = (dim == 3) ? 3 : 1;
989 }
990}
991
992void VectorFiniteElement::CalcShape(
993 const IntegrationPoint &ip, Vector &shape) const
994{
995 mfem_error("Error: Cannot use scalar CalcShape(...) function with\n"
996 " VectorFiniteElements!");
997}
998
999void VectorFiniteElement::CalcDShape(
1000 const IntegrationPoint &ip, DenseMatrix &dshape) const
1001{
1002 mfem_error("Error: Cannot use scalar CalcDShape(...) function with\n"
1003 " VectorFiniteElements!");
1004}
1005
1007{
1008 switch (map_type)
1009 {
1010 case H_DIV:
1011 deriv_type = DIV;
1014 break;
1015 case H_CURL:
1016 switch (dim)
1017 {
1018 case 3: // curl: 3D H_CURL -> 3D H_DIV
1019 deriv_type = CURL;
1022 break;
1023 case 2:
1024 // curl: 2D H_CURL -> INTEGRAL
1025 deriv_type = CURL;
1028 break;
1029 case 1:
1030 deriv_type = NONE;
1033 break;
1034 default:
1035 MFEM_ABORT("Invalid dimension, Dim = " << dim);
1036 }
1037 break;
1038 default:
1039 MFEM_ABORT("Invalid MapType = " << map_type);
1040 }
1041}
1042
1044 ElementTransformation &Trans, DenseMatrix &shape) const
1045{
1046 MFEM_ASSERT(map_type == H_DIV, "");
1049#endif
1050 CalcVShape(Trans.GetIntPoint(), vshape);
1051 MultABt(vshape, Trans.Jacobian(), shape);
1052 shape *= (1.0 / Trans.Weight());
1053}
1054
1056 ElementTransformation &Trans, DenseMatrix &shape) const
1057{
1058 MFEM_ASSERT(map_type == H_CURL, "");
1061#endif
1062 CalcVShape(Trans.GetIntPoint(), vshape);
1063 Mult(vshape, Trans.InverseJacobian(), shape);
1064}
1065
1067 const real_t *nk, const Array<int> &d2n,
1068 VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
1069{
1071 const int sdim = Trans.GetSpaceDim();
1072 MFEM_ASSERT(vc.GetVDim() == sdim, "");
1073 Vector xk(vk, sdim);
1074 const bool square_J = (dim == sdim);
1075
1076 for (int k = 0; k < dof; k++)
1077 {
1078 Trans.SetIntPoint(&Nodes.IntPoint(k));
1079 vc.Eval(xk, Trans, Nodes.IntPoint(k));
1080 // dof_k = nk^t adj(J) xk
1081 dofs(k) = Trans.AdjugateJacobian().InnerProduct(vk, nk + d2n[k]*dim);
1082 if (!square_J) { dofs(k) /= Trans.Weight(); }
1083 }
1084}
1085
1087 const real_t *nk, const Array<int> &d2n,
1088 Vector &vc, ElementTransformation &Trans, Vector &dofs) const
1089{
1090 const int sdim = Trans.GetSpaceDim();
1091 const bool square_J = (dim == sdim);
1092
1093 for (int k = 0; k < dof; k++)
1094 {
1095 Trans.SetIntPoint(&Nodes.IntPoint(k));
1096 // dof_k = nk^t adj(J) xk
1098 &vc[k*sdim], nk + d2n[k]*dim);
1099 if (!square_J) { dofs(k) /= Trans.Weight(); }
1100 }
1101}
1102
1104 const real_t *nk, const Array<int> &d2n,
1105 MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
1106{
1107 // project the rows of the matrix coefficient in an RT space
1108
1109 const int sdim = T.GetSpaceDim();
1110 MFEM_ASSERT(mc.GetWidth() == sdim, "");
1111 const bool square_J = (dim == sdim);
1112 DenseMatrix MQ(mc.GetHeight(), mc.GetWidth());
1113 Vector nk_phys(sdim), dofs_k(MQ.Height());
1114 MFEM_ASSERT(dofs.Size() == dof*MQ.Height(), "");
1115
1116 for (int k = 0; k < dof; k++)
1117 {
1118 T.SetIntPoint(&Nodes.IntPoint(k));
1119 mc.Eval(MQ, T, Nodes.IntPoint(k));
1120 // nk_phys = adj(J)^t nk
1122 if (!square_J) { nk_phys /= T.Weight(); }
1123 MQ.Mult(nk_phys, dofs_k);
1124 for (int r = 0; r < MQ.Height(); r++)
1125 {
1126 dofs(k+dof*r) = dofs_k(r);
1127 }
1128 }
1129}
1130
1132 const real_t *nk, const Array<int> &d2n, const FiniteElement &fe,
1133 ElementTransformation &Trans, DenseMatrix &I) const
1134{
1135 if (fe.GetRangeType() == SCALAR)
1136 {
1138 Vector shape(fe.GetDof());
1139 int sdim = Trans.GetSpaceDim();
1140
1141 I.SetSize(dof, sdim*fe.GetDof());
1142 for (int k = 0; k < dof; k++)
1143 {
1144 const IntegrationPoint &ip = Nodes.IntPoint(k);
1145
1146 fe.CalcShape(ip, shape);
1147 Trans.SetIntPoint(&ip);
1148 // Transform RT face normals from reference to physical space
1149 // vk = adj(J)^T nk
1151 if (fe.GetMapType() == INTEGRAL)
1152 {
1153 real_t w = 1.0/Trans.Weight();
1154 for (int d = 0; d < dim; d++)
1155 {
1156 vk[d] *= w;
1157 }
1158 }
1159
1160 for (int j = 0; j < shape.Size(); j++)
1161 {
1162 real_t s = shape(j);
1163 if (fabs(s) < 1e-12)
1164 {
1165 s = 0.0;
1166 }
1167 // Project scalar basis function multiplied by each coordinate
1168 // direction onto the transformed face normals
1169 for (int d = 0; d < sdim; d++)
1170 {
1171 I(k,j+d*shape.Size()) = s*vk[d];
1172 }
1173 }
1174 }
1175 }
1176 else
1177 {
1178 int sdim = Trans.GetSpaceDim();
1180 DenseMatrix vshape(fe.GetDof(), sdim);
1181 Vector vshapenk(fe.GetDof());
1182 const bool square_J = (dim == sdim);
1183
1184 I.SetSize(dof, fe.GetDof());
1185 for (int k = 0; k < dof; k++)
1186 {
1187 const IntegrationPoint &ip = Nodes.IntPoint(k);
1188
1189 Trans.SetIntPoint(&ip);
1190 // Transform RT face normals from reference to physical space
1191 // vk = adj(J)^T nk
1193 // Compute fe basis functions in physical space
1194 fe.CalcVShape(Trans, vshape);
1195 // Project fe basis functions onto transformed face normals
1196 vshape.Mult(vk, vshapenk);
1197 if (!square_J) { vshapenk /= Trans.Weight(); }
1198 for (int j=0; j<vshapenk.Size(); j++)
1199 {
1200 I(k,j) = vshapenk(j);
1201 }
1202 }
1203 }
1204}
1205
1207 const real_t *nk, const Array<int> &d2n, const FiniteElement &fe,
1208 ElementTransformation &Trans, DenseMatrix &grad) const
1209{
1210 if (dim != 2)
1211 {
1212 mfem_error("VectorFiniteElement::ProjectGrad_RT works only in 2D!");
1213 }
1214
1215 DenseMatrix dshape(fe.GetDof(), fe.GetDim());
1217 real_t tk[2];
1218
1220 for (int k = 0; k < dof; k++)
1221 {
1222 fe.CalcDShape(Nodes.IntPoint(k), dshape);
1223 tk[0] = nk[d2n[k]*dim+1];
1224 tk[1] = -nk[d2n[k]*dim];
1226 for (int j = 0; j < grad_k.Size(); j++)
1227 {
1229 }
1230 }
1231}
1232
1234 const real_t *tk, const Array<int> &d2t, const FiniteElement &fe,
1235 ElementTransformation &Trans, DenseMatrix &curl) const
1236{
1241#else
1242 curlshape.SetSize(fe.GetDof(), dim);
1244 JtJ.SetSize(dim, dim);
1245#endif
1246
1247 Vector curl_k(fe.GetDof());
1248
1249 curl.SetSize(dof, fe.GetDof());
1250 for (int k = 0; k < dof; k++)
1251 {
1252 const IntegrationPoint &ip = Nodes.IntPoint(k);
1253
1254 // calculate J^t * J / |J|
1255 Trans.SetIntPoint(&ip);
1256 MultAtB(Trans.Jacobian(), Trans.Jacobian(), JtJ);
1257 JtJ *= 1.0 / Trans.Weight();
1258
1259 // transform curl of shapes (rows) by J^t * J / |J|
1260 fe.CalcCurlShape(ip, curlshape);
1262
1263 curlshape_J.Mult(tk + d2t[k]*dim, curl_k);
1264 for (int j = 0; j < curl_k.Size(); j++)
1265 {
1266 curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
1267 }
1268 }
1269}
1270
1272 const real_t *nk, const Array<int> &d2n, const FiniteElement &fe,
1273 ElementTransformation &Trans, DenseMatrix &curl) const
1274{
1275 DenseMatrix curl_shape(fe.GetDof(), dim);
1276 Vector curl_k(fe.GetDof());
1277
1278 curl.SetSize(dof, fe.GetDof());
1279 for (int k = 0; k < dof; k++)
1280 {
1281 fe.CalcCurlShape(Nodes.IntPoint(k), curl_shape);
1282 curl_shape.Mult(nk + d2n[k]*dim, curl_k);
1283 for (int j = 0; j < curl_k.Size(); j++)
1284 {
1285 curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
1286 }
1287 }
1288}
1289
1291 const real_t *tk, const Array<int> &d2t,
1292 VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
1293{
1295 Vector xk(vk, vc.GetVDim());
1296
1297 for (int k = 0; k < dof; k++)
1298 {
1299 Trans.SetIntPoint(&Nodes.IntPoint(k));
1300
1301 vc.Eval(xk, Trans, Nodes.IntPoint(k));
1302 // dof_k = xk^t J tk
1303 dofs(k) = Trans.Jacobian().InnerProduct(tk + d2t[k]*dim, vk);
1304 }
1305}
1306
1308 const real_t *tk, const Array<int> &d2t,
1309 Vector &vc, ElementTransformation &Trans, Vector &dofs) const
1310{
1311 for (int k = 0; k < dof; k++)
1312 {
1313 Trans.SetIntPoint(&Nodes.IntPoint(k));
1314 // dof_k = xk^t J tk
1315 dofs(k) = Trans.Jacobian().InnerProduct(tk + d2t[k]*dim, &vc[k*dim]);
1316 }
1317}
1318
1320 const real_t *tk, const Array<int> &d2t,
1321 MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
1322{
1323 // project the rows of the matrix coefficient in an ND space
1324
1325 const int sdim = T.GetSpaceDim();
1326 MFEM_ASSERT(mc.GetWidth() == sdim, "");
1327 DenseMatrix MQ(mc.GetHeight(), mc.GetWidth());
1328 Vector tk_phys(sdim), dofs_k(MQ.Height());
1329 MFEM_ASSERT(dofs.Size() == dof*MQ.Height(), "");
1330
1331 for (int k = 0; k < dof; k++)
1332 {
1333 T.SetIntPoint(&Nodes.IntPoint(k));
1334 mc.Eval(MQ, T, Nodes.IntPoint(k));
1335 // tk_phys = J tk
1336 T.Jacobian().Mult(tk + d2t[k]*dim, tk_phys);
1337 MQ.Mult(tk_phys, dofs_k);
1338 for (int r = 0; r < MQ.Height(); r++)
1339 {
1340 dofs(k+dof*r) = dofs_k(r);
1341 }
1342 }
1343}
1344
1346 const real_t *tk, const Array<int> &d2t, const FiniteElement &fe,
1347 ElementTransformation &Trans, DenseMatrix &I) const
1348{
1349 if (fe.GetRangeType() == SCALAR)
1350 {
1351 int sdim = Trans.GetSpaceDim();
1353 Vector shape(fe.GetDof());
1354
1355 I.SetSize(dof, sdim*fe.GetDof());
1356 for (int k = 0; k < dof; k++)
1357 {
1358 const IntegrationPoint &ip = Nodes.IntPoint(k);
1359
1360 fe.CalcShape(ip, shape);
1361 Trans.SetIntPoint(&ip);
1362 // Transform ND edge tengents from reference to physical space
1363 // vk = J tk
1364 Trans.Jacobian().Mult(tk + d2t[k]*dim, vk);
1365 if (fe.GetMapType() == INTEGRAL)
1366 {
1367 real_t w = 1.0/Trans.Weight();
1368 for (int d = 0; d < sdim; d++)
1369 {
1370 vk[d] *= w;
1371 }
1372 }
1373
1374 for (int j = 0; j < shape.Size(); j++)
1375 {
1376 real_t s = shape(j);
1377 if (fabs(s) < 1e-12)
1378 {
1379 s = 0.0;
1380 }
1381 // Project scalar basis function multiplied by each coordinate
1382 // direction onto the transformed edge tangents
1383 for (int d = 0; d < sdim; d++)
1384 {
1385 I(k, j + d*shape.Size()) = s*vk[d];
1386 }
1387 }
1388 }
1389 }
1390 else
1391 {
1392 int sdim = Trans.GetSpaceDim();
1394 DenseMatrix vshape(fe.GetDof(), sdim);
1395 Vector vshapetk(fe.GetDof());
1396
1397 I.SetSize(dof, fe.GetDof());
1398 for (int k = 0; k < dof; k++)
1399 {
1400 const IntegrationPoint &ip = Nodes.IntPoint(k);
1401
1402 Trans.SetIntPoint(&ip);
1403 // Transform ND edge tangents from reference to physical space
1404 // vk = J tk
1405 Trans.Jacobian().Mult(tk + d2t[k]*dim, vk);
1406 // Compute fe basis functions in physical space
1407 fe.CalcVShape(Trans, vshape);
1408 // Project fe basis functions onto transformed edge tangents
1409 vshape.Mult(vk, vshapetk);
1410 for (int j=0; j<vshapetk.Size(); j++)
1411 {
1412 I(k, j) = vshapetk(j);
1413 }
1414 }
1415 }
1416}
1417
1419 const real_t *tk, const Array<int> &d2t, const FiniteElement &fe,
1420 ElementTransformation &Trans, DenseMatrix &grad) const
1421{
1422 MFEM_ASSERT(fe.GetMapType() == VALUE, "");
1423
1424 DenseMatrix dshape(fe.GetDof(), fe.GetDim());
1426
1428 for (int k = 0; k < dof; k++)
1429 {
1430 fe.CalcDShape(Nodes.IntPoint(k), dshape);
1432 for (int j = 0; j < grad_k.Size(); j++)
1433 {
1435 }
1436 }
1437}
1438
1440 const VectorFiniteElement &cfe, ElementTransformation &Trans,
1441 DenseMatrix &I) const
1442{
1443 Vector v(dim);
1444 IntegrationPoint tr_ip;
1445
1446 const int fs = dof, cs = cfe.GetDof();
1447 I.SetSize(fs, cs);
1448 DenseMatrix fine_shape(fs, dim), coarse_shape(cs, cfe.GetDim());
1449 DenseMatrix fine_mass(fs), fine_coarse_mass(fs, cs); // initialized with 0
1450 const int ir_order =
1451 std::max(GetOrder(), this->GetOrder()) + this->GetOrder();
1452 const IntegrationRule &ir = IntRules.Get(this->GetGeomType(), ir_order);
1453
1456 for (int i = 0; i < ir.GetNPoints(); i++)
1457 {
1458 const IntegrationPoint &ip = ir.IntPoint(i);
1459 real_t w = ip.weight;
1460 this->CalcVShape(ip, fine_shape);
1461 Trans.Transform(ip, v);
1462 tr_ip.Set(v.GetData(), dim);
1463 cfe.CalcVShape(tr_ip, coarse_shape);
1464
1466 for (int k=0; k<fs; ++k)
1467 {
1468 for (int j=0; j<cs; ++j)
1469 {
1470 real_t Mkj = 0.0;
1471 for (int d1=0; d1<dim; ++d1)
1472 {
1473 for (int d2=0; d2<dim; ++d2)
1474 {
1476 }
1477 }
1478 fine_coarse_mass(k,j) += (fabs(Mkj) < 1e-12) ? 0.0 : Mkj;
1479 }
1480 }
1481 }
1482 DenseMatrixInverse fine_mass_inv(fine_mass);
1483 fine_mass_inv.Mult(fine_coarse_mass, I);
1484}
1485
1487 const VectorFiniteElement &cfe, const real_t *nk, const Array<int> &d2n,
1488 ElementTransformation &Trans, DenseMatrix &I) const
1489{
1490 MFEM_ASSERT(map_type == cfe.GetMapType(), "");
1491
1492 if (!is_nodal) { return LocalL2Projection_RT(cfe, Trans, I); }
1493
1495 Vector xk(vk, dim);
1498 DenseMatrix vshape(cfe.GetDof(), cfe.GetDim());
1499#else
1500 DenseMatrix vshape(cfe.vshape.Data(), cfe.GetDof(), cfe.GetDim());
1501#endif
1502 I.SetSize(dof, vshape.Height());
1503
1504 // assuming Trans is linear; this should be ok for all refinement types
1507 for (int k = 0; k < dof; k++)
1508 {
1509 Trans.Transform(Nodes.IntPoint(k), xk);
1510 ip.Set3(vk);
1511 cfe.CalcVShape(ip, vshape);
1512 // xk = |J| J^{-t} n_k
1514 // I_k = vshape_k.adj(J)^t.n_k, k=1,...,dof
1515 for (int j = 0; j < vshape.Height(); j++)
1516 {
1517 real_t Ikj = 0.;
1518 for (int i = 0; i < dim; i++)
1519 {
1520 Ikj += vshape(j, i) * vk[i];
1521 }
1522 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
1523 }
1524 }
1525}
1526
1528 const VectorFiniteElement &cfe,
1529 ElementTransformation &Trans, DenseMatrix &I) const
1530{
1531 Vector v(dim);
1532 IntegrationPoint tr_ip;
1533
1534 const int fs = dof, cs = cfe.GetDof();
1535 I.SetSize(fs, cs);
1536 DenseMatrix fine_shape(fs, dim), coarse_shape(cs, cfe.GetDim());
1537 DenseMatrix fine_mass(fs), fine_coarse_mass(fs, cs); // initialized with 0
1538 const int ir_order =
1539 std::max(GetOrder(), this->GetOrder()) + this->GetOrder();
1540 const IntegrationRule &ir = IntRules.Get(this->GetGeomType(), ir_order);
1541
1543 const DenseMatrix &J = Trans.Jacobian();
1544 for (int i = 0; i < ir.GetNPoints(); i++)
1545 {
1546 const IntegrationPoint &ip = ir.IntPoint(i);
1547 this->CalcVShape(ip, fine_shape);
1548 Trans.Transform(ip, v);
1549 tr_ip.Set(v.GetData(), dim);
1550 cfe.CalcVShape(tr_ip, coarse_shape);
1551
1553 for (int k=0; k<fs; ++k)
1554 {
1555 for (int j=0; j<cs; ++j)
1556 {
1557 real_t Mkj = 0.0;
1558 for (int d1=0; d1<dim; ++d1)
1559 {
1560 for (int d2=0; d2<dim; ++d2)
1561 {
1562 Mkj += ip.weight*fine_shape(k,d1)*J(d1,d2)*coarse_shape(j,d2);
1563 }
1564 }
1565 fine_coarse_mass(k,j) += (fabs(Mkj) < 1e-12) ? 0.0 : Mkj;
1566 }
1567 }
1568 }
1569 DenseMatrixInverse fine_mass_inv(fine_mass);
1570 fine_mass_inv.Mult(fine_coarse_mass, I);
1571}
1572
1574 const VectorFiniteElement &cfe, const real_t *tk, const Array<int> &d2t,
1575 ElementTransformation &Trans, DenseMatrix &I) const
1576{
1577 if (!is_nodal) { return LocalL2Projection_ND(cfe, Trans, I); }
1578
1580 Vector xk(vk, dim);
1583 DenseMatrix vshape(cfe.GetDof(), cfe.GetDim());
1584#else
1585 DenseMatrix vshape(cfe.vshape.Data(), cfe.GetDof(), cfe.GetDim());
1586#endif
1587 I.SetSize(dof, vshape.Height());
1588
1589 // assuming Trans is linear; this should be ok for all refinement types
1591 const DenseMatrix &J = Trans.Jacobian();
1592 for (int k = 0; k < dof; k++)
1593 {
1594 Trans.Transform(Nodes.IntPoint(k), xk);
1595 ip.Set3(vk);
1596 cfe.CalcVShape(ip, vshape);
1597 // xk = J t_k
1598 J.Mult(tk + d2t[k]*dim, vk);
1599 // I_k = vshape_k.J.t_k, k=1,...,Dof
1600 for (int j = 0; j < vshape.Height(); j++)
1601 {
1602 real_t Ikj = 0.;
1603 for (int i = 0; i < dim; i++)
1604 {
1605 Ikj += vshape(j, i) * vk[i];
1606 }
1607 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
1608 }
1609 }
1610}
1611
1613 const real_t *nk, const Array<int> &d2n, ElementTransformation &Trans,
1614 DenseMatrix &R) const
1615{
1616 real_t pt_data[Geometry::MaxDim];
1618 Vector pt(pt_data, dim);
1619
1622#endif
1623
1625 const DenseMatrix &J = Trans.Jacobian();
1626 const real_t weight = Trans.Weight();
1627 for (int j = 0; j < dof; j++)
1628 {
1629 InvertLinearTrans(Trans, Nodes.IntPoint(j), pt);
1630 ip.Set(pt_data, dim);
1631 if (Geometries.CheckPoint(geom_type, ip)) // do we need an epsilon here?
1632 {
1633 CalcVShape(ip, vshape);
1634 J.MultTranspose(nk+dim*d2n[j], pt_data);
1635 pt /= weight;
1636 for (int k = 0; k < dof; k++)
1637 {
1638 real_t R_jk = 0.0;
1639 for (int d = 0; d < dim; d++)
1640 {
1641 R_jk += vshape(k,d)*pt_data[d];
1642 }
1643 R(j,k) = R_jk;
1644 }
1645 }
1646 else
1647 {
1648 // Set the whole row to avoid valgrind warnings in R.Threshold().
1649 R.SetRow(j, infinity());
1650 }
1651 }
1652 R.Threshold(1e-12);
1653}
1654
1656 const real_t *tk, const Array<int> &d2t, ElementTransformation &Trans,
1657 DenseMatrix &R) const
1658{
1659 real_t pt_data[Geometry::MaxDim];
1661 Vector pt(pt_data, dim);
1662
1665#endif
1666
1668 const DenseMatrix &Jinv = Trans.InverseJacobian();
1669 for (int j = 0; j < dof; j++)
1670 {
1671 InvertLinearTrans(Trans, Nodes.IntPoint(j), pt);
1672 ip.Set(pt_data, dim);
1673 if (Geometries.CheckPoint(geom_type, ip)) // do we need an epsilon here?
1674 {
1675 CalcVShape(ip, vshape);
1676 Jinv.Mult(tk+dim*d2t[j], pt_data);
1677 for (int k = 0; k < dof; k++)
1678 {
1679 real_t R_jk = 0.0;
1680 for (int d = 0; d < dim; d++)
1681 {
1682 R_jk += vshape(k,d)*pt_data[d];
1683 }
1684 R(j,k) = R_jk;
1685 }
1686 }
1687 else
1688 {
1689 // Set the whole row to avoid valgrind warnings in R.Threshold().
1690 R.SetRow(j, infinity());
1691 }
1692 }
1693 R.Threshold(1e-12);
1694}
1695
1696
1697Poly_1D::Basis::Basis(const int p, const real_t *nodes, EvalType etype)
1698 : etype(etype), auxiliary_basis(NULL), scale_integrated(false)
1699{
1700 switch (etype)
1701 {
1702 case ChangeOfBasis:
1703 {
1704 x.SetSize(p + 1);
1705 w.SetSize(p + 1);
1706 DenseMatrix A(p + 1);
1707 for (int i = 0; i <= p; i++)
1708 {
1709 CalcBasis(p, nodes[i], A.GetColumn(i));
1710 }
1711 Ai.Factor(A);
1712 // mfem::out << "Poly_1D::Basis(" << p << ",...) : "; Ai.TestInversion();
1713 break;
1714 }
1715 case Barycentric:
1716 {
1717 x.SetSize(p + 1);
1718 w.SetSize(p + 1);
1719 x = nodes;
1720 w = 1.0;
1721 for (int i = 0; i <= p; i++)
1722 {
1723 for (int j = 0; j < i; j++)
1724 {
1725 real_t xij = x(i) - x(j);
1726 w(i) *= xij;
1727 w(j) *= -xij;
1728 }
1729 }
1730 for (int i = 0; i <= p; i++)
1731 {
1732 w(i) = 1.0/w(i);
1733 }
1734
1735#ifdef MFEM_DEBUG
1736 // Make sure the nodes are increasing
1737 for (int i = 0; i < p; i++)
1738 {
1739 if (x(i) >= x(i+1))
1740 {
1741 mfem_error("Poly_1D::Basis::Basis : nodes are not increasing!");
1742 }
1743 }
1744#endif
1745 break;
1746 }
1747 case Positive:
1748 x.SetDataAndSize(NULL, p + 1); // use x to store (p + 1)
1749 break;
1750 case Integrated:
1751 auxiliary_basis = new Basis(
1753 u_aux.SetSize(p+2);
1754 d_aux.SetSize(p+2);
1755 d2_aux.SetSize(p+2);
1756 break;
1757 default: break;
1758 }
1759}
1760
1762{
1763 switch (etype)
1764 {
1765 case ChangeOfBasis:
1766 {
1767 CalcBasis(Ai.Width() - 1, y, x);
1768 Ai.Mult(x, u);
1769 break;
1770 }
1771 case Barycentric:
1772 {
1773 int i, k, p = x.Size() - 1;
1774 real_t l, lk;
1775
1776 if (p == 0)
1777 {
1778 u(0) = 1.0;
1779 return;
1780 }
1781
1782 lk = 1.0;
1783 for (k = 0; k < p; k++)
1784 {
1785 if (y >= (x(k) + x(k+1))/2)
1786 {
1787 lk *= y - x(k);
1788 }
1789 else
1790 {
1791 for (i = k+1; i <= p; i++)
1792 {
1793 lk *= y - x(i);
1794 }
1795 break;
1796 }
1797 }
1798 l = lk * (y - x(k));
1799
1800 for (i = 0; i < k; i++)
1801 {
1802 u(i) = l * w(i) / (y - x(i));
1803 }
1804 u(k) = lk * w(k);
1805 for (i++; i <= p; i++)
1806 {
1807 u(i) = l * w(i) / (y - x(i));
1808 }
1809 break;
1810 }
1811 case Positive:
1812 CalcBernstein(x.Size() - 1, y, u);
1813 break;
1814 case Integrated:
1815 auxiliary_basis->Eval(y, u_aux, d_aux);
1816 EvalIntegrated(d_aux, u);
1817 break;
1818 default: break;
1819 }
1820}
1821
1822void Poly_1D::Basis::Eval(const real_t y, Vector &u, Vector &d) const
1823{
1824 switch (etype)
1825 {
1826 case ChangeOfBasis:
1827 {
1828 CalcBasis(Ai.Width() - 1, y, x, w);
1829 Ai.Mult(x, u);
1830 Ai.Mult(w, d);
1831 break;
1832 }
1833 case Barycentric:
1834 {
1835 int i, k, p = x.Size() - 1;
1836 real_t l, lp, lk, sk, si;
1837
1838 if (p == 0)
1839 {
1840 u(0) = 1.0;
1841 d(0) = 0.0;
1842 return;
1843 }
1844
1845 lk = 1.0;
1846 for (k = 0; k < p; k++)
1847 {
1848 if (y >= (x(k) + x(k+1))/2)
1849 {
1850 lk *= y - x(k);
1851 }
1852 else
1853 {
1854 for (i = k+1; i <= p; i++)
1855 {
1856 lk *= y - x(i);
1857 }
1858 break;
1859 }
1860 }
1861 l = lk * (y - x(k));
1862
1863 sk = 0.0;
1864 for (i = 0; i < k; i++)
1865 {
1866 si = 1.0/(y - x(i));
1867 sk += si;
1868 u(i) = l * si * w(i);
1869 }
1870 u(k) = lk * w(k);
1871 for (i++; i <= p; i++)
1872 {
1873 si = 1.0/(y - x(i));
1874 sk += si;
1875 u(i) = l * si * w(i);
1876 }
1877 lp = l * sk + lk;
1878
1879 for (i = 0; i < k; i++)
1880 {
1881 d(i) = (lp * w(i) - u(i))/(y - x(i));
1882 }
1883 d(k) = sk * u(k);
1884 for (i++; i <= p; i++)
1885 {
1886 d(i) = (lp * w(i) - u(i))/(y - x(i));
1887 }
1888 break;
1889 }
1890 case Positive:
1891 CalcBernstein(x.Size() - 1, y, u, d);
1892 break;
1893 case Integrated:
1894 auxiliary_basis->Eval(y, u_aux, d_aux, d2_aux);
1895 EvalIntegrated(d_aux,u);
1896 EvalIntegrated(d2_aux,d);
1897 break;
1898 default: break;
1899 }
1900}
1901
1903 Vector &d2) const
1904{
1905 MFEM_VERIFY(etype == Barycentric,
1906 "Basis::Eval with second order derivatives not implemented for"
1907 " etype = " << etype);
1908 switch (etype)
1909 {
1910 case ChangeOfBasis:
1911 {
1912 CalcBasis(Ai.Width() - 1, y, x, w);
1913 Ai.Mult(x, u);
1914 Ai.Mult(w, d);
1915 // set d2 (not implemented yet)
1916 break;
1917 }
1918 case Barycentric:
1919 {
1920 int i, k, p = x.Size() - 1;
1921 real_t l, lp, lp2, lk, sk, si, sk2;
1922
1923 if (p == 0)
1924 {
1925 u(0) = 1.0;
1926 d(0) = 0.0;
1927 d2(0) = 0.0;
1928 return;
1929 }
1930
1931 lk = 1.0;
1932 for (k = 0; k < p; k++)
1933 {
1934 if (y >= (x(k) + x(k+1))/2)
1935 {
1936 lk *= y - x(k);
1937 }
1938 else
1939 {
1940 for (i = k+1; i <= p; i++)
1941 {
1942 lk *= y - x(i);
1943 }
1944 break;
1945 }
1946 }
1947 l = lk * (y - x(k));
1948
1949 sk = 0.0;
1950 sk2 = 0.0;
1951 for (i = 0; i < k; i++)
1952 {
1953 si = 1.0/(y - x(i));
1954 sk += si;
1955 sk2 -= si * si;
1956 u(i) = l * si * w(i);
1957 }
1958 u(k) = lk * w(k);
1959 for (i++; i <= p; i++)
1960 {
1961 si = 1.0/(y - x(i));
1962 sk += si;
1963 sk2 -= si * si;
1964 u(i) = l * si * w(i);
1965 }
1966 lp = l * sk + lk;
1967 lp2 = lp * sk + l * sk2 + sk * lk;
1968
1969 for (i = 0; i < k; i++)
1970 {
1971 d(i) = (lp * w(i) - u(i))/(y - x(i));
1972 d2(i) = (lp2 * w(i) - 2 * d(i))/(y - x(i));
1973 }
1974 d(k) = sk * u(k);
1975 d2(k) = sk2 * u(k) + sk * d(k);
1976 for (i++; i <= p; i++)
1977 {
1978 d(i) = (lp * w(i) - u(i))/(y - x(i));
1979 d2(i) = (lp2 * w(i) - 2 * d(i))/(y - x(i));
1980 }
1981 break;
1982 }
1983 case Positive:
1984 CalcBernstein(x.Size() - 1, y, u, d);
1985 break;
1986 case Integrated:
1987 MFEM_ABORT("Integrated basis must be evaluated with EvalIntegrated");
1988 break;
1989 default: break;
1990 }
1991}
1992
1994{
1995 MFEM_VERIFY(etype == Integrated,
1996 "EvalIntegrated is only valid for Integrated basis type");
1997 int p = d_aux_.Size() - 1;
1998 // See Gerritsma, M. (2010). "Edge functions for spectral element methods",
1999 // in Lecture Notes in Computational Science and Engineering, 199--207.
2000 u[0] = -d_aux_[0];
2001 for (int j=1; j<p; ++j)
2002 {
2003 u[j] = u[j-1] - d_aux_[j];
2004 }
2005 // If scale_integrated is true, the degrees of freedom represent mean values,
2006 // otherwise they represent subcell integrals. Generally, scale_integrated
2007 // should be true for MapType::VALUE, and false for other map types.
2008 if (scale_integrated)
2009 {
2010 Vector &aux_nodes = auxiliary_basis->x;
2011 for (int j=0; j<aux_nodes.Size()-1; ++j)
2012 {
2013 u[j] *= aux_nodes[j+1] - aux_nodes[j];
2014 }
2015 }
2016}
2017
2018void Poly_1D::Basis::ScaleIntegrated(bool scale_integrated_)
2019{
2020 scale_integrated = scale_integrated_;
2021}
2022
2024{
2025 delete auxiliary_basis;
2026}
2027
2028const int *Poly_1D::Binom(const int p)
2029{
2030 if (binom.NumCols() <= p)
2031 {
2032 binom.SetSize(p + 1, p + 1);
2033 for (int i = 0; i <= p; i++)
2034 {
2035 binom(i,0) = binom(i,i) = 1;
2036 for (int j = 1; j < i; j++)
2037 {
2038 binom(i,j) = binom(i-1,j) + binom(i-1,j-1);
2039 }
2040 }
2041 }
2042 return binom[p];
2043}
2044
2046{
2047 for (int i = 0; i <= p; i++)
2048 {
2049 // x[i] = 0.5*(1. + cos(M_PI*(p - i + 0.5)/(p + 1)));
2050 real_t s = sin(M_PI_2*(i + 0.5)/(p + 1));
2051 x[i] = s*s;
2052 }
2053}
2054
2055void Poly_1D::CalcMono(const int p, const real_t x, real_t *u)
2056{
2057 real_t xn;
2058 u[0] = xn = 1.;
2059 for (int n = 1; n <= p; n++)
2060 {
2061 u[n] = (xn *= x);
2062 }
2063}
2064
2065void Poly_1D::CalcMono(const int p, const real_t x, real_t *u, real_t *d)
2066{
2067 real_t xn;
2068 u[0] = xn = 1.;
2069 d[0] = 0.;
2070 for (int n = 1; n <= p; n++)
2071 {
2072 d[n] = n * xn;
2073 u[n] = (xn *= x);
2074 }
2075}
2076
2077void Poly_1D::CalcBinomTerms(const int p, const real_t x, const real_t y,
2078 real_t *u)
2079{
2080 if (p == 0)
2081 {
2082 u[0] = 1.;
2083 }
2084 else
2085 {
2086 int i;
2087 const int *b = Binom(p);
2088 real_t z = x;
2089
2090 for (i = 1; i < p; i++)
2091 {
2092 u[i] = b[i]*z;
2093 z *= x;
2094 }
2095 u[p] = z;
2096 z = y;
2097 for (i--; i > 0; i--)
2098 {
2099 u[i] *= z;
2100 z *= y;
2101 }
2102 u[0] = z;
2103 }
2104}
2105
2106void Poly_1D::CalcBinomTerms(const int p, const real_t x, const real_t y,
2107 real_t *u, real_t *d)
2108{
2109 if (p == 0)
2110 {
2111 u[0] = 1.;
2112 d[0] = 0.;
2113 }
2114 else
2115 {
2116 int i;
2117 const int *b = Binom(p);
2118 const real_t xpy = x + y, ptx = p*x;
2119 real_t z = 1.;
2120
2121 for (i = 1; i < p; i++)
2122 {
2123 d[i] = b[i]*z*(i*xpy - ptx);
2124 z *= x;
2125 u[i] = b[i]*z;
2126 }
2127 d[p] = p*z;
2128 u[p] = z*x;
2129 z = 1.;
2130 for (i--; i > 0; i--)
2131 {
2132 d[i] *= z;
2133 z *= y;
2134 u[i] *= z;
2135 }
2136 d[0] = -p*z;
2137 u[0] = z*y;
2138 }
2139}
2140
2141void Poly_1D::CalcDBinomTerms(const int p, const real_t x, const real_t y,
2142 real_t *d)
2143{
2144 if (p == 0)
2145 {
2146 d[0] = 0.;
2147 }
2148 else
2149 {
2150 int i;
2151 const int *b = Binom(p);
2152 const real_t xpy = x + y, ptx = p*x;
2153 real_t z = 1.;
2154
2155 for (i = 1; i < p; i++)
2156 {
2157 d[i] = b[i]*z*(i*xpy - ptx);
2158 z *= x;
2159 }
2160 d[p] = p*z;
2161 z = 1.;
2162 for (i--; i > 0; i--)
2163 {
2164 d[i] *= z;
2165 z *= y;
2166 }
2167 d[0] = -p*z;
2168 }
2169}
2170
2171void Poly_1D::CalcLegendre(const int p, const real_t x, real_t *u)
2172{
2173 // use the recursive definition for [-1,1]:
2174 // (n+1)*P_{n+1}(z) = (2*n+1)*z*P_n(z)-n*P_{n-1}(z)
2175 real_t z;
2176 u[0] = 1.;
2177 if (p == 0) { return; }
2178 u[1] = z = 2.*x - 1.;
2179 for (int n = 1; n < p; n++)
2180 {
2181 u[n+1] = ((2*n + 1)*z*u[n] - n*u[n-1])/(n + 1);
2182 }
2183}
2184
2185void Poly_1D::CalcLegendre(const int p, const real_t x, real_t *u, real_t *d)
2186{
2187 // use the recursive definition for [-1,1]:
2188 // (n+1)*P_{n+1}(z) = (2*n+1)*z*P_n(z)-n*P_{n-1}(z)
2189 // for the derivative use, z in [-1,1]:
2190 // P'_{n+1}(z) = (2*n+1)*P_n(z)+P'_{n-1}(z)
2191 real_t z;
2192 u[0] = 1.;
2193 d[0] = 0.;
2194 if (p == 0) { return; }
2195 u[1] = z = 2.*x - 1.;
2196 d[1] = 2.;
2197 for (int n = 1; n < p; n++)
2198 {
2199 u[n+1] = ((2*n + 1)*z*u[n] - n*u[n-1])/(n + 1);
2200 d[n+1] = (4*n + 2)*u[n] + d[n-1];
2201 }
2202}
2203
2204void Poly_1D::CalcChebyshev(const int p, const real_t x, real_t *u)
2205{
2206 // recursive definition, z in [-1,1]
2207 // T_0(z) = 1, T_1(z) = z
2208 // T_{n+1}(z) = 2*z*T_n(z) - T_{n-1}(z)
2209 real_t z;
2210 u[0] = 1.;
2211 if (p == 0) { return; }
2212 u[1] = z = 2.*x - 1.;
2213 for (int n = 1; n < p; n++)
2214 {
2215 u[n+1] = 2*z*u[n] - u[n-1];
2216 }
2217}
2218
2219void Poly_1D::CalcChebyshev(const int p, const real_t x, real_t *u, real_t *d)
2220{
2221 // recursive definition, z in [-1,1]
2222 // T_0(z) = 1, T_1(z) = z
2223 // T_{n+1}(z) = 2*z*T_n(z) - T_{n-1}(z)
2224 // T'_n(z) = n*U_{n-1}(z)
2225 // U_0(z) = 1 U_1(z) = 2*z
2226 // U_{n+1}(z) = 2*z*U_n(z) - U_{n-1}(z)
2227 // U_n(z) = z*U_{n-1}(z) + T_n(z) = z*T'_n(z)/n + T_n(z)
2228 // T'_{n+1}(z) = (n + 1)*(z*T'_n(z)/n + T_n(z))
2229 real_t z;
2230 u[0] = 1.;
2231 d[0] = 0.;
2232 if (p == 0) { return; }
2233 u[1] = z = 2.*x - 1.;
2234 d[1] = 2.;
2235 for (int n = 1; n < p; n++)
2236 {
2237 u[n+1] = 2*z*u[n] - u[n-1];
2238 d[n+1] = (n + 1)*(z*d[n]/n + 2*u[n]);
2239 }
2240}
2241
2242void Poly_1D::CalcChebyshev(const int p, const real_t x, real_t *u, real_t *d,
2243 real_t *dd)
2244{
2245 // recursive definition, z in [-1,1]
2246 // T_0(z) = 1, T_1(z) = z
2247 // T_{n+1}(z) = 2*z*T_n(z) - T_{n-1}(z)
2248 // T'_n(z) = n*U_{n-1}(z)
2249 // U_0(z) = 1 U_1(z) = 2*z
2250 // U_{n+1}(z) = 2*z*U_n(z) - U_{n-1}(z)
2251 // U_n(z) = z*U_{n-1}(z) + T_n(z) = z*T'_n(z)/n + T_n(z)
2252 // T'_{n+1}(z) = (n + 1)*(z*T'_n(z)/n + T_n(z))
2253 // T''_{n+1}(z) = (n + 1)*(2*(n + 1)*T'_n(z) + z*T''_n(z)) / n
2254 real_t z;
2255 u[0] = 1.;
2256 d[0] = 0.;
2257 dd[0]= 0.;
2258 if (p == 0) { return; }
2259 u[1] = z = 2.*x - 1.;
2260 d[1] = 2.;
2261 dd[1] = 0;
2262 for (int n = 1; n < p; n++)
2263 {
2264 u[n+1] = 2*z*u[n] - u[n-1];
2265 d[n+1] = (n + 1)*(z*d[n]/n + 2*u[n]);
2266 dd[n+1] = (n + 1)*(2.*(n + 1)*d[n] + z*dd[n])/n;
2267 }
2268}
2269
2270const real_t *Poly_1D::GetPoints(const int p, const int btype)
2271{
2273 BasisType::Check(btype);
2274 const int qtype = BasisType::GetQuadrature1D(btype);
2275 if (qtype == Quadrature1D::Invalid) { return NULL; }
2276
2278 #pragma omp critical (Poly1DGetPoints)
2279#endif
2280 {
2281 auto it = points_container.find(btype);
2282 if (it != points_container.end())
2283 {
2284 pts = it->second;
2285 }
2286 else
2287 {
2288 pts = new Array<real_t*>(h_mt);
2289 points_container[btype] = pts;
2290 }
2291 if (pts->Size() <= p)
2292 {
2293 pts->SetSize(p + 1, NULL);
2294 }
2295 if ((*pts)[p] == NULL)
2296 {
2297 (*pts)[p] = new real_t[p + 1];
2298 quad_func.GivePolyPoints(p + 1, (*pts)[p], qtype);
2299 }
2300 }
2301 return (*pts)[p];
2302}
2303
2304Poly_1D::Basis &Poly_1D::GetBasis(const int p, const int btype)
2305{
2306 Array<Basis*> *bases;
2307 BasisType::Check(btype);
2308
2310 #pragma omp critical (Poly1DGetBasis)
2311#endif
2312 {
2313 auto it = bases_container.find(btype);
2314 if (it != bases_container.end())
2315 {
2316 bases = it->second;
2317 }
2318 else
2319 {
2320 // we haven't been asked for basis or points of this type yet
2321 bases = new Array<Basis*>(h_mt);
2322 bases_container[btype] = bases;
2323 }
2324 if (bases->Size() <= p)
2325 {
2326 bases->SetSize(p + 1, NULL);
2327 }
2328 if ((*bases)[p] == NULL)
2329 {
2330 EvalType etype;
2331 if (btype == BasisType::Positive) { etype = Positive; }
2332 else if (btype == BasisType::IntegratedGLL) { etype = Integrated; }
2333 else { etype = Barycentric; }
2334 (*bases)[p] = new Basis(p, GetPoints(p, btype), etype);
2335 }
2336 }
2337 return *(*bases)[p];
2338}
2339
2341{
2342 for (PointsMap::iterator it = points_container.begin();
2343 it != points_container.end() ; ++it)
2344 {
2345 Array<real_t*>& pts = *it->second;
2346 for (int i = 0; i < pts.Size(); ++i)
2347 {
2348 delete [] pts[i];
2349 }
2350 delete it->second;
2351 }
2352
2353 for (BasisMap::iterator it = bases_container.begin();
2354 it != bases_container.end() ; ++it)
2355 {
2356 Array<Basis*>& bases = *it->second;
2357 for (int i = 0; i < bases.Size(); ++i)
2358 {
2359 delete bases[i];
2360 }
2361 delete it->second;
2362 }
2363}
2364
2365
2367 const int btype, const DofMapType dmtype)
2368 : b_type(btype),
2369 basis1d(poly1d.GetBasis(p, b_type))
2370{
2371 if (dmtype == H1_DOF_MAP || dmtype == Sr_DOF_MAP)
2372 {
2373 switch (dims)
2374 {
2375 case 1:
2376 {
2377 dof_map.SetSize(p + 1);
2378 dof_map[0] = 0;
2379 dof_map[p] = 1;
2380 for (int i = 1; i < p; i++)
2381 {
2382 dof_map[i] = i+1;
2383 }
2384 break;
2385 }
2386 case 2:
2387 {
2388 const int p1 = p + 1;
2389 dof_map.SetSize(p1*p1);
2390
2391 // vertices
2392 dof_map[0 + 0*p1] = 0;
2393 dof_map[p + 0*p1] = 1;
2394 dof_map[p + p*p1] = 2;
2395 dof_map[0 + p*p1] = 3;
2396
2397 // edges
2398 int o = 4;
2399 for (int i = 1; i < p; i++)
2400 {
2401 dof_map[i + 0*p1] = o++;
2402 }
2403 for (int i = 1; i < p; i++)
2404 {
2405 dof_map[p + i*p1] = o++;
2406 }
2407 for (int i = 1; i < p; i++)
2408 {
2409 dof_map[(p-i) + p*p1] = o++;
2410 }
2411 for (int i = 1; i < p; i++)
2412 {
2413 dof_map[0 + (p-i)*p1] = o++;
2414 }
2415
2416 // interior
2417 for (int j = 1; j < p; j++)
2418 {
2419 for (int i = 1; i < p; i++)
2420 {
2421 dof_map[i + j*p1] = o++;
2422 }
2423 }
2424 break;
2425 }
2426 case 3:
2427 {
2428 const int p1 = p + 1;
2429 dof_map.SetSize(p1*p1*p1);
2430
2431 // vertices
2432 dof_map[0 + (0 + 0*p1)*p1] = 0;
2433 dof_map[p + (0 + 0*p1)*p1] = 1;
2434 dof_map[p + (p + 0*p1)*p1] = 2;
2435 dof_map[0 + (p + 0*p1)*p1] = 3;
2436 dof_map[0 + (0 + p*p1)*p1] = 4;
2437 dof_map[p + (0 + p*p1)*p1] = 5;
2438 dof_map[p + (p + p*p1)*p1] = 6;
2439 dof_map[0 + (p + p*p1)*p1] = 7;
2440
2441 // edges (see Hexahedron::edges in mesh/hexahedron.cpp).
2442 // edges (see Constants<Geometry::CUBE>::Edges in fem/geom.cpp).
2443 int o = 8;
2444 for (int i = 1; i < p; i++)
2445 {
2446 dof_map[i + (0 + 0*p1)*p1] = o++; // (0,1)
2447 }
2448 for (int i = 1; i < p; i++)
2449 {
2450 dof_map[p + (i + 0*p1)*p1] = o++; // (1,2)
2451 }
2452 for (int i = 1; i < p; i++)
2453 {
2454 dof_map[i + (p + 0*p1)*p1] = o++; // (3,2)
2455 }
2456 for (int i = 1; i < p; i++)
2457 {
2458 dof_map[0 + (i + 0*p1)*p1] = o++; // (0,3)
2459 }
2460 for (int i = 1; i < p; i++)
2461 {
2462 dof_map[i + (0 + p*p1)*p1] = o++; // (4,5)
2463 }
2464 for (int i = 1; i < p; i++)
2465 {
2466 dof_map[p + (i + p*p1)*p1] = o++; // (5,6)
2467 }
2468 for (int i = 1; i < p; i++)
2469 {
2470 dof_map[i + (p + p*p1)*p1] = o++; // (7,6)
2471 }
2472 for (int i = 1; i < p; i++)
2473 {
2474 dof_map[0 + (i + p*p1)*p1] = o++; // (4,7)
2475 }
2476 for (int i = 1; i < p; i++)
2477 {
2478 dof_map[0 + (0 + i*p1)*p1] = o++; // (0,4)
2479 }
2480 for (int i = 1; i < p; i++)
2481 {
2482 dof_map[p + (0 + i*p1)*p1] = o++; // (1,5)
2483 }
2484 for (int i = 1; i < p; i++)
2485 {
2486 dof_map[p + (p + i*p1)*p1] = o++; // (2,6)
2487 }
2488 for (int i = 1; i < p; i++)
2489 {
2490 dof_map[0 + (p + i*p1)*p1] = o++; // (3,7)
2491 }
2492
2493 // faces (see Mesh::GenerateFaces in mesh/mesh.cpp)
2494 for (int j = 1; j < p; j++)
2495 {
2496 for (int i = 1; i < p; i++)
2497 {
2498 dof_map[i + ((p-j) + 0*p1)*p1] = o++; // (3,2,1,0)
2499 }
2500 }
2501 for (int j = 1; j < p; j++)
2502 {
2503 for (int i = 1; i < p; i++)
2504 {
2505 dof_map[i + (0 + j*p1)*p1] = o++; // (0,1,5,4)
2506 }
2507 }
2508 for (int j = 1; j < p; j++)
2509 {
2510 for (int i = 1; i < p; i++)
2511 {
2512 dof_map[p + (i + j*p1)*p1] = o++; // (1,2,6,5)
2513 }
2514 }
2515 for (int j = 1; j < p; j++)
2516 {
2517 for (int i = 1; i < p; i++)
2518 {
2519 dof_map[(p-i) + (p + j*p1)*p1] = o++; // (2,3,7,6)
2520 }
2521 }
2522 for (int j = 1; j < p; j++)
2523 {
2524 for (int i = 1; i < p; i++)
2525 {
2526 dof_map[0 + ((p-i) + j*p1)*p1] = o++; // (3,0,4,7)
2527 }
2528 }
2529 for (int j = 1; j < p; j++)
2530 {
2531 for (int i = 1; i < p; i++)
2532 {
2533 dof_map[i + (j + p*p1)*p1] = o++; // (4,5,6,7)
2534 }
2535 }
2536
2537 // interior
2538 for (int k = 1; k < p; k++)
2539 {
2540 for (int j = 1; j < p; j++)
2541 {
2542 for (int i = 1; i < p; i++)
2543 {
2544 dof_map[i + (j + k*p1)*p1] = o++;
2545 }
2546 }
2547 }
2548 break;
2549 }
2550 default:
2551 MFEM_ABORT("invalid dimension: " << dims);
2552 break;
2553 }
2554 }
2555 else if (dmtype == L2_DOF_MAP)
2556 {
2557 // leave dof_map empty, indicating that the dofs are ordered
2558 // lexicographically, i.e. the dof_map is identity
2559 }
2560 else
2561 {
2562 MFEM_ABORT("invalid DofMapType: " << dmtype);
2563 }
2564}
2565
2567 const FiniteElement &fe, const IntegrationRule &ir,
2568 DofToQuad::Mode mode, const Poly_1D::Basis &basis, bool closed,
2570{
2572 MFEM_VERIFY(mode == DofToQuad::TENSOR, "invalid mode requested");
2573
2576#endif
2577 {
2578 for (int i = 0; i < dof2quad_array.Size(); i++)
2579 {
2581 if (d2q->IntRule != &ir || d2q->mode != mode) { d2q = nullptr; }
2582 }
2583 if (!d2q)
2584 {
2586 const int ndof = closed ? fe.GetOrder() + 1 : fe.GetOrder();
2587 const int nqpt = (int)floor(pow(ir.GetNPoints(), 1.0/fe.GetDim()) + 0.5);
2588 d2q->FE = &fe;
2589 d2q->IntRule = &ir;
2590 d2q->mode = mode;
2591 d2q->ndof = ndof;
2592 d2q->nqpt = nqpt;
2593 d2q->B.SetSize(nqpt*ndof);
2594 d2q->Bt.SetSize(ndof*nqpt);
2595 d2q->G.SetSize(nqpt*ndof);
2596 d2q->Gt.SetSize(ndof*nqpt);
2598 for (int i = 0; i < nqpt; i++)
2599 {
2600 // The first 'nqpt' points in 'ir' have the same x-coordinates as those
2601 // of the 1D rule.
2603 for (int j = 0; j < ndof; j++)
2604 {
2605 d2q->B[i+nqpt*j] = d2q->Bt[j+ndof*i] = val(j);
2606 d2q->G[i+nqpt*j] = d2q->Gt[j+ndof*i] = grad(j);
2607 }
2608 }
2610 }
2611 }
2612 return *d2q;
2613}
2614
2616 const int p,
2617 const int btype,
2618 const DofMapType dmtype)
2619 : NodalFiniteElement(dims, GetTensorProductGeometry(dims), Pow(p + 1, dims),
2620 p, dims > 1 ? FunctionSpace::Qk : FunctionSpace::Pk),
2621 TensorBasisElement(dims, p, btype, dmtype)
2622{
2624}
2625
2627{
2629 // If we are using the "integrated" basis, the basis functions should be
2630 // scaled for MapType::VALUE, and not scaled for MapType::INTEGRAL. This
2631 // ensures spectral equivalence of the mass matrix with its low-order-refined
2632 // counterpart (cf. LORDiscretization)
2634 {
2636 }
2637}
2638
2640 const IntegrationRule &ir,
2642{
2644 {
2646 }
2647 else
2648 {
2650 }
2651}
2652
2654 Array<int> &face_map) const
2655{
2656 internal::GetTensorFaceMap(dim, order, face_id, face_map);
2657}
2658
2660 const int d,
2661 const int p,
2662 const int cbtype,
2663 const int obtype,
2664 const int M,
2665 const DofMapType dmtype)
2666 : VectorFiniteElement(dims, GetTensorProductGeometry(dims), d,
2667 p, M, FunctionSpace::Qk),
2668 TensorBasisElement(dims, p, VerifyNodal(VerifyClosed(cbtype)), dmtype),
2669 obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(obtype)))
2670{
2671 MFEM_VERIFY(dims > 1, "Constructor for VectorTensorFiniteElement with both "
2672 "open and closed bases is not valid for 1D elements.");
2673}
2674
2676 const int d,
2677 const int p,
2678 const int obtype,
2679 const int M,
2680 const DofMapType dmtype)
2681 : VectorFiniteElement(dims, GetTensorProductGeometry(dims), d,
2682 p, M, FunctionSpace::Pk),
2683 TensorBasisElement(dims, p, VerifyOpen(obtype), dmtype),
2684 obasis1d(poly1d.GetBasis(p, VerifyOpen(obtype)))
2685{
2686 MFEM_VERIFY(dims == 1, "Constructor for VectorTensorFiniteElement without "
2687 "closed basis is only valid for 1D elements.");
2688}
2689
2691{
2692 for (int i = 0; i < dof2quad_array_open.Size(); i++)
2693 {
2695 }
2696}
2697
2698}
int NumCols() const
Definition array.hpp:388
void SetSize(int m, int n)
Definition array.hpp:385
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:769
Get the corresponding Quadrature1D constant, when that makes sense; otherwise return Quadrature1D::In...
Definition fe_base.hpp:61
static int Check(int b_type)
If the input does not represents a valid BasisType, abort with an error; otherwise return the input.
Definition fe_base.hpp:45
@ GaussLobatto
Closed type.
Definition fe_base.hpp:32
@ Positive
Bernstein polynomials.
Definition fe_base.hpp:33
@ IntegratedGLL
Integrated GLL indicator functions.
Definition fe_base.hpp:39
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
void Factor()
Factor the current DenseMatrix, *a.
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication with the inverse of dense matrix.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
Definition densemat.cpp:220
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
Definition densemat.cpp:262
void Threshold(real_t eps)
Replace small entries, abs(a_ij) <= eps, with zero.
void SetRow(int r, const real_t *row)
void GetColumnReference(int c, Vector &col)
Definition densemat.hpp:312
real_t InnerProduct(const real_t *x, const real_t *y) const
Compute y^t A x.
Definition densemat.cpp:383
real_t * Data() const
Returns the matrix data array.
Definition densemat.hpp:111
real_t * GetData() const
Returns the matrix data array.
Definition densemat.hpp:115
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:105
void Invert()
Replaces the current matrix with its inverse.
Definition densemat.cpp:726
void InvRightScaling(const Vector &s)
InvRightScaling: this = this * diag(1./s);.
Definition densemat.cpp:442
void GetColumn(int c, Vector &col) const
real_t FNorm2() const
Compute the square of the Frobenius norm of the matrix.
Definition densemat.hpp:269
Structure representing the matrices/tensors needed to evaluate (in reference space) the values,...
Definition fe_base.hpp:137
Mode mode
Describes the contents of the B, Bt, G, and Gt arrays, see Mode.
Definition fe_base.hpp:170
const IntegrationRule * IntRule
IntegrationRule that defines the quadrature points at which the basis functions of the FE are evaluat...
Definition fe_base.hpp:146
Array< real_t > G
Definition fe_base.hpp:210
Mode
Type of data stored in the arrays B, Bt, G, and Gt.
Definition fe_base.hpp:150
@ FULL
Full multidimensional representation which does not use tensor product structure. The ordering of the...
Definition fe_base.hpp:154
@ LEXICOGRAPHIC_FULL
Full multidimensional representation which does not use tensor product structure. The ordering of the...
Definition fe_base.hpp:166
@ TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition fe_base.hpp:161
Array< real_t > B
Basis functions evaluated at quadrature points.
Definition fe_base.hpp:189
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:174
Array< real_t > Gt
Transpose of G.
Definition fe_base.hpp:217
const class FiniteElement * FE
The FiniteElement that created and owns this object.
Definition fe_base.hpp:141
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:178
Array< real_t > Bt
Transpose of B.
Definition fe_base.hpp:195
const DenseMatrix & Hessian()
Return the Hessian matrix of the transformation at the currently set IntegrationPoint,...
Definition eltrans.hpp:125
const DenseMatrix & InverseJacobian()
Return the inverse of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition eltrans.hpp:145
Return the adjugate of the Jacobian matrix of the transformation at the currently set IntegrationPoin...
Definition eltrans.hpp:135
const IntegrationPoint & GetIntPoint()
Get a const reference to the currently set integration point. This will return NULL if no integration...
Definition eltrans.hpp:98
virtual int GetSpaceDim() const =0
Get the dimension of the target (physical) space.
real_t Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint....
Definition eltrans.hpp:131
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint,...
Definition eltrans.hpp:119
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition eltrans.hpp:93
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
Abstract class for all finite elements.
Definition fe_base.hpp:239
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_base.cpp:40
int dof
Number of degrees of freedom.
Definition fe_base.hpp:248
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &Hessian) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition fe_base.cpp:101
virtual void ProjectMatrixCoefficient(MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
Given a matrix coefficient and a transformation, compute an approximation ("projection") in the local...
Definition fe_base.cpp:144
virtual ~FiniteElement()
Deconstruct the FiniteElement.
Definition fe_base.cpp:513
virtual void ProjectDiv(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &div) const
Compute the discrete divergence matrix from the given FiniteElement onto 'this' FiniteElement....
Definition fe_base.cpp:175
virtual void ProjectFromNodes(Vector &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector of values at the finite element nodes and a transformation, compute its projection (ap...
Definition fe_base.cpp:138
virtual void GetFaceDofs(int face, int **dofs, int *ndofs) const
Get the dofs associated with the given face. *dofs is set to an internal array of the local dofc on t...
Definition fe_base.cpp:96
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition fe_base.cpp:365
int GetRangeDim() const
Returns the vector dimension for vector-valued finite elements, which is also the dimension of the in...
Definition fe_base.hpp:320
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition fe_base.hpp:333
void CalcPhysDShape(ElementTransformation &Trans, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in physical space at the poi...
Definition fe_base.cpp:192
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition fe_base.cpp:65
Compute the discrete gradient matrix from the given FiniteElement onto 'this' FiniteElement....
Definition fe_base.cpp:161
IntegrationRule Nodes
Definition fe_base.hpp:251
virtual void GetFaceMap(const int face_id, Array< int > &face_map) const
Return the mapping from lexicographic face DOFs to lexicographic element DOFs for the given local fac...
Definition fe_base.cpp:507
int GetDim() const
Returns the reference space dimension for the finite element.
Definition fe_base.hpp:316
int vdim
Vector dimension of vector-valued basis functions.
Definition fe_base.hpp:242
virtual void CalcPhysHessian(ElementTransformation &Trans, DenseMatrix &Hessian) const
Evaluate the Hessian of all shape functions of a scalar finite element in reference space at the give...
Definition fe_base.cpp:288
FiniteElement(int D, Geometry::Type G, int Do, int O, int F=FunctionSpace::Pk)
Construct FiniteElement with given.
Definition fe_base.cpp:23
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition fe_base.cpp:150
int orders[Geometry::MaxDim]
Anisotropic orders.
Definition fe_base.hpp:250
virtual void GetTransferMatrix(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &I) const
Return interpolation matrix, I, which maps dofs from a coarse element, fe, to the fine dofs on this f...
Definition fe_base.cpp:119
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition fe_base.hpp:355
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
Definition fe_base.hpp:346
virtual void GetLocalRestriction(ElementTransformation &Trans, DenseMatrix &R) const
Return a local restriction matrix R (Dof x Dof) mapping fine dofs to coarse dofs.
Definition fe_base.cpp:113
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition fe_base.hpp:326
int cdim
Dimension of curl for vector-valued basis functions.
Definition fe_base.hpp:243
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
Definition fe_base.cpp:52
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
Definition fe_base.cpp:107
@ DIV
Implements CalcDivShape methods.
Definition fe_base.hpp:301
@ NONE
No derivatives implemented.
Definition fe_base.hpp:299
@ CURL
Implements CalcCurlShape methods.
Definition fe_base.hpp:302
Implements CalcDShape methods.
Definition fe_base.hpp:300
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
Definition fe_base.cpp:126
Geometry::Type geom_type
Geometry::Type of the reference element.
Definition fe_base.hpp:244
void CalcPhysDivShape(ElementTransformation &Trans, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in physical space at the po...
Definition fe_base.cpp:58
Container for all DofToQuad objects created by the FiniteElement.
Definition fe_base.hpp:258
virtual void CalcPhysLaplacian(ElementTransformation &Trans, Vector &Laplacian) const
Evaluate the Laplacian of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_base.cpp:203
DenseMatrix vshape
Definition fe_base.hpp:253
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void ProjectCurl(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Compute the discrete curl matrix from the given FiniteElement onto 'this' FiniteElement....
Definition fe_base.cpp:168
virtual void CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in physical space at the point de...
Definition fe_base.cpp:71
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition fe_base.hpp:329
int order
Order/degree of the shape functions.
Definition fe_base.hpp:249
virtual void CalcPhysLinLaplacian(ElementTransformation &Trans, Vector &Laplacian) const
Definition fe_base.cpp:244
int dim
Dimension of reference space.
Definition fe_base.hpp:241
void CalcPhysShape(ElementTransformation &Trans, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in physical space at the point ...
Definition fe_base.cpp:182
Describes the function space on each element.
Definition fe_base.hpp:222
static const int MaxDim
Definition geom.hpp:43
const IntegrationPoint & GetCenter(int GeomType) const
Return the center of the given Geometry::Type, GeomType.
Definition geom.hpp:71
static bool CheckPoint(int GeomType, const IntegrationPoint &ip)
Check if the given point is inside the given reference element.
Definition geom.cpp:435
Class for integration point with weight.
Definition intrules.hpp:35
void Get(real_t *p, const int dim) const
Definition intrules.hpp:60
void Set(const real_t *p, const int dim)
Definition intrules.hpp:46
void Set3(const real_t x1, const real_t x2, const real_t x3)
Definition intrules.hpp:79
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition intrules.hpp:256
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition intrules.hpp:259
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Base class for Matrix Coefficients that optionally depend on time and space.
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the matrix coefficient in the element described by T at the point ip, storing the result in ...
int GetWidth() const
Get the width of the matrix.
int GetHeight() const
Get the height of the matrix.
Class for standard nodal finite elements.
Definition fe_base.hpp:715
void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const override
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
Definition fe_base.cpp:796
void ProjectCurl_2D(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Definition fe_base.cpp:724
void ProjectDiv(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &div) const override
Compute the discrete divergence matrix from the given FiniteElement onto 'this' FiniteElement....
Definition fe_base.cpp:945
void GetLocalRestriction(ElementTransformation &Trans, DenseMatrix &R) const override
Return a local restriction matrix R (Dof x Dof) mapping fine dofs to coarse dofs.
Definition fe_base.cpp:764
Compute the discrete gradient matrix from the given FiniteElement onto 'this' FiniteElement....
Definition fe_base.cpp:916
void ProjectMatrixCoefficient(MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const override
Given a matrix coefficient and a transformation, compute an approximation ("projection") in the local...
Definition fe_base.cpp:835
Array< int > lex_ordering
Definition fe_base.hpp:720
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition fe_base.cpp:703
void GetFaceMap(const int face_id, Array< int > &face_map) const override
Return the mapping from lexicographic face DOFs to lexicographic element DOFs for the given local fac...
Definition fe_base.cpp:2653
NodalTensorFiniteElement(const int dims, const int p, const int btype, const DofMapType dmtype)
Definition fe_base.cpp:2615
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition fe_base.cpp:2639
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:2626
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
Class for evaluating 1D nodal, positive (Bernstein), or integrated (Gerritsma) bases.
Definition fe_base.hpp:991
void Eval(const real_t x, Vector &u) const
Evaluate the basis functions at point x in [0,1].
Definition fe_base.cpp:1761
bool IsIntegratedType() const
Returns true if the basis is "integrated", false otherwise.
Definition fe_base.hpp:1035
void ScaleIntegrated(bool scale_integrated_)
Set whether the "integrated" basis should be scaled by the subcell sizes. Has no effect for non-integ...
Definition fe_base.cpp:2018
void EvalIntegrated(const Vector &d, Vector &i) const
Evaluate the "integrated" basis type using pre-computed closed basis derivatives.
Definition fe_base.cpp:1993
Basis(const int p, const real_t *nodes, EvalType etype=Barycentric)
Create a nodal or positive (Bernstein) basis of degree p.
Definition fe_base.cpp:1697
static void CalcDBinomTerms(const int p, const real_t x, const real_t y, real_t *d)
Compute the derivatives (w.r.t. x) of the terms in the expansion of the binomial (x + y)^p assuming t...
Definition fe_base.cpp:2141
static const int * Binom(const int p)
Get a pointer to an array containing the binomial coefficients "pchoose k" for k=0,...
Definition fe_base.cpp:2028
Basis & GetBasis(const int p, const int btype)
Get a Poly_1D::Basis object of the given degree and BasisType, btype.
Definition fe_base.cpp:2304
static void CalcLegendre(const int p, const real_t x, real_t *u)
Definition fe_base.cpp:2171
static void CalcBernstein(const int p, const real_t x, real_t *u)
Compute the values of the Bernstein basis functions of order p at coordinate x and store the results ...
Definition fe_base.hpp:1164
EvalType
One-dimensional basis evaluation type.
Definition fe_base.hpp:980
@ ChangeOfBasis
Use change of basis, O(p^2) Evals.
Definition fe_base.hpp:981
@ Integrated
Integrated indicator functions (cf. Gerritsma)
Definition fe_base.hpp:984
@ Positive
Fast evaluation of Bernstein polynomials.
Definition fe_base.hpp:983
@ Barycentric
Use barycentric Lagrangian interpolation, O(p) Evals.
Definition fe_base.hpp:982
static void ChebyshevPoints(const int p, real_t *x)
Compute the points for the Chebyshev polynomials of order p and place them in the already allocated x...
Definition fe_base.cpp:2045
const real_t * GetPoints(const int p, const int btype)
Get the coordinates of the points of the given BasisType, btype.
Definition fe_base.cpp:2270
static void CalcBasis(const int p, const real_t x, real_t *u)
Evaluate the values of a hierarchical 1D basis at point x hierarchical = k-th basis function is degre...
Definition fe_base.hpp:1099
static void CalcBinomTerms(const int p, const real_t x, const real_t y, real_t *u)
Compute the p terms in the expansion of the binomial (x + y)^p and store them in the already allocate...
Definition fe_base.cpp:2077
static void GivePolyPoints(const int np, real_t *pts, const int type)
Definition intrules.cpp:717
Class for finite elements with basis functions that return scalar values.
Definition fe_base.hpp:656
void NodalLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I, const ScalarFiniteElement &fine_fe) const
Get the matrix I that defines nodal interpolation between this element and the refined element fine_f...
Definition fe_base.cpp:522
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:681
void ScalarLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I, const ScalarFiniteElement &fine_fe) const
Get matrix I "Interpolation" defined through local L2-projection in the space defined by the fine_fe.
Definition fe_base.cpp:561
void ScalarLocalL2Restriction(ElementTransformation &Trans, DenseMatrix &R, const ScalarFiniteElement &coarse_fe) const
Get restriction matrix R defined through local L2-projection in the space defined by the coarse_fe.
Definition fe_base.cpp:602
Poly_1D::Basis & basis1d
Definition fe_base.hpp:1200
Definition fe_base.cpp:2566
TensorBasisElement(const int dims, const int p, const int btype, const DofMapType dmtype)
Definition fe_base.cpp:2366
Base class for vector Coefficients that optionally depend on time and space.
int GetVDim()
Returns dimension of the vector.
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
Intermediate class for finite elements whose basis functions return vector values.
Definition fe_base.hpp:801
void ProjectGrad_RT(const real_t *nk, const Array< int > &d2n, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Definition fe_base.cpp:1206
void LocalRestriction_RT(const real_t *nk, const Array< int > &d2n, ElementTransformation &Trans, DenseMatrix &R) const
Definition fe_base.cpp:1612
void Project_ND(const real_t *tk, const Array< int > &d2t, VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Project a vector coefficient onto the ND basis functions.
Definition fe_base.cpp:1290
void Project_RT(const real_t *nk, const Array< int > &d2n, VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Project a vector coefficient onto the RT basis functions.
Definition fe_base.cpp:1066
void LocalInterpolation_RT(const VectorFiniteElement &cfe, const real_t *nk, const Array< int > &d2n, ElementTransformation &Trans, DenseMatrix &I) const
Definition fe_base.cpp:1486
void ProjectMatrixCoefficient_ND(const real_t *tk, const Array< int > &d2t, MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
Project the rows of the matrix coefficient in an ND space.
Definition fe_base.cpp:1319
void LocalRestriction_ND(const real_t *tk, const Array< int > &d2t, ElementTransformation &Trans, DenseMatrix &R) const
Definition fe_base.cpp:1655
void LocalL2Projection_RT(const VectorFiniteElement &cfe, ElementTransformation &Trans, DenseMatrix &I) const
Definition fe_base.cpp:1439
void ProjectGrad_ND(const real_t *tk, const Array< int > &d2t, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Definition fe_base.cpp:1418
void CalcVShape_ND(ElementTransformation &Trans, DenseMatrix &shape) const
Definition fe_base.cpp:1055
void CalcVShape_RT(ElementTransformation &Trans, DenseMatrix &shape) const
Definition fe_base.cpp:1043
VectorFiniteElement(int D, Geometry::Type G, int Do, int O, int M, int F=FunctionSpace::Pk)
Definition fe_base.cpp:977
void LocalInterpolation_ND(const VectorFiniteElement &cfe, const real_t *tk, const Array< int > &d2t, ElementTransformation &Trans, DenseMatrix &I) const
Definition fe_base.cpp:1573
void ProjectCurl_RT(const real_t *nk, const Array< int > &d2n, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Definition fe_base.cpp:1271
void LocalL2Projection_ND(const VectorFiniteElement &cfe, ElementTransformation &Trans, DenseMatrix &I) const
Definition fe_base.cpp:1527
void ProjectCurl_ND(const real_t *tk, const Array< int > &d2t, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Definition fe_base.cpp:1233
void ProjectMatrixCoefficient_RT(const real_t *nk, const Array< int > &d2n, MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
Project the rows of the matrix coefficient in an RT space.
Definition fe_base.cpp:1103
VectorTensorFiniteElement(const int dims, const int d, const int p, const int cbtype, const int obtype, const int M, const DofMapType dmtype)
Definition fe_base.cpp:2659
Vector data type.
Definition vector.hpp:80
void SetDataAndSize(real_t *d, int s)
Set the Vector data and size.
Definition vector.hpp:175
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition vector.hpp:227
void trans(const Vector &u, Vector &x)
Definition ex27.cpp:412
real_t b
Definition lissajous.cpp:42
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition vector.hpp:45
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
void mfem_error(const char *msg)
Definition error.cpp:154
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition table.cpp:476
Geometry Geometries
Definition fe.cpp:49
void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
Multiply a matrix A with the transpose of a matrix B: A*Bt.
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
void AddMult_a_VWt(const real_t a, const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += a * v w^t.
void AddMult_a_VVt(const real_t a, const Vector &v, DenseMatrix &VVt)
VVt += a * v v^t.
void AddMult_a_AAt(real_t a, const DenseMatrix &A, DenseMatrix &AAt)
AAt += a * A * A^t.
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
float real_t
Definition config.hpp:43
Poly_1D poly1d
Definition fe.cpp:28
void InvertLinearTrans(ElementTransformation &trans, const IntegrationPoint &pt, Vector &x)
Definition fe_base.cpp:748