MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
fe_base.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// 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; }
35#ifndef MFEM_THREAD_SAFE
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 {
78#ifdef MFEM_THREAD_SAFE
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{
98 MFEM_ABORT("method is not overloaded");
99}
100
102 DenseMatrix &h) const
103{
104 MFEM_ABORT("method is not overloaded");
105}
106
108 DenseMatrix &I) const
109{
110 MFEM_ABORT("method is not overloaded");
111}
112
114 DenseMatrix &) const
115{
116 MFEM_ABORT("method is not overloaded");
117}
118
121 DenseMatrix &I) const
122{
123 MFEM_ABORT("method is not overloaded");
124}
125
127 Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
128{
129 MFEM_ABORT("method is not overloaded");
130}
131
133 VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
134{
135 MFEM_ABORT("method is not overloaded");
136}
137
139 Vector &dofs) const
140{
141 mfem_error("FiniteElement::ProjectFromNodes() (vector) is not overloaded!");
142}
143
146{
147 MFEM_ABORT("method is not overloaded");
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,
163 DenseMatrix &grad) const
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, "");
196#ifdef MFEM_THREAD_SAFE
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 {
330 DenseMatrix grad(dof, dim);
331 CalcPhysDShape(Trans, grad);
332 DenseMatrix gmap(dof, size);
333 Mult(grad,Trans.Hessian(),gmap);
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
366 DofToQuad::Mode mode) const
367{
368 DofToQuad *d2q = nullptr;
369 MFEM_VERIFY(mode == DofToQuad::FULL, "invalid mode requested");
370
371#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
372 #pragma omp critical (DofToQuad)
373#endif
374 {
375 for (int i = 0; i < dof2quad_array.Size(); i++)
376 {
377 d2q = dof2quad_array[i];
378 if (d2q->IntRule != &ir || d2q->mode != mode) { d2q = nullptr; }
379 }
380 if (!d2q)
381 {
382#ifdef MFEM_THREAD_SAFE
384#endif
385 d2q = new DofToQuad;
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 {
438 case GRAD:
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 }
501 dof2quad_array.Append(d2q);
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 {
517 delete dof2quad_array[i];
518 }
519}
520
521
524 const ScalarFiniteElement &fine_fe) const
525{
527 Vector vv(v, dim);
528 IntegrationPoint f_ip;
530#ifdef MFEM_THREAD_SAFE
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(), "");
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
587 AddMult_a_VVt(ip.weight, fine_shape, fine_mass);
588 AddMult_a_VWt(ip.weight, fine_shape, coarse_shape, fine_coarse_mass);
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);
622 AddMult_a_VVt(c_ip.weight, coarse_shape, coarse_mass);
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);
636 AddMult_a_VWt(f_ip.weight*Trans.Weight(), coarse_shape, fine_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.
655 auto &d2q = GetDofToQuad(ir, DofToQuad::FULL);
656 //Undo the native ordering which is what FiniteElement::GetDofToQuad returns.
657 auto *d2q_new = new DofToQuad(d2q);
658 d2q_new->mode = DofToQuad::LEXICOGRAPHIC_FULL;
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 {
680 case GRAD: return dim;
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
700 dof2quad_array.Append(d2q_new);
701}
702
704 DofToQuad::Mode mode) const
705{
706 //Should make this loop a function of FiniteElement
707 for (int i = 0; i < dof2quad_array.Size(); i++)
708 {
709 const DofToQuad &d2q = *dof2quad_array[i];
710 if (d2q.IntRule == &ir && d2q.mode == mode) { return d2q; }
711 }
712
714 {
715 return FiniteElement::GetDofToQuad(ir, mode);
716 }
717 else
718 {
719 CreateLexicographicFullMap(ir);
720 return NodalFiniteElement::GetDofToQuad(ir, mode);
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
770#ifdef MFEM_THREAD_SAFE
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,
918 DenseMatrix &grad) const
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
925 grad.SetSize(dim*dof, fe.GetDof());
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);
932 Mult(dshape, Jinv, grad_k);
933 if (map_type == INTEGRAL)
934 {
935 grad_k *= Trans.Weight();
936 }
937 for (int j = 0; j < grad_k.Height(); j++)
938 for (int d = 0; d < dim; d++)
939 {
940 grad(k+d*dof,j) = grad_k(j,d);
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
977 Vector &dofs) const
978{
979 MFEM_ASSERT(lex_ordering.Size() == dof, "Permutation is not defined by FE.");
980 MFEM_ASSERT(dofs.Size() == ncomp * dof, "Wrong input size.");
981
982 Vector dofs_native(ncomp * dof);
983 for (int i = 0; i < dof; i++)
984 {
985 for (int c = 0; c < ncomp; c++)
986 {
987 dofs_native(c*dof + lex_ordering[i]) = dofs(c*dof + i);
988 }
989 }
990 dofs = dofs_native;
991}
992
994 int Do, int O, int M, int F)
995 : FiniteElement(D, G, Do, O, F)
996{
998 map_type = M;
1000 is_nodal = true;
1001 vdim = dim;
1002 if (map_type == H_CURL)
1003 {
1004 cdim = (dim == 3) ? 3 : 1;
1005 }
1006}
1007
1008void VectorFiniteElement::CalcShape(
1009 const IntegrationPoint &ip, Vector &shape) const
1010{
1011 mfem_error("Error: Cannot use scalar CalcShape(...) function with\n"
1012 " VectorFiniteElements!");
1013}
1014
1015void VectorFiniteElement::CalcDShape(
1016 const IntegrationPoint &ip, DenseMatrix &dshape) const
1017{
1018 mfem_error("Error: Cannot use scalar CalcDShape(...) function with\n"
1019 " VectorFiniteElements!");
1020}
1021
1023{
1024 switch (map_type)
1025 {
1026 case H_DIV:
1027 deriv_type = DIV;
1030 break;
1031 case H_CURL:
1032 switch (dim)
1033 {
1034 case 3: // curl: 3D H_CURL -> 3D H_DIV
1035 deriv_type = CURL;
1038 break;
1039 case 2:
1040 // curl: 2D H_CURL -> INTEGRAL
1041 deriv_type = CURL;
1044 break;
1045 case 1:
1046 deriv_type = NONE;
1049 break;
1050 default:
1051 MFEM_ABORT("Invalid dimension, Dim = " << dim);
1052 }
1053 break;
1054 default:
1055 MFEM_ABORT("Invalid MapType = " << map_type);
1056 }
1057}
1058
1060 ElementTransformation &Trans, DenseMatrix &shape) const
1061{
1062 MFEM_ASSERT(map_type == H_DIV, "");
1063#ifdef MFEM_THREAD_SAFE
1065#endif
1066 CalcVShape(Trans.GetIntPoint(), vshape);
1067 MultABt(vshape, Trans.Jacobian(), shape);
1068 shape *= (1.0 / Trans.Weight());
1069}
1070
1072 ElementTransformation &Trans, DenseMatrix &shape) const
1073{
1074 MFEM_ASSERT(map_type == H_CURL, "");
1075#ifdef MFEM_THREAD_SAFE
1077#endif
1078 CalcVShape(Trans.GetIntPoint(), vshape);
1079 Mult(vshape, Trans.InverseJacobian(), shape);
1080}
1081
1083 const real_t *nk, const Array<int> &d2n,
1084 VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
1085{
1087 const int sdim = Trans.GetSpaceDim();
1088 MFEM_ASSERT(vc.GetVDim() == sdim, "");
1089 Vector xk(vk, sdim);
1090 const bool square_J = (dim == sdim);
1091
1092 for (int k = 0; k < dof; k++)
1093 {
1094 Trans.SetIntPoint(&Nodes.IntPoint(k));
1095 vc.Eval(xk, Trans, Nodes.IntPoint(k));
1096 // dof_k = nk^t adj(J) xk
1097 dofs(k) = Trans.AdjugateJacobian().InnerProduct(vk, nk + d2n[k]*dim);
1098 if (!square_J) { dofs(k) /= Trans.Weight(); }
1099 }
1100}
1101
1103 const real_t *nk, const Array<int> &d2n,
1104 Vector &vc, ElementTransformation &Trans, Vector &dofs) const
1105{
1106 const int sdim = Trans.GetSpaceDim();
1107 const bool square_J = (dim == sdim);
1108
1109 for (int k = 0; k < dof; k++)
1110 {
1111 Trans.SetIntPoint(&Nodes.IntPoint(k));
1112 // dof_k = nk^t adj(J) xk
1113 dofs(k) = Trans.AdjugateJacobian().InnerProduct(
1114 &vc[k*sdim], nk + d2n[k]*dim);
1115 if (!square_J) { dofs(k) /= Trans.Weight(); }
1116 }
1117}
1118
1120 const real_t *nk, const Array<int> &d2n,
1121 MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
1122{
1123 // project the rows of the matrix coefficient in an RT space
1124
1125 const int sdim = T.GetSpaceDim();
1126 MFEM_ASSERT(mc.GetWidth() == sdim, "");
1127 const bool square_J = (dim == sdim);
1128 DenseMatrix MQ(mc.GetHeight(), mc.GetWidth());
1129 Vector nk_phys(sdim), dofs_k(MQ.Height());
1130 MFEM_ASSERT(dofs.Size() == dof*MQ.Height(), "");
1131
1132 for (int k = 0; k < dof; k++)
1133 {
1134 T.SetIntPoint(&Nodes.IntPoint(k));
1135 mc.Eval(MQ, T, Nodes.IntPoint(k));
1136 // nk_phys = adj(J)^t nk
1137 T.AdjugateJacobian().MultTranspose(nk + d2n[k]*dim, nk_phys);
1138 if (!square_J) { nk_phys /= T.Weight(); }
1139 MQ.Mult(nk_phys, dofs_k);
1140 for (int r = 0; r < MQ.Height(); r++)
1141 {
1142 dofs(k+dof*r) = dofs_k(r);
1143 }
1144 }
1145}
1146
1148 const real_t *nk, const Array<int> &d2n, const FiniteElement &fe,
1149 ElementTransformation &Trans, DenseMatrix &I) const
1150{
1151 if (fe.GetRangeType() == SCALAR)
1152 {
1154 Vector shape(fe.GetDof());
1155 int sdim = Trans.GetSpaceDim();
1156
1157 I.SetSize(dof, sdim*fe.GetDof());
1158 for (int k = 0; k < dof; k++)
1159 {
1160 const IntegrationPoint &ip = Nodes.IntPoint(k);
1161
1162 fe.CalcShape(ip, shape);
1163 Trans.SetIntPoint(&ip);
1164 // Transform RT face normals from reference to physical space
1165 // vk = adj(J)^T nk
1166 Trans.AdjugateJacobian().MultTranspose(nk + d2n[k]*dim, vk);
1167 if (fe.GetMapType() == INTEGRAL)
1168 {
1169 real_t w = 1.0/Trans.Weight();
1170 for (int d = 0; d < dim; d++)
1171 {
1172 vk[d] *= w;
1173 }
1174 }
1175
1176 for (int j = 0; j < shape.Size(); j++)
1177 {
1178 real_t s = shape(j);
1179 if (fabs(s) < 1e-12)
1180 {
1181 s = 0.0;
1182 }
1183 // Project scalar basis function multiplied by each coordinate
1184 // direction onto the transformed face normals
1185 for (int d = 0; d < sdim; d++)
1186 {
1187 I(k,j+d*shape.Size()) = s*vk[d];
1188 }
1189 }
1190 }
1191 }
1192 else
1193 {
1194 int sdim = Trans.GetSpaceDim();
1196 DenseMatrix vshape(fe.GetDof(), sdim);
1197 Vector vshapenk(fe.GetDof());
1198 const bool square_J = (dim == sdim);
1199
1200 I.SetSize(dof, fe.GetDof());
1201 for (int k = 0; k < dof; k++)
1202 {
1203 const IntegrationPoint &ip = Nodes.IntPoint(k);
1204
1205 Trans.SetIntPoint(&ip);
1206 // Transform RT face normals from reference to physical space
1207 // vk = adj(J)^T nk
1208 Trans.AdjugateJacobian().MultTranspose(nk + d2n[k]*dim, vk);
1209 // Compute fe basis functions in physical space
1210 fe.CalcVShape(Trans, vshape);
1211 // Project fe basis functions onto transformed face normals
1212 vshape.Mult(vk, vshapenk);
1213 if (!square_J) { vshapenk /= Trans.Weight(); }
1214 for (int j=0; j<vshapenk.Size(); j++)
1215 {
1216 I(k,j) = vshapenk(j);
1217 }
1218 }
1219 }
1220}
1221
1223 const real_t *nk, const Array<int> &d2n, const FiniteElement &fe,
1224 ElementTransformation &Trans, DenseMatrix &grad) const
1225{
1226 if (dim != 2)
1227 {
1228 mfem_error("VectorFiniteElement::ProjectGrad_RT works only in 2D!");
1229 }
1230
1231 DenseMatrix dshape(fe.GetDof(), fe.GetDim());
1232 Vector grad_k(fe.GetDof());
1233 real_t tk[2];
1234
1235 grad.SetSize(dof, fe.GetDof());
1236 for (int k = 0; k < dof; k++)
1237 {
1238 fe.CalcDShape(Nodes.IntPoint(k), dshape);
1239 tk[0] = nk[d2n[k]*dim+1];
1240 tk[1] = -nk[d2n[k]*dim];
1241 dshape.Mult(tk, grad_k);
1242 for (int j = 0; j < grad_k.Size(); j++)
1243 {
1244 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
1245 }
1246 }
1247}
1248
1250 const real_t *tk, const Array<int> &d2t, const FiniteElement &fe,
1251 ElementTransformation &Trans, DenseMatrix &curl) const
1252{
1253#ifdef MFEM_THREAD_SAFE
1257#else
1258 curlshape.SetSize(fe.GetDof(), dim);
1260 JtJ.SetSize(dim, dim);
1261#endif
1262
1263 Vector curl_k(fe.GetDof());
1264
1265 curl.SetSize(dof, fe.GetDof());
1266 for (int k = 0; k < dof; k++)
1267 {
1268 const IntegrationPoint &ip = Nodes.IntPoint(k);
1269
1270 // calculate J^t * J / |J|
1271 Trans.SetIntPoint(&ip);
1272 MultAtB(Trans.Jacobian(), Trans.Jacobian(), JtJ);
1273 JtJ *= 1.0 / Trans.Weight();
1274
1275 // transform curl of shapes (rows) by J^t * J / |J|
1276 fe.CalcCurlShape(ip, curlshape);
1278
1279 curlshape_J.Mult(tk + d2t[k]*dim, curl_k);
1280 for (int j = 0; j < curl_k.Size(); j++)
1281 {
1282 curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
1283 }
1284 }
1285}
1286
1288 const real_t *nk, const Array<int> &d2n, const FiniteElement &fe,
1289 ElementTransformation &Trans, DenseMatrix &curl) const
1290{
1291 DenseMatrix curl_shape(fe.GetDof(), dim);
1292 Vector curl_k(fe.GetDof());
1293
1294 curl.SetSize(dof, fe.GetDof());
1295 for (int k = 0; k < dof; k++)
1296 {
1297 fe.CalcCurlShape(Nodes.IntPoint(k), curl_shape);
1298 curl_shape.Mult(nk + d2n[k]*dim, curl_k);
1299 for (int j = 0; j < curl_k.Size(); j++)
1300 {
1301 curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
1302 }
1303 }
1304}
1305
1307 const real_t *tk, const Array<int> &d2t,
1308 VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
1309{
1311 Vector xk(vk, vc.GetVDim());
1312
1313 for (int k = 0; k < dof; k++)
1314 {
1315 Trans.SetIntPoint(&Nodes.IntPoint(k));
1316
1317 vc.Eval(xk, Trans, Nodes.IntPoint(k));
1318 // dof_k = xk^t J tk
1319 dofs(k) = Trans.Jacobian().InnerProduct(tk + d2t[k]*dim, vk);
1320 }
1321}
1322
1324 const real_t *tk, const Array<int> &d2t,
1325 Vector &vc, ElementTransformation &Trans, Vector &dofs) const
1326{
1327 for (int k = 0; k < dof; k++)
1328 {
1329 Trans.SetIntPoint(&Nodes.IntPoint(k));
1330 // dof_k = xk^t J tk
1331 dofs(k) = Trans.Jacobian().InnerProduct(tk + d2t[k]*dim, &vc[k*dim]);
1332 }
1333}
1334
1336 const real_t *tk, const Array<int> &d2t,
1337 MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
1338{
1339 // project the rows of the matrix coefficient in an ND space
1340
1341 const int sdim = T.GetSpaceDim();
1342 MFEM_ASSERT(mc.GetWidth() == sdim, "");
1343 DenseMatrix MQ(mc.GetHeight(), mc.GetWidth());
1344 Vector tk_phys(sdim), dofs_k(MQ.Height());
1345 MFEM_ASSERT(dofs.Size() == dof*MQ.Height(), "");
1346
1347 for (int k = 0; k < dof; k++)
1348 {
1349 T.SetIntPoint(&Nodes.IntPoint(k));
1350 mc.Eval(MQ, T, Nodes.IntPoint(k));
1351 // tk_phys = J tk
1352 T.Jacobian().Mult(tk + d2t[k]*dim, tk_phys);
1353 MQ.Mult(tk_phys, dofs_k);
1354 for (int r = 0; r < MQ.Height(); r++)
1355 {
1356 dofs(k+dof*r) = dofs_k(r);
1357 }
1358 }
1359}
1360
1362 const real_t *tk, const Array<int> &d2t, const FiniteElement &fe,
1363 ElementTransformation &Trans, DenseMatrix &I) const
1364{
1365 if (fe.GetRangeType() == SCALAR)
1366 {
1367 int sdim = Trans.GetSpaceDim();
1369 Vector shape(fe.GetDof());
1370
1371 I.SetSize(dof, sdim*fe.GetDof());
1372 for (int k = 0; k < dof; k++)
1373 {
1374 const IntegrationPoint &ip = Nodes.IntPoint(k);
1375
1376 fe.CalcShape(ip, shape);
1377 Trans.SetIntPoint(&ip);
1378 // Transform ND edge tengents from reference to physical space
1379 // vk = J tk
1380 Trans.Jacobian().Mult(tk + d2t[k]*dim, vk);
1381 if (fe.GetMapType() == INTEGRAL)
1382 {
1383 real_t w = 1.0/Trans.Weight();
1384 for (int d = 0; d < sdim; d++)
1385 {
1386 vk[d] *= w;
1387 }
1388 }
1389
1390 for (int j = 0; j < shape.Size(); j++)
1391 {
1392 real_t s = shape(j);
1393 if (fabs(s) < 1e-12)
1394 {
1395 s = 0.0;
1396 }
1397 // Project scalar basis function multiplied by each coordinate
1398 // direction onto the transformed edge tangents
1399 for (int d = 0; d < sdim; d++)
1400 {
1401 I(k, j + d*shape.Size()) = s*vk[d];
1402 }
1403 }
1404 }
1405 }
1406 else
1407 {
1408 int sdim = Trans.GetSpaceDim();
1410 DenseMatrix vshape(fe.GetDof(), sdim);
1411 Vector vshapetk(fe.GetDof());
1412
1413 I.SetSize(dof, fe.GetDof());
1414 for (int k = 0; k < dof; k++)
1415 {
1416 const IntegrationPoint &ip = Nodes.IntPoint(k);
1417
1418 Trans.SetIntPoint(&ip);
1419 // Transform ND edge tangents from reference to physical space
1420 // vk = J tk
1421 Trans.Jacobian().Mult(tk + d2t[k]*dim, vk);
1422 // Compute fe basis functions in physical space
1423 fe.CalcVShape(Trans, vshape);
1424 // Project fe basis functions onto transformed edge tangents
1425 vshape.Mult(vk, vshapetk);
1426 for (int j=0; j<vshapetk.Size(); j++)
1427 {
1428 I(k, j) = vshapetk(j);
1429 }
1430 }
1431 }
1432}
1433
1435 const real_t *tk, const Array<int> &d2t, const FiniteElement &fe,
1436 ElementTransformation &Trans, DenseMatrix &grad) const
1437{
1438 MFEM_ASSERT(fe.GetMapType() == VALUE, "");
1439
1440 DenseMatrix dshape(fe.GetDof(), fe.GetDim());
1441 Vector grad_k(fe.GetDof());
1442
1443 grad.SetSize(dof, fe.GetDof());
1444 for (int k = 0; k < dof; k++)
1445 {
1446 fe.CalcDShape(Nodes.IntPoint(k), dshape);
1447 dshape.Mult(tk + d2t[k]*dim, grad_k);
1448 for (int j = 0; j < grad_k.Size(); j++)
1449 {
1450 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
1451 }
1452 }
1453}
1454
1456 const VectorFiniteElement &cfe, ElementTransformation &Trans,
1457 DenseMatrix &I) const
1458{
1459 Vector v(dim);
1460 IntegrationPoint tr_ip;
1461
1462 const int fs = dof, cs = cfe.GetDof();
1463 I.SetSize(fs, cs);
1464 DenseMatrix fine_shape(fs, dim), coarse_shape(cs, cfe.GetDim());
1465 DenseMatrix fine_mass(fs), fine_coarse_mass(fs, cs); // initialized with 0
1466 const int ir_order =
1467 std::max(GetOrder(), this->GetOrder()) + this->GetOrder();
1468 const IntegrationRule &ir = IntRules.Get(this->GetGeomType(), ir_order);
1469
1471 const DenseMatrix &adjJ = Trans.AdjugateJacobian();
1472 for (int i = 0; i < ir.GetNPoints(); i++)
1473 {
1474 const IntegrationPoint &ip = ir.IntPoint(i);
1475 real_t w = ip.weight;
1476 this->CalcVShape(ip, fine_shape);
1477 Trans.Transform(ip, v);
1478 tr_ip.Set(v.GetData(), dim);
1479 cfe.CalcVShape(tr_ip, coarse_shape);
1480
1481 AddMult_a_AAt(w, fine_shape, fine_mass);
1482 for (int k=0; k<fs; ++k)
1483 {
1484 for (int j=0; j<cs; ++j)
1485 {
1486 real_t Mkj = 0.0;
1487 for (int d1=0; d1<dim; ++d1)
1488 {
1489 for (int d2=0; d2<dim; ++d2)
1490 {
1491 Mkj += w*fine_shape(k,d1)*adjJ(d2,d1)*coarse_shape(j,d2);
1492 }
1493 }
1494 fine_coarse_mass(k,j) += (fabs(Mkj) < 1e-12) ? 0.0 : Mkj;
1495 }
1496 }
1497 }
1498 DenseMatrixInverse fine_mass_inv(fine_mass);
1499 fine_mass_inv.Mult(fine_coarse_mass, I);
1500}
1501
1503 const VectorFiniteElement &cfe, const real_t *nk, const Array<int> &d2n,
1504 ElementTransformation &Trans, DenseMatrix &I) const
1505{
1506 MFEM_ASSERT(map_type == cfe.GetMapType(), "");
1507
1508 if (!is_nodal) { return LocalL2Projection_RT(cfe, Trans, I); }
1509
1511 Vector xk(vk, dim);
1513#ifdef MFEM_THREAD_SAFE
1514 DenseMatrix vshape(cfe.GetDof(), cfe.GetDim());
1515#else
1516 DenseMatrix vshape(cfe.vshape.Data(), cfe.GetDof(), cfe.GetDim());
1517#endif
1518 I.SetSize(dof, vshape.Height());
1519
1520 // assuming Trans is linear; this should be ok for all refinement types
1522 const DenseMatrix &adjJ = Trans.AdjugateJacobian();
1523 for (int k = 0; k < dof; k++)
1524 {
1525 Trans.Transform(Nodes.IntPoint(k), xk);
1526 ip.Set3(vk);
1527 cfe.CalcVShape(ip, vshape);
1528 // xk = |J| J^{-t} n_k
1529 adjJ.MultTranspose(nk + d2n[k]*dim, vk);
1530 // I_k = vshape_k.adj(J)^t.n_k, k=1,...,dof
1531 for (int j = 0; j < vshape.Height(); j++)
1532 {
1533 real_t Ikj = 0.;
1534 for (int i = 0; i < dim; i++)
1535 {
1536 Ikj += vshape(j, i) * vk[i];
1537 }
1538 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
1539 }
1540 }
1541}
1542
1544 const VectorFiniteElement &cfe,
1545 ElementTransformation &Trans, DenseMatrix &I) const
1546{
1547 Vector v(dim);
1548 IntegrationPoint tr_ip;
1549
1550 const int fs = dof, cs = cfe.GetDof();
1551 I.SetSize(fs, cs);
1552 DenseMatrix fine_shape(fs, dim), coarse_shape(cs, cfe.GetDim());
1553 DenseMatrix fine_mass(fs), fine_coarse_mass(fs, cs); // initialized with 0
1554 const int ir_order =
1555 std::max(GetOrder(), this->GetOrder()) + this->GetOrder();
1556 const IntegrationRule &ir = IntRules.Get(this->GetGeomType(), ir_order);
1557
1559 const DenseMatrix &J = Trans.Jacobian();
1560 for (int i = 0; i < ir.GetNPoints(); i++)
1561 {
1562 const IntegrationPoint &ip = ir.IntPoint(i);
1563 this->CalcVShape(ip, fine_shape);
1564 Trans.Transform(ip, v);
1565 tr_ip.Set(v.GetData(), dim);
1566 cfe.CalcVShape(tr_ip, coarse_shape);
1567
1568 AddMult_a_AAt(ip.weight, fine_shape, fine_mass);
1569 for (int k=0; k<fs; ++k)
1570 {
1571 for (int j=0; j<cs; ++j)
1572 {
1573 real_t Mkj = 0.0;
1574 for (int d1=0; d1<dim; ++d1)
1575 {
1576 for (int d2=0; d2<dim; ++d2)
1577 {
1578 Mkj += ip.weight*fine_shape(k,d1)*J(d1,d2)*coarse_shape(j,d2);
1579 }
1580 }
1581 fine_coarse_mass(k,j) += (fabs(Mkj) < 1e-12) ? 0.0 : Mkj;
1582 }
1583 }
1584 }
1585 DenseMatrixInverse fine_mass_inv(fine_mass);
1586 fine_mass_inv.Mult(fine_coarse_mass, I);
1587}
1588
1590 const VectorFiniteElement &cfe, const real_t *tk, const Array<int> &d2t,
1591 ElementTransformation &Trans, DenseMatrix &I) const
1592{
1593 if (!is_nodal) { return LocalL2Projection_ND(cfe, Trans, I); }
1594
1596 Vector xk(vk, dim);
1598#ifdef MFEM_THREAD_SAFE
1599 DenseMatrix vshape(cfe.GetDof(), cfe.GetDim());
1600#else
1601 DenseMatrix vshape(cfe.vshape.Data(), cfe.GetDof(), cfe.GetDim());
1602#endif
1603 I.SetSize(dof, vshape.Height());
1604
1605 // assuming Trans is linear; this should be ok for all refinement types
1607 const DenseMatrix &J = Trans.Jacobian();
1608 for (int k = 0; k < dof; k++)
1609 {
1610 Trans.Transform(Nodes.IntPoint(k), xk);
1611 ip.Set3(vk);
1612 cfe.CalcVShape(ip, vshape);
1613 // xk = J t_k
1614 J.Mult(tk + d2t[k]*dim, vk);
1615 // I_k = vshape_k.J.t_k, k=1,...,Dof
1616 for (int j = 0; j < vshape.Height(); j++)
1617 {
1618 real_t Ikj = 0.;
1619 for (int i = 0; i < dim; i++)
1620 {
1621 Ikj += vshape(j, i) * vk[i];
1622 }
1623 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
1624 }
1625 }
1626}
1627
1629 const real_t *nk, const Array<int> &d2n, ElementTransformation &Trans,
1630 DenseMatrix &R) const
1631{
1632 real_t pt_data[Geometry::MaxDim];
1634 Vector pt(pt_data, dim);
1635
1636#ifdef MFEM_THREAD_SAFE
1638#endif
1639
1641 const DenseMatrix &J = Trans.Jacobian();
1642 const real_t weight = Trans.Weight();
1643 for (int j = 0; j < dof; j++)
1644 {
1645 InvertLinearTrans(Trans, Nodes.IntPoint(j), pt);
1646 ip.Set(pt_data, dim);
1647 if (Geometries.CheckPoint(geom_type, ip)) // do we need an epsilon here?
1648 {
1649 CalcVShape(ip, vshape);
1650 J.MultTranspose(nk+dim*d2n[j], pt_data);
1651 pt /= weight;
1652 for (int k = 0; k < dof; k++)
1653 {
1654 real_t R_jk = 0.0;
1655 for (int d = 0; d < dim; d++)
1656 {
1657 R_jk += vshape(k,d)*pt_data[d];
1658 }
1659 R(j,k) = R_jk;
1660 }
1661 }
1662 else
1663 {
1664 // Set the whole row to avoid valgrind warnings in R.Threshold().
1665 R.SetRow(j, infinity());
1666 }
1667 }
1668 R.Threshold(1e-12);
1669}
1670
1672 const real_t *tk, const Array<int> &d2t, ElementTransformation &Trans,
1673 DenseMatrix &R) const
1674{
1675 real_t pt_data[Geometry::MaxDim];
1677 Vector pt(pt_data, dim);
1678
1679#ifdef MFEM_THREAD_SAFE
1681#endif
1682
1684 const DenseMatrix &Jinv = Trans.InverseJacobian();
1685 for (int j = 0; j < dof; j++)
1686 {
1687 InvertLinearTrans(Trans, Nodes.IntPoint(j), pt);
1688 ip.Set(pt_data, dim);
1689 if (Geometries.CheckPoint(geom_type, ip)) // do we need an epsilon here?
1690 {
1691 CalcVShape(ip, vshape);
1692 Jinv.Mult(tk+dim*d2t[j], pt_data);
1693 for (int k = 0; k < dof; k++)
1694 {
1695 real_t R_jk = 0.0;
1696 for (int d = 0; d < dim; d++)
1697 {
1698 R_jk += vshape(k,d)*pt_data[d];
1699 }
1700 R(j,k) = R_jk;
1701 }
1702 }
1703 else
1704 {
1705 // Set the whole row to avoid valgrind warnings in R.Threshold().
1706 R.SetRow(j, infinity());
1707 }
1708 }
1709 R.Threshold(1e-12);
1710}
1711
1712
1714 : etype(etype), auxiliary_basis(NULL), scale_integrated(false)
1715{
1716 switch (etype)
1717 {
1718 case ChangeOfBasis:
1719 {
1720 x.SetSize(p + 1);
1721 w.SetSize(p + 1);
1722 DenseMatrix A(p + 1);
1723 for (int i = 0; i <= p; i++)
1724 {
1725 CalcBasis(p, nodes[i], A.GetColumn(i));
1726 }
1727 Ai.Factor(A);
1728 // mfem::out << "Poly_1D::Basis(" << p << ",...) : "; Ai.TestInversion();
1729 break;
1730 }
1731 case Barycentric:
1732 {
1733 x.SetSize(p + 1);
1734 w.SetSize(p + 1);
1735 x = nodes;
1736 w = 1.0;
1737 for (int i = 0; i <= p; i++)
1738 {
1739 for (int j = 0; j < i; j++)
1740 {
1741 real_t xij = x(i) - x(j);
1742 w(i) *= xij;
1743 w(j) *= -xij;
1744 }
1745 }
1746 for (int i = 0; i <= p; i++)
1747 {
1748 w(i) = 1.0/w(i);
1749 }
1750
1751#ifdef MFEM_DEBUG
1752 // Make sure the nodes are increasing
1753 for (int i = 0; i < p; i++)
1754 {
1755 if (x(i) >= x(i+1))
1756 {
1757 mfem_error("Poly_1D::Basis::Basis : nodes are not increasing!");
1758 }
1759 }
1760#endif
1761 break;
1762 }
1763 case Positive:
1764 x.SetDataAndSize(NULL, p + 1); // use x to store (p + 1)
1765 break;
1766 case Integrated:
1767 auxiliary_basis = new Basis(
1769 u_aux.SetSize(p+2);
1770 d_aux.SetSize(p+2);
1771 d2_aux.SetSize(p+2);
1772 break;
1773 default: break;
1774 }
1775}
1776
1778{
1779 switch (etype)
1780 {
1781 case ChangeOfBasis:
1782 {
1783 CalcBasis(Ai.Width() - 1, y, x);
1784 Ai.Mult(x, u);
1785 break;
1786 }
1787 case Barycentric:
1788 {
1789 int i, k, p = x.Size() - 1;
1790 real_t l, lk;
1791
1792 if (p == 0)
1793 {
1794 u(0) = 1.0;
1795 return;
1796 }
1797
1798 lk = 1.0;
1799 for (k = 0; k < p; k++)
1800 {
1801 if (y >= (x(k) + x(k+1))/2)
1802 {
1803 lk *= y - x(k);
1804 }
1805 else
1806 {
1807 for (i = k+1; i <= p; i++)
1808 {
1809 lk *= y - x(i);
1810 }
1811 break;
1812 }
1813 }
1814 l = lk * (y - x(k));
1815
1816 for (i = 0; i < k; i++)
1817 {
1818 u(i) = l * w(i) / (y - x(i));
1819 }
1820 u(k) = lk * w(k);
1821 for (i++; i <= p; i++)
1822 {
1823 u(i) = l * w(i) / (y - x(i));
1824 }
1825 break;
1826 }
1827 case Positive:
1828 CalcBernstein(x.Size() - 1, y, u);
1829 break;
1830 case Integrated:
1831 auxiliary_basis->Eval(y, u_aux, d_aux);
1832 EvalIntegrated(d_aux, u);
1833 break;
1834 default: break;
1835 }
1836}
1837
1838void Poly_1D::Basis::Eval(const real_t y, Vector &u, Vector &d) const
1839{
1840 switch (etype)
1841 {
1842 case ChangeOfBasis:
1843 {
1844 CalcBasis(Ai.Width() - 1, y, x, w);
1845 Ai.Mult(x, u);
1846 Ai.Mult(w, d);
1847 break;
1848 }
1849 case Barycentric:
1850 {
1851 int i, k, p = x.Size() - 1;
1852 real_t l, lp, lk, sk, si;
1853
1854 if (p == 0)
1855 {
1856 u(0) = 1.0;
1857 d(0) = 0.0;
1858 return;
1859 }
1860
1861 lk = 1.0;
1862 for (k = 0; k < p; k++)
1863 {
1864 if (y >= (x(k) + x(k+1))/2)
1865 {
1866 lk *= y - x(k);
1867 }
1868 else
1869 {
1870 for (i = k+1; i <= p; i++)
1871 {
1872 lk *= y - x(i);
1873 }
1874 break;
1875 }
1876 }
1877 l = lk * (y - x(k));
1878
1879 sk = 0.0;
1880 for (i = 0; i < k; i++)
1881 {
1882 si = 1.0/(y - x(i));
1883 sk += si;
1884 u(i) = l * si * w(i);
1885 }
1886 u(k) = lk * w(k);
1887 for (i++; i <= p; i++)
1888 {
1889 si = 1.0/(y - x(i));
1890 sk += si;
1891 u(i) = l * si * w(i);
1892 }
1893 lp = l * sk + lk;
1894
1895 for (i = 0; i < k; i++)
1896 {
1897 d(i) = (lp * w(i) - u(i))/(y - x(i));
1898 }
1899 d(k) = sk * u(k);
1900 for (i++; i <= p; i++)
1901 {
1902 d(i) = (lp * w(i) - u(i))/(y - x(i));
1903 }
1904 break;
1905 }
1906 case Positive:
1907 CalcBernstein(x.Size() - 1, y, u, d);
1908 break;
1909 case Integrated:
1910 auxiliary_basis->Eval(y, u_aux, d_aux, d2_aux);
1911 EvalIntegrated(d_aux,u);
1912 EvalIntegrated(d2_aux,d);
1913 break;
1914 default: break;
1915 }
1916}
1917
1919 Vector &d2) const
1920{
1921 MFEM_VERIFY(etype == Barycentric,
1922 "Basis::Eval with second order derivatives not implemented for"
1923 " etype = " << etype);
1924 switch (etype)
1925 {
1926 case ChangeOfBasis:
1927 {
1928 CalcBasis(Ai.Width() - 1, y, x, w);
1929 Ai.Mult(x, u);
1930 Ai.Mult(w, d);
1931 // set d2 (not implemented yet)
1932 break;
1933 }
1934 case Barycentric:
1935 {
1936 int i, k, p = x.Size() - 1;
1937 real_t l, lp, lp2, lk, sk, si, sk2;
1938
1939 if (p == 0)
1940 {
1941 u(0) = 1.0;
1942 d(0) = 0.0;
1943 d2(0) = 0.0;
1944 return;
1945 }
1946
1947 lk = 1.0;
1948 for (k = 0; k < p; k++)
1949 {
1950 if (y >= (x(k) + x(k+1))/2)
1951 {
1952 lk *= y - x(k);
1953 }
1954 else
1955 {
1956 for (i = k+1; i <= p; i++)
1957 {
1958 lk *= y - x(i);
1959 }
1960 break;
1961 }
1962 }
1963 l = lk * (y - x(k));
1964
1965 sk = 0.0;
1966 sk2 = 0.0;
1967 for (i = 0; i < k; i++)
1968 {
1969 si = 1.0/(y - x(i));
1970 sk += si;
1971 sk2 -= si * si;
1972 u(i) = l * si * w(i);
1973 }
1974 u(k) = lk * w(k);
1975 for (i++; i <= p; i++)
1976 {
1977 si = 1.0/(y - x(i));
1978 sk += si;
1979 sk2 -= si * si;
1980 u(i) = l * si * w(i);
1981 }
1982 lp = l * sk + lk;
1983 lp2 = lp * sk + l * sk2 + sk * lk;
1984
1985 for (i = 0; i < k; i++)
1986 {
1987 d(i) = (lp * w(i) - u(i))/(y - x(i));
1988 d2(i) = (lp2 * w(i) - 2 * d(i))/(y - x(i));
1989 }
1990 d(k) = sk * u(k);
1991 d2(k) = sk2 * u(k) + sk * d(k);
1992 for (i++; i <= p; i++)
1993 {
1994 d(i) = (lp * w(i) - u(i))/(y - x(i));
1995 d2(i) = (lp2 * w(i) - 2 * d(i))/(y - x(i));
1996 }
1997 break;
1998 }
1999 case Positive:
2000 CalcBernstein(x.Size() - 1, y, u, d);
2001 break;
2002 case Integrated:
2003 MFEM_ABORT("Integrated basis must be evaluated with EvalIntegrated");
2004 break;
2005 default: break;
2006 }
2007}
2008
2010{
2011 MFEM_VERIFY(etype == Integrated,
2012 "EvalIntegrated is only valid for Integrated basis type");
2013 int p = d_aux_.Size() - 1;
2014 // See Gerritsma, M. (2010). "Edge functions for spectral element methods",
2015 // in Lecture Notes in Computational Science and Engineering, 199--207.
2016 u[0] = -d_aux_[0];
2017 for (int j=1; j<p; ++j)
2018 {
2019 u[j] = u[j-1] - d_aux_[j];
2020 }
2021 // If scale_integrated is true, the degrees of freedom represent mean values,
2022 // otherwise they represent subcell integrals. Generally, scale_integrated
2023 // should be true for MapType::VALUE, and false for other map types.
2024 if (scale_integrated)
2025 {
2026 Vector &aux_nodes = auxiliary_basis->x;
2027 for (int j=0; j<aux_nodes.Size()-1; ++j)
2028 {
2029 u[j] *= aux_nodes[j+1] - aux_nodes[j];
2030 }
2031 }
2032}
2033
2034void Poly_1D::Basis::ScaleIntegrated(bool scale_integrated_)
2035{
2036 scale_integrated = scale_integrated_;
2037}
2038
2040{
2041 delete auxiliary_basis;
2042}
2043
2044const int *Poly_1D::Binom(const int p)
2045{
2046 if (binom.NumCols() <= p)
2047 {
2048 binom.SetSize(p + 1, p + 1);
2049 for (int i = 0; i <= p; i++)
2050 {
2051 binom(i,0) = binom(i,i) = 1;
2052 for (int j = 1; j < i; j++)
2053 {
2054 binom(i,j) = binom(i-1,j) + binom(i-1,j-1);
2055 }
2056 }
2057 }
2058 return binom[p];
2059}
2060
2062{
2063 for (int i = 0; i <= p; i++)
2064 {
2065 // x[i] = 0.5*(1. + cos(M_PI*(p - i + 0.5)/(p + 1)));
2066 real_t s = sin(M_PI_2*(i + 0.5)/(p + 1));
2067 x[i] = s*s;
2068 }
2069}
2070
2071void Poly_1D::CalcMono(const int p, const real_t x, real_t *u)
2072{
2073 real_t xn;
2074 u[0] = xn = 1.;
2075 for (int n = 1; n <= p; n++)
2076 {
2077 u[n] = (xn *= x);
2078 }
2079}
2080
2081void Poly_1D::CalcMono(const int p, const real_t x, real_t *u, real_t *d)
2082{
2083 real_t xn;
2084 u[0] = xn = 1.;
2085 d[0] = 0.;
2086 for (int n = 1; n <= p; n++)
2087 {
2088 d[n] = n * xn;
2089 u[n] = (xn *= x);
2090 }
2091}
2092
2093void Poly_1D::CalcBinomTerms(const int p, const real_t x, const real_t y,
2094 real_t *u)
2095{
2096 if (p == 0)
2097 {
2098 u[0] = 1.;
2099 }
2100 else
2101 {
2102 int i;
2103 const int *b = Binom(p);
2104 real_t z = x;
2105
2106 for (i = 1; i < p; i++)
2107 {
2108 u[i] = b[i]*z;
2109 z *= x;
2110 }
2111 u[p] = z;
2112 z = y;
2113 for (i--; i > 0; i--)
2114 {
2115 u[i] *= z;
2116 z *= y;
2117 }
2118 u[0] = z;
2119 }
2120}
2121
2122void Poly_1D::CalcBinomTerms(const int p, const real_t x, const real_t y,
2123 real_t *u, real_t *d)
2124{
2125 if (p == 0)
2126 {
2127 u[0] = 1.;
2128 d[0] = 0.;
2129 }
2130 else
2131 {
2132 int i;
2133 const int *b = Binom(p);
2134 const real_t xpy = x + y, ptx = p*x;
2135 real_t z = 1.;
2136
2137 for (i = 1; i < p; i++)
2138 {
2139 d[i] = b[i]*z*(i*xpy - ptx);
2140 z *= x;
2141 u[i] = b[i]*z;
2142 }
2143 d[p] = p*z;
2144 u[p] = z*x;
2145 z = 1.;
2146 for (i--; i > 0; i--)
2147 {
2148 d[i] *= z;
2149 z *= y;
2150 u[i] *= z;
2151 }
2152 d[0] = -p*z;
2153 u[0] = z*y;
2154 }
2155}
2156
2157void Poly_1D::CalcDBinomTerms(const int p, const real_t x, const real_t y,
2158 real_t *d)
2159{
2160 if (p == 0)
2161 {
2162 d[0] = 0.;
2163 }
2164 else
2165 {
2166 int i;
2167 const int *b = Binom(p);
2168 const real_t xpy = x + y, ptx = p*x;
2169 real_t z = 1.;
2170
2171 for (i = 1; i < p; i++)
2172 {
2173 d[i] = b[i]*z*(i*xpy - ptx);
2174 z *= x;
2175 }
2176 d[p] = p*z;
2177 z = 1.;
2178 for (i--; i > 0; i--)
2179 {
2180 d[i] *= z;
2181 z *= y;
2182 }
2183 d[0] = -p*z;
2184 }
2185}
2186
2187void Poly_1D::CalcDxBinomTerms(const int p, const real_t x, const real_t y,
2188 real_t *u)
2189{
2190 if (p == 0)
2191 {
2192 u[0] = 0.;
2193 }
2194 else
2195 {
2196 int i;
2197 const int *b = Binom(p);
2198 real_t z = 1.;
2199
2200 for (i = 1; i < p; i++)
2201 {
2202 u[i] = i * b[i]*z;
2203 z *= x;
2204 }
2205 u[p] = i * z;
2206 z = y;
2207 for (i--; i > 0; i--)
2208 {
2209 u[i] *= z;
2210 z *= y;
2211 }
2212 u[0] = 0;
2213 }
2214}
2215
2216void Poly_1D::CalcDyBinomTerms(const int p, const real_t x, const real_t y,
2217 real_t *u)
2218{
2219 if (p == 0)
2220 {
2221 u[0] = 0.;
2222 }
2223 else
2224 {
2225 int i;
2226 const int *b = Binom(p);
2227 real_t z = x;
2228
2229 for (i = 1; i < p; i++)
2230 {
2231 u[i] = b[i]*z;
2232 z *= x;
2233 }
2234 u[p] = 0.;
2235 z = 1.;
2236 for (i--; i > 0; i--)
2237 {
2238 u[i] *= (p - i) * z;
2239 z *= y;
2240 }
2241 u[0] = p * z;
2242 }
2243}
2244
2245void Poly_1D::CalcLegendre(const int p, const real_t x, real_t *u)
2246{
2247 // use the recursive definition for [-1,1]:
2248 // (n+1)*P_{n+1}(z) = (2*n+1)*z*P_n(z)-n*P_{n-1}(z)
2249 real_t z;
2250 u[0] = 1.;
2251 if (p == 0) { return; }
2252 u[1] = z = 2.*x - 1.;
2253 for (int n = 1; n < p; n++)
2254 {
2255 u[n+1] = ((2*n + 1)*z*u[n] - n*u[n-1])/(n + 1);
2256 }
2257}
2258
2259void Poly_1D::CalcLegendre(const int p, const real_t x, real_t *u, real_t *d)
2260{
2261 // use the recursive definition for [-1,1]:
2262 // (n+1)*P_{n+1}(z) = (2*n+1)*z*P_n(z)-n*P_{n-1}(z)
2263 // for the derivative use, z in [-1,1]:
2264 // P'_{n+1}(z) = (2*n+1)*P_n(z)+P'_{n-1}(z)
2265 real_t z;
2266 u[0] = 1.;
2267 d[0] = 0.;
2268 if (p == 0) { return; }
2269 u[1] = z = 2.*x - 1.;
2270 d[1] = 2.;
2271 for (int n = 1; n < p; n++)
2272 {
2273 u[n+1] = ((2*n + 1)*z*u[n] - n*u[n-1])/(n + 1);
2274 d[n+1] = (4*n + 2)*u[n] + d[n-1];
2275 }
2276}
2277
2278void Poly_1D::CalcChebyshev(const int p, const real_t x, real_t *u)
2279{
2280 // recursive definition, z in [-1,1]
2281 // T_0(z) = 1, T_1(z) = z
2282 // T_{n+1}(z) = 2*z*T_n(z) - T_{n-1}(z)
2283 real_t z;
2284 u[0] = 1.;
2285 if (p == 0) { return; }
2286 u[1] = z = 2.*x - 1.;
2287 for (int n = 1; n < p; n++)
2288 {
2289 u[n+1] = 2*z*u[n] - u[n-1];
2290 }
2291}
2292
2293void Poly_1D::CalcChebyshev(const int p, const real_t x, real_t *u, real_t *d)
2294{
2295 // recursive definition, z in [-1,1]
2296 // T_0(z) = 1, T_1(z) = z
2297 // T_{n+1}(z) = 2*z*T_n(z) - T_{n-1}(z)
2298 // T'_n(z) = n*U_{n-1}(z)
2299 // U_0(z) = 1 U_1(z) = 2*z
2300 // U_{n+1}(z) = 2*z*U_n(z) - U_{n-1}(z)
2301 // U_n(z) = z*U_{n-1}(z) + T_n(z) = z*T'_n(z)/n + T_n(z)
2302 // T'_{n+1}(z) = (n + 1)*(z*T'_n(z)/n + T_n(z))
2303 real_t z;
2304 u[0] = 1.;
2305 d[0] = 0.;
2306 if (p == 0) { return; }
2307 u[1] = z = 2.*x - 1.;
2308 d[1] = 2.;
2309 for (int n = 1; n < p; n++)
2310 {
2311 u[n+1] = 2*z*u[n] - u[n-1];
2312 d[n+1] = (n + 1)*(z*d[n]/n + 2*u[n]);
2313 }
2314}
2315
2316void Poly_1D::CalcChebyshev(const int p, const real_t x, real_t *u, real_t *d,
2317 real_t *dd)
2318{
2319 // recursive definition, z in [-1,1]
2320 // T_0(z) = 1, T_1(z) = z
2321 // T_{n+1}(z) = 2*z*T_n(z) - T_{n-1}(z)
2322 // T'_n(z) = n*U_{n-1}(z)
2323 // U_0(z) = 1 U_1(z) = 2*z
2324 // U_{n+1}(z) = 2*z*U_n(z) - U_{n-1}(z)
2325 // U_n(z) = z*U_{n-1}(z) + T_n(z) = z*T'_n(z)/n + T_n(z)
2326 // T'_{n+1}(z) = (n + 1)*(z*T'_n(z)/n + T_n(z))
2327 // T''_{n+1}(z) = (n + 1)*(2*(n + 1)*T'_n(z) + z*T''_n(z)) / n
2328 real_t z;
2329 u[0] = 1.;
2330 d[0] = 0.;
2331 dd[0]= 0.;
2332 if (p == 0) { return; }
2333 u[1] = z = 2.*x - 1.;
2334 d[1] = 2.;
2335 dd[1] = 0;
2336 for (int n = 1; n < p; n++)
2337 {
2338 u[n+1] = 2*z*u[n] - u[n-1];
2339 d[n+1] = (n + 1)*(z*d[n]/n + 2*u[n]);
2340 dd[n+1] = (n + 1)*(2.*(n + 1)*d[n] + z*dd[n])/n;
2341 }
2342}
2343
2344const Array<real_t>* Poly_1D::GetPointsArray(const int p, const int btype)
2345{
2346 Array<real_t> *val;
2347 BasisType::Check(btype);
2348 const int qtype = BasisType::GetQuadrature1D(btype);
2349 if (qtype == Quadrature1D::Invalid) { return nullptr; }
2350
2351#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
2352 #pragma omp critical (Poly1DGetPoints)
2353#endif
2354 {
2355 std::pair<int, int> key(btype, p);
2356 auto it = points_container.find(key);
2357 if (it == points_container.end())
2358 {
2359 it = points_container.emplace(key, new Array<real_t>(p + 1, h_mt)).first;
2360 val = it->second.get();
2361 real_t* hptr = val->HostWrite();
2362 quad_func.GivePolyPoints(p + 1, hptr, qtype);
2363 }
2364 else
2365 {
2366 val = it->second.get();
2367 }
2368 }
2369 return val;
2370}
2371
2372Poly_1D::Basis &Poly_1D::GetBasis(const int p, const int btype)
2373{
2374 BasisType::Check(btype);
2375 Basis* val;
2376
2377#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
2378 #pragma omp critical (Poly1DGetBasis)
2379#endif
2380 {
2381 std::pair<int, int> key(btype, p);
2382 auto it = bases_container.find(key);
2383 if (it == bases_container.end())
2384 {
2385 EvalType etype;
2386 if (btype == BasisType::Positive) { etype = Positive; }
2387 else if (btype == BasisType::IntegratedGLL) { etype = Integrated; }
2388 else { etype = Barycentric; }
2389 it = bases_container
2390 .emplace(key, new Basis(p, GetPoints(p, btype), etype))
2391 .first;
2392 }
2393 val = it->second.get();
2394 }
2395 return *val;
2396}
2397
2398
2400 const int btype, const DofMapType dmtype)
2401 : b_type(btype),
2402 basis1d(poly1d.GetBasis(p, b_type))
2403{
2404 if (dmtype == H1_DOF_MAP || dmtype == Sr_DOF_MAP)
2405 {
2406 switch (dims)
2407 {
2408 case 1:
2409 {
2410 dof_map.SetSize(p + 1);
2411 dof_map[0] = 0;
2412 dof_map[p] = 1;
2413 for (int i = 1; i < p; i++)
2414 {
2415 dof_map[i] = i+1;
2416 }
2417 break;
2418 }
2419 case 2:
2420 {
2421 const int p1 = p + 1;
2422 dof_map.SetSize(p1*p1);
2423
2424 // vertices
2425 dof_map[0 + 0*p1] = 0;
2426 dof_map[p + 0*p1] = 1;
2427 dof_map[p + p*p1] = 2;
2428 dof_map[0 + p*p1] = 3;
2429
2430 // edges
2431 int o = 4;
2432 for (int i = 1; i < p; i++)
2433 {
2434 dof_map[i + 0*p1] = o++;
2435 }
2436 for (int i = 1; i < p; i++)
2437 {
2438 dof_map[p + i*p1] = o++;
2439 }
2440 for (int i = 1; i < p; i++)
2441 {
2442 dof_map[(p-i) + p*p1] = o++;
2443 }
2444 for (int i = 1; i < p; i++)
2445 {
2446 dof_map[0 + (p-i)*p1] = o++;
2447 }
2448
2449 // interior
2450 for (int j = 1; j < p; j++)
2451 {
2452 for (int i = 1; i < p; i++)
2453 {
2454 dof_map[i + j*p1] = o++;
2455 }
2456 }
2457 break;
2458 }
2459 case 3:
2460 {
2461 const int p1 = p + 1;
2462 dof_map.SetSize(p1*p1*p1);
2463
2464 // vertices
2465 dof_map[0 + (0 + 0*p1)*p1] = 0;
2466 dof_map[p + (0 + 0*p1)*p1] = 1;
2467 dof_map[p + (p + 0*p1)*p1] = 2;
2468 dof_map[0 + (p + 0*p1)*p1] = 3;
2469 dof_map[0 + (0 + p*p1)*p1] = 4;
2470 dof_map[p + (0 + p*p1)*p1] = 5;
2471 dof_map[p + (p + p*p1)*p1] = 6;
2472 dof_map[0 + (p + p*p1)*p1] = 7;
2473
2474 // edges (see Hexahedron::edges in mesh/hexahedron.cpp).
2475 // edges (see Constants<Geometry::CUBE>::Edges in fem/geom.cpp).
2476 int o = 8;
2477 for (int i = 1; i < p; i++)
2478 {
2479 dof_map[i + (0 + 0*p1)*p1] = o++; // (0,1)
2480 }
2481 for (int i = 1; i < p; i++)
2482 {
2483 dof_map[p + (i + 0*p1)*p1] = o++; // (1,2)
2484 }
2485 for (int i = 1; i < p; i++)
2486 {
2487 dof_map[i + (p + 0*p1)*p1] = o++; // (3,2)
2488 }
2489 for (int i = 1; i < p; i++)
2490 {
2491 dof_map[0 + (i + 0*p1)*p1] = o++; // (0,3)
2492 }
2493 for (int i = 1; i < p; i++)
2494 {
2495 dof_map[i + (0 + p*p1)*p1] = o++; // (4,5)
2496 }
2497 for (int i = 1; i < p; i++)
2498 {
2499 dof_map[p + (i + p*p1)*p1] = o++; // (5,6)
2500 }
2501 for (int i = 1; i < p; i++)
2502 {
2503 dof_map[i + (p + p*p1)*p1] = o++; // (7,6)
2504 }
2505 for (int i = 1; i < p; i++)
2506 {
2507 dof_map[0 + (i + p*p1)*p1] = o++; // (4,7)
2508 }
2509 for (int i = 1; i < p; i++)
2510 {
2511 dof_map[0 + (0 + i*p1)*p1] = o++; // (0,4)
2512 }
2513 for (int i = 1; i < p; i++)
2514 {
2515 dof_map[p + (0 + i*p1)*p1] = o++; // (1,5)
2516 }
2517 for (int i = 1; i < p; i++)
2518 {
2519 dof_map[p + (p + i*p1)*p1] = o++; // (2,6)
2520 }
2521 for (int i = 1; i < p; i++)
2522 {
2523 dof_map[0 + (p + i*p1)*p1] = o++; // (3,7)
2524 }
2525
2526 // faces (see Mesh::GenerateFaces in mesh/mesh.cpp)
2527 for (int j = 1; j < p; j++)
2528 {
2529 for (int i = 1; i < p; i++)
2530 {
2531 dof_map[i + ((p-j) + 0*p1)*p1] = o++; // (3,2,1,0)
2532 }
2533 }
2534 for (int j = 1; j < p; j++)
2535 {
2536 for (int i = 1; i < p; i++)
2537 {
2538 dof_map[i + (0 + j*p1)*p1] = o++; // (0,1,5,4)
2539 }
2540 }
2541 for (int j = 1; j < p; j++)
2542 {
2543 for (int i = 1; i < p; i++)
2544 {
2545 dof_map[p + (i + j*p1)*p1] = o++; // (1,2,6,5)
2546 }
2547 }
2548 for (int j = 1; j < p; j++)
2549 {
2550 for (int i = 1; i < p; i++)
2551 {
2552 dof_map[(p-i) + (p + j*p1)*p1] = o++; // (2,3,7,6)
2553 }
2554 }
2555 for (int j = 1; j < p; j++)
2556 {
2557 for (int i = 1; i < p; i++)
2558 {
2559 dof_map[0 + ((p-i) + j*p1)*p1] = o++; // (3,0,4,7)
2560 }
2561 }
2562 for (int j = 1; j < p; j++)
2563 {
2564 for (int i = 1; i < p; i++)
2565 {
2566 dof_map[i + (j + p*p1)*p1] = o++; // (4,5,6,7)
2567 }
2568 }
2569
2570 // interior
2571 for (int k = 1; k < p; k++)
2572 {
2573 for (int j = 1; j < p; j++)
2574 {
2575 for (int i = 1; i < p; i++)
2576 {
2577 dof_map[i + (j + k*p1)*p1] = o++;
2578 }
2579 }
2580 }
2581 break;
2582 }
2583 default:
2584 MFEM_ABORT("invalid dimension: " << dims);
2585 break;
2586 }
2587 }
2588 else if (dmtype == L2_DOF_MAP)
2589 {
2590 // leave dof_map empty, indicating that the dofs are ordered
2591 // lexicographically, i.e. the dof_map is identity
2592 }
2593 else
2594 {
2595 MFEM_ABORT("invalid DofMapType: " << dmtype);
2596 }
2597}
2598
2600 const FiniteElement &fe, const IntegrationRule &ir,
2601 DofToQuad::Mode mode, const Poly_1D::Basis &basis, bool closed,
2602 Array<DofToQuad*> &dof2quad_array)
2603{
2604 DofToQuad *d2q = nullptr;
2605 MFEM_VERIFY(mode == DofToQuad::TENSOR, "invalid mode requested");
2606
2607#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
2608 #pragma omp critical (DofToQuad)
2609#endif
2610 {
2611 for (int i = 0; i < dof2quad_array.Size(); i++)
2612 {
2613 d2q = dof2quad_array[i];
2614 if (d2q->IntRule != &ir || d2q->mode != mode) { d2q = nullptr; }
2615 }
2616 if (!d2q)
2617 {
2618 d2q = new DofToQuad;
2619 const int ndof = closed ? fe.GetOrder() + 1 : fe.GetOrder();
2620 const int nqpt = (int)floor(pow(ir.GetNPoints(), 1.0/fe.GetDim()) + 0.5);
2621 d2q->FE = &fe;
2622 d2q->IntRule = &ir;
2623 d2q->mode = mode;
2624 d2q->ndof = ndof;
2625 d2q->nqpt = nqpt;
2626 d2q->B.SetSize(nqpt*ndof);
2627 d2q->Bt.SetSize(ndof*nqpt);
2628 d2q->G.SetSize(nqpt*ndof);
2629 d2q->Gt.SetSize(ndof*nqpt);
2630 Vector val(ndof), grad(ndof);
2631 for (int i = 0; i < nqpt; i++)
2632 {
2633 // The first 'nqpt' points in 'ir' have the same x-coordinates as those
2634 // of the 1D rule.
2635 basis.Eval(ir.IntPoint(i).x, val, grad);
2636 for (int j = 0; j < ndof; j++)
2637 {
2638 d2q->B[i+nqpt*j] = d2q->Bt[j+ndof*i] = val(j);
2639 d2q->G[i+nqpt*j] = d2q->Gt[j+ndof*i] = grad(j);
2640 }
2641 }
2642 dof2quad_array.Append(d2q);
2643 }
2644 }
2645 return *d2q;
2646}
2647
2649 const int p,
2650 const int btype,
2651 const DofMapType dmtype)
2652 : NodalFiniteElement(dims, GetTensorProductGeometry(dims), Pow(p + 1, dims),
2653 p, dims > 1 ? FunctionSpace::Qk : FunctionSpace::Pk),
2654 TensorBasisElement(dims, p, btype, dmtype)
2655{
2657}
2658
2660{
2662 // If we are using the "integrated" basis, the basis functions should be
2663 // scaled for MapType::VALUE, and not scaled for MapType::INTEGRAL. This
2664 // ensures spectral equivalence of the mass matrix with its low-order-refined
2665 // counterpart (cf. LORDiscretization)
2667 {
2669 }
2670}
2671
2673 const IntegrationRule &ir,
2674 DofToQuad::Mode mode) const
2675{
2676 if (mode != DofToQuad::TENSOR)
2677 {
2678 return NodalFiniteElement::GetDofToQuad(ir, mode);
2679 }
2680 else
2681 {
2682 return GetTensorDofToQuad(*this, ir, mode, basis1d, true, dof2quad_array);
2683 }
2684}
2685
2687 Array<int> &face_map) const
2688{
2689 internal::GetTensorFaceMap(dim, order, face_id, face_map);
2690}
2691
2693 const int d,
2694 const int p,
2695 const int cbtype,
2696 const int obtype,
2697 const int M,
2698 const DofMapType dmtype)
2699 : VectorFiniteElement(dims, GetTensorProductGeometry(dims), d,
2700 p, M, FunctionSpace::Qk),
2701 TensorBasisElement(dims, p, VerifyNodal(VerifyClosed(cbtype)), dmtype),
2702 obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(obtype)))
2703{
2704 MFEM_VERIFY(dims > 1, "Constructor for VectorTensorFiniteElement with both "
2705 "open and closed bases is not valid for 1D elements.");
2706}
2707
2709 const int d,
2710 const int p,
2711 const int obtype,
2712 const int M,
2713 const DofMapType dmtype)
2714 : VectorFiniteElement(dims, GetTensorProductGeometry(dims), d,
2715 p, M, FunctionSpace::Pk),
2716 TensorBasisElement(dims, p, VerifyOpen(obtype), dmtype),
2717 obasis1d(poly1d.GetBasis(p, VerifyOpen(obtype)))
2718{
2719 MFEM_VERIFY(dims == 1, "Constructor for VectorTensorFiniteElement without "
2720 "closed basis is only valid for 1D elements.");
2721}
2722
2724{
2725 for (int i = 0; i < dof2quad_array_open.Size(); i++)
2726 {
2727 delete dof2quad_array_open[i];
2728 }
2729}
2730
2731}
int NumCols() const
Definition array.hpp:408
void SetSize(int m, int n)
Definition array.hpp:405
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
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:830
T * HostWrite()
Shortcut for mfem::Write(a.GetMemory(), a.Size(), false).
Definition array.hpp:349
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
@ GaussLobatto
Closed type.
Definition fe_base.hpp:36
@ Positive
Bernstein polynomials.
Definition fe_base.hpp:37
@ IntegratedGLL
Integrated GLL indicator functions.
Definition fe_base.hpp:43
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:125
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
Definition densemat.cpp:172
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:319
real_t InnerProduct(const real_t *x, const real_t *y) const
Compute y^t A x.
Definition densemat.cpp:301
real_t * Data() const
Returns the matrix data array.
Definition densemat.hpp:114
real_t * GetData() const
Returns the matrix data array.
Definition densemat.hpp:118
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:108
void Invert()
Replaces the current matrix with its inverse.
Definition densemat.cpp:707
void InvRightScaling(const Vector &s)
InvRightScaling: this = this * diag(1./s);.
Definition densemat.cpp:360
void GetColumn(int c, Vector &col) const
real_t FNorm2() const
Compute the square of the Frobenius norm of the matrix.
Definition densemat.hpp:276
Structure representing the matrices/tensors needed to evaluate (in reference space) the values,...
Definition fe_base.hpp:141
Mode mode
Describes the contents of the B, Bt, G, and Gt arrays, see Mode.
Definition fe_base.hpp:174
const IntegrationRule * IntRule
IntegrationRule that defines the quadrature points at which the basis functions of the FE are evaluat...
Definition fe_base.hpp:150
Array< real_t > G
Gradients/divergences/curls of basis functions evaluated at quadrature points.
Definition fe_base.hpp:214
Mode
Type of data stored in the arrays B, Bt, G, and Gt.
Definition fe_base.hpp:154
@ FULL
Full multidimensional representation which does not use tensor product structure. The ordering of the...
Definition fe_base.hpp:158
@ LEXICOGRAPHIC_FULL
Full multidimensional representation which does not use tensor product structure. The ordering of the...
Definition fe_base.hpp:170
@ TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition fe_base.hpp:165
Array< real_t > B
Basis functions evaluated at quadrature points.
Definition fe_base.hpp:193
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:178
Array< real_t > Gt
Transpose of G.
Definition fe_base.hpp:221
const class FiniteElement * FE
The FiniteElement that created and owns this object.
Definition fe_base.hpp:145
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:182
Array< real_t > Bt
Transpose of B.
Definition fe_base.hpp:199
const DenseMatrix & Hessian()
Return the Hessian matrix of the transformation at the currently set IntegrationPoint,...
Definition eltrans.hpp:138
const DenseMatrix & InverseJacobian()
Return the inverse of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition eltrans.hpp:158
const DenseMatrix & AdjugateJacobian()
Return the adjugate of the Jacobian matrix of the transformation at the currently set IntegrationPoin...
Definition eltrans.hpp:148
const IntegrationPoint & GetIntPoint()
Get a const reference to the currently set integration point. This will return NULL if no integration...
Definition eltrans.hpp:111
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:144
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint,...
Definition eltrans.hpp:132
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition eltrans.hpp:106
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:244
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:253
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
virtual const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const
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:325
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition fe_base.hpp:338
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
virtual void ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Compute the discrete gradient matrix from the given FiniteElement onto 'this' FiniteElement....
Definition fe_base.cpp:161
IntegrationRule Nodes
Definition fe_base.hpp:256
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:321
int vdim
Vector dimension of vector-valued basis functions.
Definition fe_base.hpp:247
void CalcPhysHessian(ElementTransformation &Trans, DenseMatrix &Hessian) const
Evaluate the Hessian of all shape functions of a scalar finite element in physical space at the given...
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:255
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:360
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
Definition fe_base.hpp:351
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:331
int cdim
Dimension of curl for vector-valued basis functions.
Definition fe_base.hpp:248
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:306
@ NONE
No derivatives implemented.
Definition fe_base.hpp:304
@ CURL
Implements CalcCurlShape methods.
Definition fe_base.hpp:307
@ GRAD
Implements CalcDShape methods.
Definition fe_base.hpp:305
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:249
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
Array< DofToQuad * > dof2quad_array
Container for all DofToQuad objects created by the FiniteElement.
Definition fe_base.hpp:263
void CalcPhysLaplacian(ElementTransformation &Trans, Vector &Laplacian) const
Evaluate the Laplacian of all shape functions of a scalar finite element in physical space at the giv...
Definition fe_base.cpp:203
DenseMatrix vshape
Definition fe_base.hpp:258
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:334
int order
Order/degree of the shape functions.
Definition fe_base.hpp:254
void CalcPhysLinLaplacian(ElementTransformation &Trans, Vector &Laplacian) const
Evaluate the Laplacian of all shape functions of a scalar finite element in physical space at the giv...
Definition fe_base.cpp:244
int dim
Dimension of reference space.
Definition fe_base.hpp:246
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:226
static const int MaxDim
Definition geom.hpp:47
const IntegrationPoint & GetCenter(int GeomType) const
Return the center of the given Geometry::Type, GeomType.
Definition geom.hpp:75
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:721
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 ReorderLexToNative(int ncomp, Vector &dofs) const
Definition fe_base.cpp:976
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
void ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const override
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:726
const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const override
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:2686
NodalTensorFiniteElement(const int dims, const int p, const int btype, const DofMapType dmtype)
Definition fe_base.cpp:2648
const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const override
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition fe_base.cpp:2672
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
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:1003
void Eval(const real_t x, Vector &u) const
Evaluate the basis functions at point x in [0,1].
Definition fe_base.cpp:1777
bool IsIntegratedType() const
Returns true if the basis is "integrated", false otherwise.
Definition fe_base.hpp:1047
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:2034
void EvalIntegrated(const Vector &d, Vector &i) const
Evaluate the "integrated" basis type using pre-computed closed basis derivatives.
Definition fe_base.cpp:2009
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:1713
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:2157
const real_t * GetPoints(const int p, const int btype, bool on_device=false)
Get the coordinates of the points of the given BasisType, btype.
Definition fe_base.hpp:1106
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:2044
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:2372
const Array< real_t > * GetPointsArray(const int p, const int btype)
Get the coordinates of the points of the given BasisType, btype.
Definition fe_base.cpp:2344
static void CalcDyBinomTerms(const int p, const real_t x, const real_t y, real_t *d)
Compute the derivatives (w.r.t. y) of the terms in the expansion of the binomial (x + y)^p....
Definition fe_base.cpp:2216
static void CalcLegendre(const int p, const real_t x, real_t *u)
Definition fe_base.cpp:2245
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:1215
EvalType
One-dimensional basis evaluation type.
Definition fe_base.hpp:992
@ ChangeOfBasis
Use change of basis, O(p^2) Evals.
Definition fe_base.hpp:993
@ Integrated
Integrated indicator functions (cf. Gerritsma)
Definition fe_base.hpp:996
@ Positive
Fast evaluation of Bernstein polynomials.
Definition fe_base.hpp:995
@ Barycentric
Use barycentric Lagrangian interpolation, O(p) Evals.
Definition fe_base.hpp:994
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:2061
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:1140
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:2093
static void CalcDxBinomTerms(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....
Definition fe_base.cpp:2187
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:662
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:687
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:1251
static const DofToQuad & GetTensorDofToQuad(const FiniteElement &fe, const IntegrationRule &ir, DofToQuad::Mode mode, const Poly_1D::Basis &basis, bool closed, Array< DofToQuad * > &dof2quad_array)
Definition fe_base.cpp:2599
TensorBasisElement(const int dims, const int p, const int btype, const DofMapType dmtype)
Definition fe_base.cpp:2399
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:813
void ProjectGrad_RT(const real_t *nk, const Array< int > &d2n, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Definition fe_base.cpp:1222
void LocalRestriction_RT(const real_t *nk, const Array< int > &d2n, ElementTransformation &Trans, DenseMatrix &R) const
Definition fe_base.cpp:1628
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:1306
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:1082
void LocalInterpolation_RT(const VectorFiniteElement &cfe, const real_t *nk, const Array< int > &d2n, ElementTransformation &Trans, DenseMatrix &I) const
Definition fe_base.cpp:1502
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:1335
void LocalRestriction_ND(const real_t *tk, const Array< int > &d2t, ElementTransformation &Trans, DenseMatrix &R) const
Definition fe_base.cpp:1671
void LocalL2Projection_RT(const VectorFiniteElement &cfe, ElementTransformation &Trans, DenseMatrix &I) const
Definition fe_base.cpp:1455
void ProjectGrad_ND(const real_t *tk, const Array< int > &d2t, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Definition fe_base.cpp:1434
void CalcVShape_ND(ElementTransformation &Trans, DenseMatrix &shape) const
Definition fe_base.cpp:1071
void CalcVShape_RT(ElementTransformation &Trans, DenseMatrix &shape) const
Definition fe_base.cpp:1059
VectorFiniteElement(int D, Geometry::Type G, int Do, int O, int M, int F=FunctionSpace::Pk)
Definition fe_base.cpp:993
void LocalInterpolation_ND(const VectorFiniteElement &cfe, const real_t *tk, const Array< int > &d2t, ElementTransformation &Trans, DenseMatrix &I) const
Definition fe_base.cpp:1589
void ProjectCurl_RT(const real_t *nk, const Array< int > &d2n, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Definition fe_base.cpp:1287
void LocalL2Projection_ND(const VectorFiniteElement &cfe, ElementTransformation &Trans, DenseMatrix &I) const
Definition fe_base.cpp:1543
void ProjectCurl_ND(const real_t *tk, const Array< int > &d2t, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Definition fe_base.cpp:1249
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:1119
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:2692
Vector data type.
Definition vector.hpp:82
void SetDataAndSize(real_t *d, int s)
Set the Vector data and size.
Definition vector.hpp:183
int Size() const
Returns the size of the vector.
Definition vector.hpp:226
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:558
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition vector.hpp:235
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:47
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:548
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
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition intrules.hpp:492
real_t p(const Vector &x, real_t t)
std::array< int, NCMesh::MaxFaceNodes > nodes