MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
fe_nd.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// Nedelec Finite Element classes
13
14#include "fe_nd.hpp"
15#include "face_map_utils.hpp"
16#include "../coefficient.hpp"
17
18namespace mfem
19{
20
21using namespace std;
22
23const real_t ND_HexahedronElement::tk[18] =
24{ 1.,0.,0., 0.,1.,0., 0.,0.,1., -1.,0.,0., 0.,-1.,0., 0.,0.,-1. };
25
27 const int cb_type, const int ob_type)
28 : VectorTensorFiniteElement(3, 3*p*(p + 1)*(p + 1), p, cb_type, ob_type,
29 H_CURL, DofMapType::L2_DOF_MAP),
30 dof2tk(dof), cp(poly1d.ClosedPoints(p, cb_type))
31{
32 if (obasis1d.IsIntegratedType()) { is_nodal = false; }
33
35
36 const real_t *op = poly1d.OpenPoints(p - 1, ob_type);
37 const int dof3 = dof/3;
38
39#ifndef MFEM_THREAD_SAFE
40 shape_cx.SetSize(p + 1);
41 shape_ox.SetSize(p);
42 shape_cy.SetSize(p + 1);
43 shape_oy.SetSize(p);
44 shape_cz.SetSize(p + 1);
45 shape_oz.SetSize(p);
46 dshape_cx.SetSize(p + 1);
47 dshape_cy.SetSize(p + 1);
48 dshape_cz.SetSize(p + 1);
49#endif
50
51 // edges
52 int o = 0;
53 for (int i = 0; i < p; i++) // (0,1)
54 {
55 dof_map[0*dof3 + i + (0 + 0*(p + 1))*p] = o++;
56 }
57 for (int i = 0; i < p; i++) // (1,2)
58 {
59 dof_map[1*dof3 + p + (i + 0*p)*(p + 1)] = o++;
60 }
61 for (int i = 0; i < p; i++) // (3,2)
62 {
63 dof_map[0*dof3 + i + (p + 0*(p + 1))*p] = o++;
64 }
65 for (int i = 0; i < p; i++) // (0,3)
66 {
67 dof_map[1*dof3 + 0 + (i + 0*p)*(p + 1)] = o++;
68 }
69 for (int i = 0; i < p; i++) // (4,5)
70 {
71 dof_map[0*dof3 + i + (0 + p*(p + 1))*p] = o++;
72 }
73 for (int i = 0; i < p; i++) // (5,6)
74 {
75 dof_map[1*dof3 + p + (i + p*p)*(p + 1)] = o++;
76 }
77 for (int i = 0; i < p; i++) // (7,6)
78 {
79 dof_map[0*dof3 + i + (p + p*(p + 1))*p] = o++;
80 }
81 for (int i = 0; i < p; i++) // (4,7)
82 {
83 dof_map[1*dof3 + 0 + (i + p*p)*(p + 1)] = o++;
84 }
85 for (int i = 0; i < p; i++) // (0,4)
86 {
87 dof_map[2*dof3 + 0 + (0 + i*(p + 1))*(p + 1)] = o++;
88 }
89 for (int i = 0; i < p; i++) // (1,5)
90 {
91 dof_map[2*dof3 + p + (0 + i*(p + 1))*(p + 1)] = o++;
92 }
93 for (int i = 0; i < p; i++) // (2,6)
94 {
95 dof_map[2*dof3 + p + (p + i*(p + 1))*(p + 1)] = o++;
96 }
97 for (int i = 0; i < p; i++) // (3,7)
98 {
99 dof_map[2*dof3 + 0 + (p + i*(p + 1))*(p + 1)] = o++;
100 }
101
102 // faces
103 // (3,2,1,0) -- bottom
104 for (int j = 1; j < p; j++) // x - components
105 for (int i = 0; i < p; i++)
106 {
107 dof_map[0*dof3 + i + ((p - j) + 0*(p + 1))*p] = o++;
108 }
109 for (int j = 0; j < p; j++) // y - components
110 for (int i = 1; i < p; i++)
111 {
112 dof_map[1*dof3 + i + ((p - 1 - j) + 0*p)*(p + 1)] = -1 - (o++);
113 }
114 // (0,1,5,4) -- front
115 for (int k = 1; k < p; k++) // x - components
116 for (int i = 0; i < p; i++)
117 {
118 dof_map[0*dof3 + i + (0 + k*(p + 1))*p] = o++;
119 }
120 for (int k = 0; k < p; k++) // z - components
121 for (int i = 1; i < p; i++ )
122 {
123 dof_map[2*dof3 + i + (0 + k*(p + 1))*(p + 1)] = o++;
124 }
125 // (1,2,6,5) -- right
126 for (int k = 1; k < p; k++) // y - components
127 for (int j = 0; j < p; j++)
128 {
129 dof_map[1*dof3 + p + (j + k*p)*(p + 1)] = o++;
130 }
131 for (int k = 0; k < p; k++) // z - components
132 for (int j = 1; j < p; j++)
133 {
134 dof_map[2*dof3 + p + (j + k*(p + 1))*(p + 1)] = o++;
135 }
136 // (2,3,7,6) -- back
137 for (int k = 1; k < p; k++) // x - components
138 for (int i = 0; i < p; i++)
139 {
140 dof_map[0*dof3 + (p - 1 - i) + (p + k*(p + 1))*p] = -1 - (o++);
141 }
142 for (int k = 0; k < p; k++) // z - components
143 for (int i = 1; i < p; i++)
144 {
145 dof_map[2*dof3 + (p - i) + (p + k*(p + 1))*(p + 1)] = o++;
146 }
147 // (3,0,4,7) -- left
148 for (int k = 1; k < p; k++) // y - components
149 for (int j = 0; j < p; j++)
150 {
151 dof_map[1*dof3 + 0 + ((p - 1 - j) + k*p)*(p + 1)] = -1 - (o++);
152 }
153 for (int k = 0; k < p; k++) // z - components
154 for (int j = 1; j < p; j++)
155 {
156 dof_map[2*dof3 + 0 + ((p - j) + k*(p + 1))*(p + 1)] = o++;
157 }
158 // (4,5,6,7) -- top
159 for (int j = 1; j < p; j++) // x - components
160 for (int i = 0; i < p; i++)
161 {
162 dof_map[0*dof3 + i + (j + p*(p + 1))*p] = o++;
163 }
164 for (int j = 0; j < p; j++) // y - components
165 for (int i = 1; i < p; i++)
166 {
167 dof_map[1*dof3 + i + (j + p*p)*(p + 1)] = o++;
168 }
169
170 // interior
171 // x-components
172 for (int k = 1; k < p; k++)
173 for (int j = 1; j < p; j++)
174 for (int i = 0; i < p; i++)
175 {
176 dof_map[0*dof3 + i + (j + k*(p + 1))*p] = o++;
177 }
178 // y-components
179 for (int k = 1; k < p; k++)
180 for (int j = 0; j < p; j++)
181 for (int i = 1; i < p; i++)
182 {
183 dof_map[1*dof3 + i + (j + k*p)*(p + 1)] = o++;
184 }
185 // z-components
186 for (int k = 0; k < p; k++)
187 for (int j = 1; j < p; j++)
188 for (int i = 1; i < p; i++)
189 {
190 dof_map[2*dof3 + i + (j + k*(p + 1))*(p + 1)] = o++;
191 }
192
193 // set dof2tk and Nodes
194 o = 0;
195 // x-components
196 for (int k = 0; k <= p; k++)
197 for (int j = 0; j <= p; j++)
198 for (int i = 0; i < p; i++)
199 {
200 int idx;
201 if ((idx = dof_map[o++]) < 0)
202 {
203 dof2tk[idx = -1 - idx] = 3;
204 }
205 else
206 {
207 dof2tk[idx] = 0;
208 }
209 Nodes.IntPoint(idx).Set3(op[i], cp[j], cp[k]);
210 }
211 // y-components
212 for (int k = 0; k <= p; k++)
213 for (int j = 0; j < p; j++)
214 for (int i = 0; i <= p; i++)
215 {
216 int idx;
217 if ((idx = dof_map[o++]) < 0)
218 {
219 dof2tk[idx = -1 - idx] = 4;
220 }
221 else
222 {
223 dof2tk[idx] = 1;
224 }
225 Nodes.IntPoint(idx).Set3(cp[i], op[j], cp[k]);
226 }
227 // z-components
228 for (int k = 0; k < p; k++)
229 for (int j = 0; j <= p; j++)
230 for (int i = 0; i <= p; i++)
231 {
232 int idx;
233 if ((idx = dof_map[o++]) < 0)
234 {
235 dof2tk[idx = -1 - idx] = 5;
236 }
237 else
238 {
239 dof2tk[idx] = 2;
240 }
241 Nodes.IntPoint(idx).Set3(cp[i], cp[j], op[k]);
242 }
243}
244
247 Vector &dofs) const
248{
249 MFEM_ASSERT(obasis1d.IsIntegratedType(), "Not integrated type");
251 Vector xk(vk, vc.GetVDim());
252
254 const int nqpt = ir.GetNPoints();
255
256 IntegrationPoint ip3d;
257
258 int o = 0;
259 for (int c = 0; c < 3; ++c) // loop over x, y, z components
260 {
261 const int im = c == 0 ? order - 1 : order;
262 const int jm = c == 1 ? order - 1 : order;
263 const int km = c == 2 ? order - 1 : order;
264
265 for (int k = 0; k <= km; k++)
266 for (int j = 0; j <= jm; j++)
267 for (int i = 0; i <= im; i++)
268 {
269 int idx;
270 if ((idx = dof_map[o++]) < 0)
271 {
272 idx = -1 - idx;
273 }
274
275 const int id1 = c == 0 ? i : (c == 1 ? j : k);
276 const real_t h = cp[id1+1] - cp[id1];
277
278 real_t val = 0.0;
279
280 for (int q = 0; q < nqpt; q++)
281 {
282 const IntegrationPoint &ip1d = ir.IntPoint(q);
283
284 if (c == 0)
285 {
286 ip3d.Set3(cp[i] + (h*ip1d.x), cp[j], cp[k]);
287 }
288 else if (c == 1)
289 {
290 ip3d.Set3(cp[i], cp[j] + (h*ip1d.x), cp[k]);
291 }
292 else
293 {
294 ip3d.Set3(cp[i], cp[j], cp[k] + (h*ip1d.x));
295 }
296
297 Trans.SetIntPoint(&ip3d);
298 vc.Eval(xk, Trans, ip3d);
299
300 // xk^t J tk
301 const real_t ipval = Trans.Jacobian().InnerProduct(tk + dof2tk[idx]*dim, vk);
302 val += ip1d.weight * ipval;
303 }
304
305 dofs(idx) = val*h;
306 }
307 }
308}
309
311 DenseMatrix &shape) const
312{
313 const int p = order;
314
315#ifdef MFEM_THREAD_SAFE
316 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
317 Vector shape_cz(p + 1), shape_oz(p);
318#endif
319
321 {
322#ifdef MFEM_THREAD_SAFE
323 Vector dshape_cx(p + 1), dshape_cy(p + 1), dshape_cz(p + 1);
324#endif
325 basis1d.Eval(ip.x, shape_cx, dshape_cx);
326 basis1d.Eval(ip.y, shape_cy, dshape_cy);
327 basis1d.Eval(ip.z, shape_cz, dshape_cz);
329 obasis1d.EvalIntegrated(dshape_cx, shape_ox);
330 obasis1d.EvalIntegrated(dshape_cy, shape_oy);
331 obasis1d.EvalIntegrated(dshape_cz, shape_oz);
332 }
333 else
334 {
335 basis1d.Eval(ip.x, shape_cx);
336 basis1d.Eval(ip.y, shape_cy);
337 basis1d.Eval(ip.z, shape_cz);
338 obasis1d.Eval(ip.x, shape_ox);
339 obasis1d.Eval(ip.y, shape_oy);
340 obasis1d.Eval(ip.z, shape_oz);
341 }
342
343 int o = 0;
344 // x-components
345 for (int k = 0; k <= p; k++)
346 for (int j = 0; j <= p; j++)
347 for (int i = 0; i < p; i++)
348 {
349 int idx, s;
350 if ((idx = dof_map[o++]) < 0)
351 {
352 idx = -1 - idx, s = -1;
353 }
354 else
355 {
356 s = +1;
357 }
358 shape(idx,0) = s*shape_ox(i)*shape_cy(j)*shape_cz(k);
359 shape(idx,1) = 0.;
360 shape(idx,2) = 0.;
361 }
362 // y-components
363 for (int k = 0; k <= p; k++)
364 for (int j = 0; j < p; j++)
365 for (int i = 0; i <= p; i++)
366 {
367 int idx, s;
368 if ((idx = dof_map[o++]) < 0)
369 {
370 idx = -1 - idx, s = -1;
371 }
372 else
373 {
374 s = +1;
375 }
376 shape(idx,0) = 0.;
377 shape(idx,1) = s*shape_cx(i)*shape_oy(j)*shape_cz(k);
378 shape(idx,2) = 0.;
379 }
380 // z-components
381 for (int k = 0; k < p; k++)
382 for (int j = 0; j <= p; j++)
383 for (int i = 0; i <= p; i++)
384 {
385 int idx, s;
386 if ((idx = dof_map[o++]) < 0)
387 {
388 idx = -1 - idx, s = -1;
389 }
390 else
391 {
392 s = +1;
393 }
394 shape(idx,0) = 0.;
395 shape(idx,1) = 0.;
396 shape(idx,2) = s*shape_cx(i)*shape_cy(j)*shape_oz(k);
397 }
398}
399
401 DenseMatrix &curl_shape) const
402{
403 const int p = order;
404
405#ifdef MFEM_THREAD_SAFE
406 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
407 Vector shape_cz(p + 1), shape_oz(p);
408 Vector dshape_cx(p + 1), dshape_cy(p + 1), dshape_cz(p + 1);
409#endif
410
411 basis1d.Eval(ip.x, shape_cx, dshape_cx);
412 basis1d.Eval(ip.y, shape_cy, dshape_cy);
413 basis1d.Eval(ip.z, shape_cz, dshape_cz);
415 {
417 obasis1d.EvalIntegrated(dshape_cx, shape_ox);
418 obasis1d.EvalIntegrated(dshape_cy, shape_oy);
419 obasis1d.EvalIntegrated(dshape_cz, shape_oz);
420 }
421 else
422 {
423 obasis1d.Eval(ip.x, shape_ox);
424 obasis1d.Eval(ip.y, shape_oy);
425 obasis1d.Eval(ip.z, shape_oz);
426 }
427
428 int o = 0;
429 // x-components
430 for (int k = 0; k <= p; k++)
431 for (int j = 0; j <= p; j++)
432 for (int i = 0; i < p; i++)
433 {
434 int idx, s;
435 if ((idx = dof_map[o++]) < 0)
436 {
437 idx = -1 - idx, s = -1;
438 }
439 else
440 {
441 s = +1;
442 }
443 curl_shape(idx,0) = 0.;
444 curl_shape(idx,1) = s*shape_ox(i)* shape_cy(j)*dshape_cz(k);
445 curl_shape(idx,2) = -s*shape_ox(i)*dshape_cy(j)* shape_cz(k);
446 }
447 // y-components
448 for (int k = 0; k <= p; k++)
449 for (int j = 0; j < p; j++)
450 for (int i = 0; i <= p; i++)
451 {
452 int idx, s;
453 if ((idx = dof_map[o++]) < 0)
454 {
455 idx = -1 - idx, s = -1;
456 }
457 else
458 {
459 s = +1;
460 }
461 curl_shape(idx,0) = -s* shape_cx(i)*shape_oy(j)*dshape_cz(k);
462 curl_shape(idx,1) = 0.;
463 curl_shape(idx,2) = s*dshape_cx(i)*shape_oy(j)* shape_cz(k);
464 }
465 // z-components
466 for (int k = 0; k < p; k++)
467 for (int j = 0; j <= p; j++)
468 for (int i = 0; i <= p; i++)
469 {
470 int idx, s;
471 if ((idx = dof_map[o++]) < 0)
472 {
473 idx = -1 - idx, s = -1;
474 }
475 else
476 {
477 s = +1;
478 }
479 curl_shape(idx,0) = s* shape_cx(i)*dshape_cy(j)*shape_oz(k);
480 curl_shape(idx,1) = -s*dshape_cx(i)* shape_cy(j)*shape_oz(k);
481 curl_shape(idx,2) = 0.;
482 }
483}
484
485void ND_HexahedronElement::GetFaceMap(const int face_id,
486 Array<int> &face_map) const
487{
488 const int p = order;
489 const int pp1 = p + 1;
490 const int n_face_dofs_per_component = p*pp1;
491 const int n_dof_per_dim = p*pp1*pp1;
492
493 std::vector<int> n_dofs = {p, pp1, pp1, p};
494 std::vector<int> offsets, strides;
495
496 const auto f = internal::GetFaceNormal3D(face_id);
497 const int face_normal = f.first, level = f.second;
498 if (face_normal == 0) // x-normal
499 {
500 offsets =
501 {
502 n_dof_per_dim + (level ? pp1 - 1 : 0),
503 2*n_dof_per_dim + (level ? pp1 - 1 : 0)
504 };
505 strides = {pp1, p*pp1, pp1, pp1*pp1};
506 }
507 else if (face_normal == 1) // y-normal
508 {
509 offsets =
510 {
511 level ? p*(pp1 - 1) : 0,
512 2*n_dof_per_dim + (level ? pp1*(pp1 - 1) : 0)
513 };
514 strides = {1, p*pp1, 1, pp1*pp1};
515 }
516 else if (face_normal == 2) // z-normal
517 {
518 offsets =
519 {
520 level ? p*pp1*(pp1 - 1) : 0,
521 n_dof_per_dim + (level ? p*pp1*(pp1 - 1) : 0)
522 };
523 strides = {1, p, 1, pp1};
524 }
525
526 internal::FillFaceMap(n_face_dofs_per_component, offsets, strides, n_dofs,
527 face_map);
528}
529
530const real_t ND_QuadrilateralElement::tk[8] =
531{ 1.,0., 0.,1., -1.,0., 0.,-1. };
532
534 const int cb_type,
535 const int ob_type)
536 : VectorTensorFiniteElement(2, 2*p*(p + 1), p, cb_type, ob_type,
537 H_CURL, DofMapType::L2_DOF_MAP),
538 dof2tk(dof),
539 cp(poly1d.ClosedPoints(p, cb_type))
540{
541 if (obasis1d.IsIntegratedType()) { is_nodal = false; }
542
544
545 const real_t *op = poly1d.OpenPoints(p - 1, ob_type);
546 const int dof2 = dof/2;
547
548#ifndef MFEM_THREAD_SAFE
549 shape_cx.SetSize(p + 1);
550 shape_ox.SetSize(p);
551 shape_cy.SetSize(p + 1);
552 shape_oy.SetSize(p);
553 dshape_cx.SetSize(p + 1);
554 dshape_cy.SetSize(p + 1);
555#endif
556
557 // edges
558 int o = 0;
559 for (int i = 0; i < p; i++) // (0,1)
560 {
561 dof_map[0*dof2 + i + 0*p] = o++;
562 }
563 for (int j = 0; j < p; j++) // (1,2)
564 {
565 dof_map[1*dof2 + p + j*(p + 1)] = o++;
566 }
567 for (int i = 0; i < p; i++) // (2,3)
568 {
569 dof_map[0*dof2 + (p - 1 - i) + p*p] = -1 - (o++);
570 }
571 for (int j = 0; j < p; j++) // (3,0)
572 {
573 dof_map[1*dof2 + 0 + (p - 1 - j)*(p + 1)] = -1 - (o++);
574 }
575
576 // interior
577 // x-components
578 for (int j = 1; j < p; j++)
579 for (int i = 0; i < p; i++)
580 {
581 dof_map[0*dof2 + i + j*p] = o++;
582 }
583 // y-components
584 for (int j = 0; j < p; j++)
585 for (int i = 1; i < p; i++)
586 {
587 dof_map[1*dof2 + i + j*(p + 1)] = o++;
588 }
589
590 // set dof2tk and Nodes
591 o = 0;
592 // x-components
593 for (int j = 0; j <= p; j++)
594 for (int i = 0; i < p; i++)
595 {
596 int idx;
597 if ((idx = dof_map[o++]) < 0)
598 {
599 dof2tk[idx = -1 - idx] = 2;
600 }
601 else
602 {
603 dof2tk[idx] = 0;
604 }
605 Nodes.IntPoint(idx).Set2(op[i], cp[j]);
606 }
607 // y-components
608 for (int j = 0; j < p; j++)
609 for (int i = 0; i <= p; i++)
610 {
611 int idx;
612 if ((idx = dof_map[o++]) < 0)
613 {
614 dof2tk[idx = -1 - idx] = 3;
615 }
616 else
617 {
618 dof2tk[idx] = 1;
619 }
620 Nodes.IntPoint(idx).Set2(cp[i], op[j]);
621 }
622}
623
626 Vector &dofs) const
627{
628 MFEM_ASSERT(obasis1d.IsIntegratedType(), "Not integrated type");
630 Vector xk(vk, vc.GetVDim());
631
633 const int nqpt = ir.GetNPoints();
634
635 IntegrationPoint ip2d;
636
637 int o = 0;
638 // x-components
639 for (int j = 0; j <= order; j++)
640 for (int i = 0; i < order; i++)
641 {
642 int idx;
643 if ((idx = dof_map[o++]) < 0)
644 {
645 idx = -1 - idx;
646 }
647
648 const real_t h = cp[i+1] - cp[i];
649
650 real_t val = 0.0;
651
652 for (int k = 0; k < nqpt; k++)
653 {
654 const IntegrationPoint &ip1d = ir.IntPoint(k);
655
656 ip2d.Set2(cp[i] + (h*ip1d.x), cp[j]);
657
658 Trans.SetIntPoint(&ip2d);
659 vc.Eval(xk, Trans, ip2d);
660
661 // xk^t J tk
662 const real_t ipval = Trans.Jacobian().InnerProduct(tk + dof2tk[idx]*dim, vk);
663 val += ip1d.weight * ipval;
664 }
665
666 dofs(idx) = val*h;
667 }
668 // y-components
669 for (int j = 0; j < order; j++)
670 for (int i = 0; i <= order; i++)
671 {
672 int idx;
673 if ((idx = dof_map[o++]) < 0)
674 {
675 idx = -1 - idx;
676 }
677
678 const real_t h = cp[j+1] - cp[j];
679
680 real_t val = 0.0;
681
682 for (int k = 0; k < nqpt; k++)
683 {
684 const IntegrationPoint &ip1d = ir.IntPoint(k);
685
686 ip2d.Set2(cp[i], cp[j] + (h*ip1d.x));
687
688 Trans.SetIntPoint(&ip2d);
689 vc.Eval(xk, Trans, ip2d);
690
691 // xk^t J tk
692 const real_t ipval = Trans.Jacobian().InnerProduct(tk + dof2tk[idx]*dim, vk);
693 val += ip1d.weight * ipval;
694 }
695
696 dofs(idx) = val*h;
697 }
698}
699
701 DenseMatrix &shape) const
702{
703 const int p = order;
704
705#ifdef MFEM_THREAD_SAFE
706 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
707#endif
708
710 {
711#ifdef MFEM_THREAD_SAFE
712 Vector dshape_cx(p + 1), dshape_cy(p + 1);
713#endif
714 basis1d.Eval(ip.x, shape_cx, dshape_cx);
715 basis1d.Eval(ip.y, shape_cy, dshape_cy);
717 obasis1d.EvalIntegrated(dshape_cx, shape_ox);
718 obasis1d.EvalIntegrated(dshape_cy, shape_oy);
719 }
720 else
721 {
722 basis1d.Eval(ip.x, shape_cx);
723 basis1d.Eval(ip.y, shape_cy);
724 obasis1d.Eval(ip.x, shape_ox);
725 obasis1d.Eval(ip.y, shape_oy);
726 }
727
728 int o = 0;
729 // x-components
730 for (int j = 0; j <= p; j++)
731 for (int i = 0; i < p; i++)
732 {
733 int idx, s;
734 if ((idx = dof_map[o++]) < 0)
735 {
736 idx = -1 - idx, s = -1;
737 }
738 else
739 {
740 s = +1;
741 }
742 shape(idx,0) = s*shape_ox(i)*shape_cy(j);
743 shape(idx,1) = 0.;
744 }
745 // y-components
746 for (int j = 0; j < p; j++)
747 for (int i = 0; i <= p; i++)
748 {
749 int idx, s;
750 if ((idx = dof_map[o++]) < 0)
751 {
752 idx = -1 - idx, s = -1;
753 }
754 else
755 {
756 s = +1;
757 }
758 shape(idx,0) = 0.;
759 shape(idx,1) = s*shape_cx(i)*shape_oy(j);
760 }
761}
762
764 DenseMatrix &curl_shape) const
765{
766 const int p = order;
767
768#ifdef MFEM_THREAD_SAFE
769 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
770 Vector dshape_cx(p + 1), dshape_cy(p + 1);
771#endif
772
773 basis1d.Eval(ip.x, shape_cx, dshape_cx);
774 basis1d.Eval(ip.y, shape_cy, dshape_cy);
776 {
778 obasis1d.EvalIntegrated(dshape_cx, shape_ox);
779 obasis1d.EvalIntegrated(dshape_cy, shape_oy);
780 }
781 else
782 {
783 obasis1d.Eval(ip.x, shape_ox);
784 obasis1d.Eval(ip.y, shape_oy);
785 }
786
787 int o = 0;
788 // x-components
789 for (int j = 0; j <= p; j++)
790 for (int i = 0; i < p; i++)
791 {
792 int idx, s;
793 if ((idx = dof_map[o++]) < 0)
794 {
795 idx = -1 - idx, s = -1;
796 }
797 else
798 {
799 s = +1;
800 }
801 curl_shape(idx,0) = -s*shape_ox(i)*dshape_cy(j);
802 }
803 // y-components
804 for (int j = 0; j < p; j++)
805 for (int i = 0; i <= p; i++)
806 {
807 int idx, s;
808 if ((idx = dof_map[o++]) < 0)
809 {
810 idx = -1 - idx, s = -1;
811 }
812 else
813 {
814 s = +1;
815 }
816 curl_shape(idx,0) = s*dshape_cx(i)*shape_oy(j);
817 }
818}
819
821 Array<int> &face_map) const
822{
823 const int p = order;
824 const int pp1 = order + 1;
825 const int n_face_dofs_per_component = p;
826 std::vector<int> strides = {(face_id == 0 || face_id == 2) ? 1 : pp1};
827 std::vector<int> n_dofs = {p};
828 std::vector<int> offsets;
829 switch (face_id)
830 {
831 case 0: offsets = {0}; break; // y = 0
832 case 1: offsets = {p*pp1 + pp1 - 1}; break; // x = 1
833 case 2: offsets = {p*(pp1 - 1)}; break; // y = 1
834 case 3: offsets = {p*pp1}; break; // x = 0
835 }
836 internal::FillFaceMap(n_face_dofs_per_component, offsets, strides, n_dofs,
837 face_map);
838}
839
840
841const real_t ND_TetrahedronElement::tk[18] =
842{ 1.,0.,0., 0.,1.,0., 0.,0.,1., -1.,1.,0., -1.,0.,1., 0.,-1.,1. };
843
844const real_t ND_TetrahedronElement::c = 1./4.;
845
847 : VectorFiniteElement(3, Geometry::TETRAHEDRON, p*(p + 2)*(p + 3)/2, p,
848 H_CURL, FunctionSpace::Pk), dof2tk(dof), doftrans(p)
849{
850 const real_t *eop = poly1d.OpenPoints(p - 1);
851 const real_t *fop = (p > 1) ? poly1d.OpenPoints(p - 2) : NULL;
852 const real_t *iop = (p > 2) ? poly1d.OpenPoints(p - 3) : NULL;
853
854 const int pm1 = p - 1, pm2 = p - 2, pm3 = p - 3;
855
856#ifndef MFEM_THREAD_SAFE
857 shape_x.SetSize(p);
858 shape_y.SetSize(p);
859 shape_z.SetSize(p);
860 shape_l.SetSize(p);
861 dshape_x.SetSize(p);
862 dshape_y.SetSize(p);
863 dshape_z.SetSize(p);
864 dshape_l.SetSize(p);
865 u.SetSize(dof, dim);
866#else
867 Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
868#endif
869
870 int o = 0;
871 // edges
872 for (int i = 0; i < p; i++) // (0,1)
873 {
874 Nodes.IntPoint(o).Set3(eop[i], 0., 0.);
875 dof2tk[o++] = 0;
876 }
877 for (int i = 0; i < p; i++) // (0,2)
878 {
879 Nodes.IntPoint(o).Set3(0., eop[i], 0.);
880 dof2tk[o++] = 1;
881 }
882 for (int i = 0; i < p; i++) // (0,3)
883 {
884 Nodes.IntPoint(o).Set3(0., 0., eop[i]);
885 dof2tk[o++] = 2;
886 }
887 for (int i = 0; i < p; i++) // (1,2)
888 {
889 Nodes.IntPoint(o).Set3(eop[pm1-i], eop[i], 0.);
890 dof2tk[o++] = 3;
891 }
892 for (int i = 0; i < p; i++) // (1,3)
893 {
894 Nodes.IntPoint(o).Set3(eop[pm1-i], 0., eop[i]);
895 dof2tk[o++] = 4;
896 }
897 for (int i = 0; i < p; i++) // (2,3)
898 {
899 Nodes.IntPoint(o).Set3(0., eop[pm1-i], eop[i]);
900 dof2tk[o++] = 5;
901 }
902
903 // faces
904 for (int j = 0; j <= pm2; j++) // (1,2,3)
905 for (int i = 0; i + j <= pm2; i++)
906 {
907 real_t w = fop[i] + fop[j] + fop[pm2-i-j];
908 Nodes.IntPoint(o).Set3(fop[pm2-i-j]/w, fop[i]/w, fop[j]/w);
909 dof2tk[o++] = 3;
910 Nodes.IntPoint(o).Set3(fop[pm2-i-j]/w, fop[i]/w, fop[j]/w);
911 dof2tk[o++] = 4;
912 }
913 for (int j = 0; j <= pm2; j++) // (0,3,2)
914 for (int i = 0; i + j <= pm2; i++)
915 {
916 real_t w = fop[i] + fop[j] + fop[pm2-i-j];
917 Nodes.IntPoint(o).Set3(0., fop[j]/w, fop[i]/w);
918 dof2tk[o++] = 2;
919 Nodes.IntPoint(o).Set3(0., fop[j]/w, fop[i]/w);
920 dof2tk[o++] = 1;
921 }
922 for (int j = 0; j <= pm2; j++) // (0,1,3)
923 for (int i = 0; i + j <= pm2; i++)
924 {
925 real_t w = fop[i] + fop[j] + fop[pm2-i-j];
926 Nodes.IntPoint(o).Set3(fop[i]/w, 0., fop[j]/w);
927 dof2tk[o++] = 0;
928 Nodes.IntPoint(o).Set3(fop[i]/w, 0., fop[j]/w);
929 dof2tk[o++] = 2;
930 }
931 for (int j = 0; j <= pm2; j++) // (0,2,1)
932 for (int i = 0; i + j <= pm2; i++)
933 {
934 real_t w = fop[i] + fop[j] + fop[pm2-i-j];
935 Nodes.IntPoint(o).Set3(fop[j]/w, fop[i]/w, 0.);
936 dof2tk[o++] = 1;
937 Nodes.IntPoint(o).Set3(fop[j]/w, fop[i]/w, 0.);
938 dof2tk[o++] = 0;
939 }
940
941 // interior
942 for (int k = 0; k <= pm3; k++)
943 for (int j = 0; j + k <= pm3; j++)
944 for (int i = 0; i + j + k <= pm3; i++)
945 {
946 real_t w = iop[i] + iop[j] + iop[k] + iop[pm3-i-j-k];
947 Nodes.IntPoint(o).Set3(iop[i]/w, iop[j]/w, iop[k]/w);
948 dof2tk[o++] = 0;
949 Nodes.IntPoint(o).Set3(iop[i]/w, iop[j]/w, iop[k]/w);
950 dof2tk[o++] = 1;
951 Nodes.IntPoint(o).Set3(iop[i]/w, iop[j]/w, iop[k]/w);
952 dof2tk[o++] = 2;
953 }
954
955 DenseMatrix T(dof);
956 for (int m = 0; m < dof; m++)
957 {
958 const IntegrationPoint &ip = Nodes.IntPoint(m);
959 const real_t *tm = tk + 3*dof2tk[m];
960 o = 0;
961
962 poly1d.CalcBasis(pm1, ip.x, shape_x);
963 poly1d.CalcBasis(pm1, ip.y, shape_y);
964 poly1d.CalcBasis(pm1, ip.z, shape_z);
965 poly1d.CalcBasis(pm1, 1. - ip.x - ip.y - ip.z, shape_l);
966
967 for (int k = 0; k <= pm1; k++)
968 for (int j = 0; j + k <= pm1; j++)
969 for (int i = 0; i + j + k <= pm1; i++)
970 {
971 real_t s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(pm1-i-j-k);
972 T(o++, m) = s * tm[0];
973 T(o++, m) = s * tm[1];
974 T(o++, m) = s * tm[2];
975 }
976 for (int k = 0; k <= pm1; k++)
977 for (int j = 0; j + k <= pm1; j++)
978 {
979 real_t s = shape_x(pm1-j-k)*shape_y(j)*shape_z(k);
980 T(o++, m) = s*((ip.y - c)*tm[0] - (ip.x - c)*tm[1]);
981 T(o++, m) = s*((ip.z - c)*tm[0] - (ip.x - c)*tm[2]);
982 }
983 for (int k = 0; k <= pm1; k++)
984 {
985 T(o++, m) =
986 shape_y(pm1-k)*shape_z(k)*((ip.z - c)*tm[1] - (ip.y - c)*tm[2]);
987 }
988 }
989
990 Ti.Factor(T);
991 // mfem::out << "ND_TetrahedronElement(" << p << ") : "; Ti.TestInversion();
992}
993
995 DenseMatrix &shape) const
996{
997 const int pm1 = order - 1;
998
999#ifdef MFEM_THREAD_SAFE
1000 const int p = order;
1001 Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
1002 DenseMatrix u(dof, dim);
1003#endif
1004
1005 poly1d.CalcBasis(pm1, ip.x, shape_x);
1006 poly1d.CalcBasis(pm1, ip.y, shape_y);
1007 poly1d.CalcBasis(pm1, ip.z, shape_z);
1008 poly1d.CalcBasis(pm1, 1. - ip.x - ip.y - ip.z, shape_l);
1009
1010 int n = 0;
1011 for (int k = 0; k <= pm1; k++)
1012 for (int j = 0; j + k <= pm1; j++)
1013 for (int i = 0; i + j + k <= pm1; i++)
1014 {
1015 real_t s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(pm1-i-j-k);
1016 u(n,0) = s; u(n,1) = 0.; u(n,2) = 0.; n++;
1017 u(n,0) = 0.; u(n,1) = s; u(n,2) = 0.; n++;
1018 u(n,0) = 0.; u(n,1) = 0.; u(n,2) = s; n++;
1019 }
1020 for (int k = 0; k <= pm1; k++)
1021 for (int j = 0; j + k <= pm1; j++)
1022 {
1023 real_t s = shape_x(pm1-j-k)*shape_y(j)*shape_z(k);
1024 u(n,0) = s*(ip.y - c); u(n,1) = -s*(ip.x - c); u(n,2) = 0.; n++;
1025 u(n,0) = s*(ip.z - c); u(n,1) = 0.; u(n,2) = -s*(ip.x - c); n++;
1026 }
1027 for (int k = 0; k <= pm1; k++)
1028 {
1029 real_t s = shape_y(pm1-k)*shape_z(k);
1030 u(n,0) = 0.; u(n,1) = s*(ip.z - c); u(n,2) = -s*(ip.y - c); n++;
1031 }
1032
1033 Ti.Mult(u, shape);
1034}
1035
1037 DenseMatrix &curl_shape) const
1038{
1039 const int pm1 = order - 1;
1040
1041#ifdef MFEM_THREAD_SAFE
1042 const int p = order;
1043 Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
1044 Vector dshape_x(p), dshape_y(p), dshape_z(p), dshape_l(p);
1045 DenseMatrix u(dof, dim);
1046#endif
1047
1048 poly1d.CalcBasis(pm1, ip.x, shape_x, dshape_x);
1049 poly1d.CalcBasis(pm1, ip.y, shape_y, dshape_y);
1050 poly1d.CalcBasis(pm1, ip.z, shape_z, dshape_z);
1051 poly1d.CalcBasis(pm1, 1. - ip.x - ip.y - ip.z, shape_l, dshape_l);
1052
1053 int n = 0;
1054 for (int k = 0; k <= pm1; k++)
1055 for (int j = 0; j + k <= pm1; j++)
1056 for (int i = 0; i + j + k <= pm1; i++)
1057 {
1058 int l = pm1-i-j-k;
1059 const real_t dx = (dshape_x(i)*shape_l(l) -
1060 shape_x(i)*dshape_l(l))*shape_y(j)*shape_z(k);
1061 const real_t dy = (dshape_y(j)*shape_l(l) -
1062 shape_y(j)*dshape_l(l))*shape_x(i)*shape_z(k);
1063 const real_t dz = (dshape_z(k)*shape_l(l) -
1064 shape_z(k)*dshape_l(l))*shape_x(i)*shape_y(j);
1065
1066 u(n,0) = 0.; u(n,1) = dz; u(n,2) = -dy; n++;
1067 u(n,0) = -dz; u(n,1) = 0.; u(n,2) = dx; n++;
1068 u(n,0) = dy; u(n,1) = -dx; u(n,2) = 0.; n++;
1069 }
1070 for (int k = 0; k <= pm1; k++)
1071 for (int j = 0; j + k <= pm1; j++)
1072 {
1073 int i = pm1 - j - k;
1074 // s = shape_x(i)*shape_y(j)*shape_z(k);
1075 // curl of s*(ip.y - c, -(ip.x - c), 0):
1076 u(n,0) = shape_x(i)*(ip.x - c)*shape_y(j)*dshape_z(k);
1077 u(n,1) = shape_x(i)*shape_y(j)*(ip.y - c)*dshape_z(k);
1078 u(n,2) =
1079 -((dshape_x(i)*(ip.x - c) + shape_x(i))*shape_y(j)*shape_z(k) +
1080 (dshape_y(j)*(ip.y - c) + shape_y(j))*shape_x(i)*shape_z(k));
1081 n++;
1082 // curl of s*(ip.z - c, 0, -(ip.x - c)):
1083 u(n,0) = -shape_x(i)*(ip.x - c)*dshape_y(j)*shape_z(k);
1084 u(n,1) = (shape_x(i)*shape_y(j)*(dshape_z(k)*(ip.z - c) + shape_z(k)) +
1085 (dshape_x(i)*(ip.x - c) + shape_x(i))*shape_y(j)*shape_z(k));
1086 u(n,2) = -shape_x(i)*dshape_y(j)*shape_z(k)*(ip.z - c);
1087 n++;
1088 }
1089 for (int k = 0; k <= pm1; k++)
1090 {
1091 int j = pm1 - k;
1092 // curl of shape_y(j)*shape_z(k)*(0, ip.z - c, -(ip.y - c)):
1093 u(n,0) = -((dshape_y(j)*(ip.y - c) + shape_y(j))*shape_z(k) +
1094 shape_y(j)*(dshape_z(k)*(ip.z - c) + shape_z(k)));
1095 u(n,1) = 0.;
1096 u(n,2) = 0.; n++;
1097 }
1098
1099 Ti.Mult(u, curl_shape);
1100}
1101
1102
1103const real_t ND_TriangleElement::tk[8] =
1104{ 1.,0., -1.,1., 0.,-1., 0.,1. };
1105
1106const real_t ND_TriangleElement::c = 1./3.;
1107
1109 : VectorFiniteElement(2, Geometry::TRIANGLE, p*(p + 2), p,
1110 H_CURL, FunctionSpace::Pk),
1111 dof2tk(dof), doftrans(p)
1112{
1113 const real_t *eop = poly1d.OpenPoints(p - 1);
1114 const real_t *iop = (p > 1) ? poly1d.OpenPoints(p - 2) : NULL;
1115
1116 const int pm1 = p - 1, pm2 = p - 2;
1117
1118#ifndef MFEM_THREAD_SAFE
1119 shape_x.SetSize(p);
1120 shape_y.SetSize(p);
1121 shape_l.SetSize(p);
1122 dshape_x.SetSize(p);
1123 dshape_y.SetSize(p);
1124 dshape_l.SetSize(p);
1125 u.SetSize(dof, dim);
1126 curlu.SetSize(dof);
1127#else
1128 Vector shape_x(p), shape_y(p), shape_l(p);
1129#endif
1130
1131 int n = 0;
1132 // edges
1133 for (int i = 0; i < p; i++) // (0,1)
1134 {
1135 Nodes.IntPoint(n).Set2(eop[i], 0.);
1136 dof2tk[n++] = 0;
1137 }
1138 for (int i = 0; i < p; i++) // (1,2)
1139 {
1140 Nodes.IntPoint(n).Set2(eop[pm1-i], eop[i]);
1141 dof2tk[n++] = 1;
1142 }
1143 for (int i = 0; i < p; i++) // (2,0)
1144 {
1145 Nodes.IntPoint(n).Set2(0., eop[pm1-i]);
1146 dof2tk[n++] = 2;
1147 }
1148
1149 // interior
1150 for (int j = 0; j <= pm2; j++)
1151 for (int i = 0; i + j <= pm2; i++)
1152 {
1153 real_t w = iop[i] + iop[j] + iop[pm2-i-j];
1154 Nodes.IntPoint(n).Set2(iop[i]/w, iop[j]/w);
1155 dof2tk[n++] = 0;
1156 Nodes.IntPoint(n).Set2(iop[i]/w, iop[j]/w);
1157 dof2tk[n++] = 3;
1158 }
1159
1160 DenseMatrix T(dof);
1161 for (int m = 0; m < dof; m++)
1162 {
1163 const IntegrationPoint &ip = Nodes.IntPoint(m);
1164 const real_t *tm = tk + 2*dof2tk[m];
1165 n = 0;
1166
1167 poly1d.CalcBasis(pm1, ip.x, shape_x);
1168 poly1d.CalcBasis(pm1, ip.y, shape_y);
1169 poly1d.CalcBasis(pm1, 1. - ip.x - ip.y, shape_l);
1170
1171 for (int j = 0; j <= pm1; j++)
1172 for (int i = 0; i + j <= pm1; i++)
1173 {
1174 real_t s = shape_x(i)*shape_y(j)*shape_l(pm1-i-j);
1175 T(n++, m) = s * tm[0];
1176 T(n++, m) = s * tm[1];
1177 }
1178 for (int j = 0; j <= pm1; j++)
1179 {
1180 T(n++, m) =
1181 shape_x(pm1-j)*shape_y(j)*((ip.y - c)*tm[0] - (ip.x - c)*tm[1]);
1182 }
1183 }
1184
1185 Ti.Factor(T);
1186 // mfem::out << "ND_TriangleElement(" << p << ") : "; Ti.TestInversion();
1187}
1188
1190 DenseMatrix &shape) const
1191{
1192 const int pm1 = order - 1;
1193
1194#ifdef MFEM_THREAD_SAFE
1195 const int p = order;
1196 Vector shape_x(p), shape_y(p), shape_l(p);
1197 DenseMatrix u(dof, dim);
1198#endif
1199
1200 poly1d.CalcBasis(pm1, ip.x, shape_x);
1201 poly1d.CalcBasis(pm1, ip.y, shape_y);
1202 poly1d.CalcBasis(pm1, 1. - ip.x - ip.y, shape_l);
1203
1204 int n = 0;
1205 for (int j = 0; j <= pm1; j++)
1206 for (int i = 0; i + j <= pm1; i++)
1207 {
1208 real_t s = shape_x(i)*shape_y(j)*shape_l(pm1-i-j);
1209 u(n,0) = s; u(n,1) = 0; n++;
1210 u(n,0) = 0; u(n,1) = s; n++;
1211 }
1212 for (int j = 0; j <= pm1; j++)
1213 {
1214 real_t s = shape_x(pm1-j)*shape_y(j);
1215 u(n,0) = s*(ip.y - c);
1216 u(n,1) = -s*(ip.x - c);
1217 n++;
1218 }
1219
1220 Ti.Mult(u, shape);
1221}
1222
1224 DenseMatrix &curl_shape) const
1225{
1226 const int pm1 = order - 1;
1227
1228#ifdef MFEM_THREAD_SAFE
1229 const int p = order;
1230 Vector shape_x(p), shape_y(p), shape_l(p);
1231 Vector dshape_x(p), dshape_y(p), dshape_l(p);
1232 Vector curlu(dof);
1233#endif
1234
1235 poly1d.CalcBasis(pm1, ip.x, shape_x, dshape_x);
1236 poly1d.CalcBasis(pm1, ip.y, shape_y, dshape_y);
1237 poly1d.CalcBasis(pm1, 1. - ip.x - ip.y, shape_l, dshape_l);
1238
1239 int n = 0;
1240 for (int j = 0; j <= pm1; j++)
1241 for (int i = 0; i + j <= pm1; i++)
1242 {
1243 int l = pm1-i-j;
1244 const real_t dx = (dshape_x(i)*shape_l(l) -
1245 shape_x(i)*dshape_l(l)) * shape_y(j);
1246 const real_t dy = (dshape_y(j)*shape_l(l) -
1247 shape_y(j)*dshape_l(l)) * shape_x(i);
1248
1249 curlu(n++) = -dy;
1250 curlu(n++) = dx;
1251 }
1252
1253 for (int j = 0; j <= pm1; j++)
1254 {
1255 int i = pm1 - j;
1256 // curl of shape_x(i)*shape_y(j) * (ip.y - c, -(ip.x - c), 0):
1257 curlu(n++) = -((dshape_x(i)*(ip.x - c) + shape_x(i)) * shape_y(j) +
1258 (dshape_y(j)*(ip.y - c) + shape_y(j)) * shape_x(i));
1259 }
1260
1261 Vector curl2d(curl_shape.Data(),dof);
1262 Ti.Mult(curlu, curl2d);
1263}
1264
1265
1266const real_t ND_SegmentElement::tk[1] = { 1. };
1267
1268ND_SegmentElement::ND_SegmentElement(const int p, const int ob_type)
1269 : VectorTensorFiniteElement(1, p, p - 1, ob_type, H_CURL,
1270 DofMapType::L2_DOF_MAP),
1271 dof2tk(dof)
1272{
1273 if (obasis1d.IsIntegratedType()) { is_nodal = false; }
1274
1275 const real_t *op = poly1d.OpenPoints(p - 1, ob_type);
1276
1277 // set dof2tk and Nodes
1278 for (int i = 0; i < p; i++)
1279 {
1280 dof2tk[i] = 0;
1281 Nodes.IntPoint(i).x = op[i];
1282 }
1283}
1284
1286 DenseMatrix &shape) const
1287{
1288 Vector vshape(shape.Data(), dof);
1289
1290 obasis1d.Eval(ip.x, vshape);
1291}
1292
1293const real_t ND_WedgeElement::tk[15] =
1294{ 1.,0.,0., -1.,1.,0., 0.,-1.,0., 0.,0.,1., 0.,1.,0. };
1295
1297 const int cb_type,
1298 const int ob_type)
1299 : VectorFiniteElement(3, Geometry::PRISM,
1300 3 * p * ((p + 1) * (p + 2))/2, p,
1301 H_CURL, FunctionSpace::Qk),
1302 dof2tk(dof),
1303 t_dof(dof),
1304 s_dof(dof),
1305 doftrans(p),
1306 H1TriangleFE(p, cb_type),
1307 NDTriangleFE(p),
1308 H1SegmentFE(p, cb_type),
1309 NDSegmentFE(p, ob_type)
1310{
1311 MFEM_ASSERT(H1TriangleFE.GetDof() * NDSegmentFE.GetDof() +
1312 NDTriangleFE.GetDof() * H1SegmentFE.GetDof() == dof,
1313 "Mismatch in number of degrees of freedom "
1314 "when building ND_WedgeElement!");
1315
1316#ifndef MFEM_THREAD_SAFE
1317 t1_shape.SetSize(H1TriangleFE.GetDof());
1318 s1_shape.SetSize(H1SegmentFE.GetDof());
1319 tn_shape.SetSize(NDTriangleFE.GetDof(), 2);
1320 sn_shape.SetSize(NDSegmentFE.GetDof(), 1);
1321 t1_dshape.SetSize(H1TriangleFE.GetDof(), 2);
1322 s1_dshape.SetSize(H1SegmentFE.GetDof(), 1);
1323 tn_dshape.SetSize(NDTriangleFE.GetDof(), 1);
1324#endif
1325
1326 const int pm1 = p - 1, pm2 = p - 2;
1327
1328 const IntegrationRule &t1_n = H1TriangleFE.GetNodes();
1329 const IntegrationRule &tn_n = NDTriangleFE.GetNodes();
1330 const IntegrationRule &s1_n = H1SegmentFE.GetNodes();
1331 const IntegrationRule &sn_n = NDSegmentFE.GetNodes();
1332
1333 // edges
1334 int o = 0;
1335 for (int i = 0; i < p; i++) // (0,1)
1336 {
1337 t_dof[o] = i; s_dof[o] = 0; dof2tk[o] = 0;
1338 const IntegrationPoint & t_ip = tn_n.IntPoint(t_dof[o]);
1339 Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, s1_n.IntPoint(s_dof[o]).x);
1340 o++;
1341 }
1342 for (int i = 0; i < p; i++) // (1,2)
1343 {
1344 t_dof[o] = p + i; s_dof[o] = 0; dof2tk[o] = 1;
1345 const IntegrationPoint & t_ip = tn_n.IntPoint(t_dof[o]);
1346 Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, s1_n.IntPoint(s_dof[o]).x);
1347 o++;
1348 }
1349 for (int i = 0; i < p; i++) // (2,0)
1350 {
1351 t_dof[o] = 2 * p + i; s_dof[o] = 0; dof2tk[o] = 2;
1352 const IntegrationPoint & t_ip = tn_n.IntPoint(t_dof[o]);
1353 Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, s1_n.IntPoint(s_dof[o]).x);
1354 o++;
1355 }
1356 for (int i = 0; i < p; i++) // (3,4)
1357 {
1358 t_dof[o] = i; s_dof[o] = 1; dof2tk[o] = 0;
1359 const IntegrationPoint & t_ip = tn_n.IntPoint(t_dof[o]);
1360 Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, s1_n.IntPoint(s_dof[o]).x);
1361 o++;
1362 }
1363 for (int i = 0; i < p; i++) // (4,5)
1364 {
1365 t_dof[o] = p + i; s_dof[o] = 1; dof2tk[o] = 1;
1366 const IntegrationPoint & t_ip = tn_n.IntPoint(t_dof[o]);
1367 Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, s1_n.IntPoint(s_dof[o]).x);
1368 o++;
1369 }
1370 for (int i = 0; i < p; i++) // (5,3)
1371 {
1372 t_dof[o] = 2 * p + i; s_dof[o] = 1; dof2tk[o] = 2;
1373 const IntegrationPoint & t_ip = tn_n.IntPoint(t_dof[o]);
1374 Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, s1_n.IntPoint(s_dof[o]).x);
1375 o++;
1376 }
1377 for (int i = 0; i < p; i++) // (0,3)
1378 {
1379 t_dof[o] = 0; s_dof[o] = i; dof2tk[o] = 3;
1380 const IntegrationPoint & t_ip = t1_n.IntPoint(t_dof[o]);
1381 Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, sn_n.IntPoint(s_dof[o]).x);
1382 o++;
1383 }
1384 for (int i = 0; i < p; i++) // (1,4)
1385 {
1386 t_dof[o] = 1; s_dof[o] = i; dof2tk[o] = 3;
1387 const IntegrationPoint & t_ip = t1_n.IntPoint(t_dof[o]);
1388 Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, sn_n.IntPoint(s_dof[o]).x);
1389 o++;
1390 }
1391 for (int i = 0; i < p; i++) // (2,5)
1392 {
1393 t_dof[o] = 2; s_dof[o] = i; dof2tk[o] = 3;
1394 const IntegrationPoint & t_ip = t1_n.IntPoint(t_dof[o]);
1395 Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, sn_n.IntPoint(s_dof[o]).x);
1396 o++;
1397 }
1398
1399 // faces
1400 // (0,2,1) -- bottom
1401 int l = 0;
1402 for (int j = 0; j <= pm2; j++)
1403 for (int i = 0; i + j <= pm2; i++)
1404 {
1405 l = j + ( 2 * p - 1 - i) * i / 2;
1406 t_dof[o] = 3 * p + 2*l+1; s_dof[o] = 0; dof2tk[o] = 4;
1407 const IntegrationPoint & t_ip0 = tn_n.IntPoint(t_dof[o]);
1408 Nodes.IntPoint(o).Set3(t_ip0.x, t_ip0.y, s1_n.IntPoint(s_dof[o]).x);
1409 o++;
1410 t_dof[o] = 3 * p + 2*l; s_dof[o] = 0; dof2tk[o] = 0;
1411 const IntegrationPoint & t_ip1 = tn_n.IntPoint(t_dof[o]);
1412 Nodes.IntPoint(o).Set3(t_ip1.x, t_ip1.y, s1_n.IntPoint(s_dof[o]).x);
1413 o++;
1414 }
1415 // (3,4,5) -- top
1416 int m = 0;
1417 for (int j = 0; j <= pm2; j++)
1418 for (int i = 0; i + j <= pm2; i++)
1419 {
1420 t_dof[o] = 3 * p + m; s_dof[o] = 1; dof2tk[o] = 0; m++;
1421 const IntegrationPoint & t_ip0 = tn_n.IntPoint(t_dof[o]);
1422 Nodes.IntPoint(o).Set3(t_ip0.x, t_ip0.y, s1_n.IntPoint(s_dof[o]).x);
1423 o++;
1424 t_dof[o] = 3 * p + m; s_dof[o] = 1; dof2tk[o] = 4; m++;
1425 const IntegrationPoint & t_ip1 = tn_n.IntPoint(t_dof[o]);
1426 Nodes.IntPoint(o).Set3(t_ip1.x, t_ip1.y, s1_n.IntPoint(s_dof[o]).x);
1427 o++;
1428 }
1429 // (0, 1, 4, 3) -- xz plane
1430 for (int j = 2; j <= p; j++)
1431 for (int i = 0; i < p; i++)
1432 {
1433 t_dof[o] = i; s_dof[o] = j; dof2tk[o] = 0;
1434 const IntegrationPoint & t_ip = tn_n.IntPoint(t_dof[o]);
1435 Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, s1_n.IntPoint(s_dof[o]).x);
1436 o++;
1437 }
1438 for (int j = 0; j < p; j++)
1439 for (int i = 0; i < pm1; i++)
1440 {
1441 t_dof[o] = 3 + i; s_dof[o] = j; dof2tk[o] = 3;
1442 const IntegrationPoint & t_ip = t1_n.IntPoint(t_dof[o]);
1443 Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, sn_n.IntPoint(s_dof[o]).x);
1444 o++;
1445 }
1446 // (1, 2, 5, 4) -- (y-x)z plane
1447 for (int j = 2; j <= p; j++)
1448 for (int i = 0; i < p; i++)
1449 {
1450 t_dof[o] = p + i; s_dof[o] = j; dof2tk[o] = 1;
1451 const IntegrationPoint & t_ip = tn_n.IntPoint(t_dof[o]);
1452 Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, s1_n.IntPoint(s_dof[o]).x);
1453 o++;
1454 }
1455 for (int j = 0; j < p; j++)
1456 for (int i = 0; i < pm1; i++)
1457 {
1458 t_dof[o] = p + 2 + i; s_dof[o] = j; dof2tk[o] = 3;
1459 const IntegrationPoint & t_ip = t1_n.IntPoint(t_dof[o]);
1460 Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, sn_n.IntPoint(s_dof[o]).x);
1461 o++;
1462 }
1463 // (2, 0, 3, 5) -- yz plane
1464 for (int j = 2; j <= p; j++)
1465 for (int i = 0; i < p; i++)
1466 {
1467 t_dof[o] = 2 * p + i; s_dof[o] = j; dof2tk[o] = 2;
1468 const IntegrationPoint & t_ip = tn_n.IntPoint(t_dof[o]);
1469 Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, s1_n.IntPoint(s_dof[o]).x);
1470 o++;
1471 }
1472 for (int j = 0; j < p; j++)
1473 for (int i = 0; i < pm1; i++)
1474 {
1475 t_dof[o] = 2 * p + 1 + i; s_dof[o] = j; dof2tk[o] = 3;
1476 const IntegrationPoint & t_ip = t1_n.IntPoint(t_dof[o]);
1477 Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, sn_n.IntPoint(s_dof[o]).x);
1478 o++;
1479 }
1480
1481 // interior
1482 for (int k = 2; k <= p; k++)
1483 {
1484 l = 0;
1485 for (int j = 0; j <= pm2; j++)
1486 for (int i = 0; i + j <= pm2; i++)
1487 {
1488 t_dof[o] = 3 * p + l; s_dof[o] = k; dof2tk[o] = 0; l++;
1489 const IntegrationPoint & t_ip0 = tn_n.IntPoint(t_dof[o]);
1490 Nodes.IntPoint(o).Set3(t_ip0.x, t_ip0.y, s1_n.IntPoint(s_dof[o]).x);
1491 o++;
1492 t_dof[o] = 3 * p + l; s_dof[o] = k; dof2tk[o] = 4; l++;
1493 const IntegrationPoint & t_ip1 = tn_n.IntPoint(t_dof[o]);
1494 Nodes.IntPoint(o).Set3(t_ip1.x, t_ip1.y, s1_n.IntPoint(s_dof[o]).x);
1495 o++;
1496 }
1497 }
1498 for (int k = 0; k < p; k++)
1499 {
1500 l = 0;
1501 for (int j = 0; j < pm2; j++)
1502 for (int i = 0; i + j < pm2; i++)
1503 {
1504 t_dof[o] = 3 * p + l; s_dof[o] = k; dof2tk[o] = 3; l++;
1505 const IntegrationPoint & t_ip = t1_n.IntPoint(t_dof[o]);
1506 Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, sn_n.IntPoint(s_dof[o]).x);
1507 o++;
1508 }
1509 }
1510}
1511
1513 DenseMatrix &shape) const
1514{
1515#ifdef MFEM_THREAD_SAFE
1516 Vector t1_shape(H1TriangleFE.GetDof());
1517 Vector s1_shape(H1SegmentFE.GetDof());
1518 DenseMatrix tn_shape(NDTriangleFE.GetDof(), 2);
1519 DenseMatrix sn_shape(NDSegmentFE.GetDof(), 1);
1520#endif
1521
1522 IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
1523
1524 H1TriangleFE.CalcShape(ip, t1_shape);
1525 NDTriangleFE.CalcVShape(ip, tn_shape);
1526 H1SegmentFE.CalcShape(ipz, s1_shape);
1527 NDSegmentFE.CalcVShape(ipz, sn_shape);
1528
1529 for (int i=0; i<dof; i++)
1530 {
1531 if ( dof2tk[i] != 3 )
1532 {
1533 shape(i, 0) = tn_shape(t_dof[i], 0) * s1_shape[s_dof[i]];
1534 shape(i, 1) = tn_shape(t_dof[i], 1) * s1_shape[s_dof[i]];
1535 shape(i, 2) = 0.0;
1536 }
1537 else
1538 {
1539 shape(i, 0) = 0.0;
1540 shape(i, 1) = 0.0;
1541 shape(i, 2) = t1_shape[t_dof[i]] * sn_shape(s_dof[i], 0);
1542 }
1543 }
1544}
1545
1547 DenseMatrix &curl_shape) const
1548{
1549#ifdef MFEM_THREAD_SAFE
1550 Vector s1_shape(H1SegmentFE.GetDof());
1551 DenseMatrix t1_dshape(H1TriangleFE.GetDof(), 2);
1552 DenseMatrix s1_dshape(H1SegmentFE.GetDof(), 1);
1553 DenseMatrix tn_shape(NDTriangleFE.GetDof(), 2);
1554 DenseMatrix sn_shape(NDSegmentFE.GetDof(), 1);
1555 DenseMatrix tn_dshape(NDTriangleFE.GetDof(), 1);
1556#endif
1557
1558 IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
1559
1560 H1TriangleFE.CalcDShape(ip, t1_dshape);
1561 H1SegmentFE.CalcShape(ipz, s1_shape);
1562 H1SegmentFE.CalcDShape(ipz, s1_dshape);
1563 NDTriangleFE.CalcVShape(ip, tn_shape);
1564 NDTriangleFE.CalcCurlShape(ip, tn_dshape);
1565 NDSegmentFE.CalcVShape(ipz, sn_shape);
1566
1567 for (int i=0; i<dof; i++)
1568 {
1569 if ( dof2tk[i] != 3 )
1570 {
1571 curl_shape(i, 0) = -tn_shape(t_dof[i], 1) * s1_dshape(s_dof[i], 0);
1572 curl_shape(i, 1) = tn_shape(t_dof[i], 0) * s1_dshape(s_dof[i], 0);
1573 curl_shape(i, 2) = tn_dshape(t_dof[i], 0) * s1_shape[s_dof[i]];
1574 }
1575 else
1576 {
1577 curl_shape(i, 0) = t1_dshape(t_dof[i], 1) * sn_shape(s_dof[i], 0);
1578 curl_shape(i, 1) = -t1_dshape(t_dof[i], 0) * sn_shape(s_dof[i], 0);
1579 curl_shape(i, 2) = 0.0;
1580 }
1581 }
1582}
1583
1584const real_t ND_FuentesPyramidElement::tk[27] =
1585{
1586 1., 0., 0., 0., 1., 0., 0., 0., 1.,
1587 -1., 0., 1., -1.,-1., 1., 0.,-1., 1.,
1588 -1., 0., 0., 0.,-1., 0., -M_SQRT1_2,-M_SQRT1_2,M_SQRT2
1589};
1590
1592 const int cb_type,
1593 const int ob_type)
1594 : VectorFiniteElement(3, Geometry::PYRAMID, p * (3 * p * p + 5), p,
1595 H_CURL, FunctionSpace::Uk),
1596 dof2tk(dof), doftrans(p)
1597{
1598 zmax = 0.0;
1599
1600 const real_t *eop = poly1d.OpenPoints(p - 1);
1601 const real_t *top = (p > 1) ? poly1d.OpenPoints(p - 2) : NULL;
1602 const real_t *qop = poly1d.OpenPoints(p - 1, ob_type);
1603 const real_t *qcp = poly1d.ClosedPoints(p, cb_type);
1604
1605 const int pm2 = p - 2;
1606
1607#ifndef MFEM_THREAD_SAFE
1608 tmp_E_E_ij.SetSize(p, dim);
1609 tmp_dE_E_ij.SetSize(p, dim);
1610 tmp_E_Q1_ijk.SetSize(p, p + 1, dim);
1611 tmp_dE_Q1_ijk.SetSize(p, p + 1, dim);
1612 tmp_E_Q2_ijk.SetSize(p, p + 1, dim);
1613 tmp_dE_Q2_ijk.SetSize(p, p + 1, dim);
1614 tmp_E_T_ijk.SetSize(p - 1, p, dim);
1615 tmp_dE_T_ijk.SetSize(p - 1, p, dim);
1616 tmp_phi_Q1_ij.SetSize(p + 1, p + 1);
1617 tmp_dphi_Q1_ij.SetSize(p + 1, p + 1, dim);
1618 tmp_phi_Q2_ij.SetSize(p + 1, p + 1);
1619 tmp_dphi_Q2_ij.SetSize(p + 1, p + 1, dim);
1620 tmp_phi_E_i.SetSize(p + 1);
1621 tmp_dphi_E_i.SetSize(p + 1, dim);
1622 u.SetSize(dof, dim);
1623 curlu.SetSize(dof, dim);
1624#else
1625 DenseMatrix tmp_E_E_ij(p, dim);
1626 DenseTensor tmp_E_Q1_ijk(p, p + 1, dim);
1627 DenseTensor tmp_dE_Q1_ijk(p, p + 1, dim);
1628 DenseTensor tmp_E_Q2_ijk(p, p + 1, dim);
1629 DenseTensor tmp_dE_Q2_ijk(p, p + 1, dim);
1630 DenseTensor tmp_E_T_ijk(p - 1, p, dim);
1631 DenseTensor tmp_dE_T_ijk(p - 1, p, dim);
1632 DenseMatrix tmp_phi_Q1_ij(p + 1, p + 1);
1633 DenseTensor tmp_dphi_Q1_ij(p + 1, p + 1, dim);
1634 DenseMatrix tmp_phi_Q2_ij(p + 1, p + 1);
1635 Vector tmp_phi_E_i(p + 1);
1636 DenseMatrix tmp_dphi_E_i(p + 1, dim);
1637 DenseMatrix u(dof, dim);
1638#endif
1639
1640 int o = 0;
1641
1642 // edges
1643 for (int i = 0; i < p; i++) // (0, 1)
1644 {
1645 Nodes.IntPoint(o).Set3(eop[i], 0., 0.);
1646 dof2tk[o++] = 0;
1647 }
1648 for (int i = 0; i < p; i++) // (1, 2)
1649 {
1650 Nodes.IntPoint(o).Set3(1., eop[i], 0.);
1651 dof2tk[o++] = 1;
1652 }
1653 for (int i = 0; i < p; i++) // (3, 2)
1654 {
1655 Nodes.IntPoint(o).Set3(eop[i], 1., 0.);
1656 dof2tk[o++] = 0;
1657 }
1658 for (int i = 0; i < p; i++) // (0, 3)
1659 {
1660 Nodes.IntPoint(o).Set3(0., eop[i], 0.);
1661 dof2tk[o++] = 1;
1662 }
1663 for (int i = 0; i < p; i++) // (0, 4)
1664 {
1665 Nodes.IntPoint(o).Set3(0., 0., eop[i]);
1666 dof2tk[o++] = 2;
1667 }
1668 for (int i = 0; i < p; i++) // (1, 4)
1669 {
1670 Nodes.IntPoint(o).Set3(1. - eop[i], 0., eop[i]);
1671 dof2tk[o++] = 3;
1672 }
1673 for (int i = 0; i < p; i++) // (2, 4)
1674 {
1675 Nodes.IntPoint(o).Set3(1. - eop[i], 1. - eop[i], eop[i]);
1676 dof2tk[o++] = 4;
1677 }
1678 for (int i = 0; i < p; i++) // (3, 4)
1679 {
1680 Nodes.IntPoint(o).Set3(0., 1. - eop[i], eop[i]);
1681 dof2tk[o++] = 5;
1682 }
1683
1684 // quadrilateral face (3, 2, 1, 0)
1685 // x-components
1686 for (int j = 1; j < p; j++)
1687 for (int i = 0; i < p; i++)
1688 {
1689 Nodes.IntPoint(o).Set3(qop[i], qcp[p-j], 0.);
1690 dof2tk[o++] = 0; // (1 0 0)
1691 }
1692
1693 // y-components
1694 for (int j = 0; j < p; j++)
1695 for (int i = 1; i < p; i++)
1696 {
1697 Nodes.IntPoint(o).Set3(qcp[i], qop[p-1-j], 0.);
1698 dof2tk[o++] = 7; // (0 -1 0)
1699 }
1700
1701 // triangular faces
1702 for (int j = 0; j <= pm2; j++) // (0, 1, 4)
1703 for (int i = 0; i + j <= pm2; i++)
1704 {
1705 real_t w = top[i] + top[j] + top[pm2-i-j];
1706 Nodes.IntPoint(o).Set3(top[i]/w, 0., top[j]/w);
1707 dof2tk[o++] = 0;
1708 Nodes.IntPoint(o).Set3(top[i]/w, 0., top[j]/w);
1709 dof2tk[o++] = 2;
1710 }
1711 for (int j = 0; j <= pm2; j++) // (1, 2, 4)
1712 for (int i = 0; i + j <= pm2; i++)
1713 {
1714 real_t w = top[i] + top[j] + top[pm2-i-j];
1715 Nodes.IntPoint(o).Set3((top[i] + top[pm2-i-j])/w, top[i]/w, top[j]/w);
1716 dof2tk[o++] = 1;
1717 Nodes.IntPoint(o).Set3((top[i] + top[pm2-i-j])/w, top[i]/w, top[j]/w);
1718 dof2tk[o++] = 3;
1719 }
1720 for (int j = 0; j <= pm2; j++) // (2, 3, 4)
1721 for (int i = 0; i + j <= pm2; i++)
1722 {
1723 real_t w = top[i] + top[j] + top[pm2-i-j];
1724 Nodes.IntPoint(o).Set3(top[pm2-i-j]/w, (top[i] + top[pm2-i-j])/w,
1725 top[j]/w);
1726 dof2tk[o++] = 6;
1727 Nodes.IntPoint(o).Set3(top[pm2-i-j]/w, (top[i] + top[pm2-i-j])/w,
1728 top[j]/w);
1729 dof2tk[o++] = 4;
1730 }
1731 for (int j = 0; j <= pm2; j++) // (3, 0, 4)
1732 for (int i = 0; i + j <= pm2; i++)
1733 {
1734 real_t w = top[i] + top[j] + top[pm2-i-j];
1735 Nodes.IntPoint(o).Set3(0., top[pm2-i-j]/w, top[j]/w);
1736 dof2tk[o++] = 7;
1737 Nodes.IntPoint(o).Set3(0., top[pm2-i-j]/w, top[j]/w);
1738 dof2tk[o++] = 5;
1739 }
1740
1741 // interior
1742 // x-components
1743 for (int k = 1; k < p; k++)
1744 for (int j = 1; j < p; j++)
1745 for (int i = 0; i < p; i++)
1746 {
1747 real_t w = 1.0 - qcp[k];
1748 Nodes.IntPoint(o).Set3(qop[i]*w, qcp[j]*w, qcp[k]);
1749 dof2tk[o++] = 0;
1750 }
1751 // y-components
1752 for (int k = 1; k < p; k++)
1753 for (int j = 0; j < p; j++)
1754 for (int i = 1; i < p; i++)
1755 {
1756 real_t w = 1.0 - qcp[k];
1757 Nodes.IntPoint(o).Set3(qcp[i]*w, qop[j]*w, qcp[k]);
1758 dof2tk[o++] = 1;
1759 }
1760 // z-components
1761 for (int k = 0; k < p; k++)
1762 for (int j = 1; j < p; j++)
1763 for (int i = 1; i < p; i++)
1764 {
1765 real_t w = 1.0 - qop[k];
1766 Nodes.IntPoint(o).Set3(qcp[i]*w, qcp[j]*w, qop[k]);
1767 dof2tk[o++] = 8;
1768 }
1769
1770 DenseMatrix T(dof);
1771
1772 for (int m = 0; m < dof; m++)
1773 {
1774 const IntegrationPoint &ip = Nodes.IntPoint(m);
1775 calcBasis(p, ip, tmp_E_E_ij, tmp_E_Q1_ijk, tmp_E_Q2_ijk, tmp_E_T_ijk,
1776 tmp_phi_Q1_ij, tmp_dphi_Q1_ij, tmp_phi_Q2_ij,
1777 tmp_phi_E_i, tmp_dphi_E_i, u);
1778
1779 const Vector tm({tk[3*dof2tk[m]], tk[3*dof2tk[m]+1], tk[3*dof2tk[m]+2]});
1780 u.Mult(tm, T.GetColumn(m));
1781 }
1782
1783 Ti.Factor(T);
1784}
1785
1787 DenseMatrix &shape) const
1788{
1789 const int p = order;
1790
1791#ifdef MFEM_THREAD_SAFE
1792 DenseMatrix tmp_E_E_ij(p, dim);
1793 DenseTensor tmp_E_Q1_ijk(p, p + 1, dim);
1794 DenseTensor tmp_E_Q2_ijk(p, p + 1, dim);
1795 DenseTensor tmp_E_T_ijk(p - 1, p, dim);
1796 DenseMatrix tmp_phi_Q1_ij(p + 1, p + 1);
1797 DenseTensor tmp_dphi_Q1_ij(p + 1, p + 1, dim);
1798 DenseMatrix tmp_phi_Q2_ij(p + 1, p + 1);
1799 Vector tmp_phi_E_i(p + 1);
1800 DenseMatrix tmp_dphi_E_i(p + 1, dim);
1801 DenseMatrix u(dof, dim);
1802#endif
1803
1804 calcBasis(p, ip, tmp_E_E_ij, tmp_E_Q1_ijk, tmp_E_Q2_ijk, tmp_E_T_ijk,
1805 tmp_phi_Q1_ij, tmp_dphi_Q1_ij, tmp_phi_Q2_ij,
1806 tmp_phi_E_i, tmp_dphi_E_i, u);
1807
1808 Ti.Mult(u, shape);
1809}
1810
1812 DenseMatrix &curl_shape) const
1813{
1814 const int p = order;
1815
1816#ifdef MFEM_THREAD_SAFE
1817 DenseMatrix tmp_E_E_ij(p, dim);
1818 DenseMatrix tmp_dE_E_ij(p, dim);
1819 DenseTensor tmp_E_Q1_ijk(p, p + 1, dim);
1820 DenseTensor tmp_dE_Q1_ijk(p, p + 1, dim);
1821 DenseTensor tmp_E_Q2_ijk(p, p + 1, dim);
1822 DenseTensor tmp_dE_Q2_ijk(p, p + 1, dim);
1823 DenseTensor tmp_E_T_ijk(p - 1, p, dim);
1824 DenseTensor tmp_dE_T_ijk(p - 1, p, dim);
1825 DenseMatrix tmp_phi_Q2_ij(p + 1, p + 1);
1826 DenseTensor tmp_dphi_Q2_ij(p + 1, p + 1, dim);
1827 Vector tmp_phi_E_i(p + 1);
1828 DenseMatrix tmp_dphi_E_i(p + 1, dim);
1829 DenseMatrix curlu(dof, dim);
1830#endif
1831
1832 calcCurlBasis(p, ip, tmp_E_E_ij, tmp_dE_E_ij, tmp_E_Q1_ijk, tmp_dE_Q1_ijk,
1833 tmp_E_Q2_ijk, tmp_dE_Q2_ijk, tmp_E_T_ijk, tmp_dE_T_ijk,
1834 tmp_phi_Q2_ij, tmp_dphi_Q2_ij, tmp_phi_E_i, tmp_dphi_E_i,
1835 curlu);
1836
1837 Ti.Mult(curlu, curl_shape);
1838}
1839
1841 DenseMatrix &shape) const
1842{
1843 const int p = order;
1844
1845#ifdef MFEM_THREAD_SAFE
1846 DenseMatrix tmp_E_E_ij(p, dim);
1847 DenseTensor tmp_E_Q1_ijk(p, p + 1, dim);
1848 DenseTensor tmp_E_Q2_ijk(p, p + 1, dim);
1849 DenseTensor tmp_E_T_ijk(p - 1, p, dim);
1850 DenseMatrix tmp_phi_Q1_ij(p + 1, p + 1);
1851 DenseTensor tmp_dphi_Q1_ij(p + 1, p + 1, dim);
1852 DenseMatrix tmp_phi_Q2_ij(p + 1, p + 1);
1853 Vector tmp_phi_E_i(p + 1);
1854 DenseMatrix tmp_dphi_E_i(p + 1, dim);
1855#endif
1856
1857 calcBasis(p, ip, tmp_E_E_ij, tmp_E_Q1_ijk, tmp_E_Q2_ijk, tmp_E_T_ijk,
1858 tmp_phi_Q1_ij, tmp_dphi_Q1_ij, tmp_phi_Q2_ij,
1859 tmp_phi_E_i, tmp_dphi_E_i, shape);
1860}
1861
1863 DenseMatrix &dshape) const
1864{
1865 const int p = order;
1866
1867#ifdef MFEM_THREAD_SAFE
1868 DenseMatrix tmp_E_E_ij(p, dim);
1869 DenseMatrix tmp_dE_E_ij(p, dim);
1870 DenseTensor tmp_E_Q1_ijk(p, p + 1, dim);
1871 DenseTensor tmp_dE_Q1_ijk(p, p + 1, dim);
1872 DenseTensor tmp_E_Q2_ijk(p, p + 1, dim);
1873 DenseTensor tmp_dE_Q2_ijk(p, p + 1, dim);
1874 DenseTensor tmp_E_T_ijk(p - 1, p, dim);
1875 DenseTensor tmp_dE_T_ijk(p - 1, p, dim);
1876 DenseMatrix tmp_phi_Q2_ij(p + 1, p + 1);
1877 DenseTensor tmp_dphi_Q2_ij(p + 1, p + 1, dim);
1878 Vector tmp_phi_E_i(p + 1);
1879 DenseMatrix tmp_dphi_E_i(p + 1, dim);
1880#endif
1881
1882 calcCurlBasis(p, ip, tmp_E_E_ij, tmp_dE_E_ij, tmp_E_Q1_ijk, tmp_dE_Q1_ijk,
1883 tmp_E_Q2_ijk, tmp_dE_Q2_ijk, tmp_E_T_ijk, tmp_dE_T_ijk,
1884 tmp_phi_Q2_ij, tmp_dphi_Q2_ij, tmp_phi_E_i, tmp_dphi_E_i,
1885 dshape);
1886}
1887
1888void ND_FuentesPyramidElement::calcBasis(const int p,
1889 const IntegrationPoint &ip,
1890 DenseMatrix & E_E_ik,
1891 DenseTensor & E_Q1_ijk,
1892 DenseTensor & E_Q2_ijk,
1893 DenseTensor & E_T_ijk,
1894 DenseMatrix & phi_Q1_ij,
1895 DenseTensor & dphi_Q1_ij,
1896 DenseMatrix & phi_Q2_ij,
1897 Vector & phi_E_k,
1898 DenseMatrix & dphi_E_k,
1899 DenseMatrix &W) const
1900{
1901 real_t x = ip.x;
1902 real_t y = ip.y;
1903 real_t z = ip.z;
1904 Vector xy({x,y}), dmu(3);
1905 real_t mu, mu2;
1906
1907 if (std::fabs(1.0 - z) < apex_tol)
1908 {
1909 z = 1.0 - apex_tol;
1910 y = 0.5 * (1.0 - z);
1911 x = 0.5 * (1.0 - z);
1912 xy(0) = x; xy(1) = y;
1913 }
1914 zmax = std::max(z, zmax);
1915
1916 W = 0.0;
1917
1918 int o = 0;
1919
1920 // Mixed Edges
1921 if (z < 1.0)
1922 {
1923 // (a, b) = (1, 2), c = 0
1924 mu = mu0(z, xy, 2);
1925 E_E(p, nu01(z, xy, 1), nu01_grad_nu01(z, xy, 1), E_E_ik);
1926 for (int i=0; i<p; i++, o++)
1927 for (int k=0; k<3; k++)
1928 {
1929 W(o, k) = mu * E_E_ik(i, k);
1930 }
1931
1932 // (a, b) = (1, 2), c = 1
1933 mu = mu1(z, xy, 2);
1934 for (int i=0; i<p; i++, o++)
1935 for (int k=0; k<3; k++)
1936 {
1937 W(o, k) = mu * E_E_ik(i, k);
1938 }
1939
1940 // (a, b) = (2, 1), c = 0
1941 mu = mu0(z, xy, 1);
1942 E_E(p, nu01(z, xy, 2), nu01_grad_nu01(z, xy, 2), E_E_ik);
1943 for (int i=0; i<p; i++, o++)
1944 for (int k=0; k<3; k++)
1945 {
1946 W(o, k) = mu * E_E_ik(i, k);
1947 }
1948
1949 // (a, b) = (2, 1), c = 1
1950 mu = mu1(z, xy, 1);
1951 for (int i=0; i<p; i++, o++)
1952 for (int k=0; k<3; k++)
1953 {
1954 W(o, k) = mu * E_E_ik(i, k);
1955 }
1956 }
1957
1958 // Triangle Edges
1959 if (z < 1.0)
1960 {
1961 E_E(p, lam15(x, y, z), lam15_grad_lam15(x, y, z), E_E_ik);
1962 for (int i=0; i<p; i++, o++)
1963 for (int k=0; k<3; k++)
1964 {
1965 W(o, k) = E_E_ik(i, k);
1966 }
1967
1968 E_E(p, lam25(x, y, z), lam25_grad_lam25(x, y, z), E_E_ik);
1969 for (int i=0; i<p; i++, o++)
1970 for (int k=0; k<3; k++)
1971 {
1972 W(o, k) = E_E_ik(i, k);
1973 }
1974
1975 E_E(p, lam35(x, y, z), lam35_grad_lam35(x, y, z), E_E_ik);
1976 for (int i=0; i<p; i++, o++)
1977 for (int k=0; k<3; k++)
1978 {
1979 W(o, k) = E_E_ik(i, k);
1980 }
1981
1982 E_E(p, lam45(x, y, z), lam45_grad_lam45(x, y, z), E_E_ik);
1983 for (int i=0; i<p; i++, o++)
1984 for (int k=0; k<3; k++)
1985 {
1986 W(o, k) = E_E_ik(i, k);
1987 }
1988 }
1989
1990 // Quadrilateral Face
1991 if (z < 1.0 && p >= 2)
1992 {
1993 mu = mu0(z);
1994 mu2 = mu * mu;
1995
1996 // Family I
1997 E_Q(p, mu01(z, xy, 1), mu01_grad_mu01(z, xy, 1), mu01(z, xy, 2),
1998 E_Q1_ijk);
1999 for (int j=2; j<=p; j++)
2000 for (int i=0; i<p; i++, o++)
2001 for (int k=0; k<3; k++)
2002 {
2003 W(o, k) = mu2 * E_Q1_ijk(i, j, k);
2004 }
2005
2006 // Family II
2007 E_Q(p, mu01(z, xy, 2), mu01_grad_mu01(z, xy, 2), mu01(z, xy, 1),
2008 E_Q2_ijk);
2009 for (int j=2; j<=p; j++)
2010 for (int i=0; i<p; i++, o++)
2011 for (int k=0; k<3; k++)
2012 {
2013 W(o, k) = mu2 * E_Q2_ijk(i, j, k);
2014 }
2015 }
2016
2017 // Triangular Faces
2018 if (z < 1.0 && p >= 2)
2019 {
2020 // Family I
2021 // (a, b) = (1, 2), c = 0
2022 mu = mu0(z, xy, 2);
2023 E_T(p, nu012(z, xy, 1), nu01_grad_nu01(z, xy, 1), E_T_ijk);
2024 for (int j=1; j<p; j++)
2025 for (int i=0; i+j<p; i++, o++)
2026 for (int k=0; k<3; k++)
2027 {
2028 W(o, k) = mu * E_T_ijk(i, j, k);
2029 }
2030
2031 // (a, b) = (1, 2), c = 1
2032 mu = mu1(z, xy, 2);
2033 for (int j=1; j<p; j++)
2034 for (int i=0; i+j<p; i++, o++)
2035 for (int k=0; k<3; k++)
2036 {
2037 W(o, k) = mu * E_T_ijk(i, j, k);
2038 }
2039
2040 // (a, b) = (2, 1), c = 0
2041 mu = mu0(z, xy, 1);
2042 E_T(p, nu012(z, xy, 2), nu01_grad_nu01(z, xy, 2), E_T_ijk);
2043 for (int j=1; j<p; j++)
2044 for (int i=0; i+j<p; i++, o++)
2045 for (int k=0; k<3; k++)
2046 {
2047 W(o, k) = mu * E_T_ijk(i, j, k);
2048 }
2049
2050 // (a, b) = (2, 1), c = 1
2051 mu = mu1(z, xy, 1);
2052 for (int j=1; j<p; j++)
2053 for (int i=0; i+j<p; i++, o++)
2054 for (int k=0; k<3; k++)
2055 {
2056 W(o, k) = mu * E_T_ijk(i, j, k);
2057 }
2058
2059 // Family II
2060 // (a, b) = (1, 2), c = 0
2061 mu = mu0(z, xy, 2);
2062 E_T(p, nu120(z, xy, 1), nu12_grad_nu12(z, xy, 1), E_T_ijk);
2063 for (int j=1; j<p; j++)
2064 for (int i=0; i+j<p; i++, o++)
2065 for (int k=0; k<3; k++)
2066 {
2067 W(o, k) = mu * E_T_ijk(i, j, k);
2068 }
2069
2070 // (a, b) = (1, 2), c = 1
2071 mu = mu1(z, xy, 2);
2072 for (int j=1; j<p; j++)
2073 for (int i=0; i+j<p; i++, o++)
2074 for (int k=0; k<3; k++)
2075 {
2076 W(o, k) = mu * E_T_ijk(i, j, k);
2077 }
2078
2079 // (a, b) = (2, 1), c = 0
2080 mu = mu0(z, xy, 1);
2081 E_T(p, nu120(z, xy, 2), nu12_grad_nu12(z, xy, 2), E_T_ijk);
2082 for (int j=1; j<p; j++)
2083 for (int i=0; i+j<p; i++, o++)
2084 for (int k=0; k<3; k++)
2085 {
2086 W(o, k) = mu * E_T_ijk(i, j, k);
2087 }
2088
2089 // (a, b) = (2, 1), c = 1
2090 mu = mu1(z, xy, 1);
2091 for (int j=1; j<p; j++)
2092 for (int i=0; i+j<p; i++, o++)
2093 for (int k=0; k<3; k++)
2094 {
2095 W(o, k) = mu * E_T_ijk(i, j, k);
2096 }
2097 }
2098
2099 // Interior
2100 if (z < 1.0 && p >= 2)
2101 {
2102 // Family I
2103 phi_Q(p, mu01(z, xy, 1), grad_mu01(z, xy, 1), mu01(z, xy, 2),
2104 grad_mu01(z, xy, 2), phi_Q1_ij, dphi_Q1_ij);
2105 phi_E(p, mu01(z), grad_mu01(z), phi_E_k, dphi_E_k);
2106 for (int k=2; k<=p; k++)
2107 for (int j=2; j<=p; j++)
2108 for (int i=2; i<=p; i++, o++)
2109 for (int l=0; l<3; l++)
2110 W(o, l) = dphi_Q1_ij(i, j, l) * phi_E_k(k) +
2111 phi_Q1_ij(i, j) * dphi_E_k(k, l);
2112
2113 // Family II
2114 mu = mu0(z);
2115 for (int k=2; k<=p; k++)
2116 for (int j=2; j<=p; j++)
2117 for (int i=0; i<p; i++, o++)
2118 for (int l=0; l<3; l++)
2119 {
2120 W(o, l) = mu * E_Q1_ijk(i, j, l) * phi_E_k(k);
2121 }
2122
2123 // Family III
2124 for (int k=2; k<=p; k++)
2125 for (int j=2; j<=p; j++)
2126 for (int i=0; i<p; i++, o++)
2127 for (int l=0; l<3; l++)
2128 {
2129 W(o, l) = mu * E_Q2_ijk(i, j, l) * phi_E_k(k);
2130 }
2131
2132 // Family IV
2133 // Re-using mu from Family II
2134 dmu = grad_mu0(z);
2135 phi_Q(p, mu01(z, xy, 2), mu01(z, xy, 1), phi_Q2_ij);
2136 for (int j=2; j<=p; j++)
2137 for (int i=2; i<=p; i++, o++)
2138 {
2139 const int n = std::max(i,j);
2140 const real_t nmu = n * pow(mu, n-1);
2141 for (int l=0; l<3; l++)
2142 {
2143 W(o, l) = nmu * phi_Q2_ij(i, j) * dmu(l);
2144 }
2145 }
2146 }
2147}
2148
2149void ND_FuentesPyramidElement::calcCurlBasis(const int p,
2150 const IntegrationPoint &ip,
2151 DenseMatrix & E_E_ik,
2152 DenseMatrix & dE_E_ik,
2153 DenseTensor & E_Q1_ijk,
2154 DenseTensor & dE_Q1_ijk,
2155 DenseTensor & E_Q2_ijk,
2156 DenseTensor & dE_Q2_ijk,
2157 DenseTensor & E_T_ijk,
2158 DenseTensor & dE_T_ijk,
2159 DenseMatrix & phi_Q2_ij,
2160 DenseTensor & dphi_Q2_ij,
2161 Vector & phi_E_k,
2162 DenseMatrix & dphi_E_k,
2163 DenseMatrix & dW) const
2164{
2165 real_t x = ip.x;
2166 real_t y = ip.y;
2167 real_t z = ip.z;
2168 Vector xy({x,y}), dmu(3);
2169 Vector dmuxE(3), E(3), dphi(3), muphi(3);
2170
2171 real_t mu, mu2;
2172
2173 if (std::fabs(1.0 - z) < apex_tol)
2174 {
2175 z = 1.0 - apex_tol;
2176 y = 0.5 * (1.0 - z);
2177 x = 0.5 * (1.0 - z);
2178 xy(0) = x; xy(1) = y;
2179 }
2180 zmax = std::max(z, zmax);
2181
2182 dW = 0.0;
2183
2184 int o = 0;
2185
2186 // Mixed Edges
2187 if (z < 1.0)
2188 {
2189 // (a, b) = (1, 2), c = 0
2190 mu = mu0(z, xy, 2);
2191 dmu = grad_mu0(z, xy, 2);
2192 E_E(p, nu01(z, xy, 1), grad_nu01(z, xy, 1), E_E_ik, dE_E_ik);
2193 for (int i=0; i<p; i++, o++)
2194 {
2195 E(0) = E_E_ik(i, 0); E(1) = E_E_ik(i, 1); E(2) = E_E_ik(i, 2);
2196 dmu.cross3D(E, dmuxE);
2197 for (int k=0; k<3; k++)
2198 {
2199 dW(o, k) = mu * dE_E_ik(i, k) + dmuxE(k);
2200 }
2201 }
2202
2203 // (a, b) = (1, 2), c = 1
2204 mu = mu1(z, xy, 2);
2205 dmu = grad_mu1(z, xy, 2);
2206 for (int i=0; i<p; i++, o++)
2207 {
2208 E(0) = E_E_ik(i, 0); E(1) = E_E_ik(i, 1); E(2) = E_E_ik(i, 2);
2209 dmu.cross3D(E, dmuxE);
2210 for (int k=0; k<3; k++)
2211 {
2212 dW(o, k) = mu * dE_E_ik(i, k) + dmuxE(k);
2213 }
2214 }
2215
2216 // (a, b) = (2, 1), c = 0
2217 mu = mu0(z, xy, 1);
2218 dmu = grad_mu0(z, xy, 1);
2219 E_E(p, nu01(z, xy, 2), grad_nu01(z, xy, 2), E_E_ik, dE_E_ik);
2220 for (int i=0; i<p; i++, o++)
2221 {
2222 E(0) = E_E_ik(i, 0); E(1) = E_E_ik(i, 1); E(2) = E_E_ik(i, 2);
2223 dmu.cross3D(E, dmuxE);
2224 for (int k=0; k<3; k++)
2225 {
2226 dW(o, k) = mu * dE_E_ik(i, k) + dmuxE(k);
2227 }
2228 }
2229
2230 // (a, b) = (2, 1), c = 1
2231 mu = mu1(z, xy, 1);
2232 dmu = grad_mu1(z, xy, 1);
2233 for (int i=0; i<p; i++, o++)
2234 {
2235 E(0) = E_E_ik(i, 0); E(1) = E_E_ik(i, 1); E(2) = E_E_ik(i, 2);
2236 dmu.cross3D(E, dmuxE);
2237 for (int k=0; k<3; k++)
2238 {
2239 dW(o, k) = mu * dE_E_ik(i, k) + dmuxE(k);
2240 }
2241 }
2242 }
2243
2244 // Triangle Edges
2245 if (z < 1.0)
2246 {
2247 E_E(p, lam15(x, y, z), grad_lam15(x, y, z), E_E_ik, dE_E_ik);
2248 for (int i=0; i<p; i++, o++)
2249 for (int k=0; k<3; k++)
2250 {
2251 dW(o, k) = dE_E_ik(i, k);
2252 }
2253
2254 E_E(p, lam25(x, y, z), grad_lam25(x, y, z), E_E_ik, dE_E_ik);
2255 for (int i=0; i<p; i++, o++)
2256 for (int k=0; k<3; k++)
2257 {
2258 dW(o, k) = dE_E_ik(i, k);
2259 }
2260
2261 E_E(p, lam35(x, y, z), grad_lam35(x, y, z), E_E_ik, dE_E_ik);
2262 for (int i=0; i<p; i++, o++)
2263 for (int k=0; k<3; k++)
2264 {
2265 dW(o, k) = dE_E_ik(i, k);
2266 }
2267
2268 E_E(p, lam45(x, y, z), grad_lam45(x, y, z), E_E_ik, dE_E_ik);
2269 for (int i=0; i<p; i++, o++)
2270 for (int k=0; k<3; k++)
2271 {
2272 dW(o, k) = dE_E_ik(i, k);
2273 }
2274 }
2275
2276 // Quadrilateral Face
2277 if (z < 1.0 && p >= 2)
2278 {
2279 mu = mu0(z);
2280 mu2 = mu * mu;
2281 dmu = grad_mu0(z);
2282
2283 // Family I
2284 E_Q(p, mu01(z, xy, 1), grad_mu01(z, xy, 1),
2285 mu01(z, xy, 2), grad_mu01(z, xy, 2), E_Q1_ijk, dE_Q1_ijk);
2286 for (int j=2; j<=p; j++)
2287 for (int i=0; i<p; i++, o++)
2288 {
2289 E(0) = E_Q1_ijk(i, j, 0);
2290 E(1) = E_Q1_ijk(i, j, 1);
2291 E(2) = E_Q1_ijk(i, j, 2);
2292 dmu.cross3D(E, dmuxE);
2293 for (int k=0; k<3; k++)
2294 {
2295 dW(o, k) = mu2 * dE_Q1_ijk(i, j, k) + 2.0 * mu * dmuxE(k);
2296 }
2297 }
2298
2299 // Family II
2300 E_Q(p, mu01(z, xy, 2), grad_mu01(z, xy, 2),
2301 mu01(z, xy, 1), grad_mu01(z, xy, 1), E_Q2_ijk, dE_Q2_ijk);
2302 for (int j=2; j<=p; j++)
2303 for (int i=0; i<p; i++, o++)
2304 {
2305 E(0) = E_Q2_ijk(i, j, 0);
2306 E(1) = E_Q2_ijk(i, j, 1);
2307 E(2) = E_Q2_ijk(i, j, 2);
2308 dmu.cross3D(E, dmuxE);
2309 for (int k=0; k<3; k++)
2310 {
2311 dW(o, k) = mu2 * dE_Q2_ijk(i, j, k) + 2.0 * mu * dmuxE(k);
2312 }
2313 }
2314 }
2315
2316 // Triangular Faces
2317 if (z < 1.0 && p >= 2)
2318 {
2319 // Family I
2320 // (a, b) = (1, 2), c = 0
2321 mu = mu0(z, xy, 2);
2322 dmu = grad_mu0(z, xy, 2);
2323 E_T(p, nu012(z, xy, 1), grad_nu012(z, xy, 1), E_T_ijk, dE_T_ijk);
2324 for (int j=1; j<p; j++)
2325 for (int i=0; i+j<p; i++, o++)
2326 {
2327 E(0) = E_T_ijk(i, j, 0);
2328 E(1) = E_T_ijk(i, j, 1);
2329 E(2) = E_T_ijk(i, j, 2);
2330 dmu.cross3D(E, dmuxE);
2331 for (int k=0; k<3; k++)
2332 {
2333 dW(o, k) = mu * dE_T_ijk(i, j, k) + dmuxE(k);
2334 }
2335 }
2336
2337 // (a, b) = (1, 2), c = 1
2338 mu = mu1(z, xy, 2);
2339 dmu = grad_mu1(z, xy, 2);
2340 for (int j=1; j<p; j++)
2341 for (int i=0; i+j<p; i++, o++)
2342 {
2343 E(0) = E_T_ijk(i, j, 0);
2344 E(1) = E_T_ijk(i, j, 1);
2345 E(2) = E_T_ijk(i, j, 2);
2346 dmu.cross3D(E, dmuxE);
2347 for (int k=0; k<3; k++)
2348 {
2349 dW(o, k) = mu * dE_T_ijk(i, j, k) + dmuxE(k);
2350 }
2351 }
2352
2353 // (a, b) = (2, 1), c = 0
2354 mu = mu0(z, xy, 1);
2355 dmu = grad_mu0(z, xy, 1);
2356 E_T(p, nu012(z, xy, 2), grad_nu012(z, xy, 2), E_T_ijk, dE_T_ijk);
2357 for (int j=1; j<p; j++)
2358 for (int i=0; i+j<p; i++, o++)
2359 {
2360 E(0) = E_T_ijk(i, j, 0);
2361 E(1) = E_T_ijk(i, j, 1);
2362 E(2) = E_T_ijk(i, j, 2);
2363 dmu.cross3D(E, dmuxE);
2364 for (int k=0; k<3; k++)
2365 {
2366 dW(o, k) = mu * dE_T_ijk(i, j, k) + dmuxE(k);
2367 }
2368 }
2369
2370 // (a, b) = (2, 1), c = 1
2371 mu = mu1(z, xy, 1);
2372 dmu = grad_mu1(z, xy, 1);
2373 for (int j=1; j<p; j++)
2374 for (int i=0; i+j<p; i++, o++)
2375 {
2376 E(0) = E_T_ijk(i, j, 0);
2377 E(1) = E_T_ijk(i, j, 1);
2378 E(2) = E_T_ijk(i, j, 2);
2379 dmu.cross3D(E, dmuxE);
2380 for (int k=0; k<3; k++)
2381 {
2382 dW(o, k) = mu * dE_T_ijk(i, j, k) + dmuxE(k);
2383 }
2384 }
2385
2386 // Family II
2387 // (a, b) = (1, 2), c = 0
2388 mu = mu0(z, xy, 2);
2389 dmu = grad_mu0(z, xy, 2);
2390 E_T(p, nu120(z, xy, 1), grad_nu120(z, xy, 1), E_T_ijk, dE_T_ijk);
2391 for (int j=1; j<p; j++)
2392 for (int i=0; i+j<p; i++, o++)
2393 {
2394 E(0) = E_T_ijk(i, j, 0);
2395 E(1) = E_T_ijk(i, j, 1);
2396 E(2) = E_T_ijk(i, j, 2);
2397 dmu.cross3D(E, dmuxE);
2398 for (int k=0; k<3; k++)
2399 {
2400 dW(o, k) = mu * dE_T_ijk(i, j, k) + dmuxE(k);
2401 }
2402 }
2403
2404 // (a, b) = (1, 2), c = 1
2405 mu = mu1(z, xy, 2);
2406 dmu = grad_mu1(z, xy, 2);
2407 for (int j=1; j<p; j++)
2408 for (int i=0; i+j<p; i++, o++)
2409 {
2410 E(0) = E_T_ijk(i, j, 0);
2411 E(1) = E_T_ijk(i, j, 1);
2412 E(2) = E_T_ijk(i, j, 2);
2413 dmu.cross3D(E, dmuxE);
2414 for (int k=0; k<3; k++)
2415 {
2416 dW(o, k) = mu * dE_T_ijk(i, j, k) + dmuxE(k);
2417 }
2418 }
2419
2420 // (a, b) = (2, 1), c = 0
2421 mu = mu0(z, xy, 1);
2422 dmu = grad_mu0(z, xy, 1);
2423 E_T(p, nu120(z, xy, 2), grad_nu120(z, xy, 2), E_T_ijk, dE_T_ijk);
2424 for (int j=1; j<p; j++)
2425 for (int i=0; i+j<p; i++, o++)
2426 {
2427 E(0) = E_T_ijk(i, j, 0);
2428 E(1) = E_T_ijk(i, j, 1);
2429 E(2) = E_T_ijk(i, j, 2);
2430 dmu.cross3D(E, dmuxE);
2431 for (int k=0; k<3; k++)
2432 {
2433 dW(o, k) = mu * dE_T_ijk(i, j, k) + dmuxE(k);
2434 }
2435 }
2436
2437 // (a, b) = (2, 1), c = 1
2438 mu = mu1(z, xy, 1);
2439 dmu = grad_mu1(z, xy, 1);
2440 for (int j=1; j<p; j++)
2441 for (int i=0; i+j<p; i++, o++)
2442 {
2443 E(0) = E_T_ijk(i, j, 0);
2444 E(1) = E_T_ijk(i, j, 1);
2445 E(2) = E_T_ijk(i, j, 2);
2446 dmu.cross3D(E, dmuxE);
2447 for (int k=0; k<3; k++)
2448 {
2449 dW(o, k) = mu * dE_T_ijk(i, j, k) + dmuxE(k);
2450 }
2451 }
2452 }
2453
2454 // Interior
2455 if (z < 1.0 && p >= 2)
2456 {
2457 // Family I
2458 // Curl is zero so skip these functions
2459 o += (p - 1) * (p - 1) * (p - 1);
2460
2461 // Family II
2462 mu = mu0(z);
2463 dmu = grad_mu0(z);
2464 phi_E(p, mu01(z), grad_mu01(z), phi_E_k, dphi_E_k);
2465 for (int k=2; k<=p; k++)
2466 {
2467 dphi(0) = dphi_E_k(k, 0);
2468 dphi(1) = dphi_E_k(k, 1);
2469 dphi(2) = dphi_E_k(k, 2);
2470 add(mu, dphi, phi_E_k(k), dmu, muphi);
2471
2472 for (int j=2; j<=p; j++)
2473 for (int i=0; i<p; i++, o++)
2474 {
2475 E(0) = E_Q1_ijk(i, j, 0);
2476 E(1) = E_Q1_ijk(i, j, 1);
2477 E(2) = E_Q1_ijk(i, j, 2);
2478 muphi.cross3D(E, dmuxE);
2479 for (int l=0; l<3; l++)
2480 {
2481 dW(o, l) = mu * dE_Q1_ijk(i, j, l) * phi_E_k(k) + dmuxE(l);
2482 }
2483 }
2484 }
2485
2486 // Family III
2487 for (int k=2; k<=p; k++)
2488 {
2489 dphi(0) = dphi_E_k(k, 0);
2490 dphi(1) = dphi_E_k(k, 1);
2491 dphi(2) = dphi_E_k(k, 2);
2492 add(mu, dphi, phi_E_k(k), dmu, muphi);
2493
2494 for (int j=2; j<=p; j++)
2495 for (int i=0; i<p; i++, o++)
2496 {
2497 E(0) = E_Q2_ijk(i, j, 0);
2498 E(1) = E_Q2_ijk(i, j, 1);
2499 E(2) = E_Q2_ijk(i, j, 2);
2500 muphi.cross3D(E, dmuxE);
2501 for (int l=0; l<3; l++)
2502 {
2503 dW(o, l) = mu * dE_Q2_ijk(i, j, l) * phi_E_k(k) + dmuxE(l);
2504 }
2505 }
2506 }
2507
2508 // Family IV
2509 // Re-using mu from Family II
2510 dmu = grad_mu0(z);
2511 phi_Q(p, mu01(z, xy, 2), grad_mu01(z, xy, 2), mu01(z, xy, 1),
2512 grad_mu01(z, xy, 1), phi_Q2_ij, dphi_Q2_ij);
2513 for (int j=2; j<=p; j++)
2514 for (int i=2; i<=p; i++, o++)
2515 {
2516 const int n = std::max(i,j);
2517 const real_t nmu = n * pow(mu, n-1);
2518
2519 dphi(0) = dphi_Q2_ij(i, j, 0);
2520 dphi(1) = dphi_Q2_ij(i, j, 1);
2521 dphi(2) = dphi_Q2_ij(i, j, 2);
2522 dphi.cross3D(dmu, muphi);
2523
2524 for (int l=0; l<3; l++)
2525 {
2526 dW(o, l) = nmu * muphi(l);
2527 }
2528 }
2529 }
2530}
2531
2533 : VectorFiniteElement(1, Geometry::POINT, 2, p,
2534 H_CURL, FunctionSpace::Pk)
2535{
2536 // VectorFiniteElement::SetDerivMembers doesn't support 0D H_CURL elements
2537 // so we mimic a 1D element and then correct the dimension here.
2538 dim = 0;
2539 vdim = 2;
2540 cdim = 0;
2541}
2542
2544 DenseMatrix &shape) const
2545{
2546 shape(0,0) = 1.0;
2547 shape(0,1) = 0.0;
2548
2549 shape(1,0) = 0.0;
2550 shape(1,1) = 1.0;
2551}
2552
2554 DenseMatrix &shape) const
2555{
2556 CalcVShape(Trans.GetIntPoint(), shape);
2557}
2558
2559const real_t ND_R1D_SegmentElement::tk[9] = { 1.,0.,0., 0.,1.,0., 0.,0.,1. };
2560
2562 const int cb_type,
2563 const int ob_type)
2564 : VectorFiniteElement(1, Geometry::SEGMENT, 3 * p + 2, p,
2565 H_CURL, FunctionSpace::Pk),
2566 dof2tk(dof),
2567 cbasis1d(poly1d.GetBasis(p, VerifyClosed(cb_type))),
2568 obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(ob_type)))
2569{
2570 // Override default types for VectorFiniteElements
2571 deriv_type = CURL;
2574
2575 // Override default dimensions for VectorFiniteElements
2576 vdim = 3;
2577 cdim = 3;
2578
2579 const real_t *cp = poly1d.ClosedPoints(p, cb_type);
2580 const real_t *op = poly1d.OpenPoints(p - 1, ob_type);
2581
2582#ifndef MFEM_THREAD_SAFE
2583 shape_cx.SetSize(p + 1);
2584 shape_ox.SetSize(p);
2585 dshape_cx.SetSize(p + 1);
2586#endif
2587
2588 dof_map.SetSize(dof);
2589
2590 int o = 0;
2591 // nodes
2592 // (0)
2593 Nodes.IntPoint(o).x = cp[0]; // y-directed
2594 dof_map[p] = o; dof2tk[o++] = 1;
2595 Nodes.IntPoint(o).x = cp[0]; // z-directed
2596 dof_map[2*p+1] = o; dof2tk[o++] = 2;
2597
2598 // (1)
2599 Nodes.IntPoint(o).x = cp[p]; // y-directed
2600 dof_map[2*p] = o; dof2tk[o++] = 1;
2601 Nodes.IntPoint(o).x = cp[p]; // z-directed
2602 dof_map[3*p+1] = o; dof2tk[o++] = 2;
2603
2604 // interior
2605 // x-components
2606 for (int i = 0; i < p; i++)
2607 {
2608 Nodes.IntPoint(o).x = op[i];
2609 dof_map[i] = o; dof2tk[o++] = 0;
2610 }
2611 // y-components
2612 for (int i = 1; i < p; i++)
2613 {
2614 Nodes.IntPoint(o).x = cp[i];
2615 dof_map[p+i] = o; dof2tk[o++] = 1;
2616 }
2617 // z-components
2618 for (int i = 1; i < p; i++)
2619 {
2620 Nodes.IntPoint(o).x = cp[i];
2621 dof_map[2*p+1+i] = o; dof2tk[o++] = 2;
2622 }
2623}
2624
2626 DenseMatrix &shape) const
2627{
2628 const int p = order;
2629
2630#ifdef MFEM_THREAD_SAFE
2631 Vector shape_cx(p + 1), shape_ox(p);
2632#endif
2633
2634 cbasis1d.Eval(ip.x, shape_cx);
2635 obasis1d.Eval(ip.x, shape_ox);
2636
2637 int o = 0;
2638 // x-components
2639 for (int i = 0; i < p; i++)
2640 {
2641 int idx = dof_map[o++];
2642 shape(idx,0) = shape_ox(i);
2643 shape(idx,1) = 0.;
2644 shape(idx,2) = 0.;
2645 }
2646 // y-components
2647 for (int i = 0; i <= p; i++)
2648 {
2649 int idx = dof_map[o++];
2650 shape(idx,0) = 0.;
2651 shape(idx,1) = shape_cx(i);
2652 shape(idx,2) = 0.;
2653 }
2654 // z-components
2655 for (int i = 0; i <= p; i++)
2656 {
2657 int idx = dof_map[o++];
2658 shape(idx,0) = 0.;
2659 shape(idx,1) = 0.;
2660 shape(idx,2) = shape_cx(i);
2661 }
2662}
2663
2665 DenseMatrix &shape) const
2666{
2667 CalcVShape(Trans.GetIntPoint(), shape);
2668 const DenseMatrix & JI = Trans.InverseJacobian();
2669 MFEM_ASSERT(JI.Width() == 1 && JI.Height() == 1,
2670 "ND_R1D_SegmentElement cannot be embedded in "
2671 "2 or 3 dimensional spaces");
2672 for (int i=0; i<dof; i++)
2673 {
2674 shape(i, 0) *= JI(0,0);
2675 }
2676}
2677
2679 DenseMatrix &curl_shape) const
2680{
2681 const int p = order;
2682
2683#ifdef MFEM_THREAD_SAFE
2684 Vector shape_cx(p + 1), shape_ox(p);
2685 Vector dshape_cx(p + 1);
2686#endif
2687
2688 cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
2689 obasis1d.Eval(ip.x, shape_ox);
2690
2691 int o = 0;
2692 // x-components
2693 for (int i = 0; i < p; i++)
2694 {
2695 int idx = dof_map[o++];
2696 curl_shape(idx,0) = 0.;
2697 curl_shape(idx,1) = 0.;
2698 curl_shape(idx,2) = 0.;
2699 }
2700 // y-components
2701 for (int i = 0; i <= p; i++)
2702 {
2703 int idx = dof_map[o++];
2704 curl_shape(idx,0) = 0.;
2705 curl_shape(idx,1) = 0.;
2706 curl_shape(idx,2) = dshape_cx(i);
2707 }
2708 // z-components
2709 for (int i = 0; i <= p; i++)
2710 {
2711 int idx = dof_map[o++];
2712 curl_shape(idx,0) = 0.;
2713 curl_shape(idx,1) = -dshape_cx(i);
2714 curl_shape(idx,2) = 0.;
2715 }
2716}
2717
2719 DenseMatrix &curl_shape) const
2720{
2721 CalcCurlShape(Trans.GetIntPoint(), curl_shape);
2722 const DenseMatrix & J = Trans.Jacobian();
2723 MFEM_ASSERT(J.Width() == 1 && J.Height() == 1,
2724 "ND_R1D_SegmentElement cannot be embedded in "
2725 "2 or 3 dimensional spaces");
2726 curl_shape *= (1.0 / J.Weight());
2727}
2728
2730 ElementTransformation &Trans,
2731 Vector &dofs) const
2732{
2733 real_t data[3];
2734 Vector vk(data, 3);
2735
2736 for (int k = 0; k < dof; k++)
2737 {
2738 Trans.SetIntPoint(&Nodes.IntPoint(k));
2739
2740 vc.Eval(vk, Trans, Nodes.IntPoint(k));
2741 // dof_k = vk^t J tk
2742 Vector t(const_cast<real_t*>(&tk[dof2tk[k] * 3]), 3);
2743 dofs(k) = Trans.Jacobian()(0,0) * t(0) * vk(0) +
2744 t(1) * vk(1) + t(2) * vk(2);
2745 }
2746
2747}
2748
2750 ElementTransformation &Trans,
2751 DenseMatrix &I) const
2752{
2753 if (fe.GetRangeType() == SCALAR)
2754 {
2756 Vector shape(fe.GetDof());
2757
2758 real_t * tk_ptr = const_cast<real_t*>(tk);
2759
2760 I.SetSize(dof, vdim*fe.GetDof());
2761 for (int k = 0; k < dof; k++)
2762 {
2763 const IntegrationPoint &ip = Nodes.IntPoint(k);
2764
2765 Vector t1(&tk_ptr[dof2tk[k] * 3], 1);
2766 Vector t3(&tk_ptr[dof2tk[k] * 3], 3);
2767
2768 fe.CalcShape(ip, shape);
2769 Trans.SetIntPoint(&ip);
2770 // Transform ND edge tengents from reference to physical space
2771 // vk = J tk
2772 Trans.Jacobian().Mult(t1, vk);
2773 vk[1] = t3[1];
2774 vk[2] = t3[2];
2775 if (fe.GetMapType() == INTEGRAL)
2776 {
2777 real_t w = 1.0/Trans.Weight();
2778 for (int d = 0; d < vdim; d++)
2779 {
2780 vk[d] *= w;
2781 }
2782 }
2783
2784 for (int j = 0; j < shape.Size(); j++)
2785 {
2786 real_t s = shape(j);
2787 if (fabs(s) < 1e-12)
2788 {
2789 s = 0.0;
2790 }
2791 // Project scalar basis function multiplied by each coordinate
2792 // direction onto the transformed edge tangents
2793 for (int d = 0; d < vdim; d++)
2794 {
2795 I(k, j + d*shape.Size()) = s*vk[d];
2796 }
2797 }
2798 }
2799 }
2800 else
2801 {
2804
2805 real_t * tk_ptr = const_cast<real_t*>(tk);
2806
2807 I.SetSize(dof, fe.GetDof());
2808 for (int k = 0; k < dof; k++)
2809 {
2810 const IntegrationPoint &ip = Nodes.IntPoint(k);
2811
2812 Vector t1(&tk_ptr[dof2tk[k] * 3], 1);
2813 Vector t3(&tk_ptr[dof2tk[k] * 3], 3);
2814
2815 Trans.SetIntPoint(&ip);
2816 // Transform ND edge tangents from reference to physical space
2817 // vk = J tk
2818 Trans.Jacobian().Mult(t1, vk);
2819 // Compute fe basis functions in physical space
2820 fe.CalcVShape(Trans, vshape);
2821 // Project fe basis functions onto transformed edge tangents
2822 for (int j=0; j<vshape.Height(); j++)
2823 {
2824 I(k, j) = 0.0;
2825 I(k, j) += vshape(j, 0) * vk[0];
2826 if (vshape.Width() == 3)
2827 {
2828 I(k, j) += vshape(j, 1) * t3(1);
2829 I(k, j) += vshape(j, 2) * t3(2);
2830 }
2831 }
2832 }
2833 }
2834}
2835
2836const real_t ND_R2D_SegmentElement::tk[4] = { 1.,0., 0.,1. };
2837
2839 const int cb_type,
2840 const int ob_type)
2841 : VectorFiniteElement(1, Geometry::SEGMENT, 2 * p + 1, p,
2842 H_CURL, FunctionSpace::Pk),
2843 dof2tk(dof),
2844 cbasis1d(poly1d.GetBasis(p, VerifyClosed(cb_type))),
2845 obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(ob_type)))
2846{
2847 // Override default dimensions for VectorFiniteElements
2848 vdim = 2;
2849 cdim = 1;
2850
2851 const real_t *cp = poly1d.ClosedPoints(p, cb_type);
2852 const real_t *op = poly1d.OpenPoints(p - 1, ob_type);
2853
2854#ifndef MFEM_THREAD_SAFE
2855 shape_cx.SetSize(p + 1);
2856 shape_ox.SetSize(p);
2857 dshape_cx.SetSize(p + 1);
2858#endif
2859
2860 dof_map.SetSize(dof);
2861
2862 int o = 0;
2863 // nodes
2864 // (0)
2865 Nodes.IntPoint(o).x = cp[0]; // z-directed
2866 dof_map[p] = o; dof2tk[o++] = 1;
2867
2868 // (1)
2869 Nodes.IntPoint(o).x = cp[p]; // z-directed
2870 dof_map[2*p] = o; dof2tk[o++] = 1;
2871
2872 // interior
2873 // x-components
2874 for (int i = 0; i < p; i++)
2875 {
2876 Nodes.IntPoint(o).x = op[i];
2877 dof_map[i] = o; dof2tk[o++] = 0;
2878 }
2879 // z-components
2880 for (int i = 1; i < p; i++)
2881 {
2882 Nodes.IntPoint(o).x = cp[i];
2883 dof_map[p+i] = o; dof2tk[o++] = 1;
2884 }
2885}
2886
2888 DenseMatrix &shape) const
2889{
2890 const int p = order;
2891
2892#ifdef MFEM_THREAD_SAFE
2893 Vector shape_cx(p + 1), shape_ox(p);
2894#endif
2895
2896 cbasis1d.Eval(ip.x, shape_cx);
2897 obasis1d.Eval(ip.x, shape_ox);
2898
2899 int o = 0;
2900 // x-components
2901 for (int i = 0; i < p; i++)
2902 {
2903 int idx = dof_map[o++];
2904 shape(idx,0) = shape_ox(i);
2905 shape(idx,1) = 0.;
2906 }
2907 // z-components
2908 for (int i = 0; i <= p; i++)
2909 {
2910 int idx = dof_map[o++];
2911 shape(idx,0) = 0.;
2912 shape(idx,1) = shape_cx(i);
2913 }
2914}
2915
2917 DenseMatrix &shape) const
2918{
2919 CalcVShape(Trans.GetIntPoint(), shape);
2920 const DenseMatrix & JI = Trans.InverseJacobian();
2921 MFEM_ASSERT(JI.Width() == 1 && JI.Height() == 1,
2922 "ND_R2D_SegmentElement cannot be embedded in "
2923 "2 or 3 dimensional spaces");
2924 for (int i=0; i<dof; i++)
2925 {
2926 shape(i, 0) *= JI(0,0);
2927 }
2928}
2929
2931 DenseMatrix &curl_shape) const
2932{
2933 const int p = order;
2934
2935#ifdef MFEM_THREAD_SAFE
2936 Vector shape_cx(p + 1), shape_ox(p);
2937 Vector dshape_cx(p + 1);
2938#endif
2939
2940 cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
2941 obasis1d.Eval(ip.x, shape_ox);
2942
2943 int o = 0;
2944 // x-components
2945 for (int i = 0; i < p; i++)
2946 {
2947 int idx = dof_map[o++];
2948 curl_shape(idx,0) = 0.;
2949 }
2950 // z-components
2951 for (int i = 0; i <= p; i++)
2952 {
2953 int idx = dof_map[o++];
2954 curl_shape(idx,0) = -dshape_cx(i);
2955 }
2956}
2957
2958void ND_R2D_SegmentElement::LocalInterpolation(const VectorFiniteElement &cfe,
2959 ElementTransformation &Trans,
2960 DenseMatrix &I) const
2961{
2962 real_t vk[Geometry::MaxDim]; vk[1] = 0.0; vk[2] = 0.0;
2963 Vector xk(vk, dim);
2965 DenseMatrix vshape(cfe.GetDof(), vdim);
2966
2967 real_t * tk_ptr = const_cast<real_t*>(tk);
2968
2969 I.SetSize(dof, vshape.Height());
2970
2971 // assuming Trans is linear; this should be ok for all refinement types
2973 const DenseMatrix &J = Trans.Jacobian();
2974 for (int k = 0; k < dof; k++)
2975 {
2976 Vector t1(&tk_ptr[dof2tk[k] * 2], 1);
2977 Vector t2(&tk_ptr[dof2tk[k] * 2], 2);
2978
2979 Trans.Transform(Nodes.IntPoint(k), xk);
2980 ip.Set3(vk);
2981 cfe.CalcVShape(ip, vshape);
2982 // xk = J t_k
2983 J.Mult(t1, vk);
2984 // I_k = vshape_k.J.t_k, k=1,...,Dof
2985 for (int j = 0; j < vshape.Height(); j++)
2986 {
2987 real_t Ikj = 0.;
2988 for (int i = 0; i < dim; i++)
2989 {
2990 Ikj += vshape(j, i) * vk[i];
2991 }
2992 Ikj += vshape(j, 1) * t2(1);
2993 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
2994 }
2995 }
2996}
2997
2999 ElementTransformation &Trans,
3000 Vector &dofs) const
3001{
3002 real_t data[3];
3003 Vector vk1(data, 1);
3004 Vector vk2(data, 2);
3005 Vector vk3(data, 3);
3006
3007 real_t * tk_ptr = const_cast<real_t*>(tk);
3008
3009 for (int k = 0; k < dof; k++)
3010 {
3011 Trans.SetIntPoint(&Nodes.IntPoint(k));
3012
3013 vc.Eval(vk3, Trans, Nodes.IntPoint(k));
3014 // dof_k = vk^t J tk
3015 Vector t1(&tk_ptr[dof2tk[k] * 2], 1);
3016 Vector t2(&tk_ptr[dof2tk[k] * 2], 2);
3017
3018 dofs(k) = Trans.Jacobian().InnerProduct(t1, vk2) + t2(1) * vk3(2);
3019 }
3020
3021}
3022
3024 const real_t *tk_fe)
3025 : VectorFiniteElement(2, G, Do, p,
3026 H_CURL, FunctionSpace::Pk),
3027 tk(tk_fe),
3028 dof_map(dof),
3029 dof2tk(dof)
3030{
3031 // Override default types for VectorFiniteElements
3032 deriv_type = CURL;
3035
3036 // Override default dimensions for VectorFiniteElements
3037 vdim = 3;
3038 cdim = 3;
3039}
3040
3042 DenseMatrix &shape) const
3043{
3044 CalcVShape(Trans.GetIntPoint(), shape);
3045 const DenseMatrix & JI = Trans.InverseJacobian();
3046 MFEM_ASSERT(JI.Width() == 2 && JI.Height() == 2,
3047 "ND_R2D_FiniteElement cannot be embedded in "
3048 "3 dimensional spaces");
3049 for (int i=0; i<dof; i++)
3050 {
3051 real_t sx = shape(i, 0);
3052 real_t sy = shape(i, 1);
3053 shape(i, 0) = sx * JI(0, 0) + sy * JI(1, 0);
3054 shape(i, 1) = sx * JI(0, 1) + sy * JI(1, 1);
3055 }
3056}
3057
3059 DenseMatrix &curl_shape) const
3060{
3061 CalcCurlShape(Trans.GetIntPoint(), curl_shape);
3062 const DenseMatrix & J = Trans.Jacobian();
3063 MFEM_ASSERT(J.Width() == 2 && J.Height() == 2,
3064 "ND_R2D_FiniteElement cannot be embedded in "
3065 "3 dimensional spaces");
3066 for (int i=0; i<dof; i++)
3067 {
3068 real_t sx = curl_shape(i, 0);
3069 real_t sy = curl_shape(i, 1);
3070 curl_shape(i, 0) = sx * J(0, 0) + sy * J(0, 1);
3071 curl_shape(i, 1) = sx * J(1, 0) + sy * J(1, 1);
3072 }
3073 curl_shape *= (1.0 / Trans.Weight());
3074}
3075
3076void ND_R2D_FiniteElement::LocalInterpolation(
3077 const VectorFiniteElement &cfe,
3078 ElementTransformation &Trans,
3079 DenseMatrix &I) const
3080{
3081 real_t vk[Geometry::MaxDim]; vk[2] = 0.0;
3082 Vector xk(vk, dim);
3084#ifdef MFEM_THREAD_SAFE
3085 DenseMatrix vshape(cfe.GetDof(), vdim);
3086#else
3087 vshape.SetSize(cfe.GetDof(), vdim);
3088#endif
3089
3090 real_t * tk_ptr = const_cast<real_t*>(tk);
3091
3092 I.SetSize(dof, vshape.Height());
3093
3094 // assuming Trans is linear; this should be ok for all refinement types
3096 const DenseMatrix &J = Trans.Jacobian();
3097 for (int k = 0; k < dof; k++)
3098 {
3099 Vector t2(&tk_ptr[dof2tk[k] * 3], 2);
3100 Vector t3(&tk_ptr[dof2tk[k] * 3], 3);
3101
3102 Trans.Transform(Nodes.IntPoint(k), xk);
3103 ip.Set3(vk);
3104 cfe.CalcVShape(ip, vshape);
3105 // xk = J t_k
3106 J.Mult(t2, vk);
3107 // I_k = vshape_k.J.t_k, k=1,...,Dof
3108 for (int j = 0; j < vshape.Height(); j++)
3109 {
3110 real_t Ikj = 0.;
3111 for (int i = 0; i < dim; i++)
3112 {
3113 Ikj += vshape(j, i) * vk[i];
3114 }
3115 Ikj += vshape(j, 2) * t3(2);
3116 I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
3117 }
3118 }
3119}
3120
3122 DenseMatrix &R) const
3123{
3124 real_t pt_data[Geometry::MaxDim];
3126 Vector pt(pt_data, dim);
3127
3128#ifdef MFEM_THREAD_SAFE
3130#endif
3131
3132 real_t * tk_ptr = const_cast<real_t*>(tk);
3133
3135 const DenseMatrix &Jinv = Trans.InverseJacobian();
3136 for (int j = 0; j < dof; j++)
3137 {
3138 Vector t2(&tk_ptr[dof2tk[j] * 3], 2);
3139 Vector t3(&tk_ptr[dof2tk[j] * 3], 3);
3140
3141 InvertLinearTrans(Trans, Nodes.IntPoint(j), pt);
3142 ip.Set(pt_data, dim);
3143 if (Geometries.CheckPoint(geom_type, ip)) // do we need an epsilon here?
3144 {
3145 CalcVShape(ip, vshape);
3146 Jinv.Mult(t2, pt_data);
3147 for (int k = 0; k < dof; k++)
3148 {
3149 real_t R_jk = 0.0;
3150 for (int d = 0; d < dim; d++)
3151 {
3152 R_jk += vshape(k,d)*pt_data[d];
3153 }
3154 R_jk += vshape(k, 2) * t3(2);
3155 R(j,k) = R_jk;
3156 }
3157 }
3158 else
3159 {
3160 // Set the whole row to avoid valgrind warnings in R.Threshold().
3161 R.SetRow(j, infinity());
3162 }
3163 }
3164 R.Threshold(1e-12);
3165}
3166
3168 ElementTransformation &Trans,
3169 Vector &dofs) const
3170{
3171 real_t data[3];
3172 Vector vk2(data, 2);
3173 Vector vk3(data, 3);
3174
3175 real_t * tk_ptr = const_cast<real_t*>(tk);
3176
3177 for (int k = 0; k < dof; k++)
3178 {
3179 Trans.SetIntPoint(&Nodes.IntPoint(k));
3180
3181 vc.Eval(vk3, Trans, Nodes.IntPoint(k));
3182 // dof_k = vk^t J tk
3183 Vector t2(&tk_ptr[dof2tk[k] * 3], 2);
3184 Vector t3(&tk_ptr[dof2tk[k] * 3], 3);
3185
3186 dofs(k) = Trans.Jacobian().InnerProduct(t2, vk2) + t3(2) * vk3(2);
3187 }
3188
3189}
3190
3192 ElementTransformation &Trans,
3193 DenseMatrix &I) const
3194{
3195 if (fe.GetRangeType() == SCALAR)
3196 {
3198 Vector shape(fe.GetDof());
3199
3200 real_t * tk_ptr = const_cast<real_t*>(tk);
3201
3202 I.SetSize(dof, vdim*fe.GetDof());
3203 for (int k = 0; k < dof; k++)
3204 {
3205 const IntegrationPoint &ip = Nodes.IntPoint(k);
3206
3207 Vector t2(&tk_ptr[dof2tk[k] * 3], 2);
3208 Vector t3(&tk_ptr[dof2tk[k] * 3], 3);
3209
3210 fe.CalcShape(ip, shape);
3211 Trans.SetIntPoint(&ip);
3212 // Transform ND edge tengents from reference to physical space
3213 // vk = J tk
3214 Trans.Jacobian().Mult(t2, vk);
3215 vk[2] = t3[2];
3216 if (fe.GetMapType() == INTEGRAL)
3217 {
3218 real_t w = 1.0/Trans.Weight();
3219 for (int d = 0; d < vdim; d++)
3220 {
3221 vk[d] *= w;
3222 }
3223 }
3224
3225 for (int j = 0; j < shape.Size(); j++)
3226 {
3227 real_t s = shape(j);
3228 if (fabs(s) < 1e-12)
3229 {
3230 s = 0.0;
3231 }
3232 // Project scalar basis function multiplied by each coordinate
3233 // direction onto the transformed edge tangents
3234 for (int d = 0; d < vdim; d++)
3235 {
3236 I(k, j + d*shape.Size()) = s*vk[d];
3237 }
3238 }
3239 }
3240 }
3241 else
3242 {
3245
3246 real_t * tk_ptr = const_cast<real_t*>(tk);
3247
3248 I.SetSize(dof, fe.GetDof());
3249 for (int k = 0; k < dof; k++)
3250 {
3251 const IntegrationPoint &ip = Nodes.IntPoint(k);
3252
3253 Vector t2(&tk_ptr[dof2tk[k] * 3], 2);
3254 Vector t3(&tk_ptr[dof2tk[k] * 3], 3);
3255
3256 Trans.SetIntPoint(&ip);
3257 // Transform ND edge tangents from reference to physical space
3258 // vk = J tk
3259 Trans.Jacobian().Mult(t2, vk);
3260 // Compute fe basis functions in physical space
3261 fe.CalcVShape(Trans, vshape);
3262 // Project fe basis functions onto transformed edge tangents
3263 for (int j=0; j<vshape.Height(); j++)
3264 {
3265 I(k, j) = 0.0;
3266 for (int i=0; i<2; i++)
3267 {
3268 I(k, j) += vshape(j, i) * vk[i];
3269 }
3270 if (vshape.Width() == 3)
3271 {
3272 I(k, j) += vshape(j, 2) * t3(2);
3273 }
3274 }
3275 }
3276 }
3277}
3278
3280 ElementTransformation &Trans,
3281 DenseMatrix &grad) const
3282{
3283 MFEM_ASSERT(fe.GetMapType() == VALUE, "");
3284
3285 DenseMatrix dshape(fe.GetDof(), fe.GetDim());
3286 Vector grad_k(fe.GetDof());
3287
3288 grad.SetSize(dof, fe.GetDof());
3289 for (int k = 0; k < dof; k++)
3290 {
3291 fe.CalcDShape(Nodes.IntPoint(k), dshape);
3292 dshape.Mult(tk + dof2tk[k]*vdim, grad_k);
3293 for (int j = 0; j < grad_k.Size(); j++)
3294 {
3295 grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
3296 }
3297 }
3298}
3299
3300const real_t ND_R2D_TriangleElement::tk_t[15] =
3301{ 1.,0.,0., -1.,1.,0., 0.,-1.,0., 0.,1.,0., 0.,0.,1. };
3302
3304 const int cb_type)
3305 : ND_R2D_FiniteElement(p, Geometry::TRIANGLE, ((3*p + 1)*(p + 2))/2, tk_t),
3306 ND_FE(p),
3307 H1_FE(p, cb_type)
3308{
3309 int pm1 = p - 1, pm2 = p - 2;
3310
3311#ifndef MFEM_THREAD_SAFE
3312 nd_shape.SetSize(ND_FE.GetDof(), 2);
3313 h1_shape.SetSize(H1_FE.GetDof());
3314 nd_dshape.SetSize(ND_FE.GetDof(), 1);
3315 h1_dshape.SetSize(H1_FE.GetDof(), 2);
3316#endif
3317
3318 int o = 0;
3319 int n = 0;
3320 int h = 0;
3321 // Three nodes
3322 dof_map[o] = -1 - h++; dof2tk[o++] = 4;
3323 dof_map[o] = -1 - h++; dof2tk[o++] = 4;
3324 dof_map[o] = -1 - h++; dof2tk[o++] = 4;
3325
3326 // Three edges
3327 for (int e=0; e<3; e++)
3328 {
3329 // Dofs in the plane
3330 for (int i=0; i<p; i++)
3331 {
3332 dof_map[o] = n++; dof2tk[o++] = e;
3333 }
3334 // z-directed dofs
3335 for (int i=0; i<pm1; i++)
3336 {
3337 dof_map[o] = -1 - h++; dof2tk[o++] = 4;
3338 }
3339 }
3340
3341 // Interior dofs in the plane
3342 for (int j = 0; j <= pm2; j++)
3343 for (int i = 0; i + j <= pm2; i++)
3344 {
3345 dof_map[o] = n++; dof2tk[o++] = 0;
3346 dof_map[o] = n++; dof2tk[o++] = 3;
3347 }
3348
3349 // Interior z-directed dofs
3350 for (int j = 0; j < pm1; j++)
3351 for (int i = 0; i + j < pm2; i++)
3352 {
3353 dof_map[o] = -1 - h++; dof2tk[o++] = 4;
3354 }
3355
3356 MFEM_VERIFY(n == ND_FE.GetDof(),
3357 "ND_R2D_Triangle incorrect number of ND dofs.");
3358 MFEM_VERIFY(h == H1_FE.GetDof(),
3359 "ND_R2D_Triangle incorrect number of H1 dofs.");
3360 MFEM_VERIFY(o == GetDof(),
3361 "ND_R2D_Triangle incorrect number of dofs.");
3362
3363 const IntegrationRule & nd_Nodes = ND_FE.GetNodes();
3364 const IntegrationRule & h1_Nodes = H1_FE.GetNodes();
3365
3366 for (int i=0; i<dof; i++)
3367 {
3368 int idx = dof_map[i];
3369 if (idx >= 0)
3370 {
3371 const IntegrationPoint & ip = nd_Nodes.IntPoint(idx);
3372 Nodes.IntPoint(i).Set2(ip.x, ip.y);
3373 }
3374 else
3375 {
3376 const IntegrationPoint & ip = h1_Nodes.IntPoint(-idx-1);
3377 Nodes.IntPoint(i).Set2(ip.x, ip.y);
3378 }
3379 }
3380}
3381
3383 DenseMatrix &shape) const
3384{
3385#ifdef MFEM_THREAD_SAFE
3386 DenseMatrix nd_shape(ND_FE.GetDof(), 2);
3387 Vector h1_shape(H1_FE.GetDof());
3388#endif
3389
3390 ND_FE.CalcVShape(ip, nd_shape);
3391 H1_FE.CalcShape(ip, h1_shape);
3392
3393 for (int i=0; i<dof; i++)
3394 {
3395 int idx = dof_map[i];
3396 if (idx >= 0)
3397 {
3398 shape(i, 0) = nd_shape(idx, 0);
3399 shape(i, 1) = nd_shape(idx, 1);
3400 shape(i, 2) = 0.0;
3401 }
3402 else
3403 {
3404 shape(i, 0) = 0.0;
3405 shape(i, 1) = 0.0;
3406 shape(i, 2) = h1_shape(-idx-1);
3407 }
3408 }
3409}
3410
3412 DenseMatrix &curl_shape) const
3413{
3414#ifdef MFEM_THREAD_SAFE
3415 DenseMatrix nd_dshape(ND_FE.GetDof(), 1);
3416 DenseMatrix h1_dshape(H1_FE.GetDof(), 2);
3417#endif
3418
3419 ND_FE.CalcCurlShape(ip, nd_dshape);
3420 H1_FE.CalcDShape(ip, h1_dshape);
3421
3422 for (int i=0; i<dof; i++)
3423 {
3424 int idx = dof_map[i];
3425 if (idx >= 0)
3426 {
3427 curl_shape(i, 0) = 0.0;
3428 curl_shape(i, 1) = 0.0;
3429 curl_shape(i, 2) = nd_dshape(idx, 0);
3430 }
3431 else
3432 {
3433 curl_shape(i, 0) = h1_dshape(-idx-1, 1);
3434 curl_shape(i, 1) = -h1_dshape(-idx-1, 0);
3435 curl_shape(i, 2) = 0.0;
3436 }
3437 }
3438}
3439
3440
3441const real_t ND_R2D_QuadrilateralElement::tk_q[15] =
3442{ 1.,0.,0., 0.,1.,0., -1.,0.,0., 0.,-1.,0., 0.,0.,1. };
3443
3445 const int cb_type,
3446 const int ob_type)
3447 : ND_R2D_FiniteElement(p, Geometry::SQUARE, ((3*p + 1)*(p + 1)), tk_q),
3448 cbasis1d(poly1d.GetBasis(p, VerifyClosed(cb_type))),
3449 obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(ob_type)))
3450{
3451 const real_t *cp = poly1d.ClosedPoints(p, cb_type);
3452 const real_t *op = poly1d.OpenPoints(p - 1, ob_type);
3453 const int dofx = p*(p+1);
3454 const int dofy = p*(p+1);
3455 const int dofxy = dofx+dofy;
3456
3457#ifndef MFEM_THREAD_SAFE
3458 shape_cx.SetSize(p + 1);
3459 shape_ox.SetSize(p);
3460 shape_cy.SetSize(p + 1);
3461 shape_oy.SetSize(p);
3462 dshape_cx.SetSize(p + 1);
3463 dshape_cy.SetSize(p + 1);
3464#endif
3465
3467
3468 int o = 0;
3469 // nodes
3470 dof_map[dofxy] = o++; // (0)
3471 dof_map[dofxy+p] = o++; // (1)
3472 dof_map[dof-1] = o++; // (2)
3473 dof_map[dof-p-1] = o++; // (3)
3474
3475 // edges
3476 for (int i = 0; i < p; i++) // (0,1) - x-directed
3477 {
3478 dof_map[i + 0*p] = o++;
3479 }
3480 for (int i = 1; i < p; i++) // (0,1) - z-directed
3481 {
3482 dof_map[dofxy + i + 0*(p+1)] = o++;
3483 }
3484 for (int j = 0; j < p; j++) // (1,2) - y-directed
3485 {
3486 dof_map[dofx + p + j*(p + 1)] = o++;
3487 }
3488 for (int j = 1; j < p; j++) // (1,2) - z-directed
3489 {
3490 dof_map[dofxy + p + j*(p + 1)] = o++;
3491 }
3492 for (int i = 0; i < p; i++) // (2,3) - x-directed
3493 {
3494 dof_map[(p - 1 - i) + p*p] = -1 - (o++);
3495 }
3496 for (int i = 1; i < p; i++) // (2,3) - z-directed
3497 {
3498 dof_map[dofxy + (p - i) + p*(p + 1)] = o++;
3499 }
3500 for (int j = 0; j < p; j++) // (3,0) - y-directed
3501 {
3502 dof_map[dofx + 0 + (p - 1 - j)*(p + 1)] = -1 - (o++);
3503 }
3504 for (int j = 1; j < p; j++) // (3,0) - z-directed
3505 {
3506 dof_map[dofxy + (p - j)*(p + 1)] = o++;
3507 }
3508
3509 // interior
3510 // x-components
3511 for (int j = 1; j < p; j++)
3512 for (int i = 0; i < p; i++)
3513 {
3514 dof_map[i + j*p] = o++;
3515 }
3516 // y-components
3517 for (int j = 0; j < p; j++)
3518 for (int i = 1; i < p; i++)
3519 {
3520 dof_map[dofx + i + j*(p + 1)] = o++;
3521 }
3522 // z-components
3523 for (int j = 1; j < p; j++)
3524 for (int i = 1; i < p; i++)
3525 {
3526 dof_map[dofxy + i + j*(p + 1)] = o++;
3527 }
3528
3529 // set dof2tk and Nodes
3530 o = 0;
3531 // x-components
3532 for (int j = 0; j <= p; j++)
3533 for (int i = 0; i < p; i++)
3534 {
3535 int idx;
3536 if ((idx = dof_map[o++]) < 0)
3537 {
3538 dof2tk[idx = -1 - idx] = 2;
3539 }
3540 else
3541 {
3542 dof2tk[idx] = 0;
3543 }
3544 Nodes.IntPoint(idx).Set2(op[i], cp[j]);
3545 }
3546 // y-components
3547 for (int j = 0; j < p; j++)
3548 for (int i = 0; i <= p; i++)
3549 {
3550 int idx;
3551 if ((idx = dof_map[o++]) < 0)
3552 {
3553 dof2tk[idx = -1 - idx] = 3;
3554 }
3555 else
3556 {
3557 dof2tk[idx] = 1;
3558 }
3559 Nodes.IntPoint(idx).Set2(cp[i], op[j]);
3560 }
3561 // z-components
3562 for (int j = 0; j <= p; j++)
3563 for (int i = 0; i <= p; i++)
3564 {
3565 int idx = dof_map[o++];
3566 dof2tk[idx] = 4;
3567 Nodes.IntPoint(idx).Set2(cp[i], cp[j]);
3568 }
3569}
3570
3572 DenseMatrix &shape) const
3573{
3574 const int p = order;
3575
3576#ifdef MFEM_THREAD_SAFE
3577 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
3578#endif
3579
3580 cbasis1d.Eval(ip.x, shape_cx);
3581 obasis1d.Eval(ip.x, shape_ox);
3582 cbasis1d.Eval(ip.y, shape_cy);
3583 obasis1d.Eval(ip.y, shape_oy);
3584
3585 int o = 0;
3586 // x-components
3587 for (int j = 0; j <= p; j++)
3588 for (int i = 0; i < p; i++)
3589 {
3590 int idx, s;
3591 if ((idx = dof_map[o++]) < 0)
3592 {
3593 idx = -1 - idx, s = -1;
3594 }
3595 else
3596 {
3597 s = +1;
3598 }
3599 shape(idx,0) = s*shape_ox(i)*shape_cy(j);
3600 shape(idx,1) = 0.;
3601 shape(idx,2) = 0.;
3602 }
3603 // y-components
3604 for (int j = 0; j < p; j++)
3605 for (int i = 0; i <= p; i++)
3606 {
3607 int idx, s;
3608 if ((idx = dof_map[o++]) < 0)
3609 {
3610 idx = -1 - idx, s = -1;
3611 }
3612 else
3613 {
3614 s = +1;
3615 }
3616 shape(idx,0) = 0.;
3617 shape(idx,1) = s*shape_cx(i)*shape_oy(j);
3618 shape(idx,2) = 0.;
3619 }
3620 // z-components
3621 for (int j = 0; j <= p; j++)
3622 for (int i = 0; i <= p; i++)
3623 {
3624 int idx = dof_map[o++];
3625 shape(idx,0) = 0.;
3626 shape(idx,1) = 0.;
3627 shape(idx,2) = shape_cx(i)*shape_cy(j);
3628 }
3629}
3630
3632 DenseMatrix &curl_shape) const
3633{
3634 const int p = order;
3635
3636#ifdef MFEM_THREAD_SAFE
3637 Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
3638 Vector dshape_cx(p + 1), dshape_cy(p + 1);
3639#endif
3640
3641 cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
3642 obasis1d.Eval(ip.x, shape_ox);
3643 cbasis1d.Eval(ip.y, shape_cy, dshape_cy);
3644 obasis1d.Eval(ip.y, shape_oy);
3645
3646 int o = 0;
3647 // x-components
3648 for (int j = 0; j <= p; j++)
3649 for (int i = 0; i < p; i++)
3650 {
3651 int idx, s;
3652 if ((idx = dof_map[o++]) < 0)
3653 {
3654 idx = -1 - idx, s = -1;
3655 }
3656 else
3657 {
3658 s = +1;
3659 }
3660 curl_shape(idx,0) = 0.;
3661 curl_shape(idx,1) = 0.;
3662 curl_shape(idx,2) = -s*shape_ox(i)*dshape_cy(j);
3663 }
3664 // y-components
3665 for (int j = 0; j < p; j++)
3666 for (int i = 0; i <= p; i++)
3667 {
3668 int idx, s;
3669 if ((idx = dof_map[o++]) < 0)
3670 {
3671 idx = -1 - idx, s = -1;
3672 }
3673 else
3674 {
3675 s = +1;
3676 }
3677 curl_shape(idx,0) = 0.;
3678 curl_shape(idx,1) = 0.;
3679 curl_shape(idx,2) = s*dshape_cx(i)*shape_oy(j);
3680 }
3681 // z-components
3682 for (int j = 0; j <= p; j++)
3683 for (int i = 0; i <= p; i++)
3684 {
3685 int idx = dof_map[o++];
3686 curl_shape(idx,0) = shape_cx(i)*dshape_cy(j);
3687 curl_shape(idx,1) = -dshape_cx(i)*shape_cy(j);
3688 curl_shape(idx,2) = 0.;
3689 }
3690}
3691
3692}
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 Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
Definition densemat.cpp:125
void Threshold(real_t eps)
Replace small entries, abs(a_ij) <= eps, with zero.
void SetRow(int r, const real_t *row)
real_t InnerProduct(const real_t *x, const real_t *y) const
Compute y^t A x.
Definition densemat.cpp:301
real_t * Data() const
Returns the matrix data array.
Definition densemat.hpp:114
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:108
real_t Weight() const
Definition densemat.cpp:573
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)
const DenseMatrix & InverseJacobian()
Return the inverse of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition eltrans.hpp:158
const IntegrationPoint & GetIntPoint()
Get a const reference to the currently set integration point. This will return NULL if no integration...
Definition eltrans.hpp:111
real_t Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint....
Definition eltrans.hpp:144
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint,...
Definition eltrans.hpp:132
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition eltrans.hpp:106
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
Abstract class for all finite elements.
Definition fe_base.hpp:244
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_base.cpp:40
int dof
Number of degrees of freedom.
Definition fe_base.hpp:253
int GetRangeDim() const
Returns the vector dimension for vector-valued finite elements, which is also the dimension of the in...
Definition fe_base.hpp:325
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition fe_base.cpp:65
IntegrationRule Nodes
Definition fe_base.hpp:256
int GetDim() const
Returns the reference space dimension for the finite element.
Definition fe_base.hpp:321
int vdim
Vector dimension of vector-valued basis functions.
Definition fe_base.hpp:247
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition fe_base.hpp:360
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
Definition fe_base.hpp:351
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition fe_base.hpp:400
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
int cdim
Dimension of curl for vector-valued basis functions.
Definition fe_base.hpp:248
@ CURL
Implements CalcCurlShape methods.
Definition fe_base.hpp:307
Geometry::Type geom_type
Geometry::Type of the reference element.
Definition fe_base.hpp:249
DenseMatrix vshape
Definition fe_base.hpp:258
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
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)
static Vector lam45_grad_lam45(real_t x, real_t y, real_t z)
void phi_Q(int p, Vector s, Vector t, DenseMatrix &u) const
static constexpr real_t apex_tol
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 DenseMatrix grad_nu120(real_t z, Vector xy, unsigned int ab)
static DenseMatrix grad_lam25(real_t x, real_t y, real_t z)
void E_E(int p, Vector s, Vector sds, DenseMatrix &u) const
static Vector nu01_grad_nu01(real_t z, Vector xy, unsigned int ab)
static Vector lam35(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 lam25_grad_lam25(real_t x, real_t y, real_t z)
static Vector mu01(real_t z)
static DenseMatrix grad_nu012(real_t z, Vector xy, unsigned int ab)
static Vector grad_mu1(real_t z)
void E_T(int p, Vector s, Vector sds, DenseTensor &u) const
static real_t mu0(real_t z)
static Vector nu120(real_t z, Vector xy, unsigned int ab)
static Vector lam15_grad_lam15(real_t x, real_t y, real_t z)
Computes .
static Vector lam35_grad_lam35(real_t x, real_t y, real_t z)
static Vector grad_mu0(real_t z)
static DenseMatrix grad_lam35(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.
void E_Q(int p, Vector s, Vector ds, Vector t, DenseTensor &u) const
static Vector lam15(real_t x, real_t y, real_t z)
Two component vectors associated with edges touching the apex.
static Vector mu01_grad_mu01(real_t z, Vector xy, unsigned int ab)
static real_t mu1(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 Vector nu12_grad_nu12(real_t z, Vector xy, unsigned int ab)
Describes the function space on each element.
Definition fe_base.hpp:226
static const int MaxDim
Definition geom.hpp:47
const IntegrationPoint & GetCenter(int GeomType) const
Return the center of the given Geometry::Type, GeomType.
Definition geom.hpp:75
static bool CheckPoint(int GeomType, const IntegrationPoint &ip)
Check if the given point is inside the given reference element.
Definition geom.cpp:435
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
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 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 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
Class for integration point with weight.
Definition intrules.hpp:35
void Set(const real_t *p, const int dim)
Definition intrules.hpp:46
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
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition intrules.hpp:256
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition intrules.hpp:259
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
void CalcRawVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Definition fe_nd.cpp:1840
void CalcRawCurlShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Definition fe_nd.cpp:1862
ND_FuentesPyramidElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Definition fe_nd.cpp:1591
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_nd.cpp:1786
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition fe_nd.cpp:1811
ND_HexahedronElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Construct the ND_HexahedronElement of order p and closed and open BasisType cb_type and ob_type.
Definition fe_nd.cpp:26
void GetFaceMap(const int face_id, Array< int > &face_map) const override
Return the mapping from lexicographic face DOFs to lexicographic element DOFs for the given local fac...
Definition fe_nd.cpp:485
void ProjectIntegrated(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Definition fe_nd.cpp:245
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_nd.cpp:310
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition fe_nd.cpp:400
ND_QuadrilateralElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Construct the ND_QuadrilateralElement of order p and closed and open BasisType cb_type and ob_type.
Definition fe_nd.cpp:533
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_nd.cpp:700
void GetFaceMap(const int face_id, Array< int > &face_map) const override
Return the mapping from lexicographic face DOFs to lexicographic element DOFs for the given local fac...
Definition fe_nd.cpp:820
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition fe_nd.cpp:763
void ProjectIntegrated(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Definition fe_nd.cpp:624
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_nd.cpp:2543
ND_R1D_PointElement(int p)
Construct the ND_R1D_PointElement.
Definition fe_nd.cpp:2532
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_nd.cpp:2625
ND_R1D_SegmentElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Construct the ND_R1D_SegmentElement of order p and closed and open BasisType cb_type and ob_type.
Definition fe_nd.cpp:2561
void CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in physical space at the point de...
Definition fe_nd.cpp:2718
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition fe_nd.cpp:2678
void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const override
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
Definition fe_nd.cpp:2729
ND_R2D_FiniteElement(int p, Geometry::Type G, int Do, const real_t *tk_fe)
Definition fe_nd.cpp:3023
void CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in physical space at the point de...
Definition fe_nd.cpp:3058
void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const override
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
Definition fe_nd.cpp:3167
void CalcVShape(ElementTransformation &Trans, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in physical space at the point ...
Definition fe_nd.cpp:3041
void ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const override
Compute the discrete gradient matrix from the given FiniteElement onto 'this' FiniteElement....
Definition fe_nd.cpp:3279
void GetLocalRestriction(ElementTransformation &Trans, DenseMatrix &R) const override
Return a local restriction matrix R (Dof x Dof) mapping fine dofs to coarse dofs.
Definition fe_nd.cpp:3121
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition fe_nd.cpp:3631
ND_R2D_QuadrilateralElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Construct the ND_R2D_QuadrilateralElement of order p and closed and open BasisType cb_type and ob_typ...
Definition fe_nd.cpp:3444
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_nd.cpp:3571
void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const override
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
Definition fe_nd.cpp:2998
ND_R2D_SegmentElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Construct the ND_R2D_SegmentElement of order p and closed and open BasisType cb_type and ob_type.
Definition fe_nd.cpp:2838
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_nd.cpp:2887
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition fe_nd.cpp:2930
ND_R2D_TriangleElement(const int p, const int cb_type=BasisType::GaussLobatto)
Construct the ND_R2D_TriangleElement of order p.
Definition fe_nd.cpp:3303
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition fe_nd.cpp:3411
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_nd.cpp:3382
ND_SegmentElement(const int p, const int ob_type=BasisType::GaussLegendre)
Construct the ND_SegmentElement of order p and open BasisType ob_type.
Definition fe_nd.cpp:1268
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_nd.cpp:1285
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition fe_nd.cpp:1036
ND_TetrahedronElement(const int p)
Construct the ND_TetrahedronElement of order p.
Definition fe_nd.cpp:846
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_nd.cpp:994
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_nd.cpp:1189
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition fe_nd.cpp:1223
ND_TriangleElement(const int p)
Construct the ND_TriangleElement of order p.
Definition fe_nd.cpp:1108
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_nd.cpp:1512
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition fe_nd.cpp:1546
ND_WedgeElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Definition fe_nd.cpp:1296
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
void Eval(const real_t x, Vector &u) const
Evaluate the basis functions at point x in [0,1].
Definition fe_base.cpp:1777
bool IsIntegratedType() const
Returns true if the basis is "integrated", false otherwise.
Definition fe_base.hpp:1047
void ScaleIntegrated(bool scale_integrated_)
Set whether the "integrated" basis should be scaled by the subcell sizes. Has no effect for non-integ...
Definition fe_base.cpp:2034
void EvalIntegrated(const Vector &d, Vector &i) const
Evaluate the "integrated" basis type using pre-computed closed basis derivatives.
Definition fe_base.cpp:2009
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
const real_t * OpenPoints(const int p, const int btype=BasisType::GaussLegendre, bool on_device=false)
Get coordinates of an open (GaussLegendre) set of points if degree p.
Definition fe_base.hpp:1113
static void CalcBasis(const int p, const real_t x, real_t *u)
Evaluate the values of a hierarchical 1D basis at point x hierarchical = k-th basis function is degre...
Definition fe_base.hpp:1140
Poly_1D::Basis & basis1d
Definition fe_base.hpp:1251
Base class for vector Coefficients that optionally depend on time and space.
int GetVDim()
Returns dimension of the vector.
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
Intermediate class for finite elements whose basis functions return vector values.
Definition fe_base.hpp:813
Vector data type.
Definition vector.hpp:82
int Size() const
Returns the size of the vector.
Definition vector.hpp:226
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:558
real_t mu
Definition ex25.cpp:140
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition vector.hpp:47
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
Geometry Geometries
Definition fe.cpp:49
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition vector.cpp:391
float real_t
Definition config.hpp:43
Poly_1D poly1d
Definition fe.cpp:28
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
void InvertLinearTrans(ElementTransformation &trans, const IntegrationPoint &pt, Vector &x)
Definition fe_base.cpp:748
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition intrules.hpp:492
real_t p(const Vector &x, real_t t)