MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
fe_l2.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// L2 Finite Element classes
13
14#include "fe_l2.hpp"
15#include "fe_h1.hpp"
16#include "../coefficient.hpp"
17
18namespace mfem
19{
20
21using namespace std;
22
23L2_SegmentElement::L2_SegmentElement(const int p, const int btype)
24 : NodalTensorFiniteElement(1, p, VerifyOpen(btype), L2_DOF_MAP)
25{
26 const real_t *op = poly1d.OpenPoints(p, btype);
27
28#ifndef MFEM_THREAD_SAFE
29 shape_x.SetSize(p + 1);
30 dshape_x.SetDataAndSize(NULL, p + 1);
31#endif
32
33 for (int i = 0; i <= p; i++)
34 {
35 Nodes.IntPoint(i).x = op[i];
36 }
37}
38
40 Vector &shape) const
41{
43 basis1d.Eval(ip.x, shape);
44}
45
47 DenseMatrix &dshape) const
48{
49#ifdef MFEM_THREAD_SAFE
50 Vector shape_x(dof), dshape_x(dshape.Data(), dof);
51#else
52 dshape_x.SetData(dshape.Data());
53#endif
55 basis1d.Eval(ip.x, shape_x, dshape_x);
56}
57
58void L2_SegmentElement::ProjectDelta(int vertex, Vector &dofs) const
59{
60 const int p = order;
61 const real_t *op = poly1d.OpenPoints(p, b_type);
62
63 switch (vertex)
64 {
65 case 0:
66 for (int i = 0; i <= p; i++)
67 {
68 dofs(i) = poly1d.CalcDelta(p,(1.0 - op[i]));
69 }
70 break;
71
72 case 1:
73 for (int i = 0; i <= p; i++)
74 {
75 dofs(i) = poly1d.CalcDelta(p,op[i]);
76 }
77 break;
78 }
79}
80
81
83 : NodalTensorFiniteElement(2, p, VerifyOpen(btype), L2_DOF_MAP)
84{
85 const real_t *op = poly1d.OpenPoints(p, b_type);
86
87#ifndef MFEM_THREAD_SAFE
88 shape_x.SetSize(p + 1);
89 shape_y.SetSize(p + 1);
90 dshape_x.SetSize(p + 1);
91 dshape_y.SetSize(p + 1);
92#endif
93
94 for (int o = 0, j = 0; j <= p; j++)
95 for (int i = 0; i <= p; i++)
96 {
97 Nodes.IntPoint(o++).Set2(op[i], op[j]);
98 }
99}
100
102 Vector &shape) const
103{
104 const int p = order;
105
106#ifdef MFEM_THREAD_SAFE
107 Vector shape_x(p+1), shape_y(p+1);
108#endif
109
111 basis1d.Eval(ip.x, shape_x);
112 basis1d.Eval(ip.y, shape_y);
113
114 for (int o = 0, j = 0; j <= p; j++)
115 for (int i = 0; i <= p; i++)
116 {
117 shape(o++) = shape_x(i)*shape_y(j);
118 }
119}
120
122 DenseMatrix &dshape) const
123{
124 const int p = order;
125
126#ifdef MFEM_THREAD_SAFE
127 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
128#endif
129
131 basis1d.Eval(ip.x, shape_x, dshape_x);
132 basis1d.Eval(ip.y, shape_y, dshape_y);
133
134 for (int o = 0, j = 0; j <= p; j++)
135 for (int i = 0; i <= p; i++)
136 {
137 dshape(o,0) = dshape_x(i)* shape_y(j);
138 dshape(o,1) = shape_x(i)*dshape_y(j); o++;
139 }
140}
141
143{
144 const int p = order;
145 const real_t *op = poly1d.OpenPoints(p, b_type);
146
147#ifdef MFEM_THREAD_SAFE
148 Vector shape_x(p+1), shape_y(p+1);
149#endif
150
151 for (int i = 0; i <= p; i++)
152 {
153 shape_x(i) = poly1d.CalcDelta(p,(1.0 - op[i]));
154 shape_y(i) = poly1d.CalcDelta(p,op[i]);
155 }
156
157 switch (vertex)
158 {
159 case 0:
160 for (int o = 0, j = 0; j <= p; j++)
161 for (int i = 0; i <= p; i++)
162 {
163 dofs[o++] = shape_x(i)*shape_x(j);
164 }
165 break;
166 case 1:
167 for (int o = 0, j = 0; j <= p; j++)
168 for (int i = 0; i <= p; i++)
169 {
170 dofs[o++] = shape_y(i)*shape_x(j);
171 }
172 break;
173 case 2:
174 for (int o = 0, j = 0; j <= p; j++)
175 for (int i = 0; i <= p; i++)
176 {
177 dofs[o++] = shape_y(i)*shape_y(j);
178 }
179 break;
180 case 3:
181 for (int o = 0, j = 0; j <= p; j++)
182 for (int i = 0; i <= p; i++)
183 {
184 dofs[o++] = shape_x(i)*shape_y(j);
185 }
186 break;
187 }
188}
189
192 DenseMatrix &div) const
193{
195 {
196 // Compute subcell integrals of the divergence
197 const int fe_ndof = fe.GetDof();
198 Vector div_shape(fe_ndof);
199 div.SetSize(dof, fe_ndof);
200 div = 0.0;
201
202 const IntegrationRule &ir = IntRules.Get(geom_type, fe.GetOrder());
204
205 // Loop over subcells
206 for (int iy = 0; iy < order+1; ++iy)
207 {
208 const real_t hy = gll_pts[iy+1] - gll_pts[iy];
209 for (int ix = 0; ix < order+1; ++ix)
210 {
211 const int i = ix + iy*(order+1);
212 const real_t hx = gll_pts[ix+1] - gll_pts[ix];
213 // Loop over subcell quadrature points
214 for (int iq = 0; iq < ir.Size(); ++iq)
215 {
216 IntegrationPoint ip = ir[iq];
217 ip.x = gll_pts[ix] + hx*ip.x;
218 ip.y = gll_pts[iy] + hy*ip.y;
219 Trans.SetIntPoint(&ip);
220 fe.CalcDivShape(ip, div_shape);
221 real_t w = ip.weight;
222 if (map_type == VALUE)
223 {
224 const real_t detJ = Trans.Weight();
225 w /= detJ;
226 }
227 else if (map_type == INTEGRAL)
228 {
229 w *= hx*hy;
230 }
231 for (int j = 0; j < fe_ndof; j++)
232 {
233 const real_t div_j = div_shape(j);
234 div(i,j) += w*div_j;
235 }
236 }
237 }
238 }
239 // Filter small entries
240 for (int i = 0; i < dof; ++i)
241 {
242 for (int j = 0; j < fe_ndof; j++)
243 {
244 if (std::fabs(div(i,j)) < 1e-12) { div(i,j) = 0.0; }
245 }
246 }
247 }
248 else
249 {
250 // Fall back on standard nodal interpolation
251 NodalFiniteElement::ProjectDiv(fe, Trans, div);
252 }
253}
254
257 Vector &dofs) const
258{
260 {
263
264 dofs = 0.0;
265 // Loop over subcells
266 for (int iy = 0; iy < order+1; ++iy)
267 {
268 const real_t hy = gll_pts[iy+1] - gll_pts[iy];
269 for (int ix = 0; ix < order+1; ++ix)
270 {
271 const int i = ix + iy*(order+1);
272 const real_t hx = gll_pts[ix+1] - gll_pts[ix];
273 // Loop over subcell quadrature points
274 for (int iq = 0; iq < ir.Size(); ++iq)
275 {
276 IntegrationPoint ip = ir[iq];
277 ip.x = gll_pts[ix] + hx*ip.x;
278 ip.y = gll_pts[iy] + hy*ip.y;
279 Trans.SetIntPoint(&ip);
280 const real_t val = coeff.Eval(Trans, ip);
281 real_t w = ip.weight;
282 if (map_type == INTEGRAL)
283 {
284 w *= hx*hy*Trans.Weight();
285 }
286 dofs[i] += val*w;
287 }
288 }
289 }
290 }
291 else
292 {
293 NodalFiniteElement::Project(coeff, Trans, dofs);
294 }
295}
296
297
299 : NodalTensorFiniteElement(3, p, VerifyOpen(btype), L2_DOF_MAP)
300{
301 const real_t *op = poly1d.OpenPoints(p, btype);
302
303#ifndef MFEM_THREAD_SAFE
304 shape_x.SetSize(p + 1);
305 shape_y.SetSize(p + 1);
306 shape_z.SetSize(p + 1);
307 dshape_x.SetSize(p + 1);
308 dshape_y.SetSize(p + 1);
309 dshape_z.SetSize(p + 1);
310#endif
311
312 for (int o = 0, k = 0; k <= p; k++)
313 for (int j = 0; j <= p; j++)
314 for (int i = 0; i <= p; i++)
315 {
316 Nodes.IntPoint(o++).Set3(op[i], op[j], op[k]);
317 }
318}
319
321 Vector &shape) const
322{
323 const int p = order;
324
325#ifdef MFEM_THREAD_SAFE
326 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
327#endif
328
330 basis1d.Eval(ip.x, shape_x);
331 basis1d.Eval(ip.y, shape_y);
332 basis1d.Eval(ip.z, shape_z);
333
334 for (int o = 0, k = 0; k <= p; k++)
335 for (int j = 0; j <= p; j++)
336 for (int i = 0; i <= p; i++)
337 {
338 shape(o++) = shape_x(i)*shape_y(j)*shape_z(k);
339 }
340}
341
343 DenseMatrix &dshape) const
344{
345 const int p = order;
346
347#ifdef MFEM_THREAD_SAFE
348 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
349 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
350#endif
351
353 basis1d.Eval(ip.x, shape_x, dshape_x);
354 basis1d.Eval(ip.y, shape_y, dshape_y);
355 basis1d.Eval(ip.z, shape_z, dshape_z);
356
357 for (int o = 0, k = 0; k <= p; k++)
358 for (int j = 0; j <= p; j++)
359 for (int i = 0; i <= p; i++)
360 {
361 dshape(o,0) = dshape_x(i)* shape_y(j)* shape_z(k);
362 dshape(o,1) = shape_x(i)*dshape_y(j)* shape_z(k);
363 dshape(o,2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
364 }
365}
366
367void L2_HexahedronElement::ProjectDelta(int vertex, Vector &dofs) const
368{
369 const int p = order;
370 const real_t *op = poly1d.OpenPoints(p, b_type);
371
372#ifdef MFEM_THREAD_SAFE
373 Vector shape_x(p+1), shape_y(p+1);
374#endif
375
376 for (int i = 0; i <= p; i++)
377 {
378 shape_x(i) = poly1d.CalcDelta(p,(1.0 - op[i]));
379 shape_y(i) = poly1d.CalcDelta(p,op[i]);
380 }
381
382 switch (vertex)
383 {
384 case 0:
385 for (int o = 0, k = 0; k <= p; k++)
386 for (int j = 0; j <= p; j++)
387 for (int i = 0; i <= p; i++)
388 {
389 dofs[o++] = shape_x(i)*shape_x(j)*shape_x(k);
390 }
391 break;
392 case 1:
393 for (int o = 0, k = 0; k <= p; k++)
394 for (int j = 0; j <= p; j++)
395 for (int i = 0; i <= p; i++)
396 {
397 dofs[o++] = shape_y(i)*shape_x(j)*shape_x(k);
398 }
399 break;
400 case 2:
401 for (int o = 0, k = 0; k <= p; k++)
402 for (int j = 0; j <= p; j++)
403 for (int i = 0; i <= p; i++)
404 {
405 dofs[o++] = shape_y(i)*shape_y(j)*shape_x(k);
406 }
407 break;
408 case 3:
409 for (int o = 0, k = 0; k <= p; k++)
410 for (int j = 0; j <= p; j++)
411 for (int i = 0; i <= p; i++)
412 {
413 dofs[o++] = shape_x(i)*shape_y(j)*shape_x(k);
414 }
415 break;
416 case 4:
417 for (int o = 0, k = 0; k <= p; k++)
418 for (int j = 0; j <= p; j++)
419 for (int i = 0; i <= p; i++)
420 {
421 dofs[o++] = shape_x(i)*shape_x(j)*shape_y(k);
422 }
423 break;
424 case 5:
425 for (int o = 0, k = 0; k <= p; k++)
426 for (int j = 0; j <= p; j++)
427 for (int i = 0; i <= p; i++)
428 {
429 dofs[o++] = shape_y(i)*shape_x(j)*shape_y(k);
430 }
431 break;
432 case 6:
433 for (int o = 0, k = 0; k <= p; k++)
434 for (int j = 0; j <= p; j++)
435 for (int i = 0; i <= p; i++)
436 {
437 dofs[o++] = shape_y(i)*shape_y(j)*shape_y(k);
438 }
439 break;
440 case 7:
441 for (int o = 0, k = 0; k <= p; k++)
442 for (int j = 0; j <= p; j++)
443 for (int i = 0; i <= p; i++)
444 {
445 dofs[o++] = shape_x(i)*shape_y(j)*shape_y(k);
446 }
447 break;
448 }
449}
450
453 DenseMatrix &div) const
454{
456 {
457 // Compute subcell integrals of the divergence
458 const int fe_ndof = fe.GetDof();
459 Vector div_shape(fe_ndof);
460 div.SetSize(dof, fe_ndof);
461 div = 0.0;
462
463 const IntegrationRule &ir = IntRules.Get(geom_type, fe.GetOrder());
465
466 // Loop over subcells
467 for (int iz = 0; iz < order+1; ++iz)
468 {
469 const real_t hz = gll_pts[iz+1] - gll_pts[iz];
470 for (int iy = 0; iy < order+1; ++iy)
471 {
472 const real_t hy = gll_pts[iy+1] - gll_pts[iy];
473 for (int ix = 0; ix < order+1; ++ix)
474 {
475 const int i = ix + iy*(order+1) + iz*(order+1)*(order+1);
476 const real_t hx = gll_pts[ix+1] - gll_pts[ix];
477 // Loop over subcell quadrature points
478 for (int iq = 0; iq < ir.Size(); ++iq)
479 {
480 IntegrationPoint ip = ir[iq];
481 ip.x = gll_pts[ix] + hx*ip.x;
482 ip.y = gll_pts[iy] + hy*ip.y;
483 ip.z = gll_pts[iz] + hz*ip.z;
484 Trans.SetIntPoint(&ip);
485 fe.CalcDivShape(ip, div_shape);
486 real_t w = ip.weight;
487 if (map_type == VALUE)
488 {
489 const real_t detJ = Trans.Weight();
490 w /= detJ;
491 }
492 else if (map_type == INTEGRAL)
493 {
494 w *= hx*hy*hz;
495 }
496 for (int j = 0; j < fe_ndof; j++)
497 {
498 const real_t div_j = div_shape(j);
499 div(i,j) += w*div_j;
500 }
501 }
502 }
503 }
504 }
505 // Filter small entries
506 for (int i = 0; i < dof; ++i)
507 {
508 for (int j = 0; j < fe_ndof; j++)
509 {
510 if (std::fabs(div(i,j)) < 1e-12) { div(i,j) = 0.0; }
511 }
512 }
513 }
514 else
515 {
516 // Fall back on standard nodal interpolation
517 NodalFiniteElement::ProjectDiv(fe, Trans, div);
518 }
519}
520
523 Vector &dofs) const
524{
526 {
529
530 dofs = 0.0;
531 // Loop over subcells
532 for (int iz = 0; iz < order+1; ++iz)
533 {
534 const real_t hz = gll_pts[iz+1] - gll_pts[iz];
535 for (int iy = 0; iy < order+1; ++iy)
536 {
537 const real_t hy = gll_pts[iy+1] - gll_pts[iy];
538 for (int ix = 0; ix < order+1; ++ix)
539 {
540 const real_t hx = gll_pts[ix+1] - gll_pts[ix];
541 const int i = ix + iy*(order+1) + iz*(order+1)*(order+1);
542 // Loop over subcell quadrature points
543 for (int iq = 0; iq < ir.Size(); ++iq)
544 {
545 IntegrationPoint ip = ir[iq];
546 ip.x = gll_pts[ix] + hx*ip.x;
547 ip.y = gll_pts[iy] + hy*ip.y;
548 ip.z = gll_pts[iz] + hz*ip.z;
549 Trans.SetIntPoint(&ip);
550 const real_t val = coeff.Eval(Trans, ip);
551 real_t w = ip.weight;
552 if (map_type == INTEGRAL)
553 {
554 const real_t detJ = Trans.Weight();
555 w *= detJ*hx*hy*hz;
556 }
557 dofs[i] += val*w;
558 }
559 }
560 }
561 }
562 }
563 else
564 {
565 NodalFiniteElement::Project(coeff, Trans, dofs);
566 }
567}
568
569
570L2_TriangleElement::L2_TriangleElement(const int p, const int btype)
571 : NodalFiniteElement(2, Geometry::TRIANGLE, ((p + 1)*(p + 2))/2, p,
572 FunctionSpace::Pk)
573{
574 const real_t *op = poly1d.OpenPoints(p, VerifyOpen(btype));
575
576#ifndef MFEM_THREAD_SAFE
577 shape_x.SetSize(p + 1);
578 shape_y.SetSize(p + 1);
579 shape_l.SetSize(p + 1);
580 dshape_x.SetSize(p + 1);
581 dshape_y.SetSize(p + 1);
582 dshape_l.SetSize(p + 1);
583 u.SetSize(dof);
584 du.SetSize(dof, dim);
585#else
586 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
587#endif
588
589 for (int o = 0, j = 0; j <= p; j++)
590 for (int i = 0; i + j <= p; i++)
591 {
592 real_t w = op[i] + op[j] + op[p-i-j];
593 Nodes.IntPoint(o++).Set2(op[i]/w, op[j]/w);
594 }
595
596 DenseMatrix T(dof);
597 for (int k = 0; k < dof; k++)
598 {
600 poly1d.CalcBasis(p, ip.x, shape_x);
601 poly1d.CalcBasis(p, ip.y, shape_y);
602 poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
603
604 for (int o = 0, j = 0; j <= p; j++)
605 for (int i = 0; i + j <= p; i++)
606 {
607 T(o++, k) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
608 }
609 }
610
611 Ti.Factor(T);
612 // mfem::out << "L2_TriangleElement(" << p << ") : "; Ti.TestInversion();
613}
614
616 Vector &shape) const
617{
618 const int p = order;
619
620#ifdef MFEM_THREAD_SAFE
621 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1), u(dof);
622#endif
623
624 poly1d.CalcBasis(p, ip.x, shape_x);
625 poly1d.CalcBasis(p, ip.y, shape_y);
626 poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
627
628 for (int o = 0, j = 0; j <= p; j++)
629 for (int i = 0; i + j <= p; i++)
630 {
631 u(o++) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
632 }
633
634 Ti.Mult(u, shape);
635}
636
638 DenseMatrix &dshape) const
639{
640 const int p = order;
641
642#ifdef MFEM_THREAD_SAFE
643 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
644 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
645 DenseMatrix du(dof, dim);
646#endif
647
648 poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
649 poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
650 poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l, dshape_l);
651
652 for (int o = 0, j = 0; j <= p; j++)
653 for (int i = 0; i + j <= p; i++)
654 {
655 int k = p - i - j;
656 du(o,0) = ((dshape_x(i)* shape_l(k)) -
657 ( shape_x(i)*dshape_l(k)))*shape_y(j);
658 du(o,1) = ((dshape_y(j)* shape_l(k)) -
659 ( shape_y(j)*dshape_l(k)))*shape_x(i);
660 o++;
661 }
662
663 Ti.Mult(du, dshape);
664}
665
666void L2_TriangleElement::ProjectDelta(int vertex, Vector &dofs) const
667{
668 switch (vertex)
669 {
670 case 0:
671 for (int i = 0; i < dof; i++)
672 {
673 const IntegrationPoint &ip = Nodes.IntPoint(i);
674 dofs[i] = pow(1.0 - ip.x - ip.y, order);
675 }
676 break;
677 case 1:
678 for (int i = 0; i < dof; i++)
679 {
680 const IntegrationPoint &ip = Nodes.IntPoint(i);
681 dofs[i] = pow(ip.x, order);
682 }
683 break;
684 case 2:
685 for (int i = 0; i < dof; i++)
686 {
687 const IntegrationPoint &ip = Nodes.IntPoint(i);
688 dofs[i] = pow(ip.y, order);
689 }
690 break;
691 }
692}
693
694
696 : NodalFiniteElement(3, Geometry::TETRAHEDRON, ((p + 1)*(p + 2)*(p + 3))/6,
697 p, FunctionSpace::Pk)
698{
699 const real_t *op = poly1d.OpenPoints(p, VerifyOpen(btype));
700
701#ifndef MFEM_THREAD_SAFE
702 shape_x.SetSize(p + 1);
703 shape_y.SetSize(p + 1);
704 shape_z.SetSize(p + 1);
705 shape_l.SetSize(p + 1);
706 dshape_x.SetSize(p + 1);
707 dshape_y.SetSize(p + 1);
708 dshape_z.SetSize(p + 1);
709 dshape_l.SetSize(p + 1);
710 u.SetSize(dof);
711 du.SetSize(dof, dim);
712#else
713 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
714#endif
715
716 for (int o = 0, k = 0; k <= p; k++)
717 for (int j = 0; j + k <= p; j++)
718 for (int i = 0; i + j + k <= p; i++)
719 {
720 real_t w = op[i] + op[j] + op[k] + op[p-i-j-k];
721 Nodes.IntPoint(o++).Set3(op[i]/w, op[j]/w, op[k]/w);
722 }
723
724 DenseMatrix T(dof);
725 for (int m = 0; m < dof; m++)
726 {
728 poly1d.CalcBasis(p, ip.x, shape_x);
729 poly1d.CalcBasis(p, ip.y, shape_y);
730 poly1d.CalcBasis(p, ip.z, shape_z);
731 poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
732
733 for (int o = 0, k = 0; k <= p; k++)
734 for (int j = 0; j + k <= p; j++)
735 for (int i = 0; i + j + k <= p; i++)
736 {
737 T(o++, m) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
738 }
739 }
740
741 Ti.Factor(T);
742 // mfem::out << "L2_TetrahedronElement(" << p << ") : "; Ti.TestInversion();
743}
744
746 Vector &shape) const
747{
748 const int p = order;
749
750#ifdef MFEM_THREAD_SAFE
751 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
752 Vector u(dof);
753#endif
754
755 poly1d.CalcBasis(p, ip.x, shape_x);
756 poly1d.CalcBasis(p, ip.y, shape_y);
757 poly1d.CalcBasis(p, ip.z, shape_z);
758 poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
759
760 for (int o = 0, k = 0; k <= p; k++)
761 for (int j = 0; j + k <= p; j++)
762 for (int i = 0; i + j + k <= p; i++)
763 {
764 u(o++) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
765 }
766
767 Ti.Mult(u, shape);
768}
769
771 DenseMatrix &dshape) const
772{
773 const int p = order;
774
775#ifdef MFEM_THREAD_SAFE
776 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
777 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
778 DenseMatrix du(dof, dim);
779#endif
780
781 poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
782 poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
783 poly1d.CalcBasis(p, ip.z, shape_z, dshape_z);
784 poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l, dshape_l);
785
786 for (int o = 0, k = 0; k <= p; k++)
787 for (int j = 0; j + k <= p; j++)
788 for (int i = 0; i + j + k <= p; i++)
789 {
790 int l = p - i - j - k;
791 du(o,0) = ((dshape_x(i)* shape_l(l)) -
792 ( shape_x(i)*dshape_l(l)))*shape_y(j)*shape_z(k);
793 du(o,1) = ((dshape_y(j)* shape_l(l)) -
794 ( shape_y(j)*dshape_l(l)))*shape_x(i)*shape_z(k);
795 du(o,2) = ((dshape_z(k)* shape_l(l)) -
796 ( shape_z(k)*dshape_l(l)))*shape_x(i)*shape_y(j);
797 o++;
798 }
799
800 Ti.Mult(du, dshape);
801}
802
803void L2_TetrahedronElement::ProjectDelta(int vertex, Vector &dofs) const
804{
805 switch (vertex)
806 {
807 case 0:
808 for (int i = 0; i < dof; i++)
809 {
810 const IntegrationPoint &ip = Nodes.IntPoint(i);
811 dofs[i] = pow(1.0 - ip.x - ip.y - ip.z, order);
812 }
813 break;
814 case 1:
815 for (int i = 0; i < dof; i++)
816 {
817 const IntegrationPoint &ip = Nodes.IntPoint(i);
818 dofs[i] = pow(ip.x, order);
819 }
820 break;
821 case 2:
822 for (int i = 0; i < dof; i++)
823 {
824 const IntegrationPoint &ip = Nodes.IntPoint(i);
825 dofs[i] = pow(ip.y, order);
826 }
827 break;
828 case 3:
829 for (int i = 0; i < dof; i++)
830 {
831 const IntegrationPoint &ip = Nodes.IntPoint(i);
832 dofs[i] = pow(ip.z, order);
833 }
834 break;
835 }
836}
837
838
839L2_WedgeElement::L2_WedgeElement(const int p, const int btype)
840 : NodalFiniteElement(3, Geometry::PRISM, ((p + 1)*(p + 1)*(p + 2))/2,
841 p, FunctionSpace::Qk),
842 TriangleFE(p, btype),
843 SegmentFE(p, btype)
844{
845#ifndef MFEM_THREAD_SAFE
846 t_shape.SetSize(TriangleFE.GetDof());
847 s_shape.SetSize(SegmentFE.GetDof());
848 t_dshape.SetSize(TriangleFE.GetDof(), 2);
849 s_dshape.SetSize(SegmentFE.GetDof(), 1);
850#endif
851
852 t_dof.SetSize(dof);
853 s_dof.SetSize(dof);
854
855 // Interior DoFs
856 int m=0;
857 for (int k=0; k<=p; k++)
858 {
859 int l=0;
860 for (int j=0; j<=p; j++)
861 {
862 for (int i=0; i<=j; i++)
863 {
864 t_dof[m] = l;
865 s_dof[m] = k;
866 l++; m++;
867 }
868 }
869 }
870
871 // Define Nodes
872 const IntegrationRule & t_Nodes = TriangleFE.GetNodes();
873 const IntegrationRule & s_Nodes = SegmentFE.GetNodes();
874 for (int i=0; i<dof; i++)
875 {
876 Nodes.IntPoint(i).x = t_Nodes.IntPoint(t_dof[i]).x;
877 Nodes.IntPoint(i).y = t_Nodes.IntPoint(t_dof[i]).y;
878 Nodes.IntPoint(i).z = s_Nodes.IntPoint(s_dof[i]).x;
879 }
880}
881
883 Vector &shape) const
884{
885#ifdef MFEM_THREAD_SAFE
886 Vector t_shape(TriangleFE.GetDof());
887 Vector s_shape(SegmentFE.GetDof());
888#endif
889
890 IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
891
892 TriangleFE.CalcShape(ip, t_shape);
893 SegmentFE.CalcShape(ipz, s_shape);
894
895 for (int i=0; i<dof; i++)
896 {
897 shape[i] = t_shape[t_dof[i]] * s_shape[s_dof[i]];
898 }
899}
900
902 DenseMatrix &dshape) const
903{
904#ifdef MFEM_THREAD_SAFE
905 Vector t_shape(TriangleFE.GetDof());
906 DenseMatrix t_dshape(TriangleFE.GetDof(), 2);
907 Vector s_shape(SegmentFE.GetDof());
908 DenseMatrix s_dshape(SegmentFE.GetDof(), 1);
909#endif
910
911 IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
912
913 TriangleFE.CalcShape(ip, t_shape);
914 TriangleFE.CalcDShape(ip, t_dshape);
915 SegmentFE.CalcShape(ipz, s_shape);
916 SegmentFE.CalcDShape(ipz, s_dshape);
917
918 for (int i=0; i<dof; i++)
919 {
920 dshape(i, 0) = t_dshape(t_dof[i],0) * s_shape[s_dof[i]];
921 dshape(i, 1) = t_dshape(t_dof[i],1) * s_shape[s_dof[i]];
922 dshape(i, 2) = t_shape[t_dof[i]] * s_dshape(s_dof[i],0);
923 }
924}
925
927 : NodalFiniteElement(3, Geometry::PYRAMID, ((p + 1)*(p + 1)*(p + 1)),
928 p, FunctionSpace::Uk)
929{
930 const real_t *op = poly1d.OpenPoints(p, VerifyOpen(btype));
931
932 // These basis functions are not independent on a closed set of
933 // interpolation points when p >= 1. For this reason we force the points
934 // to be open in the z direction whenever closed points are requested.
935 // This should be regarded as a limitation of this choice of basis function.
936 // If a truly closed set of points is needed consider using
937 // L2_BergotPyramidElement instead.
938 real_t a = 1.0;
939 if (IsClosedType(btype) && p > 0)
940 {
942 }
943
944
945#ifndef MFEM_THREAD_SAFE
946 shape_x.SetSize(p + 1);
947 shape_y.SetSize(p + 1);
948 shape_z.SetSize(p + 1);
949 dshape_x.SetSize(p + 1);
950 dshape_y.SetSize(p + 1);
951 dshape_z.SetSize(p + 1);
952 u.SetSize(dof);
953 du.SetSize(dof, dim);
954#else
955 Vector shape_x(p + 1);
956 Vector shape_y(p + 1);
957 Vector shape_z(p + 1);
958#endif
959
960 int o = 0;
961 for (int k = 0; k <= p; k++)
962 for (int j = 0; j <= p; j++)
963 for (int i = 0; i <= p; i++)
964 {
965 Nodes.IntPoint(o++).Set3(op[i] * (1.0 - a * op[k]),
966 op[j] * (1.0 - a * op[k]),
967 a * op[k]);
968 }
969
970 MFEM_ASSERT(o == dof,
971 "Number of nodes does not match the "
972 "number of degrees of freedom");
973 DenseMatrix T(dof);
974
975 for (int m = 0; m < dof; m++)
976 {
977 const IntegrationPoint &ip = Nodes.IntPoint(m);
978 real_t x = ip.x;
979 real_t y = ip.y;
980 real_t z = ip.z;
981 Vector xy({x,y});
982 CalcHomogenizedScaLegendre(p, mu0(z, xy, 1), mu1(z, xy, 1), shape_x);
983 CalcHomogenizedScaLegendre(p, mu0(z, xy, 2), mu1(z, xy, 2), shape_y);
984 CalcHomogenizedScaLegendre(p, mu0(z), mu1(z), shape_z);
985
986 o = 0;
987 for (int k = 0; k <= p; k++)
988 {
989 for (int j = 0; j <= p; j++)
990 {
991 for (int i = 0; i <= p; i++, o++)
992 {
993 T(o, m) = shape_x[i] * shape_y[j] * shape_z[k];
994 }
995 }
996 }
997 }
998
999 Ti.Factor(T);
1000}
1001
1003 Vector &shape) const
1004{
1005 const int p = order;
1006
1007#ifdef MFEM_THREAD_SAFE
1008 Vector shape_x(p + 1);
1009 Vector shape_y(p + 1);
1010 Vector shape_z(p + 1);
1011 Vector u(dof);
1012#endif
1013 real_t x = ip.x;
1014 real_t y = ip.y;
1015 real_t z = ip.z;
1016 Vector xy({x,y});
1017
1018 if (z < 1.0)
1019 {
1020 CalcHomogenizedScaLegendre(p, mu0(z, xy, 1), mu1(z, xy, 1), shape_x);
1021 CalcHomogenizedScaLegendre(p, mu0(z, xy, 2), mu1(z, xy, 2), shape_y);
1022 }
1023 else
1024 {
1025 shape_x = 0.0; shape_x(0) = 1.0;
1026 shape_y = 0.0; shape_y(0) = 1.0;
1027 }
1028 CalcHomogenizedScaLegendre(p, mu0(z), mu1(z), shape_z);
1029
1030 int o = 0;
1031 for (int k = 0; k <= p; k++)
1032 for (int j = 0; j <= p; j++)
1033 for (int i = 0; i <= p; i++, o++)
1034 {
1035 u[o] = shape_x[i] * shape_y[j] * shape_z[k];
1036 }
1037
1038 Ti.Mult(u, shape);
1039}
1040
1042 DenseMatrix &dshape) const
1043{
1044 const int p = order;
1045
1046#ifdef MFEM_THREAD_SAFE
1047 Vector shape_x(p + 1);
1048 Vector shape_y(p + 1);
1049 Vector shape_z(p + 1);
1050 Vector dshape_x(p + 1);
1051 Vector dshape_y(p + 1);
1052 Vector dshape_z(p + 1);
1053 DenseMatrix du(dof, dim);
1054#endif
1055
1056 Poly_1D::CalcLegendre(p, ip.x / (1.0 - ip.z), shape_x.GetData(),
1057 dshape_x.GetData());
1058 Poly_1D::CalcLegendre(p, ip.y / (1.0 - ip.z), shape_y.GetData(),
1059 dshape_y.GetData());
1060 Poly_1D::CalcLegendre(p, ip.z, shape_z.GetData(), dshape_z.GetData());
1061
1062 int o = 0;
1063 for (int k = 0; k <= p; k++)
1064 for (int j = 0; j <= p; j++)
1065 for (int i = 0; i <= p; i++, o++)
1066 {
1067 du(o, 0) = dshape_x[i] * shape_y[j] * shape_z[k] / (1.0 - ip.z);
1068 du(o, 1) = shape_x[i] * dshape_y[j] * shape_z[k] / (1.0 - ip.z);
1069 du(o, 2) = shape_x[i] * shape_y[j] * dshape_z[k] +
1070 (ip.x * dshape_x[i] * shape_y[j] +
1071 ip.y * shape_x[i] * dshape_y[j]) *
1072 shape_z[k] / pow(1.0 - ip.z, 2);
1073 }
1074 Ti.Mult(du, dshape);
1075}
1076
1078 : NodalFiniteElement(3, Geometry::PYRAMID, (p + 1)*(p + 2)*(2*p + 3)/6,
1079 p, FunctionSpace::Pk)
1080{
1081 const real_t *op = poly1d.OpenPoints(p, VerifyOpen(btype));
1082
1083#ifndef MFEM_THREAD_SAFE
1084 shape_x.SetSize(p + 1);
1085 shape_y.SetSize(p + 1);
1086 shape_z.SetSize(p + 1);
1087 dshape_x.SetSize(p + 1);
1088 dshape_y.SetSize(p + 1);
1089 dshape_z.SetSize(p + 1);
1090 dshape_z_dt.SetSize(p + 1);
1091 u.SetSize(dof);
1092 du.SetSize(dof, dim);
1093#else
1094 Vector shape_x(p + 1);
1095 Vector shape_y(p + 1);
1096 Vector shape_z(p + 1);
1097 Vector dshape_z_dt(p + 1);
1098#endif
1099
1100 int o = 0;
1101 for (int k = 0; k <= p; k++)
1102 for (int j = 0; j <= p - k; j++)
1103 {
1104 const real_t wjk = op[j] + op[k] + op[p-j-k];
1105 for (int i = 0; i <= p - k; i++)
1106 {
1107 const real_t wik = op[i] + op[k] + op[p-i-k];
1108 const real_t w = wik * wjk * op[p-k];
1109 Nodes.IntPoint(o++).Set3(op[i] * (op[j] + op[p-j-k]) / w,
1110 op[j] * (op[j] + op[p-j-k]) / w,
1111 op[k] * op[p-k] / w);
1112 }
1113 }
1114
1115 MFEM_ASSERT(o == dof,
1116 "Number of nodes does not match the "
1117 "number of degrees of freedom");
1118 DenseMatrix T(dof);
1119
1120 for (int m = 0; m < dof; m++)
1121 {
1122 const IntegrationPoint &ip = Nodes.IntPoint(m);
1123
1124 const real_t x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
1125 const real_t y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
1126 const real_t z = ip.z;
1127
1128 poly1d.CalcLegendre(p, x, shape_x.GetData());
1129 poly1d.CalcLegendre(p, y, shape_y.GetData());
1130
1131 o = 0;
1132 for (int i = 0; i <= p; i++)
1133 {
1134 for (int j = 0; j <= p; j++)
1135 {
1136 int maxij = std::max(i, j);
1137 FuentesPyramid::CalcScaledJacobi(p-maxij, 2.0 * (maxij + 1.0),
1138 z, 1.0, shape_z);
1139
1140 for (int k = 0; k <= p - maxij; k++)
1141 {
1142 T(o++, m) = shape_x(i) * shape_y(j) * shape_z(k) *
1143 pow(1.0 - ip.z, maxij);
1144 }
1145 }
1146 }
1147 }
1148
1149 Ti.Factor(T);
1150}
1151
1153 Vector &shape) const
1154{
1155 const int p = order;
1156
1157#ifdef MFEM_THREAD_SAFE
1158 Vector shape_x(p + 1);
1159 Vector shape_y(p + 1);
1160 Vector shape_z(p + 1);
1161 Vector u(dof);
1162#endif
1163
1164 const real_t x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
1165 const real_t y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
1166 const real_t z = ip.z;
1167
1168 poly1d.CalcLegendre(p, x, shape_x.GetData());
1169 poly1d.CalcLegendre(p, y, shape_y.GetData());
1170
1171 int o = 0;
1172 for (int i = 0; i <= p; i++)
1173 {
1174 for (int j = 0; j <= p; j++)
1175 {
1176 int maxij = std::max(i, j);
1177 FuentesPyramid::CalcScaledJacobi(p-maxij, 2.0 * (maxij + 1.0), z, 1.0,
1178 shape_z);
1179
1180 for (int k = 0; k <= p - maxij; k++)
1181 {
1182 u[o++] = shape_x(i) * shape_y(j) * shape_z(k) *
1183 pow(1.0 - ip.z, maxij);
1184 }
1185 }
1186 }
1187
1188 Ti.Mult(u, shape);
1189}
1190
1192 DenseMatrix &dshape) const
1193{
1194 const int p = order;
1195
1196#ifdef MFEM_THREAD_SAFE
1197 Vector shape_x(p + 1);
1198 Vector shape_y(p + 1);
1199 Vector shape_z(p + 1);
1200 Vector dshape_x(p + 1);
1201 Vector dshape_y(p + 1);
1202 Vector dshape_z(p + 1);
1203 Vector dshape_z_dt(p + 1);
1204 DenseMatrix du(dof, dim);
1205#endif
1206
1207 const real_t x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
1208 const real_t y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
1209 const real_t z = ip.z;
1210
1211 Poly_1D::CalcLegendre(p, x, shape_x.GetData(), dshape_x.GetData());
1212 Poly_1D::CalcLegendre(p, y, shape_y.GetData(), dshape_y.GetData());
1213
1214 int o = 0;
1215 for (int i = 0; i <= p; i++)
1216 {
1217 for (int j = 0; j <= p; j++)
1218 {
1219 int maxij = std::max(i, j);
1220 FuentesPyramid::CalcScaledJacobi(p-maxij, 2.0 * (maxij + 1.0), z, 1.0,
1221 shape_z, dshape_z, dshape_z_dt);
1222
1223 for (int k = 0; k <= p - maxij; k++, o++)
1224 {
1225 du(o,0) = dshape_x(i) * shape_y(j) * shape_z(k) *
1226 pow(1.0 - ip.z, maxij - 1);
1227 du(o,1) = shape_x(i) * dshape_y(j) * shape_z(k) *
1228 pow(1.0 - ip.z, maxij - 1);
1229 du(o,2) = shape_x(i) * shape_y(j) * dshape_z(k) *
1230 pow(1.0 - ip.z, maxij) +
1231 (ip.x * dshape_x(i) * shape_y(j) +
1232 ip.y * shape_x(i) * dshape_y(j)) *
1233 shape_z(k) * pow(1.0 - ip.z, maxij - 2) -
1234 ((maxij > 0) ? (maxij * shape_x(i) * shape_y(j) * shape_z(k) *
1235 pow(1.0 - ip.z, maxij - 1)) : 0.0);
1236 }
1237 }
1238 }
1239
1240 Ti.Mult(du, dshape);
1241}
1242
1243}
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
@ GaussLobatto
Closed type.
Definition fe_base.hpp:36
@ GaussLegendre
Open type.
Definition fe_base.hpp:35
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
real_t * Data() const
Returns the matrix data array.
Definition densemat.hpp:114
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:108
real_t Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint....
Definition eltrans.hpp:144
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition eltrans.hpp:106
Abstract class for all finite elements.
Definition fe_base.hpp:244
int dof
Number of degrees of freedom.
Definition fe_base.hpp:253
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
IntegrationRule Nodes
Definition fe_base.hpp:256
static bool IsClosedType(int b_type)
Return true if the BasisType of b_type is closed (has Quadrature1D points on the boundary).
Definition fe_base.hpp:618
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition fe_base.hpp:400
static int VerifyOpen(int b_type)
Ensure that the BasisType of b_type is open (doesn't have Quadrature1D points on the boundary).
Definition fe_base.hpp:645
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
Geometry::Type geom_type
Geometry::Type of the reference element.
Definition fe_base.hpp:249
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
int dim
Dimension of reference space.
Definition fe_base.hpp:246
static real_t mu0(real_t z)
static real_t mu1(real_t z)
static void CalcScaledJacobi(int p, real_t alpha, real_t x, real_t t, real_t *u)
Shifted and Scaled Jacobi Polynomials.
static void CalcHomogenizedScaLegendre(int p, real_t s0, real_t s1, real_t *u)
Describes the function space on each element.
Definition fe_base.hpp:226
Class for integration point with weight.
Definition intrules.hpp:35
void Set2(const real_t x1, const real_t x2)
Definition intrules.hpp:89
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
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.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_l2.cpp:1152
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_l2.cpp:1191
L2_BergotPyramidElement(const int p, const int btype=BasisType::GaussLegendre)
Construct the L2_PyramidElement of order p and BasisType btype.
Definition fe_l2.cpp:1077
L2_FuentesPyramidElement(const int p, const int btype=BasisType::GaussLegendre)
Construct the L2_PyramidElement of order p and BasisType btype.
Definition fe_l2.cpp:926
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_l2.cpp:1041
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_l2.cpp:1002
void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const override
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_l2.cpp:342
L2_HexahedronElement(const int p, const int btype=BasisType::GaussLegendre)
Construct the L2_HexahedronElement of order p and BasisType btype.
Definition fe_l2.cpp:298
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_l2.cpp:451
void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_l2.cpp:320
void ProjectDelta(int vertex, Vector &dofs) const override
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition fe_l2.cpp:367
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_l2.cpp:521
void ProjectDelta(int vertex, Vector &dofs) const override
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition fe_l2.cpp:142
L2_QuadrilateralElement(const int p, const int btype=BasisType::GaussLegendre)
Construct the L2_QuadrilateralElement of order p and BasisType btype.
Definition fe_l2.cpp:82
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_l2.cpp:190
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_l2.cpp:255
void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_l2.cpp:101
void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const override
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_l2.cpp:121
void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_l2.cpp:39
void ProjectDelta(int vertex, Vector &dofs) const override
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition fe_l2.cpp:58
L2_SegmentElement(const int p, const int btype=BasisType::GaussLegendre)
Construct the L2_SegmentElement of order p and BasisType btype.
Definition fe_l2.cpp:23
void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const override
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_l2.cpp:46
void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const override
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_l2.cpp:770
void ProjectDelta(int vertex, Vector &dofs) const override
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition fe_l2.cpp:803
L2_TetrahedronElement(const int p, const int btype=BasisType::GaussLegendre)
Construct the L2_TetrahedronElement of order p and BasisType btype.
Definition fe_l2.cpp:695
void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_l2.cpp:745
L2_TriangleElement(const int p, const int btype=BasisType::GaussLegendre)
Construct the L2_TriangleElement of order p and BasisType btype.
Definition fe_l2.cpp:570
void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const override
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_l2.cpp:637
void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_l2.cpp:615
void ProjectDelta(int vertex, Vector &dofs) const override
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition fe_l2.cpp:666
L2_WedgeElement(const int p, const int btype=BasisType::GaussLegendre)
Construct the L2_WedgeElement of order p and BasisType btype.
Definition fe_l2.cpp:839
void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const override
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_l2.cpp:901
void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_l2.cpp:882
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 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 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
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 real_t CalcDelta(const int p, const real_t x)
Evaluate a representation of a Delta function at point x.
Definition fe_base.hpp:1180
static void CalcLegendre(const int p, const real_t x, real_t *u)
Definition fe_base.cpp:2245
const real_t * OpenPoints(const int p, const int btype=BasisType::GaussLegendre, bool on_device=false)
Get coordinates of an open (GaussLegendre) set of points if degree p.
Definition fe_base.hpp:1113
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
Poly_1D::Basis & basis1d
Definition fe_base.hpp:1251
Vector data type.
Definition vector.hpp:82
void SetDataAndSize(real_t *d, int s)
Set the Vector data and size.
Definition vector.hpp:183
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 SetData(real_t *d)
Definition vector.hpp:176
real_t a
Definition lissajous.cpp:41
Linear1DFiniteElement SegmentFE
Definition segment.cpp:52
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
float real_t
Definition config.hpp:43
Poly_1D poly1d
Definition fe.cpp:28
MFEM_EXPORT Linear2DFiniteElement TriangleFE
Definition fe.cpp:32
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)