MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
fe_l2.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
39void L2_SegmentElement::CalcShape(const IntegrationPoint &ip,
40 Vector &shape) const
41{
42 basis1d.ScaleIntegrated(map_type == VALUE);
43 basis1d.Eval(ip.x, shape);
44}
45
46void L2_SegmentElement::CalcDShape(const IntegrationPoint &ip,
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
54 basis1d.ScaleIntegrated(map_type == VALUE);
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
82L2_QuadrilateralElement::L2_QuadrilateralElement(const int p, const int btype)
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
101void L2_QuadrilateralElement::CalcShape(const IntegrationPoint &ip,
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
110 basis1d.ScaleIntegrated(map_type == VALUE);
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
121void L2_QuadrilateralElement::CalcDShape(const IntegrationPoint &ip,
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
130 basis1d.ScaleIntegrated(map_type == VALUE);
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
142void L2_QuadrilateralElement::ProjectDelta(int vertex, Vector &dofs) const
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
190void L2_QuadrilateralElement::ProjectDiv(const FiniteElement &fe,
191 ElementTransformation &Trans,
192 DenseMatrix &div) const
193{
194 if (basis1d.IsIntegratedType())
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());
203 const real_t *gll_pts = poly1d.GetPoints(order+1, BasisType::GaussLobatto);
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
255void L2_QuadrilateralElement::Project(Coefficient &coeff,
256 ElementTransformation &Trans,
257 Vector &dofs) const
258{
259 if (basis1d.IsIntegratedType())
260 {
261 const IntegrationRule &ir = IntRules.Get(geom_type, order);
262 const real_t *gll_pts = poly1d.GetPoints(order+1, BasisType::GaussLobatto);
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
298L2_HexahedronElement::L2_HexahedronElement(const int p, const int btype)
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
320void L2_HexahedronElement::CalcShape(const IntegrationPoint &ip,
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
329 basis1d.ScaleIntegrated(map_type == VALUE);
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
342void L2_HexahedronElement::CalcDShape(const IntegrationPoint &ip,
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
352 basis1d.ScaleIntegrated(map_type == VALUE);
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
451void L2_HexahedronElement::ProjectDiv(const FiniteElement &fe,
452 ElementTransformation &Trans,
453 DenseMatrix &div) const
454{
455 if (basis1d.IsIntegratedType())
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());
464 const real_t *gll_pts = poly1d.GetPoints(order+1, BasisType::GaussLobatto);
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
521void L2_HexahedronElement::Project(Coefficient &coeff,
522 ElementTransformation &Trans,
523 Vector &dofs) const
524{
525 if (basis1d.IsIntegratedType())
526 {
527 const IntegrationRule &ir = IntRules.Get(geom_type, order);
528 const real_t *gll_pts = poly1d.GetPoints(order+1, BasisType::GaussLobatto);
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 {
599 IntegrationPoint &ip = Nodes.IntPoint(k);
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
615void L2_TriangleElement::CalcShape(const IntegrationPoint &ip,
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
637void L2_TriangleElement::CalcDShape(const IntegrationPoint &ip,
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
695L2_TetrahedronElement::L2_TetrahedronElement(const int p, const int btype)
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 {
727 IntegrationPoint &ip = Nodes.IntPoint(m);
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
745void L2_TetrahedronElement::CalcShape(const IntegrationPoint &ip,
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
770void L2_TetrahedronElement::CalcDShape(const IntegrationPoint &ip,
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
882void L2_WedgeElement::CalcShape(const IntegrationPoint &ip,
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
901void L2_WedgeElement::CalcDShape(const IntegrationPoint &ip,
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
926}
mfem::Array::SetSize
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
mfem::Array::Size
int Size() const
Return the logical size of the array.
Definition array.hpp:144
mfem::BasisType::GaussLobatto
@ GaussLobatto
Closed type.
Definition fe_base.hpp:32
mfem::Coefficient
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition coefficient.hpp:42
mfem::Coefficient::Eval
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
mfem::DenseMatrixInverse::Factor
void Factor()
Factor the current DenseMatrix, *a.
Definition densemat.cpp:4185
mfem::DenseMatrixInverse::Mult
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition densemat.cpp:4231
mfem::DenseMatrix
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
mfem::DenseMatrix::Data
real_t * Data() const
Returns the matrix data array.
Definition densemat.hpp:111
mfem::DenseMatrix::SetSize
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:105
mfem::ElementTransformation::Weight
real_t Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint....
Definition eltrans.hpp:131
mfem::ElementTransformation::SetIntPoint
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition eltrans.hpp:93
mfem::FiniteElement
Abstract class for all finite elements.
Definition fe_base.hpp:239
mfem::FiniteElement::dof
int dof
Number of degrees of freedom.
Definition fe_base.hpp:248
mfem::FiniteElement::GetOrder
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition fe_base.hpp:333
mfem::FiniteElement::Nodes
IntegrationRule Nodes
Definition fe_base.hpp:251
mfem::FiniteElement::GetNodes
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition fe_base.hpp:395
mfem::FiniteElement::VerifyOpen
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:639
mfem::FiniteElement::CalcDivShape
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
mfem::FiniteElement::geom_type
Geometry::Type geom_type
Geometry::Type of the reference element.
Definition fe_base.hpp:244
mfem::FiniteElement::GetDof
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition fe_base.hpp:329
mfem::FiniteElement::order
int order
Order/degree of the shape functions.
Definition fe_base.hpp:249
mfem::FiniteElement::dim
int dim
Dimension of reference space.
Definition fe_base.hpp:241
mfem::FunctionSpace
Describes the function space on each element.
Definition fe_base.hpp:222
mfem::IntegrationPoint
Class for integration point with weight.
Definition intrules.hpp:35
mfem::IntegrationPoint::Set2
void Set2(const real_t x1, const real_t x2)
Definition intrules.hpp:89
mfem::IntegrationPoint::Set3
void Set3(const real_t x1, const real_t x2, const real_t x3)
Definition intrules.hpp:79
mfem::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
mfem::IntegrationRule::IntPoint
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition intrules.hpp:259
mfem::IntegrationRules::Get
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition intrules.cpp:1006
mfem::L2_HexahedronElement::CalcShape
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:320
mfem::L2_HexahedronElement::CalcDShape
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:342
mfem::L2_HexahedronElement::ProjectDelta
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_l2.cpp:367
mfem::L2_HexahedronElement::L2_HexahedronElement
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
mfem::L2_HexahedronElement::Project
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_l2.cpp:521
mfem::L2_HexahedronElement::ProjectDiv
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_l2.cpp:451
mfem::L2_QuadrilateralElement::L2_QuadrilateralElement
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
mfem::L2_QuadrilateralElement::ProjectDelta
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_l2.cpp:142
mfem::L2_QuadrilateralElement::ProjectDiv
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_l2.cpp:190
mfem::L2_QuadrilateralElement::Project
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_l2.cpp:255
mfem::L2_QuadrilateralElement::CalcDShape
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:121
mfem::L2_QuadrilateralElement::CalcShape
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:101
mfem::L2_SegmentElement::CalcShape
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:39
mfem::L2_SegmentElement::ProjectDelta
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_l2.cpp:58
mfem::L2_SegmentElement::L2_SegmentElement
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
mfem::L2_SegmentElement::CalcDShape
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:46
mfem::L2_TetrahedronElement::CalcShape
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:745
mfem::L2_TetrahedronElement::CalcDShape
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:770
mfem::L2_TetrahedronElement::L2_TetrahedronElement
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
mfem::L2_TetrahedronElement::ProjectDelta
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_l2.cpp:803
mfem::L2_TriangleElement::L2_TriangleElement
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
mfem::L2_TriangleElement::ProjectDelta
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_l2.cpp:666
mfem::L2_TriangleElement::CalcShape
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:615
mfem::L2_TriangleElement::CalcDShape
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:637
mfem::L2_WedgeElement::L2_WedgeElement
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
mfem::L2_WedgeElement::CalcDShape
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:901
mfem::L2_WedgeElement::CalcShape
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:882
mfem::NodalFiniteElement
Class for standard nodal finite elements.
Definition fe_base.hpp:715
mfem::NodalFiniteElement::Project
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
mfem::NodalFiniteElement::ProjectDiv
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
mfem::Poly_1D::Basis::Eval
void Eval(const real_t x, Vector &u) const
Evaluate the basis functions at point x in [0,1].
Definition fe_base.cpp:1761
mfem::Poly_1D::Basis::IsIntegratedType
bool IsIntegratedType() const
Returns true if the basis is "integrated", false otherwise.
Definition fe_base.hpp:1035
mfem::Poly_1D::Basis::ScaleIntegrated
void ScaleIntegrated(bool scale_integrated_)
Set whether the "integrated" basis should be scaled by the subcell sizes. Has no effect for non-integ...
Definition fe_base.cpp:2018
mfem::Poly_1D::OpenPoints
const real_t * OpenPoints(const int p, const int btype=BasisType::GaussLegendre)
Get coordinates of an open (GaussLegendre) set of points if degree p.
Definition fe_base.hpp:1078
mfem::Poly_1D::CalcDelta
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:1139
mfem::Poly_1D::GetPoints
const real_t * GetPoints(const int p, const int btype)
Get the coordinates of the points of the given BasisType, btype.
Definition fe_base.cpp:2270
mfem::Poly_1D::CalcBasis
static void CalcBasis(const int p, const real_t x, real_t *u)
Evaluate the values of a hierarchical 1D basis at point x hierarchical = k-th basis function is degre...
Definition fe_base.hpp:1099
mfem::TensorBasisElement::basis1d
Poly_1D::Basis & basis1d
Definition fe_base.hpp:1200
mfem::Vector
Vector data type.
Definition vector.hpp:80
mfem::Vector::SetDataAndSize
void SetDataAndSize(real_t *d, int s)
Set the Vector data and size.
Definition vector.hpp:175
mfem::Vector::SetSize
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
mfem::Vector::SetData
void SetData(real_t *d)
Definition vector.hpp:168
fe_h1.hpp
fe_l2.hpp
mfem::SegmentFE
Linear1DFiniteElement SegmentFE
Definition segment.cpp:52
mfem::u
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
mfem::real_t
float real_t
Definition config.hpp:43
mfem::poly1d
Poly_1D poly1d
Definition fe.cpp:28
mfem::TriangleFE
MFEM_EXPORT Linear2DFiniteElement TriangleFE
Definition fe.cpp:32
mfem::IntRules
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition intrules.hpp:486
p
real_t p(const Vector &x, real_t t)
Definition navier_mms.cpp:53
doxygen 1.11.0