MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
fe_h1.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12// H1 Finite Element classes
13
14#include "fe_h1.hpp"
15
16namespace mfem
17{
18
19using namespace std;
20
21H1_SegmentElement::H1_SegmentElement(const int p, const int btype)
22 : NodalTensorFiniteElement(1, p, VerifyClosed(btype), H1_DOF_MAP)
23{
24 const real_t *cp = poly1d.ClosedPoints(p, b_type);
25
26#ifndef MFEM_THREAD_SAFE
27 shape_x.SetSize(p+1);
28 dshape_x.SetSize(p+1);
29 d2shape_x.SetSize(p+1);
30#endif
31
32 Nodes.IntPoint(0).x = cp[0];
33 Nodes.IntPoint(1).x = cp[p];
34 for (int i = 1; i < p; i++)
35 {
36 Nodes.IntPoint(i+1).x = cp[i];
37 }
38}
39
41 Vector &shape) const
42{
43 const int p = order;
44
45#ifdef MFEM_THREAD_SAFE
46 Vector shape_x(p+1);
47#endif
48
49 basis1d.Eval(ip.x, shape_x);
50
51 shape(0) = shape_x(0);
52 shape(1) = shape_x(p);
53 for (int i = 1; i < p; i++)
54 {
55 shape(i+1) = shape_x(i);
56 }
57}
58
60 DenseMatrix &dshape) const
61{
62 const int p = order;
63
64#ifdef MFEM_THREAD_SAFE
65 Vector shape_x(p+1), dshape_x(p+1);
66#endif
67
68 basis1d.Eval(ip.x, shape_x, dshape_x);
69
70 dshape(0,0) = dshape_x(0);
71 dshape(1,0) = dshape_x(p);
72 for (int i = 1; i < p; i++)
73 {
74 dshape(i+1,0) = dshape_x(i);
75 }
76}
77
79 DenseMatrix &Hessian) const
80{
81 const int p = order;
82
83#ifdef MFEM_THREAD_SAFE
84 Vector shape_x(p+1), dshape_x(p+1), d2shape_x(p+1);
85#endif
86
87 basis1d.Eval(ip.x, shape_x, dshape_x, d2shape_x);
88
89 Hessian(0,0) = d2shape_x(0);
90 Hessian(1,0) = d2shape_x(p);
91 for (int i = 1; i < p; i++)
92 {
93 Hessian(i+1,0) = d2shape_x(i);
94 }
95}
96
97void H1_SegmentElement::ProjectDelta(int vertex, Vector &dofs) const
98{
99 const int p = order;
100 const real_t *cp = poly1d.ClosedPoints(p, b_type);
101
102 switch (vertex)
103 {
104 case 0:
105 dofs(0) = poly1d.CalcDelta(p, (1.0 - cp[0]));
106 dofs(1) = poly1d.CalcDelta(p, (1.0 - cp[p]));
107 for (int i = 1; i < p; i++)
108 {
109 dofs(i+1) = poly1d.CalcDelta(p, (1.0 - cp[i]));
110 }
111 break;
112
113 case 1:
114 dofs(0) = poly1d.CalcDelta(p, cp[0]);
115 dofs(1) = poly1d.CalcDelta(p, cp[p]);
116 for (int i = 1; i < p; i++)
117 {
118 dofs(i+1) = poly1d.CalcDelta(p, cp[i]);
119 }
120 break;
121 }
122}
123
124
126 : NodalTensorFiniteElement(2, p, VerifyClosed(btype), H1_DOF_MAP)
127{
128 const real_t *cp = poly1d.ClosedPoints(p, b_type);
129
130#ifndef MFEM_THREAD_SAFE
131 const int p1 = p + 1;
132
133 shape_x.SetSize(p1);
134 shape_y.SetSize(p1);
135 dshape_x.SetSize(p1);
136 dshape_y.SetSize(p1);
137 d2shape_x.SetSize(p1);
138 d2shape_y.SetSize(p1);
139#endif
140
141 int o = 0;
142 for (int j = 0; j <= p; j++)
143 {
144 for (int i = 0; i <= p; i++)
145 {
146 Nodes.IntPoint(dof_map[o++]).Set2(cp[i], cp[j]);
147 }
148 }
149}
150
152 Vector &shape) const
153{
154 const int p = order;
155
156#ifdef MFEM_THREAD_SAFE
157 Vector shape_x(p+1), shape_y(p+1);
158#endif
159
160 basis1d.Eval(ip.x, shape_x);
161 basis1d.Eval(ip.y, shape_y);
162
163 for (int o = 0, j = 0; j <= p; j++)
164 for (int i = 0; i <= p; i++)
165 {
166 shape(dof_map[o++]) = shape_x(i)*shape_y(j);
167 }
168}
169
171 DenseMatrix &dshape) const
172{
173 const int p = order;
174
175#ifdef MFEM_THREAD_SAFE
176 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
177#endif
178
179 basis1d.Eval(ip.x, shape_x, dshape_x);
180 basis1d.Eval(ip.y, shape_y, dshape_y);
181
182 for (int o = 0, j = 0; j <= p; j++)
183 {
184 for (int i = 0; i <= p; i++)
185 {
186 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j);
187 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j); o++;
188 }
189 }
190}
191
193 DenseMatrix &Hessian) const
194{
195 const int p = order;
196
197#ifdef MFEM_THREAD_SAFE
198 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1),
199 d2shape_x(p+1), d2shape_y(p+1);
200#endif
201
202 basis1d.Eval(ip.x, shape_x, dshape_x, d2shape_x);
203 basis1d.Eval(ip.y, shape_y, dshape_y, d2shape_y);
204
205 for (int o = 0, j = 0; j <= p; j++)
206 {
207 for (int i = 0; i <= p; i++)
208 {
209 Hessian(dof_map[o],0) = d2shape_x(i)* shape_y(j);
210 Hessian(dof_map[o],1) = dshape_x(i)* dshape_y(j);
211 Hessian(dof_map[o],2) = shape_x(i)*d2shape_y(j); o++;
212 }
213 }
214}
215
217{
218 const int p = order;
219 const real_t *cp = poly1d.ClosedPoints(p, b_type);
220
221#ifdef MFEM_THREAD_SAFE
222 Vector shape_x(p+1), shape_y(p+1);
223#endif
224
225 for (int i = 0; i <= p; i++)
226 {
227 shape_x(i) = poly1d.CalcDelta(p, (1.0 - cp[i]));
228 shape_y(i) = poly1d.CalcDelta(p, cp[i]);
229 }
230
231 switch (vertex)
232 {
233 case 0:
234 for (int o = 0, j = 0; j <= p; j++)
235 for (int i = 0; i <= p; i++)
236 {
237 dofs(dof_map[o++]) = shape_x(i)*shape_x(j);
238 }
239 break;
240 case 1:
241 for (int o = 0, j = 0; j <= p; j++)
242 for (int i = 0; i <= p; i++)
243 {
244 dofs(dof_map[o++]) = shape_y(i)*shape_x(j);
245 }
246 break;
247 case 2:
248 for (int o = 0, j = 0; j <= p; j++)
249 for (int i = 0; i <= p; i++)
250 {
251 dofs(dof_map[o++]) = shape_y(i)*shape_y(j);
252 }
253 break;
254 case 3:
255 for (int o = 0, j = 0; j <= p; j++)
256 for (int i = 0; i <= p; i++)
257 {
258 dofs(dof_map[o++]) = shape_x(i)*shape_y(j);
259 }
260 break;
261 }
262}
263
264
266 : NodalTensorFiniteElement(3, p, VerifyClosed(btype), H1_DOF_MAP)
267{
268 const real_t *cp = poly1d.ClosedPoints(p, b_type);
269
270#ifndef MFEM_THREAD_SAFE
271 const int p1 = p + 1;
272
273 shape_x.SetSize(p1);
274 shape_y.SetSize(p1);
275 shape_z.SetSize(p1);
276 dshape_x.SetSize(p1);
277 dshape_y.SetSize(p1);
278 dshape_z.SetSize(p1);
279 d2shape_x.SetSize(p1);
280 d2shape_y.SetSize(p1);
281 d2shape_z.SetSize(p1);
282#endif
283
284 int o = 0;
285 for (int k = 0; k <= p; k++)
286 for (int j = 0; j <= p; j++)
287 for (int i = 0; i <= p; i++)
288 {
289 Nodes.IntPoint(dof_map[o++]).Set3(cp[i], cp[j], cp[k]);
290 }
291}
292
294 Vector &shape) const
295{
296 const int p = order;
297
298#ifdef MFEM_THREAD_SAFE
299 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
300#endif
301
302 basis1d.Eval(ip.x, shape_x);
303 basis1d.Eval(ip.y, shape_y);
304 basis1d.Eval(ip.z, shape_z);
305
306 for (int o = 0, k = 0; k <= p; k++)
307 for (int j = 0; j <= p; j++)
308 for (int i = 0; i <= p; i++)
309 {
310 shape(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_z(k);
311 }
312}
313
315 DenseMatrix &dshape) const
316{
317 const int p = order;
318
319#ifdef MFEM_THREAD_SAFE
320 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
321 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
322#endif
323
324 basis1d.Eval(ip.x, shape_x, dshape_x);
325 basis1d.Eval(ip.y, shape_y, dshape_y);
326 basis1d.Eval(ip.z, shape_z, dshape_z);
327
328 for (int o = 0, k = 0; k <= p; k++)
329 for (int j = 0; j <= p; j++)
330 for (int i = 0; i <= p; i++)
331 {
332 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j)* shape_z(k);
333 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j)* shape_z(k);
334 dshape(dof_map[o],2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
335 }
336}
337
339 DenseMatrix &Hessian) const
340{
341 const int p = order;
342
343#ifdef MFEM_THREAD_SAFE
344 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
345 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
346 Vector d2shape_x(p+1), d2shape_y(p+1), d2shape_z(p+1);
347#endif
348
349 basis1d.Eval(ip.x, shape_x, dshape_x, d2shape_x);
350 basis1d.Eval(ip.y, shape_y, dshape_y, d2shape_y);
351 basis1d.Eval(ip.z, shape_z, dshape_z, d2shape_z);
352
353 for (int o = 0, k = 0; k <= p; k++)
354 for (int j = 0; j <= p; j++)
355 for (int i = 0; i <= p; i++)
356 {
357 Hessian(dof_map[o],0) = d2shape_x(i)* shape_y(j)* shape_z(k);
358 Hessian(dof_map[o],1) = dshape_x(i)* dshape_y(j)* shape_z(k);
359 Hessian(dof_map[o],2) = dshape_x(i)* shape_y(j)* dshape_z(k);
360 Hessian(dof_map[o],3) = shape_x(i)*d2shape_y(j)* shape_z(k);
361 Hessian(dof_map[o],4) = shape_x(i)* dshape_y(j)* dshape_z(k);
362 Hessian(dof_map[o],5) = shape_x(i)* shape_y(j)*d2shape_z(k);
363 o++;
364 }
365}
366
367void H1_HexahedronElement::ProjectDelta(int vertex, Vector &dofs) const
368{
369 const int p = order;
370 const real_t *cp = poly1d.ClosedPoints(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 - cp[i]));
379 shape_y(i) = poly1d.CalcDelta(p, cp[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(dof_map[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(dof_map[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(dof_map[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(dof_map[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(dof_map[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(dof_map[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(dof_map[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(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_y(k);
446 }
447 break;
448 }
449}
450
451H1_TriangleElement::H1_TriangleElement(const int p, const int btype)
452 : NodalFiniteElement(2, Geometry::TRIANGLE, ((p + 1)*(p + 2))/2, p,
453 FunctionSpace::Pk)
454{
455 const real_t *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
456
457#ifndef MFEM_THREAD_SAFE
458 shape_x.SetSize(p + 1);
459 shape_y.SetSize(p + 1);
460 shape_l.SetSize(p + 1);
461 dshape_x.SetSize(p + 1);
462 dshape_y.SetSize(p + 1);
463 dshape_l.SetSize(p + 1);
464 ddshape_x.SetSize(p + 1);
465 ddshape_y.SetSize(p + 1);
466 ddshape_l.SetSize(p + 1);
467 u.SetSize(dof);
468 du.SetSize(dof, dim);
469 ddu.SetSize(dof, (dim * (dim + 1)) / 2 );
470#else
471 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
472#endif
473
474 int p2p3 = 2*p + 3;
475 auto idx = [p2p3](int i, int j) { return ((p2p3-j)*j)/2+i; };
477
478 // vertices
479 lex_ordering[idx(0,0)] = 0;
480 Nodes.IntPoint(0).Set2(cp[0], cp[0]);
481 lex_ordering[idx(p,0)] = 1;
482 Nodes.IntPoint(1).Set2(cp[p], cp[0]);
483 lex_ordering[idx(0,p)] = 2;
484 Nodes.IntPoint(2).Set2(cp[0], cp[p]);
485
486 // edges
487 int o = 3;
488 for (int i = 1; i < p; i++)
489 {
490 lex_ordering[idx(i,0)] = o;
491 Nodes.IntPoint(o++).Set2(cp[i], cp[0]);
492 }
493 for (int i = 1; i < p; i++)
494 {
495 lex_ordering[idx(p-i,i)] = o;
496 Nodes.IntPoint(o++).Set2(cp[p-i], cp[i]);
497 }
498 for (int i = 1; i < p; i++)
499 {
500 lex_ordering[idx(0,p-i)] = o;
501 Nodes.IntPoint(o++).Set2(cp[0], cp[p-i]);
502 }
503
504 // interior
505 for (int j = 1; j < p; j++)
506 for (int i = 1; i + j < p; i++)
507 {
508 const real_t w = cp[i] + cp[j] + cp[p-i-j];
509 lex_ordering[idx(i,j)] = o;
510 Nodes.IntPoint(o++).Set2(cp[i]/w, cp[j]/w);
511 }
512
513 DenseMatrix T(dof);
514 for (int k = 0; k < dof; k++)
515 {
517 poly1d.CalcBasis(p, ip.x, shape_x);
518 poly1d.CalcBasis(p, ip.y, shape_y);
519 poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
520
521 o = 0;
522 for (int j = 0; j <= p; j++)
523 for (int i = 0; i + j <= p; i++)
524 {
525 T(o++, k) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
526 }
527 }
528
529 Ti.Factor(T);
530 // mfem::out << "H1_TriangleElement(" << p << ") : "; Ti.TestInversion();
531}
532
534 Vector &shape) const
535{
536 const int p = order;
537
538#ifdef MFEM_THREAD_SAFE
539 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1), u(dof);
540#endif
541
542 poly1d.CalcBasis(p, ip.x, shape_x);
543 poly1d.CalcBasis(p, ip.y, shape_y);
544 poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
545
546 for (int o = 0, j = 0; j <= p; j++)
547 for (int i = 0; i + j <= p; i++)
548 {
549 u(o++) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
550 }
551
552 Ti.Mult(u, shape);
553}
554
556 DenseMatrix &dshape) const
557{
558 const int p = order;
559
560#ifdef MFEM_THREAD_SAFE
561 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
562 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
563 DenseMatrix du(dof, dim);
564#endif
565
566 poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
567 poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
568 poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l, dshape_l);
569
570 for (int o = 0, j = 0; j <= p; j++)
571 for (int i = 0; i + j <= p; i++)
572 {
573 int k = p - i - j;
574 du(o,0) = ((dshape_x(i)* shape_l(k)) -
575 ( shape_x(i)*dshape_l(k)))*shape_y(j);
576 du(o,1) = ((dshape_y(j)* shape_l(k)) -
577 ( shape_y(j)*dshape_l(k)))*shape_x(i);
578 o++;
579 }
580
581 Ti.Mult(du, dshape);
582}
583
585 DenseMatrix &ddshape) const
586{
587 const int p = order;
588#ifdef MFEM_THREAD_SAFE
589 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
590 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
591 Vector ddshape_x(p + 1), ddshape_y(p + 1), ddshape_l(p + 1);
592 DenseMatrix ddu(dof, dim);
593#endif
594
595 poly1d.CalcBasis(p, ip.x, shape_x, dshape_x, ddshape_x);
596 poly1d.CalcBasis(p, ip.y, shape_y, dshape_y, ddshape_y);
597 poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l, dshape_l, ddshape_l);
598
599 for (int o = 0, j = 0; j <= p; j++)
600 for (int i = 0; i + j <= p; i++)
601 {
602 int k = p - i - j;
603 // u_xx, u_xy, u_yy
604 ddu(o,0) = ((ddshape_x(i) * shape_l(k)) - 2. * (dshape_x(i) * dshape_l(k)) +
605 (shape_x(i) * ddshape_l(k))) * shape_y(j);
606 ddu(o,1) = (((shape_x(i) * ddshape_l(k)) - dshape_x(i) * dshape_l(k)) * shape_y(
607 j)) + (((dshape_x(i) * shape_l(k)) - (shape_x(i) * dshape_l(k))) * dshape_y(j));
608 ddu(o,2) = ((ddshape_y(j) * shape_l(k)) - 2. * (dshape_y(j) * dshape_l(k)) +
609 (shape_y(j) * ddshape_l(k))) * shape_x(i);
610 o++;
611 }
612
613 Ti.Mult(ddu, ddshape);
614}
615
616
618 : NodalFiniteElement(3, Geometry::TETRAHEDRON, ((p + 1)*(p + 2)*(p + 3))/6,
619 p, FunctionSpace::Pk)
620{
621 const real_t *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
622
623#ifndef MFEM_THREAD_SAFE
624 shape_x.SetSize(p + 1);
625 shape_y.SetSize(p + 1);
626 shape_z.SetSize(p + 1);
627 shape_l.SetSize(p + 1);
628 dshape_x.SetSize(p + 1);
629 dshape_y.SetSize(p + 1);
630 dshape_z.SetSize(p + 1);
631 dshape_l.SetSize(p + 1);
632 ddshape_x.SetSize(p + 1);
633 ddshape_y.SetSize(p + 1);
634 ddshape_z.SetSize(p + 1);
635 ddshape_l.SetSize(p + 1);
636 u.SetSize(dof);
637 du.SetSize(dof, dim);
638 ddu.SetSize(dof, (dim * (dim + 1)) / 2);
639#else
640 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
641#endif
642
643 auto tri = [](int k) { return (k*(k + 1))/2; };
644 auto tet = [](int k) { return (k*(k + 1)*(k + 2))/6; };
645 int ndof = tet(p+1);
646 auto idx = [tri, tet, p, ndof](int i, int j, int k)
647 {
648 return ndof - tet(p - k) - tri(p + 1 - k - j) + i;
649 };
650
652
653 // vertices
654 lex_ordering[idx(0,0,0)] = 0;
655 Nodes.IntPoint(0).Set3(cp[0], cp[0], cp[0]);
656 lex_ordering[idx(p,0,0)] = 1;
657 Nodes.IntPoint(1).Set3(cp[p], cp[0], cp[0]);
658 lex_ordering[idx(0,p,0)] = 2;
659 Nodes.IntPoint(2).Set3(cp[0], cp[p], cp[0]);
660 lex_ordering[idx(0,0,p)] = 3;
661 Nodes.IntPoint(3).Set3(cp[0], cp[0], cp[p]);
662
663 // edges (see Tetrahedron::edges in mesh/tetrahedron.cpp)
664 int o = 4;
665 for (int i = 1; i < p; i++) // (0,1)
666 {
667 lex_ordering[idx(i,0,0)] = o;
668 Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[0]);
669 }
670 for (int i = 1; i < p; i++) // (0,2)
671 {
672 lex_ordering[idx(0,i,0)] = o;
673 Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[0]);
674 }
675 for (int i = 1; i < p; i++) // (0,3)
676 {
677 lex_ordering[idx(0,0,i)] = o;
678 Nodes.IntPoint(o++).Set3(cp[0], cp[0], cp[i]);
679 }
680 for (int i = 1; i < p; i++) // (1,2)
681 {
682 lex_ordering[idx(p-i,i,0)] = o;
683 Nodes.IntPoint(o++).Set3(cp[p-i], cp[i], cp[0]);
684 }
685 for (int i = 1; i < p; i++) // (1,3)
686 {
687 lex_ordering[idx(p-i,0,i)] = o;
688 Nodes.IntPoint(o++).Set3(cp[p-i], cp[0], cp[i]);
689 }
690 for (int i = 1; i < p; i++) // (2,3)
691 {
692 lex_ordering[idx(0,p-i,i)] = o;
693 Nodes.IntPoint(o++).Set3(cp[0], cp[p-i], cp[i]);
694 }
695
696 // faces (see Mesh::GenerateFaces in mesh/mesh.cpp)
697 for (int j = 1; j < p; j++)
698 for (int i = 1; i + j < p; i++) // (1,2,3)
699 {
700 lex_ordering[idx(p-i-j,i,j)] = o;
701 real_t w = cp[i] + cp[j] + cp[p-i-j];
702 Nodes.IntPoint(o++).Set3(cp[p-i-j]/w, cp[i]/w, cp[j]/w);
703 }
704 for (int j = 1; j < p; j++)
705 for (int i = 1; i + j < p; i++) // (0,3,2)
706 {
707 lex_ordering[idx(0,j,i)] = o;
708 real_t w = cp[i] + cp[j] + cp[p-i-j];
709 Nodes.IntPoint(o++).Set3(cp[0], cp[j]/w, cp[i]/w);
710 }
711 for (int j = 1; j < p; j++)
712 for (int i = 1; i + j < p; i++) // (0,1,3)
713 {
714 lex_ordering[idx(i,0,j)] = o;
715 real_t w = cp[i] + cp[j] + cp[p-i-j];
716 Nodes.IntPoint(o++).Set3(cp[i]/w, cp[0], cp[j]/w);
717 }
718 for (int j = 1; j < p; j++)
719 for (int i = 1; i + j < p; i++) // (0,2,1)
720 {
721 lex_ordering[idx(j,i,0)] = o;
722 real_t w = cp[i] + cp[j] + cp[p-i-j];
723 Nodes.IntPoint(o++).Set3(cp[j]/w, cp[i]/w, cp[0]);
724 }
725
726 // interior
727 for (int k = 1; k < p; k++)
728 for (int j = 1; j + k < p; j++)
729 for (int i = 1; i + j + k < p; i++)
730 {
731 lex_ordering[idx(i,j,k)] = o;
732 real_t w = cp[i] + cp[j] + cp[k] + cp[p-i-j-k];
733 Nodes.IntPoint(o++).Set3(cp[i]/w, cp[j]/w, cp[k]/w);
734 }
735
736 DenseMatrix T(dof);
737 for (int m = 0; m < dof; m++)
738 {
740 poly1d.CalcBasis(p, ip.x, shape_x);
741 poly1d.CalcBasis(p, ip.y, shape_y);
742 poly1d.CalcBasis(p, ip.z, shape_z);
743 poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
744
745 o = 0;
746 for (int k = 0; k <= p; k++)
747 for (int j = 0; j + k <= p; j++)
748 for (int i = 0; i + j + k <= p; i++)
749 {
750 T(o++, m) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
751 }
752 }
753
754 Ti.Factor(T);
755 // mfem::out << "H1_TetrahedronElement(" << p << ") : "; Ti.TestInversion();
756}
757
759 Vector &shape) const
760{
761 const int p = order;
762
763#ifdef MFEM_THREAD_SAFE
764 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
765 Vector u(dof);
766#endif
767
768 poly1d.CalcBasis(p, ip.x, shape_x);
769 poly1d.CalcBasis(p, ip.y, shape_y);
770 poly1d.CalcBasis(p, ip.z, shape_z);
771 poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
772
773 for (int o = 0, k = 0; k <= p; k++)
774 for (int j = 0; j + k <= p; j++)
775 for (int i = 0; i + j + k <= p; i++)
776 {
777 u(o++) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
778 }
779
780 Ti.Mult(u, shape);
781}
782
784 DenseMatrix &dshape) const
785{
786 const int p = order;
787
788#ifdef MFEM_THREAD_SAFE
789 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
790 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
791 DenseMatrix du(dof, dim);
792#endif
793
794 poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
795 poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
796 poly1d.CalcBasis(p, ip.z, shape_z, dshape_z);
797 poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l, dshape_l);
798
799 for (int o = 0, k = 0; k <= p; k++)
800 for (int j = 0; j + k <= p; j++)
801 for (int i = 0; i + j + k <= p; i++)
802 {
803 int l = p - i - j - k;
804 du(o,0) = ((dshape_x(i)* shape_l(l)) -
805 ( shape_x(i)*dshape_l(l)))*shape_y(j)*shape_z(k);
806 du(o,1) = ((dshape_y(j)* shape_l(l)) -
807 ( shape_y(j)*dshape_l(l)))*shape_x(i)*shape_z(k);
808 du(o,2) = ((dshape_z(k)* shape_l(l)) -
809 ( shape_z(k)*dshape_l(l)))*shape_x(i)*shape_y(j);
810 o++;
811 }
812
813 Ti.Mult(du, dshape);
814}
815
817 DenseMatrix &ddshape) const
818{
819 const int p = order;
820
821#ifdef MFEM_THREAD_SAFE
822 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
823 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
824 Vector ddshape_x(p + 1), ddshape_y(p + 1), ddshape_z(p + 1), ddshape_l(p + 1);
825 DenseMatrix ddu(dof, ((dim + 1) * dim) / 2);
826#endif
827
828 poly1d.CalcBasis(p, ip.x, shape_x, dshape_x, ddshape_x);
829 poly1d.CalcBasis(p, ip.y, shape_y, dshape_y, ddshape_y);
830 poly1d.CalcBasis(p, ip.z, shape_z, dshape_z, ddshape_z);
831 poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l, dshape_l, ddshape_l);
832
833 for (int o = 0, k = 0; k <= p; k++)
834 for (int j = 0; j + k <= p; j++)
835 for (int i = 0; i + j + k <= p; i++)
836 {
837 // u_xx, u_xy, u_xz, u_yy, u_yz, u_zz
838 int l = p - i - j - k;
839 ddu(o,0) = ((ddshape_x(i) * shape_l(l)) - 2. * (dshape_x(i) * dshape_l(l)) +
840 (shape_x(i) * ddshape_l(l))) * shape_y(j) * shape_z(k);
841 ddu(o,1) = ((dshape_y(j) * ((dshape_x(i) * shape_l(l)) -
842 (shape_x(i) * dshape_l(l)))) +
843 (shape_y(j) * ((ddshape_l(l) * shape_x(i)) -
844 (dshape_x(i) * dshape_l(l)))))* shape_z(k);
845 ddu(o,2) = ((dshape_z(k) * ((dshape_x(i) * shape_l(l)) -
846 (shape_x(i) * dshape_l(l)))) +
847 (shape_z(k) * ((ddshape_l(l) * shape_x(i)) -
848 (dshape_x(i) * dshape_l(l)))))* shape_y(j);
849 ddu(o,3) = ((ddshape_y(j) * shape_l(l)) - 2. * (dshape_y(j) * dshape_l(l)) +
850 (shape_y(j) * ddshape_l(l))) * shape_x(i) * shape_z(k);
851 ddu(o,4) = ((dshape_z(k) * ((dshape_y(j) * shape_l(l)) -
852 (shape_y(j)*dshape_l(l))) ) +
853 (shape_z(k)* ((ddshape_l(l)*shape_y(j)) -
854 (dshape_y(j) * dshape_l(l)) ) ) )* shape_x(i);
855 ddu(o,5) = ((ddshape_z(k) * shape_l(l)) - 2. * (dshape_z(k) * dshape_l(l)) +
856 (shape_z(k) * ddshape_l(l))) * shape_y(j) * shape_x(i);
857 o++;
858 }
859 Ti.Mult(ddu, ddshape);
860}
861
862// TODO: use a FunctionSpace specific to wedges instead of Qk.
864 const int btype)
865 : NodalFiniteElement(3, Geometry::PRISM, ((p + 1)*(p + 1)*(p + 2))/2,
866 p, FunctionSpace::Qk),
867 TriangleFE(p, btype),
868 SegmentFE(p, btype)
869{
870#ifndef MFEM_THREAD_SAFE
871 t_shape.SetSize(TriangleFE.GetDof());
872 s_shape.SetSize(SegmentFE.GetDof());
873 t_dshape.SetSize(TriangleFE.GetDof(), 2);
874 s_dshape.SetSize(SegmentFE.GetDof(), 1);
875#endif
876
877 t_dof.SetSize(dof);
878 s_dof.SetSize(dof);
879
880 int p2p3 = 2*p + 3, ntri = ((p + 1)*(p + 2))/2;
881 auto idx = [p2p3,ntri](int i, int j, int k)
882 {
883 return k*ntri + ((p2p3-j)*j)/2+i;
884 };
885
887 int o = 0;
888
889 // Nodal DoFs
890 lex_ordering[idx(0,0,0)] = o++;
891 lex_ordering[idx(p,0,0)] = o++;
892 lex_ordering[idx(0,p,0)] = o++;
893 lex_ordering[idx(0,0,p)] = o++;
894 lex_ordering[idx(p,0,p)] = o++;
895 lex_ordering[idx(0,p,p)] = o++;
896 t_dof[0] = 0; s_dof[0] = 0;
897 t_dof[1] = 1; s_dof[1] = 0;
898 t_dof[2] = 2; s_dof[2] = 0;
899 t_dof[3] = 0; s_dof[3] = 1;
900 t_dof[4] = 1; s_dof[4] = 1;
901 t_dof[5] = 2; s_dof[5] = 1;
902
903 // Edge DoFs
904 int k = 0;
905 int ne = p-1;
906 for (int i=1; i<p; i++)
907 {
908 lex_ordering[idx(i,0,0)] = o + 0*ne + k;
909 lex_ordering[idx(p-i,i,0)] = o + 1*ne + k;
910 lex_ordering[idx(0,p-i,0)] = o + 2*ne + k;
911 lex_ordering[idx(i,0,p)] = o + 3*ne + k;
912 lex_ordering[idx(p-i,i,p)] = o + 4*ne + k;
913 lex_ordering[idx(0,p-i,p)] = o + 5*ne + k;
914 lex_ordering[idx(0,0,i)] = o + 6*ne + k;
915 lex_ordering[idx(p,0,i)] = o + 7*ne + k;
916 lex_ordering[idx(0,p,i)] = o + 8*ne + k;
917 t_dof[5 + 0 * ne + i] = 2 + 0 * ne + i; s_dof[5 + 0 * ne + i] = 0;
918 t_dof[5 + 1 * ne + i] = 2 + 1 * ne + i; s_dof[5 + 1 * ne + i] = 0;
919 t_dof[5 + 2 * ne + i] = 2 + 2 * ne + i; s_dof[5 + 2 * ne + i] = 0;
920 t_dof[5 + 3 * ne + i] = 2 + 0 * ne + i; s_dof[5 + 3 * ne + i] = 1;
921 t_dof[5 + 4 * ne + i] = 2 + 1 * ne + i; s_dof[5 + 4 * ne + i] = 1;
922 t_dof[5 + 5 * ne + i] = 2 + 2 * ne + i; s_dof[5 + 5 * ne + i] = 1;
923 t_dof[5 + 6 * ne + i] = 0; s_dof[5 + 6 * ne + i] = i + 1;
924 t_dof[5 + 7 * ne + i] = 1; s_dof[5 + 7 * ne + i] = i + 1;
925 t_dof[5 + 8 * ne + i] = 2; s_dof[5 + 8 * ne + i] = i + 1;
926 ++k;
927 }
928 o += 9*ne;
929
930 // Triangular Face DoFs
931 k=0;
932 int nt = (p-1)*(p-2)/2;
933 for (int j=1; j<p; j++)
934 {
935 for (int i=1; i<p-j; i++)
936 {
937 int l = j - p + (((2 * p - 1) - i) * i) / 2;
938 lex_ordering[idx(i,j,0)] = o+l;
939 lex_ordering[idx(i,j,p)] = o+nt+k;
940 t_dof[6 + 9 * ne + k] = 3 * p + l; s_dof[6 + 9 * ne + k] = 0;
941 t_dof[6 + 9 * ne + nt + k] = 3 * p + k; s_dof[6 + 9 * ne + nt + k] = 1;
942 k++;
943 }
944 }
945 o += 2*nt;
946
947 // Quadrilateral Face DoFs
948 k=0;
949 int nq = (p-1)*(p-1);
950 for (int j=1; j<p; j++)
951 {
952 for (int i=1; i<p; i++)
953 {
954 lex_ordering[idx(i,0,j)] = o+k;
955 lex_ordering[idx(p-i,i,j)] = o+nq+k;
956 lex_ordering[idx(0,p-i,j)] = o+2*nq+k;
957
958 t_dof[6 + 9 * ne + 2 * nt + 0 * nq + k] = 2 + 0 * ne + i;
959 t_dof[6 + 9 * ne + 2 * nt + 1 * nq + k] = 2 + 1 * ne + i;
960 t_dof[6 + 9 * ne + 2 * nt + 2 * nq + k] = 2 + 2 * ne + i;
961
962 s_dof[6 + 9 * ne + 2 * nt + 0 * nq + k] = 1 + j;
963 s_dof[6 + 9 * ne + 2 * nt + 1 * nq + k] = 1 + j;
964 s_dof[6 + 9 * ne + 2 * nt + 2 * nq + k] = 1 + j;
965
966 k++;
967 }
968 }
969 o += 3*nq;
970
971 // Interior DoFs
972 int m=0;
973 for (k=1; k<p; k++)
974 {
975 int l=0;
976 for (int j=1; j<p; j++)
977 {
978 for (int i=1; i+j<p; i++)
979 {
980 lex_ordering[idx(i,j,k)] = o++;
981 t_dof[6 + 9 * ne + 2 * nt + 3 * nq + m] = 3 * p + l;
982 s_dof[6 + 9 * ne + 2 * nt + 3 * nq + m] = 1 + k;
983 l++; m++;
984 }
985 }
986 }
987
988 // Define Nodes
989 const IntegrationRule & t_Nodes = TriangleFE.GetNodes();
990 const IntegrationRule & s_Nodes = SegmentFE.GetNodes();
991 for (int i=0; i<dof; i++)
992 {
993 Nodes.IntPoint(i).x = t_Nodes.IntPoint(t_dof[i]).x;
994 Nodes.IntPoint(i).y = t_Nodes.IntPoint(t_dof[i]).y;
995 Nodes.IntPoint(i).z = s_Nodes.IntPoint(s_dof[i]).x;
996 }
997}
998
1000 Vector &shape) const
1001{
1002#ifdef MFEM_THREAD_SAFE
1003 Vector t_shape(TriangleFE.GetDof());
1004 Vector s_shape(SegmentFE.GetDof());
1005#endif
1006
1007 IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
1008
1009 TriangleFE.CalcShape(ip, t_shape);
1010 SegmentFE.CalcShape(ipz, s_shape);
1011
1012 for (int i=0; i<dof; i++)
1013 {
1014 shape[i] = t_shape[t_dof[i]] * s_shape[s_dof[i]];
1015 }
1016}
1017
1019 DenseMatrix &dshape) const
1020{
1021#ifdef MFEM_THREAD_SAFE
1022 Vector t_shape(TriangleFE.GetDof());
1023 DenseMatrix t_dshape(TriangleFE.GetDof(), 2);
1024 Vector s_shape(SegmentFE.GetDof());
1025 DenseMatrix s_dshape(SegmentFE.GetDof(), 1);
1026#endif
1027
1028 IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
1029
1030 TriangleFE.CalcShape(ip, t_shape);
1031 TriangleFE.CalcDShape(ip, t_dshape);
1032 SegmentFE.CalcShape(ipz, s_shape);
1033 SegmentFE.CalcDShape(ipz, s_dshape);
1034
1035 for (int i=0; i<dof; i++)
1036 {
1037 dshape(i, 0) = t_dshape(t_dof[i],0) * s_shape[s_dof[i]];
1038 dshape(i, 1) = t_dshape(t_dof[i],1) * s_shape[s_dof[i]];
1039 dshape(i, 2) = t_shape[t_dof[i]] * s_dshape(s_dof[i],0);
1040 }
1041}
1042
1044 : NodalFiniteElement(3, Geometry::PYRAMID,
1045 p * (p * p + 3) + 1, // Fuentes et al
1046 p, FunctionSpace::Uk)
1047{
1048 zmax = 0.0;
1049
1050 const real_t *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
1051
1052#ifndef MFEM_THREAD_SAFE
1053 tmp_i.SetSize(p + 1);
1054 tmp1_ij.SetSize(p + 1, p + 1);
1055 tmp2_ij.SetSize(p + 1, dim);
1056 tmp_ijk.SetSize(p + 1, p + 1, dim);
1057 tmp_u.SetSize(dof);
1058 tmp_du.SetSize(dof, dim);
1059#else
1060 Vector tmp_i(p + 1);
1061 DenseMatrix tmp1_ij(p + 1, p + 1);
1062#endif
1063
1064 // vertices
1065 Nodes.IntPoint(0).Set3(cp[0], cp[0], cp[0]);
1066 Nodes.IntPoint(1).Set3(cp[p], cp[0], cp[0]);
1067 Nodes.IntPoint(2).Set3(cp[p], cp[p], cp[0]);
1068 Nodes.IntPoint(3).Set3(cp[0], cp[p], cp[0]);
1069 Nodes.IntPoint(4).Set3(cp[0], cp[0], cp[p]);
1070
1071 // edges
1072 int o = 5;
1073 for (int i = 1; i < p; i++) // (0,1)
1074 {
1075 Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[0]);
1076 }
1077 for (int i = 1; i < p; i++) // (1,2)
1078 {
1079 Nodes.IntPoint(o++).Set3(cp[p], cp[i], cp[0]);
1080 }
1081 for (int i = 1; i < p; i++) // (3,2)
1082 {
1083 Nodes.IntPoint(o++).Set3(cp[i], cp[p], cp[0]);
1084 }
1085 for (int i = 1; i < p; i++) // (0,3)
1086 {
1087 Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[0]);
1088 }
1089 for (int i = 1; i < p; i++) // (0,4)
1090 {
1091 Nodes.IntPoint(o++).Set3(cp[0], cp[0], cp[i]);
1092 }
1093 for (int i = 1; i < p; i++) // (1,4)
1094 {
1095 Nodes.IntPoint(o++).Set3(cp[p-i], cp[0], cp[i]);
1096 }
1097 for (int i = 1; i < p; i++) // (2,4)
1098 {
1099 Nodes.IntPoint(o++).Set3(cp[p-i], cp[p-i], cp[i]);
1100 }
1101 for (int i = 1; i < p; i++) // (3,4)
1102 {
1103 Nodes.IntPoint(o++).Set3(cp[0], cp[p-i], cp[i]);
1104 }
1105
1106 // quadrilateral face
1107 for (int j = 1; j < p; j++)
1108 {
1109 for (int i = 1; i < p; i++)
1110 {
1111 Nodes.IntPoint(o++).Set3(cp[i], cp[p-j], cp[0]);
1112 }
1113 }
1114
1115 // triangular faces
1116 for (int j = 1; j < p; j++)
1117 for (int i = 1; i + j < p; i++) // (0,1,4)
1118 {
1119 real_t w = cp[i] + cp[j] + cp[p-i-j];
1120 Nodes.IntPoint(o++).Set3(cp[i]/w, cp[0], cp[j]/w);
1121 }
1122 for (int j = 1; j < p; j++)
1123 for (int i = 1; i + j < p; i++) // (1,2,4)
1124 {
1125 real_t w = cp[i] + cp[j] + cp[p-i-j];
1126 Nodes.IntPoint(o++).Set3((cp[i] + cp[p-i-j])/w, cp[i]/w, cp[j]/w);
1127 }
1128 for (int j = 1; j < p; j++)
1129 for (int i = 1; i + j < p; i++) // (2,3,4)
1130 {
1131 real_t w = cp[i] + cp[j] + cp[p-i-j];
1132 Nodes.IntPoint(o++).Set3(cp[p-i-j]/w, (cp[i] + cp[p-i-j])/w, cp[j]/w);
1133 }
1134 for (int j = 1; j < p; j++)
1135 for (int i = 1; i + j < p; i++) // (3,0,4)
1136 {
1137 real_t w = cp[i] + cp[j] + cp[p-i-j];
1138 Nodes.IntPoint(o++).Set3(cp[0], cp[p-i-j]/w, cp[j]/w);
1139 }
1140
1141 // Points based on Fuentes' interior bubbles
1142 for (int k = 1; k < p; k++)
1143 {
1144 for (int j = 1; j < p; j++)
1145 {
1146 for (int i = 1; i < p; i++)
1147 {
1148 Nodes.IntPoint(o++).Set3(cp[i] * (1.0 - cp[k]),
1149 cp[j] * (1.0 - cp[k]),
1150 cp[k]);
1151 }
1152 }
1153 }
1154
1155 MFEM_ASSERT(o == dof,
1156 "Number of nodes does not match the "
1157 "number of degrees of freedom");
1158 DenseMatrix T(dof);
1159
1160 for (int m = 0; m < dof; m++)
1161 {
1162 const IntegrationPoint &ip = Nodes.IntPoint(m);
1163 Vector col(T.GetColumn(m), dof);
1164 calcBasis(order, ip, tmp_i, tmp1_ij, col);
1165 }
1166
1167 Ti.Factor(T);
1168}
1169
1171 Vector &shape) const
1172{
1173 const int p = order;
1174
1175#ifdef MFEM_THREAD_SAFE
1176 Vector tmp_i(p + 1);
1177 Vector tmp_u(dof);
1178 DenseMatrix tmp1_ij(p + 1, p + 1);
1179#endif
1180
1181 calcBasis(p, ip, tmp_i, tmp1_ij, tmp_u);
1182
1183 Ti.Mult(tmp_u, shape);
1184}
1185
1187 DenseMatrix &dshape) const
1188{
1189 const int p = order;
1190
1191#ifdef MFEM_THREAD_SAFE
1192 Vector tmp_i(p + 1);
1193 DenseMatrix tmp1_ij(p + 1, p + 1);
1194 DenseMatrix tmp2_ij(p + 1, dim);
1195 DenseTensor tmp_ijk(p + 1, p + 1, dim);
1196 DenseMatrix tmp_du(dof, dim);
1197#endif
1198
1199 calcGradBasis(p, ip, tmp_i, tmp2_ij, tmp1_ij, tmp_ijk, tmp_du);
1200 Ti.Mult(tmp_du, dshape);
1201}
1202
1204 Vector &shape) const
1205{
1206 const int p = order;
1207
1208#ifdef MFEM_THREAD_SAFE
1209 Vector tmp_i(p + 1);
1210 DenseMatrix tmp1_ij(p + 1, p + 1);
1211#endif
1212
1213 calcBasis(p, ip, tmp_i, tmp1_ij, shape);
1214}
1215
1217 DenseMatrix &dshape) const
1218{
1219 const int p = order;
1220
1221#ifdef MFEM_THREAD_SAFE
1222 Vector tmp_i(p + 1);
1223 DenseMatrix tmp1_ij(p + 1, p + 1);
1224 DenseMatrix tmp2_ij(p + 1, dim);
1225 DenseTensor tmp_ijk(p + 1, p + 1, dim);
1226#endif
1227
1228 calcGradBasis(p, ip, tmp_i, tmp2_ij, tmp1_ij, tmp_ijk, dshape);
1229}
1230
1231void H1_FuentesPyramidElement::calcBasis(const int p,
1232 const IntegrationPoint &ip,
1233 Vector &phi_i, DenseMatrix &phi_ij,
1234 Vector &u) const
1235{
1236 real_t x = ip.x;
1237 real_t y = ip.y;
1238 real_t z = ip.z;
1239 Vector xy({x,y});
1240
1241 zmax = std::max(z, zmax);
1242
1243 real_t mu;
1244
1245 int o = 0;
1246
1247 // Vertices
1248 u[0] = lam1(x, y, z);
1249 u[1] = lam2(x, y, z);
1250 u[2] = lam3(x, y, z);
1251 u[3] = lam4(x, y, z);
1252 u[4] = lam5(x, y, z);
1253
1254 o += 5;
1255
1256 // Mixed edges (base edges)
1257 if (CheckZ(z) && p >= 2)
1258 {
1259 // (a,b) = (1,2), c = 0
1260 phi_E(p, nu01(z, xy, 1), phi_i);
1261 mu = mu0(z, xy, 2);
1262 for (int i = 2; i <= p; i++, o++)
1263 {
1264 u[o] = mu * phi_i[i];
1265 }
1266 // (a,b) = (1,2), c = 1
1267 mu = mu1(z, xy, 2);
1268 for (int i = 2; i <= p; i++, o++)
1269 {
1270 u[o] = mu * phi_i[i];
1271 }
1272 // (a,b) = (2,1), c = 0
1273 phi_E(p, nu01(z, xy, 2), phi_i);
1274 mu = mu0(z, xy, 1);
1275 for (int i = 2; i <= p; i++, o++)
1276 {
1277 u[o] = mu * phi_i[i];
1278 }
1279 // (a,b) = (2,1), c = 1
1280 mu = mu1(z, xy, 1);
1281 for (int i = 2; i <= p; i++, o++)
1282 {
1283 u[o] = mu * phi_i[i];
1284 }
1285 }
1286 else
1287 {
1288 for (int i = 0; i < 4 * (p - 1); i++, o++)
1289 {
1290 u[o] = 0.0;
1291 }
1292 }
1293
1294 // Triangle edges (upright edges)
1295 if (p >= 2)
1296 {
1297 phi_E(p, lam15(x, y, z), phi_i);
1298 for (int i = 2; i<= p; i++, o++)
1299 {
1300 u[o] = phi_i[i];
1301 }
1302 phi_E(p, lam25(x, y, z), phi_i);
1303 for (int i = 2; i<= p; i++, o++)
1304 {
1305 u[o] = phi_i[i];
1306 }
1307 phi_E(p, lam35(x, y, z), phi_i);
1308 for (int i = 2; i<= p; i++, o++)
1309 {
1310 u[o] = phi_i[i];
1311 }
1312 phi_E(p, lam45(x, y, z), phi_i);
1313 for (int i = 2; i<= p; i++, o++)
1314 {
1315 u[o] = phi_i[i];
1316 }
1317 }
1318
1319 // Quadrilateral face
1320 if (CheckZ(z) && p >= 2)
1321 {
1322 phi_Q(p, mu01(z, xy, 1), mu01(z, xy, 2), phi_ij);
1323 mu = mu0(z);
1324 for (int j = 2; j <= p; j++)
1325 {
1326 for (int i = 2; i <= p; i++, o++)
1327 {
1328 u[o] = mu * phi_ij(i,j);
1329 }
1330 }
1331 }
1332 else
1333 {
1334 for (int j = 2; j <= p; j++)
1335 {
1336 for (int i = 2; i <= p; i++, o++)
1337 {
1338 u[o] = 0.0;
1339 }
1340 }
1341 }
1342
1343 // Triangular faces
1344 if (CheckZ(z) && p >= 3)
1345 {
1346 // (a,b) = (1,2), c = 0
1347 phi_T(p, nu012(z, xy, 1), phi_ij);
1348 mu = mu0(z, xy, 2);
1349 for (int i = 2; i < p; i++)
1350 for (int j = 1; i + j <= p; j++, o++)
1351 {
1352 u[o] = mu * phi_ij(i,j);
1353 }
1354 // (a,b) = (1,2), c = 1
1355 mu = mu1(z, xy, 2);
1356 for (int i = 2; i < p; i++)
1357 for (int j = 1; i + j <= p; j++, o++)
1358 {
1359 u[o] = mu * phi_ij(i,j);
1360 }
1361 // (a,b) = (2,1), c = 0
1362 phi_T(p, nu012(z, xy, 2), phi_ij);
1363 mu = mu0(z, xy, 1);
1364 for (int i = 2; i < p; i++)
1365 for (int j = 1; i + j <= p; j++, o++)
1366 {
1367 u[o] = mu * phi_ij(i,j);
1368 }
1369 // (a,b) = (2,1), c = 1
1370 mu = mu1(z, xy, 1);
1371 for (int i = 2; i < p; i++)
1372 for (int j = 1; i + j <= p; j++, o++)
1373 {
1374 u[o] = mu * phi_ij(i,j);
1375 }
1376 }
1377 else
1378 {
1379 for (int i = 0; i < 2 * (p - 1) * (p - 2); i++, o++)
1380 {
1381 u[o] = 0.0;
1382 }
1383 }
1384
1385 // Interior
1386 if (CheckZ(z) && p >= 2)
1387 {
1388 phi_Q(p, mu01(z, xy, 1), mu01(z, xy, 2), phi_ij);
1389 phi_E(p, mu01(z), phi_i);
1390 for (int k = 2; k <= p; k++)
1391 {
1392 for (int j = 2; j <= p; j++)
1393 {
1394 for (int i = 2; i <= p; i++, o++)
1395 {
1396 u[o] = phi_ij(i,j) * phi_i(k);
1397 }
1398 }
1399 }
1400 }
1401 else
1402 {
1403 for (int i = 0; i < (p - 1) * (p - 1) * (p - 1); i++, o++)
1404 {
1405 u[o]= 0.0;
1406 }
1407 }
1408}
1409
1410void H1_FuentesPyramidElement::calcGradBasis(const int p,
1411 const IntegrationPoint &ip,
1412 Vector &phi_i,
1413 DenseMatrix &dphi_i,
1414 DenseMatrix &phi_ij,
1415 DenseTensor &dphi_ij,
1416 DenseMatrix &du) const
1417{
1418 real_t x = ip.x;
1419 real_t y = ip.y;
1420 real_t z = ip.z;
1421 Vector xy({x,y});
1422
1423 zmax = std::max(z, zmax);
1424
1425 real_t mu;
1426 Vector dmu(3);
1427 Vector dlam(3);
1428
1429 int o = 0;
1430
1431 // Vertices
1432 dlam = grad_lam1(x, y, z);
1433 for (int d=0; d<3; d++) { du(0, d) = dlam(d); }
1434 dlam = grad_lam2(x, y, z);
1435 for (int d=0; d<3; d++) { du(1, d) = dlam(d); }
1436 dlam = grad_lam3(x, y, z);
1437 for (int d=0; d<3; d++) { du(2, d) = dlam(d); }
1438 dlam = grad_lam4(x, y, z);
1439 for (int d=0; d<3; d++) { du(3, d) = dlam(d); }
1440 dlam = grad_lam5(x, y, z);
1441 for (int d=0; d<3; d++) { du(4, d) = dlam(d); }
1442
1443 o += 5;
1444
1445 // Mixed edges (base edges)
1446 if (CheckZ(z) && p >= 2)
1447 {
1448 // (a,b) = (1,2), c = 0
1449 phi_E(p, nu01(z, xy, 1), grad_nu01(z, xy, 1), phi_i, dphi_i);
1450 mu = mu0(z, xy, 2);
1451 dmu = grad_mu0(z, xy, 2);;
1452 for (int i = 2; i <= p; i++, o++)
1453 for (int d=0; d<3; d++)
1454 {
1455 du(o, d) = dmu(d) * phi_i[i] + mu * dphi_i(i, d);
1456 }
1457
1458 // (a,b) = (1,2), c = 1
1459 mu = mu1(z, xy, 2);
1460 dmu = grad_mu1(z, xy, 2);;
1461 for (int i = 2; i <= p; i++, o++)
1462 for (int d=0; d<3; d++)
1463 {
1464 du(o, d) = dmu(d) * phi_i[i] + mu * dphi_i(i, d);
1465 }
1466
1467 // (a,b) = (2,1), c = 0
1468 phi_E(p, nu01(z, xy, 2), grad_nu01(z, xy, 2), phi_i, dphi_i);
1469 mu = mu0(z, xy, 1);
1470 dmu = grad_mu0(z, xy, 1);;
1471 for (int i = 2; i <= p; i++, o++)
1472 for (int d=0; d<3; d++)
1473 {
1474 du(o, d) = dmu(d) * phi_i[i] + mu * dphi_i(i, d);
1475 }
1476
1477 // (a,b) = (2,1), c = 1
1478 mu = mu1(z, xy, 1);
1479 dmu = grad_mu1(z, xy, 1);;
1480 for (int i = 2; i <= p; i++, o++)
1481 for (int d=0; d<3; d++)
1482 {
1483 du(o, d) = dmu(d) * phi_i[i] + mu * dphi_i(i, d);
1484 }
1485 }
1486 else
1487 {
1488 for (int i = 0; i < 4 * (p - 1); i++, o++)
1489 for (int d=0; d<3; d++)
1490 {
1491 du(o, d) = 0.0;
1492 }
1493 }
1494
1495 // Triangle edges (upright edges)
1496 if (p >= 2)
1497 {
1498 phi_E(p, lam15(x, y, z), grad_lam15(x,y,z), phi_i, dphi_i);
1499 for (int i = 2; i<= p; i++, o++)
1500 for (int d=0; d<3; d++)
1501 {
1502 du(o, d) = dphi_i(i, d);
1503 }
1504
1505 phi_E(p, lam25(x, y, z), grad_lam25(x, y, z), phi_i, dphi_i);
1506 for (int i = 2; i<= p; i++, o++)
1507 for (int d=0; d<3; d++)
1508 {
1509 du(o, d) = dphi_i(i, d);
1510 }
1511
1512 phi_E(p, lam35(x, y, z), grad_lam35(x, y, z), phi_i, dphi_i);
1513 for (int i = 2; i<= p; i++, o++)
1514 for (int d=0; d<3; d++)
1515 {
1516 du(o, d) = dphi_i(i, d);
1517 }
1518
1519 phi_E(p, lam45(x, y, z), grad_lam45(x, y, z), phi_i, dphi_i);
1520 for (int i = 2; i<= p; i++, o++)
1521 for (int d=0; d<3; d++)
1522 {
1523 du(o, d) = dphi_i(i, d);
1524 }
1525 }
1526
1527 // Quadrilateral face
1528 if (CheckZ(z) && p >= 2)
1529 {
1530 phi_Q(p, mu01(z, xy, 1), grad_mu01(z, xy, 1),
1531 mu01(z, xy, 2), grad_mu01(z, xy, 2), phi_ij, dphi_ij);
1532 mu = mu0(z);
1533 dmu = grad_mu0(z);
1534 for (int j = 2; j <= p; j++)
1535 for (int i = 2; i <= p; i++, o++)
1536 for (int d=0; d<3; d++)
1537 {
1538 du(o, d) = dmu(d) * phi_ij(i, j) + mu * dphi_ij(i, j, d);
1539 }
1540 }
1541 else
1542 {
1543 for (int j = 2; j <= p; j++)
1544 for (int i = 2; i <= p; i++, o++)
1545 for (int d=0; d<3; d++)
1546 {
1547 du(o, d) = 0.0;
1548 }
1549 }
1550
1551 // Triangular faces
1552 if (CheckZ(z) && p >= 3)
1553 {
1554 // (a,b) = (1,2), c = 0
1555 phi_T(p, nu012(z, xy, 1), grad_nu012(z, xy, 1), phi_ij, dphi_ij);
1556 mu = mu0(z, xy, 2);
1557 dmu = grad_mu0(z, xy, 2);
1558 for (int i = 2; i < p; i++)
1559 for (int j = 1; i + j <= p; j++, o++)
1560 for (int d=0; d<3; d++)
1561 {
1562 du(o, d) = dmu(d) * phi_ij(i, j) + mu * dphi_ij(i, j, d);
1563 }
1564
1565 // (a,b) = (1,2), c = 1
1566 mu = mu1(z, xy, 2);
1567 dmu = grad_mu1(z, xy, 2);
1568 for (int i = 2; i < p; i++)
1569 for (int j = 1; i + j <= p; j++, o++)
1570 for (int d=0; d<3; d++)
1571 {
1572 du(o, d) = dmu(d) * phi_ij(i, j) + mu * dphi_ij(i, j, d);
1573 }
1574
1575 // (a,b) = (2,1), c = 0
1576 phi_T(p, nu012(z, xy, 2), grad_nu012(z, xy, 2), phi_ij, dphi_ij);
1577 mu = mu0(z, xy, 1);
1578 dmu = grad_mu0(z, xy, 1);
1579 for (int i = 2; i < p; i++)
1580 for (int j = 1; i + j <= p; j++, o++)
1581 for (int d=0; d<3; d++)
1582 {
1583 du(o, d) = dmu(d) * phi_ij(i, j) + mu * dphi_ij(i, j, d);
1584 }
1585
1586 // (a,b) = (2,1), c = 1
1587 mu = mu1(z, xy, 1);
1588 dmu = grad_mu1(z, xy, 1);
1589 for (int i = 2; i < p; i++)
1590 for (int j = 1; i + j <= p; j++, o++)
1591 for (int d=0; d<3; d++)
1592 {
1593 du(o, d) = dmu(d) * phi_ij(i, j) + mu * dphi_ij(i, j, d);
1594 }
1595 }
1596 else
1597 {
1598 for (int i = 0; i < 2 * (p - 1) * (p - 2); i++, o++)
1599 for (int d=0; d<3; d++)
1600 {
1601 du(o, d) = 0.0;
1602 }
1603 }
1604
1605 // Interior
1606 if (CheckZ(z) && p >= 2)
1607 {
1608 phi_Q(p, mu01(z, xy, 1), grad_mu01(z, xy, 1),
1609 mu01(z, xy, 2), grad_mu01(z, xy, 2), phi_ij, dphi_ij);
1610 phi_E(p, mu01(z), grad_mu01(z), phi_i, dphi_i);
1611 for (int k = 2; k <= p; k++)
1612 for (int j = 2; j <= p; j++)
1613 for (int i = 2; i <= p; i++, o++)
1614 for (int d=0; d<3; d++)
1615 du(o, d) = dphi_ij(i, j, d) * phi_i(k) +
1616 phi_ij(i, j) * dphi_i(k, d);
1617 }
1618 else
1619 {
1620 for (int i = 0; i < (p - 1) * (p - 1) * (p - 1); i++, o++)
1621 for (int d=0; d<3; d++)
1622 {
1623 du(o, d) = 0.0;
1624 }
1625 }
1626}
1627
1629 : NodalFiniteElement(3, Geometry::PYRAMID,
1630 (p + 1) * (p + 2) * (2 * p + 3) / 6, // Bergot (JSC)
1631 p, FunctionSpace::Uk)
1632{
1633 const real_t *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
1634
1635#ifndef MFEM_THREAD_SAFE
1636 shape_x.SetSize(p + 1);
1637 shape_y.SetSize(p + 1);
1638 shape_z.SetSize(p + 1);
1639 dshape_x.SetSize(p + 1);
1640 dshape_y.SetSize(p + 1);
1641 dshape_z.SetSize(p + 1);
1642 dshape_z_dt.SetSize(p + 1);
1643 ddshape_x.SetSize(p + 1);
1644 ddshape_y.SetSize(p + 1);
1645 ddshape_z.SetSize(p + 1);
1646 u.SetSize(dof);
1647 du.SetSize(dof, dim);
1648 ddu.SetSize(dof, (dim * (dim + 1)) / 2);
1649#else
1650 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1);
1651#endif
1652
1653 // vertices
1654 Nodes.IntPoint(0).Set3(cp[0], cp[0], cp[0]);
1655 Nodes.IntPoint(1).Set3(cp[p], cp[0], cp[0]);
1656 Nodes.IntPoint(2).Set3(cp[p], cp[p], cp[0]);
1657 Nodes.IntPoint(3).Set3(cp[0], cp[p], cp[0]);
1658 Nodes.IntPoint(4).Set3(cp[0], cp[0], cp[p]);
1659
1660 // edges
1661 int o = 5;
1662 for (int i = 1; i < p; i++) // (0,1)
1663 {
1664 Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[0]);
1665 }
1666 for (int i = 1; i < p; i++) // (1,2)
1667 {
1668 Nodes.IntPoint(o++).Set3(cp[p], cp[i], cp[0]);
1669 }
1670 for (int i = 1; i < p; i++) // (3,2)
1671 {
1672 Nodes.IntPoint(o++).Set3(cp[i], cp[p], cp[0]);
1673 }
1674 for (int i = 1; i < p; i++) // (0,3)
1675 {
1676 Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[0]);
1677 }
1678 for (int i = 1; i < p; i++) // (0,4)
1679 {
1680 Nodes.IntPoint(o++).Set3(cp[0], cp[0], cp[i]);
1681 }
1682 for (int i = 1; i < p; i++) // (1,4)
1683 {
1684 Nodes.IntPoint(o++).Set3(cp[p-i], cp[0], cp[i]);
1685 }
1686 for (int i = 1; i < p; i++) // (2,4)
1687 {
1688 Nodes.IntPoint(o++).Set3(cp[p-i], cp[p-i], cp[i]);
1689 }
1690 for (int i = 1; i < p; i++) // (3,4)
1691 {
1692 Nodes.IntPoint(o++).Set3(cp[0], cp[p-i], cp[i]);
1693 }
1694
1695 // quadrilateral face
1696 for (int j = 1; j < p; j++)
1697 {
1698 for (int i = 1; i < p; i++)
1699 {
1700 Nodes.IntPoint(o++).Set3(cp[i], cp[j], cp[0]);
1701 }
1702 }
1703
1704 // triangular faces
1705 for (int j = 1; j < p; j++)
1706 for (int i = 1; i + j < p; i++) // (0,1,4)
1707 {
1708 real_t w = cp[i] + cp[j] + cp[p-i-j];
1709 Nodes.IntPoint(o++).Set3(cp[i]/w, cp[0], cp[j]/w);
1710 }
1711 for (int j = 1; j < p; j++)
1712 for (int i = 1; i + j < p; i++) // (1,2,4)
1713 {
1714 real_t w = cp[i] + cp[j] + cp[p-i-j];
1715 Nodes.IntPoint(o++).Set3(1.0 - cp[j]/w, cp[i]/w, cp[j]/w);
1716 }
1717 for (int j = 1; j < p; j++)
1718 for (int i = 1; i + j < p; i++) // (3,4,2)
1719 {
1720 real_t w = cp[i] + cp[j] + cp[p-i-j];
1721 Nodes.IntPoint(o++).Set3(cp[j]/w, 1.0 - cp[i]/w, cp[i]/w);
1722 }
1723 for (int j = 1; j < p; j++)
1724 for (int i = 1; i + j < p; i++) // (0,4,3)
1725 {
1726 real_t w = cp[i] + cp[j] + cp[p-i-j];
1727 Nodes.IntPoint(o++).Set3(cp[0], cp[j]/w, cp[i]/w);
1728 }
1729
1730 // interior
1731 for (int k = 1; k < p - 1; k++)
1732 {
1733 for (int j = 1; j < p - k; j++)
1734 {
1735 real_t wjk = cp[j] + cp[k] + cp[p-j-k];
1736 for (int i = 1; i < p - k; i++)
1737 {
1738 real_t wik = cp[i] + cp[k] + cp[p-i-k];
1739 real_t w = wik * wjk * cp[p-k];
1740 Nodes.IntPoint(o++).Set3(cp[i] * (cp[j] + cp[p-j-k]) / w,
1741 cp[j] * (cp[i] + cp[p-i-k]) / w,
1742 cp[k] * cp[p-k] / w);
1743 }
1744 }
1745 }
1746
1747 MFEM_ASSERT(o == dof,
1748 "Number of nodes does not match the "
1749 "number of degrees of freedom");
1750 DenseMatrix T(dof);
1751
1752 for (int m = 0; m < dof; m++)
1753 {
1754 const IntegrationPoint &ip = Nodes.IntPoint(m);
1755
1756 real_t x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
1757 real_t y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
1758 real_t z = ip.z;
1759
1760 poly1d.CalcLegendre(p, x, shape_x.GetData());
1761 poly1d.CalcLegendre(p, y, shape_y.GetData());
1762
1763 o = 0;
1764 for (int i = 0; i <= p; i++)
1765 {
1766 for (int j = 0; j <= p; j++)
1767 {
1768 int maxij = std::max(i, j);
1769 FuentesPyramid::CalcScaledJacobi(p-maxij, 2.0 * (maxij + 1.0),
1770 z, 1.0, shape_z);
1771
1772 for (int k = 0; k <= p - maxij; k++)
1773 {
1774 T(o++, m) = shape_x(i) * shape_y(j) * shape_z(k) *
1775 pow(1.0 - ip.z, maxij);
1776 }
1777 }
1778 }
1779 }
1780
1781 Ti.Factor(T);
1782}
1783
1785 Vector &shape) const
1786{
1787 const int p = order;
1788
1789#ifdef MFEM_THREAD_SAFE
1790 Vector shape_x(order+1);
1791 Vector shape_y(order+1);
1792 Vector shape_z(order+1);
1793 Vector u(dof);
1794#endif
1795
1796 real_t x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
1797 real_t y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
1798 real_t z = ip.z;
1799
1800 poly1d.CalcLegendre(p, x, shape_x.GetData());
1801 poly1d.CalcLegendre(p, y, shape_y.GetData());
1802
1803 int o = 0;
1804 for (int i = 0; i <= p; i++)
1805 for (int j = 0; j <= p; j++)
1806 {
1807 int maxij = std::max(i, j);
1808 FuentesPyramid::CalcScaledJacobi(p-maxij, 2.0 * (maxij + 1.0), z, 1.0,
1809 shape_z);
1810 for (int k = 0; k <= p - maxij; k++)
1811 u[o++] = shape_x(i) * shape_y(j) * shape_z(k) *
1812 pow(1.0 - ip.z, maxij);
1813 }
1814
1815 Ti.Mult(u, shape);
1816}
1817
1819 DenseMatrix &dshape) const
1820{
1821 const int p = order;
1822
1823#ifdef MFEM_THREAD_SAFE
1824 DenseMatrix du(dof, dim);
1825 Vector shape_x(order+1);
1826 Vector shape_y(order+1);
1827 Vector shape_z(order+1);
1828 Vector dshape_x(order+1);
1829 Vector dshape_y(order+1);
1830 Vector dshape_z(order+1);
1831 Vector dshape_z_dt(order+1);
1832#endif
1833 real_t x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
1834 real_t y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
1835 real_t z = ip.z;
1836
1837 poly1d.CalcLegendre(p, x, shape_x.GetData(), dshape_x.GetData());
1838 poly1d.CalcLegendre(p, y, shape_y.GetData(), dshape_y.GetData());
1839
1840 int o = 0;
1841 for (int i = 0; i <= p; i++)
1842 for (int j = 0; j <= p; j++)
1843 {
1844 int maxij = std::max(i, j);
1845 FuentesPyramid::CalcScaledJacobi(p-maxij, 2.0 * (maxij + 1.0), z, 1.0,
1846 shape_z, dshape_z, dshape_z_dt);
1847
1848 for (int k = 0; k <= p - maxij; k++, o++)
1849 {
1850 du(o,0) = dshape_x(i) * shape_y(j) * shape_z(k) *
1851 pow(1.0 - ip.z, maxij - 1);
1852 du(o,1) = shape_x(i) * dshape_y(j) * shape_z(k) *
1853 pow(1.0 - ip.z, maxij - 1);
1854 du(o,2) = shape_x(i) * shape_y(j) * dshape_z(k) *
1855 pow(1.0 - ip.z, maxij) +
1856 (ip.x * dshape_x(i) * shape_y(j) +
1857 ip.y * shape_x(i) * dshape_y(j)) *
1858 shape_z(k) * pow(1.0 - ip.z, maxij - 2) -
1859 maxij * shape_x(i) * shape_y(j) * shape_z(k) *
1860 pow(1.0 - ip.z, maxij - 1);
1861 }
1862 }
1863
1864 Ti.Mult(du, dshape);
1865}
1866
1867
1868}
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:758
void Factor()
Factor the current DenseMatrix, *a.
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication with the inverse of dense matrix.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:108
void GetColumn(int c, Vector &col) const
Rank 3 tensor (array of matrices)
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
static int VerifyNodal(int b_type)
Ensure that the BasisType of b_type nodal (satisfies the interpolation property).
Definition fe_base.hpp:653
int dof
Number of degrees of freedom.
Definition fe_base.hpp:253
IntegrationRule Nodes
Definition fe_base.hpp:256
static int VerifyClosed(int b_type)
Ensure that the BasisType of b_type is closed (has Quadrature1D points on the boundary).
Definition fe_base.hpp:636
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition fe_base.hpp:400
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition fe_base.hpp:334
int order
Order/degree of the shape functions.
Definition fe_base.hpp:254
int dim
Dimension of reference space.
Definition fe_base.hpp:246
static DenseMatrix grad_mu01(real_t z)
void phi_Q(int p, Vector s, Vector t, DenseMatrix &u) const
static Vector lam45(real_t x, real_t y, real_t z)
static void phi_E(int p, real_t s0, real_t s1, real_t *u)
static Vector nu012(real_t z, Vector xy, unsigned int ab)
static real_t lam4(real_t x, real_t y, real_t z)
static DenseMatrix grad_lam25(real_t x, real_t y, real_t z)
static Vector lam35(real_t x, real_t y, real_t z)
static Vector grad_lam5(real_t x, real_t y, real_t z)
static real_t lam3(real_t x, real_t y, real_t z)
static Vector lam25(real_t x, real_t y, real_t z)
static Vector nu01(real_t z, Vector xy, unsigned int ab)
static Vector mu01(real_t z)
static DenseMatrix grad_nu012(real_t z, Vector xy, unsigned int ab)
static Vector grad_lam4(real_t x, real_t y, real_t z)
static Vector grad_mu1(real_t z)
static real_t mu0(real_t z)
static real_t lam2(real_t x, real_t y, real_t z)
static real_t lam1(real_t x, real_t y, real_t z)
Pyramid "Affine" Coordinates.
static Vector grad_mu0(real_t z)
static DenseMatrix grad_lam35(real_t x, real_t y, real_t z)
static real_t lam5(real_t x, real_t y, real_t z)
static DenseMatrix grad_lam15(real_t x, real_t y, real_t z)
Gradients of the above two component vectors.
static Vector grad_lam1(real_t x, real_t y, real_t z)
Gradients of the "Affine" Coordinates.
static Vector lam15(real_t x, real_t y, real_t z)
Two component vectors associated with edges touching the apex.
static real_t mu1(real_t z)
static Vector grad_lam2(real_t x, real_t y, real_t z)
static DenseMatrix grad_lam45(real_t x, real_t y, real_t z)
static DenseMatrix grad_nu01(real_t z, Vector xy, unsigned int ab)
static bool CheckZ(real_t z)
static void CalcScaledJacobi(int p, real_t alpha, real_t x, real_t t, real_t *u)
Shifted and Scaled Jacobi Polynomials.
void phi_T(int p, Vector nu, DenseMatrix &u) const
static Vector grad_lam3(real_t x, real_t y, real_t z)
Describes the function space on each element.
Definition fe_base.hpp:226
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_h1.cpp:1784
H1_BergotPyramidElement(const int p, const int btype=BasisType::GaussLobatto)
Definition fe_h1.cpp:1628
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_h1.cpp:1818
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_h1.cpp:1170
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_h1.cpp:1186
void CalcRawShape(const IntegrationPoint &ip, Vector &shape) const
Definition fe_h1.cpp:1203
H1_FuentesPyramidElement(const int p, const int btype=BasisType::GaussLobatto)
Definition fe_h1.cpp:1043
void CalcRawDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Definition fe_h1.cpp:1216
H1_HexahedronElement(const int p, const int btype=BasisType::GaussLobatto)
Construct the H1_HexahedronElement of order p and BasisType btype.
Definition fe_h1.cpp:265
void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_h1.cpp:293
void ProjectDelta(int vertex, Vector &dofs) const override
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition fe_h1.cpp:367
void CalcHessian(const IntegrationPoint &ip, DenseMatrix &Hessian) const override
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition fe_h1.cpp:338
void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const override
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_h1.cpp:314
void ProjectDelta(int vertex, Vector &dofs) const override
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition fe_h1.cpp:216
void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_h1.cpp:151
void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const override
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_h1.cpp:170
void CalcHessian(const IntegrationPoint &ip, DenseMatrix &Hessian) const override
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition fe_h1.cpp:192
H1_QuadrilateralElement(const int p, const int btype=BasisType::GaussLobatto)
Construct the H1_QuadrilateralElement of order p and BasisType btype.
Definition fe_h1.cpp:125
void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_h1.cpp:40
H1_SegmentElement(const int p, const int btype=BasisType::GaussLobatto)
Construct the H1_SegmentElement of order p and BasisType btype.
Definition fe_h1.cpp:21
void CalcHessian(const IntegrationPoint &ip, DenseMatrix &Hessian) const override
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition fe_h1.cpp:78
void ProjectDelta(int vertex, Vector &dofs) const override
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition fe_h1.cpp:97
void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const override
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_h1.cpp:59
void CalcHessian(const IntegrationPoint &ip, DenseMatrix &ddshape) const override
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition fe_h1.cpp:816
H1_TetrahedronElement(const int p, const int btype=BasisType::GaussLobatto)
Construct the H1_TetrahedronElement of order p and BasisType btype.
Definition fe_h1.cpp:617
void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_h1.cpp:758
void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const override
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_h1.cpp:783
void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const override
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_h1.cpp:555
void CalcHessian(const IntegrationPoint &ip, DenseMatrix &ddshape) const override
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition fe_h1.cpp:584
void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_h1.cpp:533
H1_TriangleElement(const int p, const int btype=BasisType::GaussLobatto)
Construct the H1_TriangleElement of order p and BasisType btype.
Definition fe_h1.cpp:451
void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const override
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_h1.cpp:1018
void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_h1.cpp:999
H1_WedgeElement(const int p, const int btype=BasisType::GaussLobatto)
Construct the H1_WedgeElement of order p and BasisType btype.
Definition fe_h1.cpp:863
Class for integration point with weight.
Definition intrules.hpp:35
void Set2(const real_t x1, const real_t x2)
Definition intrules.hpp:89
void Set3(const real_t x1, const real_t x2, const real_t x3)
Definition intrules.hpp:79
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition intrules.hpp:259
Class for standard nodal finite elements.
Definition fe_base.hpp:721
Array< int > lex_ordering
Definition fe_base.hpp:726
void Eval(const real_t x, Vector &u) const
Evaluate the basis functions at point x in [0,1].
Definition fe_base.cpp:1777
static real_t CalcDelta(const int p, const real_t x)
Evaluate a representation of a Delta function at point x.
Definition fe_base.hpp:1180
static void CalcLegendre(const int p, const real_t x, real_t *u)
Definition fe_base.cpp:2245
const real_t * ClosedPoints(const int p, const int btype=BasisType::GaussLobatto, bool on_device=false)
Get coordinates of a closed (GaussLegendre) set of points if degree p.
Definition fe_base.hpp:1121
static void CalcBasis(const int p, const real_t x, real_t *u)
Evaluate the values of a hierarchical 1D basis at point x hierarchical = k-th basis function is degre...
Definition fe_base.hpp:1140
Poly_1D::Basis & basis1d
Definition fe_base.hpp:1251
Vector data type.
Definition vector.hpp:82
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:558
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition vector.hpp:235
real_t mu
Definition ex25.cpp:140
Linear1DFiniteElement SegmentFE
Definition segment.cpp:52
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
float real_t
Definition config.hpp:43
Poly_1D poly1d
Definition fe.cpp:28
MFEM_EXPORT Linear2DFiniteElement TriangleFE
Definition fe.cpp:32
real_t p(const Vector &x, real_t t)