MFEM  v4.4.0
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 
17 namespace mfem
18 {
19 
20 using namespace std;
21 
22 L2_SegmentElement::L2_SegmentElement(const int p, const int btype)
23  : NodalTensorFiniteElement(1, p, VerifyOpen(btype), L2_DOF_MAP)
24 {
25  const double *op = poly1d.OpenPoints(p, btype);
26 
27 #ifndef MFEM_THREAD_SAFE
28  shape_x.SetSize(p + 1);
29  dshape_x.SetDataAndSize(NULL, p + 1);
30 #endif
31 
32  for (int i = 0; i <= p; i++)
33  {
34  Nodes.IntPoint(i).x = op[i];
35  }
36 }
37 
38 void L2_SegmentElement::CalcShape(const IntegrationPoint &ip,
39  Vector &shape) const
40 {
41  basis1d.Eval(ip.x, shape);
42 }
43 
44 void L2_SegmentElement::CalcDShape(const IntegrationPoint &ip,
45  DenseMatrix &dshape) const
46 {
47 #ifdef MFEM_THREAD_SAFE
48  Vector shape_x(dof), dshape_x(dshape.Data(), dof);
49 #else
50  dshape_x.SetData(dshape.Data());
51 #endif
52  basis1d.Eval(ip.x, shape_x, dshape_x);
53 }
54 
55 void L2_SegmentElement::ProjectDelta(int vertex, Vector &dofs) const
56 {
57  const int p = order;
58  const double *op = poly1d.OpenPoints(p, b_type);
59 
60  switch (vertex)
61  {
62  case 0:
63  for (int i = 0; i <= p; i++)
64  {
65  dofs(i) = poly1d.CalcDelta(p,(1.0 - op[i]));
66  }
67  break;
68 
69  case 1:
70  for (int i = 0; i <= p; i++)
71  {
72  dofs(i) = poly1d.CalcDelta(p,op[i]);
73  }
74  break;
75  }
76 }
77 
78 
79 L2_QuadrilateralElement::L2_QuadrilateralElement(const int p, const int btype)
80  : NodalTensorFiniteElement(2, p, VerifyOpen(btype), L2_DOF_MAP)
81 {
82  const double *op = poly1d.OpenPoints(p, b_type);
83 
84 #ifndef MFEM_THREAD_SAFE
85  shape_x.SetSize(p + 1);
86  shape_y.SetSize(p + 1);
87  dshape_x.SetSize(p + 1);
88  dshape_y.SetSize(p + 1);
89 #endif
90 
91  for (int o = 0, j = 0; j <= p; j++)
92  for (int i = 0; i <= p; i++)
93  {
94  Nodes.IntPoint(o++).Set2(op[i], op[j]);
95  }
96 }
97 
98 void L2_QuadrilateralElement::CalcShape(const IntegrationPoint &ip,
99  Vector &shape) const
100 {
101  const int p = order;
102 
103 #ifdef MFEM_THREAD_SAFE
104  Vector shape_x(p+1), shape_y(p+1);
105 #endif
106 
107  basis1d.Eval(ip.x, shape_x);
108  basis1d.Eval(ip.y, shape_y);
109 
110  for (int o = 0, j = 0; j <= p; j++)
111  for (int i = 0; i <= p; i++)
112  {
113  shape(o++) = shape_x(i)*shape_y(j);
114  }
115 }
116 
117 void L2_QuadrilateralElement::CalcDShape(const IntegrationPoint &ip,
118  DenseMatrix &dshape) const
119 {
120  const int p = order;
121 
122 #ifdef MFEM_THREAD_SAFE
123  Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
124 #endif
125 
126  basis1d.Eval(ip.x, shape_x, dshape_x);
127  basis1d.Eval(ip.y, shape_y, dshape_y);
128 
129  for (int o = 0, j = 0; j <= p; j++)
130  for (int i = 0; i <= p; i++)
131  {
132  dshape(o,0) = dshape_x(i)* shape_y(j);
133  dshape(o,1) = shape_x(i)*dshape_y(j); o++;
134  }
135 }
136 
137 void L2_QuadrilateralElement::ProjectDelta(int vertex, Vector &dofs) const
138 {
139  const int p = order;
140  const double *op = poly1d.OpenPoints(p, b_type);
141 
142 #ifdef MFEM_THREAD_SAFE
143  Vector shape_x(p+1), shape_y(p+1);
144 #endif
145 
146  for (int i = 0; i <= p; i++)
147  {
148  shape_x(i) = poly1d.CalcDelta(p,(1.0 - op[i]));
149  shape_y(i) = poly1d.CalcDelta(p,op[i]);
150  }
151 
152  switch (vertex)
153  {
154  case 0:
155  for (int o = 0, j = 0; j <= p; j++)
156  for (int i = 0; i <= p; i++)
157  {
158  dofs[o++] = shape_x(i)*shape_x(j);
159  }
160  break;
161  case 1:
162  for (int o = 0, j = 0; j <= p; j++)
163  for (int i = 0; i <= p; i++)
164  {
165  dofs[o++] = shape_y(i)*shape_x(j);
166  }
167  break;
168  case 2:
169  for (int o = 0, j = 0; j <= p; j++)
170  for (int i = 0; i <= p; i++)
171  {
172  dofs[o++] = shape_y(i)*shape_y(j);
173  }
174  break;
175  case 3:
176  for (int o = 0, j = 0; j <= p; j++)
177  for (int i = 0; i <= p; i++)
178  {
179  dofs[o++] = shape_x(i)*shape_y(j);
180  }
181  break;
182  }
183 }
184 
185 
186 L2_HexahedronElement::L2_HexahedronElement(const int p, const int btype)
187  : NodalTensorFiniteElement(3, p, VerifyOpen(btype), L2_DOF_MAP)
188 {
189  const double *op = poly1d.OpenPoints(p, btype);
190 
191 #ifndef MFEM_THREAD_SAFE
192  shape_x.SetSize(p + 1);
193  shape_y.SetSize(p + 1);
194  shape_z.SetSize(p + 1);
195  dshape_x.SetSize(p + 1);
196  dshape_y.SetSize(p + 1);
197  dshape_z.SetSize(p + 1);
198 #endif
199 
200  for (int o = 0, k = 0; k <= p; k++)
201  for (int j = 0; j <= p; j++)
202  for (int i = 0; i <= p; i++)
203  {
204  Nodes.IntPoint(o++).Set3(op[i], op[j], op[k]);
205  }
206 }
207 
208 void L2_HexahedronElement::CalcShape(const IntegrationPoint &ip,
209  Vector &shape) const
210 {
211  const int p = order;
212 
213 #ifdef MFEM_THREAD_SAFE
214  Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
215 #endif
216 
217  basis1d.Eval(ip.x, shape_x);
218  basis1d.Eval(ip.y, shape_y);
219  basis1d.Eval(ip.z, shape_z);
220 
221  for (int o = 0, k = 0; k <= p; k++)
222  for (int j = 0; j <= p; j++)
223  for (int i = 0; i <= p; i++)
224  {
225  shape(o++) = shape_x(i)*shape_y(j)*shape_z(k);
226  }
227 }
228 
229 void L2_HexahedronElement::CalcDShape(const IntegrationPoint &ip,
230  DenseMatrix &dshape) const
231 {
232  const int p = order;
233 
234 #ifdef MFEM_THREAD_SAFE
235  Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
236  Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
237 #endif
238 
239  basis1d.Eval(ip.x, shape_x, dshape_x);
240  basis1d.Eval(ip.y, shape_y, dshape_y);
241  basis1d.Eval(ip.z, shape_z, dshape_z);
242 
243  for (int o = 0, k = 0; k <= p; k++)
244  for (int j = 0; j <= p; j++)
245  for (int i = 0; i <= p; i++)
246  {
247  dshape(o,0) = dshape_x(i)* shape_y(j)* shape_z(k);
248  dshape(o,1) = shape_x(i)*dshape_y(j)* shape_z(k);
249  dshape(o,2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
250  }
251 }
252 
253 void L2_HexahedronElement::ProjectDelta(int vertex, Vector &dofs) const
254 {
255  const int p = order;
256  const double *op = poly1d.OpenPoints(p, b_type);
257 
258 #ifdef MFEM_THREAD_SAFE
259  Vector shape_x(p+1), shape_y(p+1);
260 #endif
261 
262  for (int i = 0; i <= p; i++)
263  {
264  shape_x(i) = poly1d.CalcDelta(p,(1.0 - op[i]));
265  shape_y(i) = poly1d.CalcDelta(p,op[i]);
266  }
267 
268  switch (vertex)
269  {
270  case 0:
271  for (int o = 0, k = 0; k <= p; k++)
272  for (int j = 0; j <= p; j++)
273  for (int i = 0; i <= p; i++)
274  {
275  dofs[o++] = shape_x(i)*shape_x(j)*shape_x(k);
276  }
277  break;
278  case 1:
279  for (int o = 0, k = 0; k <= p; k++)
280  for (int j = 0; j <= p; j++)
281  for (int i = 0; i <= p; i++)
282  {
283  dofs[o++] = shape_y(i)*shape_x(j)*shape_x(k);
284  }
285  break;
286  case 2:
287  for (int o = 0, k = 0; k <= p; k++)
288  for (int j = 0; j <= p; j++)
289  for (int i = 0; i <= p; i++)
290  {
291  dofs[o++] = shape_y(i)*shape_y(j)*shape_x(k);
292  }
293  break;
294  case 3:
295  for (int o = 0, k = 0; k <= p; k++)
296  for (int j = 0; j <= p; j++)
297  for (int i = 0; i <= p; i++)
298  {
299  dofs[o++] = shape_x(i)*shape_y(j)*shape_x(k);
300  }
301  break;
302  case 4:
303  for (int o = 0, k = 0; k <= p; k++)
304  for (int j = 0; j <= p; j++)
305  for (int i = 0; i <= p; i++)
306  {
307  dofs[o++] = shape_x(i)*shape_x(j)*shape_y(k);
308  }
309  break;
310  case 5:
311  for (int o = 0, k = 0; k <= p; k++)
312  for (int j = 0; j <= p; j++)
313  for (int i = 0; i <= p; i++)
314  {
315  dofs[o++] = shape_y(i)*shape_x(j)*shape_y(k);
316  }
317  break;
318  case 6:
319  for (int o = 0, k = 0; k <= p; k++)
320  for (int j = 0; j <= p; j++)
321  for (int i = 0; i <= p; i++)
322  {
323  dofs[o++] = shape_y(i)*shape_y(j)*shape_y(k);
324  }
325  break;
326  case 7:
327  for (int o = 0, k = 0; k <= p; k++)
328  for (int j = 0; j <= p; j++)
329  for (int i = 0; i <= p; i++)
330  {
331  dofs[o++] = shape_x(i)*shape_y(j)*shape_y(k);
332  }
333  break;
334  }
335 }
336 
337 
338 L2_TriangleElement::L2_TriangleElement(const int p, const int btype)
339  : NodalFiniteElement(2, Geometry::TRIANGLE, ((p + 1)*(p + 2))/2, p,
340  FunctionSpace::Pk)
341 {
342  const double *op = poly1d.OpenPoints(p, VerifyOpen(btype));
343 
344 #ifndef MFEM_THREAD_SAFE
345  shape_x.SetSize(p + 1);
346  shape_y.SetSize(p + 1);
347  shape_l.SetSize(p + 1);
348  dshape_x.SetSize(p + 1);
349  dshape_y.SetSize(p + 1);
350  dshape_l.SetSize(p + 1);
351  u.SetSize(dof);
352  du.SetSize(dof, dim);
353 #else
354  Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
355 #endif
356 
357  for (int o = 0, j = 0; j <= p; j++)
358  for (int i = 0; i + j <= p; i++)
359  {
360  double w = op[i] + op[j] + op[p-i-j];
361  Nodes.IntPoint(o++).Set2(op[i]/w, op[j]/w);
362  }
363 
364  DenseMatrix T(dof);
365  for (int k = 0; k < dof; k++)
366  {
367  IntegrationPoint &ip = Nodes.IntPoint(k);
368  poly1d.CalcBasis(p, ip.x, shape_x);
369  poly1d.CalcBasis(p, ip.y, shape_y);
370  poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
371 
372  for (int o = 0, j = 0; j <= p; j++)
373  for (int i = 0; i + j <= p; i++)
374  {
375  T(o++, k) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
376  }
377  }
378 
379  Ti.Factor(T);
380  // mfem::out << "L2_TriangleElement(" << p << ") : "; Ti.TestInversion();
381 }
382 
383 void L2_TriangleElement::CalcShape(const IntegrationPoint &ip,
384  Vector &shape) const
385 {
386  const int p = order;
387 
388 #ifdef MFEM_THREAD_SAFE
389  Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1), u(dof);
390 #endif
391 
392  poly1d.CalcBasis(p, ip.x, shape_x);
393  poly1d.CalcBasis(p, ip.y, shape_y);
394  poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
395 
396  for (int o = 0, j = 0; j <= p; j++)
397  for (int i = 0; i + j <= p; i++)
398  {
399  u(o++) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
400  }
401 
402  Ti.Mult(u, shape);
403 }
404 
405 void L2_TriangleElement::CalcDShape(const IntegrationPoint &ip,
406  DenseMatrix &dshape) const
407 {
408  const int p = order;
409 
410 #ifdef MFEM_THREAD_SAFE
411  Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
412  Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
413  DenseMatrix du(dof, dim);
414 #endif
415 
416  poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
417  poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
418  poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l, dshape_l);
419 
420  for (int o = 0, j = 0; j <= p; j++)
421  for (int i = 0; i + j <= p; i++)
422  {
423  int k = p - i - j;
424  du(o,0) = ((dshape_x(i)* shape_l(k)) -
425  ( shape_x(i)*dshape_l(k)))*shape_y(j);
426  du(o,1) = ((dshape_y(j)* shape_l(k)) -
427  ( shape_y(j)*dshape_l(k)))*shape_x(i);
428  o++;
429  }
430 
431  Ti.Mult(du, dshape);
432 }
433 
434 void L2_TriangleElement::ProjectDelta(int vertex, Vector &dofs) const
435 {
436  switch (vertex)
437  {
438  case 0:
439  for (int i = 0; i < dof; i++)
440  {
441  const IntegrationPoint &ip = Nodes.IntPoint(i);
442  dofs[i] = pow(1.0 - ip.x - ip.y, order);
443  }
444  break;
445  case 1:
446  for (int i = 0; i < dof; i++)
447  {
448  const IntegrationPoint &ip = Nodes.IntPoint(i);
449  dofs[i] = pow(ip.x, order);
450  }
451  break;
452  case 2:
453  for (int i = 0; i < dof; i++)
454  {
455  const IntegrationPoint &ip = Nodes.IntPoint(i);
456  dofs[i] = pow(ip.y, order);
457  }
458  break;
459  }
460 }
461 
462 
463 L2_TetrahedronElement::L2_TetrahedronElement(const int p, const int btype)
464  : NodalFiniteElement(3, Geometry::TETRAHEDRON, ((p + 1)*(p + 2)*(p + 3))/6,
465  p, FunctionSpace::Pk)
466 {
467  const double *op = poly1d.OpenPoints(p, VerifyNodal(VerifyOpen(btype)));
468 
469 #ifndef MFEM_THREAD_SAFE
470  shape_x.SetSize(p + 1);
471  shape_y.SetSize(p + 1);
472  shape_z.SetSize(p + 1);
473  shape_l.SetSize(p + 1);
474  dshape_x.SetSize(p + 1);
475  dshape_y.SetSize(p + 1);
476  dshape_z.SetSize(p + 1);
477  dshape_l.SetSize(p + 1);
478  u.SetSize(dof);
479  du.SetSize(dof, dim);
480 #else
481  Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
482 #endif
483 
484  for (int o = 0, k = 0; k <= p; k++)
485  for (int j = 0; j + k <= p; j++)
486  for (int i = 0; i + j + k <= p; i++)
487  {
488  double w = op[i] + op[j] + op[k] + op[p-i-j-k];
489  Nodes.IntPoint(o++).Set3(op[i]/w, op[j]/w, op[k]/w);
490  }
491 
492  DenseMatrix T(dof);
493  for (int m = 0; m < dof; m++)
494  {
495  IntegrationPoint &ip = Nodes.IntPoint(m);
496  poly1d.CalcBasis(p, ip.x, shape_x);
497  poly1d.CalcBasis(p, ip.y, shape_y);
498  poly1d.CalcBasis(p, ip.z, shape_z);
499  poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
500 
501  for (int o = 0, k = 0; k <= p; k++)
502  for (int j = 0; j + k <= p; j++)
503  for (int i = 0; i + j + k <= p; i++)
504  {
505  T(o++, m) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
506  }
507  }
508 
509  Ti.Factor(T);
510  // mfem::out << "L2_TetrahedronElement(" << p << ") : "; Ti.TestInversion();
511 }
512 
513 void L2_TetrahedronElement::CalcShape(const IntegrationPoint &ip,
514  Vector &shape) const
515 {
516  const int p = order;
517 
518 #ifdef MFEM_THREAD_SAFE
519  Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
520  Vector u(dof);
521 #endif
522 
523  poly1d.CalcBasis(p, ip.x, shape_x);
524  poly1d.CalcBasis(p, ip.y, shape_y);
525  poly1d.CalcBasis(p, ip.z, shape_z);
526  poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
527 
528  for (int o = 0, k = 0; k <= p; k++)
529  for (int j = 0; j + k <= p; j++)
530  for (int i = 0; i + j + k <= p; i++)
531  {
532  u(o++) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
533  }
534 
535  Ti.Mult(u, shape);
536 }
537 
538 void L2_TetrahedronElement::CalcDShape(const IntegrationPoint &ip,
539  DenseMatrix &dshape) const
540 {
541  const int p = order;
542 
543 #ifdef MFEM_THREAD_SAFE
544  Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
545  Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
546  DenseMatrix du(dof, dim);
547 #endif
548 
549  poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
550  poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
551  poly1d.CalcBasis(p, ip.z, shape_z, dshape_z);
552  poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l, dshape_l);
553 
554  for (int o = 0, k = 0; k <= p; k++)
555  for (int j = 0; j + k <= p; j++)
556  for (int i = 0; i + j + k <= p; i++)
557  {
558  int l = p - i - j - k;
559  du(o,0) = ((dshape_x(i)* shape_l(l)) -
560  ( shape_x(i)*dshape_l(l)))*shape_y(j)*shape_z(k);
561  du(o,1) = ((dshape_y(j)* shape_l(l)) -
562  ( shape_y(j)*dshape_l(l)))*shape_x(i)*shape_z(k);
563  du(o,2) = ((dshape_z(k)* shape_l(l)) -
564  ( shape_z(k)*dshape_l(l)))*shape_x(i)*shape_y(j);
565  o++;
566  }
567 
568  Ti.Mult(du, dshape);
569 }
570 
571 void L2_TetrahedronElement::ProjectDelta(int vertex, Vector &dofs) const
572 {
573  switch (vertex)
574  {
575  case 0:
576  for (int i = 0; i < dof; i++)
577  {
578  const IntegrationPoint &ip = Nodes.IntPoint(i);
579  dofs[i] = pow(1.0 - ip.x - ip.y - ip.z, order);
580  }
581  break;
582  case 1:
583  for (int i = 0; i < dof; i++)
584  {
585  const IntegrationPoint &ip = Nodes.IntPoint(i);
586  dofs[i] = pow(ip.x, order);
587  }
588  break;
589  case 2:
590  for (int i = 0; i < dof; i++)
591  {
592  const IntegrationPoint &ip = Nodes.IntPoint(i);
593  dofs[i] = pow(ip.y, order);
594  }
595  break;
596  case 3:
597  for (int i = 0; i < dof; i++)
598  {
599  const IntegrationPoint &ip = Nodes.IntPoint(i);
600  dofs[i] = pow(ip.z, order);
601  }
602  break;
603  }
604 }
605 
606 
607 L2_WedgeElement::L2_WedgeElement(const int p, const int btype)
608  : NodalFiniteElement(3, Geometry::PRISM, ((p + 1)*(p + 1)*(p + 2))/2,
609  p, FunctionSpace::Qk),
610  TriangleFE(p, btype),
611  SegmentFE(p, btype)
612 {
613 #ifndef MFEM_THREAD_SAFE
614  t_shape.SetSize(TriangleFE.GetDof());
615  s_shape.SetSize(SegmentFE.GetDof());
616  t_dshape.SetSize(TriangleFE.GetDof(), 2);
617  s_dshape.SetSize(SegmentFE.GetDof(), 1);
618 #endif
619 
620  t_dof.SetSize(dof);
621  s_dof.SetSize(dof);
622 
623  // Interior DoFs
624  int m=0;
625  for (int k=0; k<=p; k++)
626  {
627  int l=0;
628  for (int j=0; j<=p; j++)
629  {
630  for (int i=0; i<=j; i++)
631  {
632  t_dof[m] = l;
633  s_dof[m] = k;
634  l++; m++;
635  }
636  }
637  }
638 
639  // Define Nodes
640  const IntegrationRule & t_Nodes = TriangleFE.GetNodes();
641  const IntegrationRule & s_Nodes = SegmentFE.GetNodes();
642  for (int i=0; i<dof; i++)
643  {
644  Nodes.IntPoint(i).x = t_Nodes.IntPoint(t_dof[i]).x;
645  Nodes.IntPoint(i).y = t_Nodes.IntPoint(t_dof[i]).y;
646  Nodes.IntPoint(i).z = s_Nodes.IntPoint(s_dof[i]).x;
647  }
648 }
649 
650 void L2_WedgeElement::CalcShape(const IntegrationPoint &ip,
651  Vector &shape) const
652 {
653 #ifdef MFEM_THREAD_SAFE
654  Vector t_shape(TriangleFE.GetDof());
655  Vector s_shape(SegmentFE.GetDof());
656 #endif
657 
658  IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
659 
660  TriangleFE.CalcShape(ip, t_shape);
661  SegmentFE.CalcShape(ipz, s_shape);
662 
663  for (int i=0; i<dof; i++)
664  {
665  shape[i] = t_shape[t_dof[i]] * s_shape[s_dof[i]];
666  }
667 }
668 
669 void L2_WedgeElement::CalcDShape(const IntegrationPoint &ip,
670  DenseMatrix &dshape) const
671 {
672 #ifdef MFEM_THREAD_SAFE
673  Vector t_shape(TriangleFE.GetDof());
674  DenseMatrix t_dshape(TriangleFE.GetDof(), 2);
675  Vector s_shape(SegmentFE.GetDof());
676  DenseMatrix s_dshape(SegmentFE.GetDof(), 1);
677 #endif
678 
679  IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
680 
681  TriangleFE.CalcShape(ip, t_shape);
682  TriangleFE.CalcDShape(ip, t_dshape);
683  SegmentFE.CalcShape(ipz, s_shape);
684  SegmentFE.CalcDShape(ipz, s_dshape);
685 
686  for (int i=0; i<dof; i++)
687  {
688  dshape(i, 0) = t_dshape(t_dof[i],0) * s_shape[s_dof[i]];
689  dshape(i, 1) = t_dshape(t_dof[i],1) * s_shape[s_dof[i]];
690  dshape(i, 2) = t_shape[t_dof[i]] * s_dshape(s_dof[i],0);
691  }
692 }
693 
694 }
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:98
mfem::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
mfem::Poly_1D::CalcDelta
static double CalcDelta(const int p, const double x)
Evaluate a representation of a Delta function at point x.
Definition: fe_base.hpp:1109
mfem::Vector::SetSize
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:521
mfem::FiniteElement::dim
int dim
Dimension of reference space.
Definition: fe_base.hpp:238
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:117
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:22
mfem::DenseMatrix
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
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:38
mfem::DenseMatrixInverse::Factor
void Factor()
Factor the current DenseMatrix, *a.
Definition: densemat.cpp:3212
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:571
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:513
mfem::Poly_1D::CalcBasis
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:1085
mfem::NodalFiniteElement
Class for standard nodal finite elements.
Definition: fe_base.hpp:706
mfem::IntegrationRule::IntPoint
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
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:44
fe_l2.hpp
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:405
mfem::TensorBasisElement::basis1d
Poly_1D::Basis & basis1d
Definition: fe_base.hpp:1159
mfem::FiniteElement::GetNodes
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:390
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:253
mfem::FiniteElement::VerifyOpen
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
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:650
mfem::IntegrationPoint::Set2
void Set2(const double x1, const double x2)
Definition: intrules.hpp:80
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:383
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:79
mfem::Vector::SetData
void SetData(double *d)
Definition: vector.hpp:149
mfem::Poly_1D::OpenPoints
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:1064
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:607
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:463
mfem::ad::pow
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
mfem::IntegrationPoint::Set3
void Set3(const double x1, const double x2, const double x3)
Definition: intrules.hpp:70
mfem::FiniteElement::Nodes
IntegrationRule Nodes
Definition: fe_base.hpp:248
mfem::DenseMatrix::Data
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
p
double p(const Vector &x, double t)
Definition: navier_mms.cpp:53
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:55
mfem::FiniteElement::GetDof
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
mfem::Array::SetSize
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
mfem::DenseMatrixInverse::Mult
void Mult(const double *x, double *y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition: densemat.cpp:3252
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:538
mfem::Vector::SetDataAndSize
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:156
mfem::TriangleFE
Linear2DFiniteElement TriangleFE
Definition: fe.cpp:32
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:669
mfem::IntegrationPoint
Class for integration point with weight.
Definition: intrules.hpp:25
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:434
mfem::FiniteElement::dof
int dof
Number of degrees of freedom.
Definition: fe_base.hpp:245
mfem::FiniteElement::VerifyNodal
static int VerifyNodal(int b_type)
Ensure that the BasisType of b_type nodal (satisfies the interpolation property). ...
Definition: fe_base.hpp:620
mfem::Poly_1D::Basis::Eval
void Eval(const double x, Vector &u) const
Evaluate the basis functions at point x in [0,1].
Definition: fe_base.cpp:1622
mfem::Vector
Vector data type.
Definition: vector.hpp:60
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:229
mfem::FunctionSpace
Describes the function space on each element.
Definition: fe_base.hpp:218
fe_h1.hpp
mfem::DenseMatrix::SetSize
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
mfem::SegmentFE
Linear1DFiniteElement SegmentFE
Definition: segment.cpp:49
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:186
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:208
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:338
mfem::FiniteElement::order
int order
Order/degree of the shape functions.
Definition: fe_base.hpp:245
mfem::poly1d
Poly_1D poly1d
Definition: fe.cpp:28
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:137
Generated on Mon Mar 21 2022 16:49:49 for MFEM by   doxygen 1.8.5