MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
fe_l2.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 
18 namespace mfem
19 {
20 
21 using namespace std;
22 
23 L2_SegmentElement::L2_SegmentElement(const int p, const int btype)
24  : NodalTensorFiniteElement(1, p, VerifyOpen(btype), L2_DOF_MAP)
25 {
26  const double *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 
58 void L2_SegmentElement::ProjectDelta(int vertex, Vector &dofs) const
59 {
60  const int p = order;
61  const double *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 double *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 
142 void L2_QuadrilateralElement::ProjectDelta(int vertex, Vector &dofs) const
143 {
144  const int p = order;
145  const double *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());
203  const double *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 double 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 double 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  double w = ip.weight;
222  if (map_type == VALUE)
223  {
224  const double 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 double 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  {
262  const double *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 double 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 double 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 double val = coeff.Eval(Trans, ip);
281  double 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 double *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 
367 void L2_HexahedronElement::ProjectDelta(int vertex, Vector &dofs) const
368 {
369  const int p = order;
370  const double *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());
464  const double *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 double hz = gll_pts[iz+1] - gll_pts[iz];
470  for (int iy = 0; iy < order+1; ++iy)
471  {
472  const double 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 double 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  double w = ip.weight;
487  if (map_type == VALUE)
488  {
489  const double 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 double 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  {
528  const double *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 double hz = gll_pts[iz+1] - gll_pts[iz];
535  for (int iy = 0; iy < order+1; ++iy)
536  {
537  const double hy = gll_pts[iy+1] - gll_pts[iy];
538  for (int ix = 0; ix < order+1; ++ix)
539  {
540  const double 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 double val = coeff.Eval(Trans, ip);
551  double w = ip.weight;
552  if (map_type == INTEGRAL)
553  {
554  const double 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 
570 L2_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 double *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  double 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 
666 void 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 double *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  double 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 
803 void 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 
839 L2_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 
926 }
Abstract class for all finite elements.
Definition: fe_base.hpp:235
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
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
static double CalcDelta(const int p, const double x)
Evaluate a representation of a Delta function at point x.
Definition: fe_base.hpp:1110
void ScaleIntegrated(bool scale_integrated_)
Set whether the &quot;integrated&quot; basis should be scaled by the subcell sizes. Has no effect for non-integ...
Definition: fe_base.cpp:1879
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
bool IsIntegratedType() const
Returns true if the basis is &quot;integrated&quot;, false otherwise.
Definition: fe_base.hpp:1022
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
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:656
int dim
Dimension of reference space.
Definition: fe_base.hpp:238
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
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
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
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:327
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
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
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
void Factor()
Factor the current DenseMatrix, *a.
Definition: densemat.cpp:3878
Geometry::Type geom_type
Geometry::Type of the reference element.
Definition: fe_base.hpp:241
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
virtual void ProjectDiv(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &div) const
Compute the discrete divergence matrix from the given FiniteElement onto &#39;this&#39; FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
Definition: fe_l2.cpp:190
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
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
static void CalcBasis(const int p, const double x, double *u)
Evaluate the values of a hierarchical 1D basis at point x hierarchical = k-th basis function is degre...
Definition: fe_base.hpp:1086
const double * GetPoints(const int p, const int btype)
Get the coordinates of the points of the given BasisType, btype.
Definition: fe_base.cpp:2131
virtual void ProjectDiv(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &div) const
Compute the discrete divergence matrix from the given FiniteElement onto &#39;this&#39; FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
Definition: fe_l2.cpp:451
Class for standard nodal finite elements.
Definition: fe_base.hpp:706
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
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
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
Poly_1D::Basis & basis1d
Definition: fe_base.hpp:1160
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:390
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
static int VerifyOpen(int b_type)
Ensure that the BasisType of b_type is open (doesn&#39;t have Quadrature1D points on the boundary)...
Definition: fe_base.hpp:612
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:131
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
void Set2(const double x1, const double x2)
Definition: intrules.hpp:80
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
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 SetData(double *d)
Definition: vector.hpp:150
const double * 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:1065
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
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 Set3(const double x1, const double x2, const double x3)
Definition: intrules.hpp:70
IntegrationRule Nodes
Definition: fe_base.hpp:248
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
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
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
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:51
void Mult(const double *x, double *y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition: densemat.cpp:3924
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
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:157
Linear2DFiniteElement TriangleFE
Definition: fe.cpp:32
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
Class for integration point with weight.
Definition: intrules.hpp:25
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
int dof
Number of degrees of freedom.
Definition: fe_base.hpp:245
virtual void ProjectDiv(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &div) const
Compute the discrete divergence matrix from the given FiniteElement onto &#39;this&#39; FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
Definition: fe_base.cpp:805
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
void Eval(const double x, Vector &u) const
Evaluate the basis functions at point x in [0,1].
Definition: fe_base.cpp:1622
Vector data type.
Definition: vector.hpp:60
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
Describes the function space on each element.
Definition: fe_base.hpp:218
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
Linear1DFiniteElement SegmentFE
Definition: segment.cpp:49
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
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
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
int order
Order/degree of the shape functions.
Definition: fe_base.hpp:245
Poly_1D poly1d
Definition: fe.cpp:28
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