MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
fe_h1.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12// H1 Finite Element classes
13
14#include "fe_h1.hpp"
15
16namespace mfem
17{
18
19using namespace std;
20
21H1_SegmentElement::H1_SegmentElement(const int p, const int btype)
22 : NodalTensorFiniteElement(1, p, VerifyClosed(btype), H1_DOF_MAP)
23{
24 const real_t *cp = poly1d.ClosedPoints(p, b_type);
25
26#ifndef MFEM_THREAD_SAFE
27 shape_x.SetSize(p+1);
28 dshape_x.SetSize(p+1);
29 d2shape_x.SetSize(p+1);
30#endif
31
32 Nodes.IntPoint(0).x = cp[0];
33 Nodes.IntPoint(1).x = cp[p];
34 for (int i = 1; i < p; i++)
35 {
36 Nodes.IntPoint(i+1).x = cp[i];
37 }
38}
39
40void H1_SegmentElement::CalcShape(const IntegrationPoint &ip,
41 Vector &shape) const
42{
43 const int p = order;
44
45#ifdef MFEM_THREAD_SAFE
46 Vector shape_x(p+1);
47#endif
48
49 basis1d.Eval(ip.x, shape_x);
50
51 shape(0) = shape_x(0);
52 shape(1) = shape_x(p);
53 for (int i = 1; i < p; i++)
54 {
55 shape(i+1) = shape_x(i);
56 }
57}
58
59void H1_SegmentElement::CalcDShape(const IntegrationPoint &ip,
60 DenseMatrix &dshape) const
61{
62 const int p = order;
63
64#ifdef MFEM_THREAD_SAFE
65 Vector shape_x(p+1), dshape_x(p+1);
66#endif
67
68 basis1d.Eval(ip.x, shape_x, dshape_x);
69
70 dshape(0,0) = dshape_x(0);
71 dshape(1,0) = dshape_x(p);
72 for (int i = 1; i < p; i++)
73 {
74 dshape(i+1,0) = dshape_x(i);
75 }
76}
77
78void H1_SegmentElement::CalcHessian(const IntegrationPoint &ip,
79 DenseMatrix &Hessian) const
80{
81 const int p = order;
82
83#ifdef MFEM_THREAD_SAFE
84 Vector shape_x(p+1), dshape_x(p+1), d2shape_x(p+1);
85#endif
86
87 basis1d.Eval(ip.x, shape_x, dshape_x, d2shape_x);
88
89 Hessian(0,0) = d2shape_x(0);
90 Hessian(1,0) = d2shape_x(p);
91 for (int i = 1; i < p; i++)
92 {
93 Hessian(i+1,0) = d2shape_x(i);
94 }
95}
96
97void H1_SegmentElement::ProjectDelta(int vertex, Vector &dofs) const
98{
99 const int p = order;
100 const real_t *cp = poly1d.ClosedPoints(p, b_type);
101
102 switch (vertex)
103 {
104 case 0:
105 dofs(0) = poly1d.CalcDelta(p, (1.0 - cp[0]));
106 dofs(1) = poly1d.CalcDelta(p, (1.0 - cp[p]));
107 for (int i = 1; i < p; i++)
108 {
109 dofs(i+1) = poly1d.CalcDelta(p, (1.0 - cp[i]));
110 }
111 break;
112
113 case 1:
114 dofs(0) = poly1d.CalcDelta(p, cp[0]);
115 dofs(1) = poly1d.CalcDelta(p, cp[p]);
116 for (int i = 1; i < p; i++)
117 {
118 dofs(i+1) = poly1d.CalcDelta(p, cp[i]);
119 }
120 break;
121 }
122}
123
124
125H1_QuadrilateralElement::H1_QuadrilateralElement(const int p, const int btype)
126 : NodalTensorFiniteElement(2, p, VerifyClosed(btype), H1_DOF_MAP)
127{
128 const real_t *cp = poly1d.ClosedPoints(p, b_type);
129
130#ifndef MFEM_THREAD_SAFE
131 const int p1 = p + 1;
132
133 shape_x.SetSize(p1);
134 shape_y.SetSize(p1);
135 dshape_x.SetSize(p1);
136 dshape_y.SetSize(p1);
137 d2shape_x.SetSize(p1);
138 d2shape_y.SetSize(p1);
139#endif
140
141 int o = 0;
142 for (int j = 0; j <= p; j++)
143 {
144 for (int i = 0; i <= p; i++)
145 {
146 Nodes.IntPoint(dof_map[o++]).Set2(cp[i], cp[j]);
147 }
148 }
149}
150
151void H1_QuadrilateralElement::CalcShape(const IntegrationPoint &ip,
152 Vector &shape) const
153{
154 const int p = order;
155
156#ifdef MFEM_THREAD_SAFE
157 Vector shape_x(p+1), shape_y(p+1);
158#endif
159
160 basis1d.Eval(ip.x, shape_x);
161 basis1d.Eval(ip.y, shape_y);
162
163 for (int o = 0, j = 0; j <= p; j++)
164 for (int i = 0; i <= p; i++)
165 {
166 shape(dof_map[o++]) = shape_x(i)*shape_y(j);
167 }
168}
169
170void H1_QuadrilateralElement::CalcDShape(const IntegrationPoint &ip,
171 DenseMatrix &dshape) const
172{
173 const int p = order;
174
175#ifdef MFEM_THREAD_SAFE
176 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
177#endif
178
179 basis1d.Eval(ip.x, shape_x, dshape_x);
180 basis1d.Eval(ip.y, shape_y, dshape_y);
181
182 for (int o = 0, j = 0; j <= p; j++)
183 {
184 for (int i = 0; i <= p; i++)
185 {
186 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j);
187 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j); o++;
188 }
189 }
190}
191
192void H1_QuadrilateralElement::CalcHessian(const IntegrationPoint &ip,
193 DenseMatrix &Hessian) const
194{
195 const int p = order;
196
197#ifdef MFEM_THREAD_SAFE
198 Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1),
199 d2shape_x(p+1), d2shape_y(p+1);
200#endif
201
202 basis1d.Eval(ip.x, shape_x, dshape_x, d2shape_x);
203 basis1d.Eval(ip.y, shape_y, dshape_y, d2shape_y);
204
205 for (int o = 0, j = 0; j <= p; j++)
206 {
207 for (int i = 0; i <= p; i++)
208 {
209 Hessian(dof_map[o],0) = d2shape_x(i)* shape_y(j);
210 Hessian(dof_map[o],1) = dshape_x(i)* dshape_y(j);
211 Hessian(dof_map[o],2) = shape_x(i)*d2shape_y(j); o++;
212 }
213 }
214}
215
216void H1_QuadrilateralElement::ProjectDelta(int vertex, Vector &dofs) const
217{
218 const int p = order;
219 const real_t *cp = poly1d.ClosedPoints(p, b_type);
220
221#ifdef MFEM_THREAD_SAFE
222 Vector shape_x(p+1), shape_y(p+1);
223#endif
224
225 for (int i = 0; i <= p; i++)
226 {
227 shape_x(i) = poly1d.CalcDelta(p, (1.0 - cp[i]));
228 shape_y(i) = poly1d.CalcDelta(p, cp[i]);
229 }
230
231 switch (vertex)
232 {
233 case 0:
234 for (int o = 0, j = 0; j <= p; j++)
235 for (int i = 0; i <= p; i++)
236 {
237 dofs(dof_map[o++]) = shape_x(i)*shape_x(j);
238 }
239 break;
240 case 1:
241 for (int o = 0, j = 0; j <= p; j++)
242 for (int i = 0; i <= p; i++)
243 {
244 dofs(dof_map[o++]) = shape_y(i)*shape_x(j);
245 }
246 break;
247 case 2:
248 for (int o = 0, j = 0; j <= p; j++)
249 for (int i = 0; i <= p; i++)
250 {
251 dofs(dof_map[o++]) = shape_y(i)*shape_y(j);
252 }
253 break;
254 case 3:
255 for (int o = 0, j = 0; j <= p; j++)
256 for (int i = 0; i <= p; i++)
257 {
258 dofs(dof_map[o++]) = shape_x(i)*shape_y(j);
259 }
260 break;
261 }
262}
263
264
265H1_HexahedronElement::H1_HexahedronElement(const int p, const int btype)
266 : NodalTensorFiniteElement(3, p, VerifyClosed(btype), H1_DOF_MAP)
267{
268 const real_t *cp = poly1d.ClosedPoints(p, b_type);
269
270#ifndef MFEM_THREAD_SAFE
271 const int p1 = p + 1;
272
273 shape_x.SetSize(p1);
274 shape_y.SetSize(p1);
275 shape_z.SetSize(p1);
276 dshape_x.SetSize(p1);
277 dshape_y.SetSize(p1);
278 dshape_z.SetSize(p1);
279 d2shape_x.SetSize(p1);
280 d2shape_y.SetSize(p1);
281 d2shape_z.SetSize(p1);
282#endif
283
284 int o = 0;
285 for (int k = 0; k <= p; k++)
286 for (int j = 0; j <= p; j++)
287 for (int i = 0; i <= p; i++)
288 {
289 Nodes.IntPoint(dof_map[o++]).Set3(cp[i], cp[j], cp[k]);
290 }
291}
292
293void H1_HexahedronElement::CalcShape(const IntegrationPoint &ip,
294 Vector &shape) const
295{
296 const int p = order;
297
298#ifdef MFEM_THREAD_SAFE
299 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
300#endif
301
302 basis1d.Eval(ip.x, shape_x);
303 basis1d.Eval(ip.y, shape_y);
304 basis1d.Eval(ip.z, shape_z);
305
306 for (int o = 0, k = 0; k <= p; k++)
307 for (int j = 0; j <= p; j++)
308 for (int i = 0; i <= p; i++)
309 {
310 shape(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_z(k);
311 }
312}
313
314void H1_HexahedronElement::CalcDShape(const IntegrationPoint &ip,
315 DenseMatrix &dshape) const
316{
317 const int p = order;
318
319#ifdef MFEM_THREAD_SAFE
320 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
321 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
322#endif
323
324 basis1d.Eval(ip.x, shape_x, dshape_x);
325 basis1d.Eval(ip.y, shape_y, dshape_y);
326 basis1d.Eval(ip.z, shape_z, dshape_z);
327
328 for (int o = 0, k = 0; k <= p; k++)
329 for (int j = 0; j <= p; j++)
330 for (int i = 0; i <= p; i++)
331 {
332 dshape(dof_map[o],0) = dshape_x(i)* shape_y(j)* shape_z(k);
333 dshape(dof_map[o],1) = shape_x(i)*dshape_y(j)* shape_z(k);
334 dshape(dof_map[o],2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
335 }
336}
337
338void H1_HexahedronElement::CalcHessian(const IntegrationPoint &ip,
339 DenseMatrix &Hessian) const
340{
341 const int p = order;
342
343#ifdef MFEM_THREAD_SAFE
344 Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
345 Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
346 Vector d2shape_x(p+1), d2shape_y(p+1), d2shape_z(p+1);
347#endif
348
349 basis1d.Eval(ip.x, shape_x, dshape_x, d2shape_x);
350 basis1d.Eval(ip.y, shape_y, dshape_y, d2shape_y);
351 basis1d.Eval(ip.z, shape_z, dshape_z, d2shape_z);
352
353 for (int o = 0, k = 0; k <= p; k++)
354 for (int j = 0; j <= p; j++)
355 for (int i = 0; i <= p; i++)
356 {
357 Hessian(dof_map[o],0) = d2shape_x(i)* shape_y(j)* shape_z(k);
358 Hessian(dof_map[o],1) = dshape_x(i)* dshape_y(j)* shape_z(k);
359 Hessian(dof_map[o],2) = dshape_x(i)* shape_y(j)* dshape_z(k);
360 Hessian(dof_map[o],3) = shape_x(i)*d2shape_y(j)* shape_z(k);
361 Hessian(dof_map[o],4) = shape_x(i)* dshape_y(j)* dshape_z(k);
362 Hessian(dof_map[o],5) = shape_x(i)* shape_y(j)*d2shape_z(k);
363 o++;
364 }
365}
366
367void H1_HexahedronElement::ProjectDelta(int vertex, Vector &dofs) const
368{
369 const int p = order;
370 const real_t *cp = poly1d.ClosedPoints(p,b_type);
371
372#ifdef MFEM_THREAD_SAFE
373 Vector shape_x(p+1), shape_y(p+1);
374#endif
375
376 for (int i = 0; i <= p; i++)
377 {
378 shape_x(i) = poly1d.CalcDelta(p, (1.0 - cp[i]));
379 shape_y(i) = poly1d.CalcDelta(p, cp[i]);
380 }
381
382 switch (vertex)
383 {
384 case 0:
385 for (int o = 0, k = 0; k <= p; k++)
386 for (int j = 0; j <= p; j++)
387 for (int i = 0; i <= p; i++)
388 {
389 dofs(dof_map[o++]) = shape_x(i)*shape_x(j)*shape_x(k);
390 }
391 break;
392 case 1:
393 for (int o = 0, k = 0; k <= p; k++)
394 for (int j = 0; j <= p; j++)
395 for (int i = 0; i <= p; i++)
396 {
397 dofs(dof_map[o++]) = shape_y(i)*shape_x(j)*shape_x(k);
398 }
399 break;
400 case 2:
401 for (int o = 0, k = 0; k <= p; k++)
402 for (int j = 0; j <= p; j++)
403 for (int i = 0; i <= p; i++)
404 {
405 dofs(dof_map[o++]) = shape_y(i)*shape_y(j)*shape_x(k);
406 }
407 break;
408 case 3:
409 for (int o = 0, k = 0; k <= p; k++)
410 for (int j = 0; j <= p; j++)
411 for (int i = 0; i <= p; i++)
412 {
413 dofs(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_x(k);
414 }
415 break;
416 case 4:
417 for (int o = 0, k = 0; k <= p; k++)
418 for (int j = 0; j <= p; j++)
419 for (int i = 0; i <= p; i++)
420 {
421 dofs(dof_map[o++]) = shape_x(i)*shape_x(j)*shape_y(k);
422 }
423 break;
424 case 5:
425 for (int o = 0, k = 0; k <= p; k++)
426 for (int j = 0; j <= p; j++)
427 for (int i = 0; i <= p; i++)
428 {
429 dofs(dof_map[o++]) = shape_y(i)*shape_x(j)*shape_y(k);
430 }
431 break;
432 case 6:
433 for (int o = 0, k = 0; k <= p; k++)
434 for (int j = 0; j <= p; j++)
435 for (int i = 0; i <= p; i++)
436 {
437 dofs(dof_map[o++]) = shape_y(i)*shape_y(j)*shape_y(k);
438 }
439 break;
440 case 7:
441 for (int o = 0, k = 0; k <= p; k++)
442 for (int j = 0; j <= p; j++)
443 for (int i = 0; i <= p; i++)
444 {
445 dofs(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_y(k);
446 }
447 break;
448 }
449}
450
451H1_TriangleElement::H1_TriangleElement(const int p, const int btype)
452 : NodalFiniteElement(2, Geometry::TRIANGLE, ((p + 1)*(p + 2))/2, p,
453 FunctionSpace::Pk)
454{
455 const real_t *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
456
457#ifndef MFEM_THREAD_SAFE
458 shape_x.SetSize(p + 1);
459 shape_y.SetSize(p + 1);
460 shape_l.SetSize(p + 1);
461 dshape_x.SetSize(p + 1);
462 dshape_y.SetSize(p + 1);
463 dshape_l.SetSize(p + 1);
464 ddshape_x.SetSize(p + 1);
465 ddshape_y.SetSize(p + 1);
466 ddshape_l.SetSize(p + 1);
467 u.SetSize(dof);
468 du.SetSize(dof, dim);
469 ddu.SetSize(dof, (dim * (dim + 1)) / 2 );
470#else
471 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
472#endif
473
474 int p2p3 = 2*p + 3;
475 auto idx = [p2p3](int i, int j) { return ((p2p3-j)*j)/2+i; };
476 lex_ordering.SetSize(dof);
477
478 // vertices
479 lex_ordering[idx(0,0)] = 0;
480 Nodes.IntPoint(0).Set2(cp[0], cp[0]);
481 lex_ordering[idx(p,0)] = 1;
482 Nodes.IntPoint(1).Set2(cp[p], cp[0]);
483 lex_ordering[idx(0,p)] = 2;
484 Nodes.IntPoint(2).Set2(cp[0], cp[p]);
485
486 // edges
487 int o = 3;
488 for (int i = 1; i < p; i++)
489 {
490 lex_ordering[idx(i,0)] = o;
491 Nodes.IntPoint(o++).Set2(cp[i], cp[0]);
492 }
493 for (int i = 1; i < p; i++)
494 {
495 lex_ordering[idx(p-i,i)] = o;
496 Nodes.IntPoint(o++).Set2(cp[p-i], cp[i]);
497 }
498 for (int i = 1; i < p; i++)
499 {
500 lex_ordering[idx(0,p-i)] = o;
501 Nodes.IntPoint(o++).Set2(cp[0], cp[p-i]);
502 }
503
504 // interior
505 for (int j = 1; j < p; j++)
506 for (int i = 1; i + j < p; i++)
507 {
508 const real_t w = cp[i] + cp[j] + cp[p-i-j];
509 lex_ordering[idx(i,j)] = o;
510 Nodes.IntPoint(o++).Set2(cp[i]/w, cp[j]/w);
511 }
512
513 DenseMatrix T(dof);
514 for (int k = 0; k < dof; k++)
515 {
516 IntegrationPoint &ip = Nodes.IntPoint(k);
517 poly1d.CalcBasis(p, ip.x, shape_x);
518 poly1d.CalcBasis(p, ip.y, shape_y);
519 poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
520
521 o = 0;
522 for (int j = 0; j <= p; j++)
523 for (int i = 0; i + j <= p; i++)
524 {
525 T(o++, k) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
526 }
527 }
528
529 Ti.Factor(T);
530 // mfem::out << "H1_TriangleElement(" << p << ") : "; Ti.TestInversion();
531}
532
533void H1_TriangleElement::CalcShape(const IntegrationPoint &ip,
534 Vector &shape) const
535{
536 const int p = order;
537
538#ifdef MFEM_THREAD_SAFE
539 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1), u(dof);
540#endif
541
542 poly1d.CalcBasis(p, ip.x, shape_x);
543 poly1d.CalcBasis(p, ip.y, shape_y);
544 poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
545
546 for (int o = 0, j = 0; j <= p; j++)
547 for (int i = 0; i + j <= p; i++)
548 {
549 u(o++) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
550 }
551
552 Ti.Mult(u, shape);
553}
554
555void H1_TriangleElement::CalcDShape(const IntegrationPoint &ip,
556 DenseMatrix &dshape) const
557{
558 const int p = order;
559
560#ifdef MFEM_THREAD_SAFE
561 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
562 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
563 DenseMatrix du(dof, dim);
564#endif
565
566 poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
567 poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
568 poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l, dshape_l);
569
570 for (int o = 0, j = 0; j <= p; j++)
571 for (int i = 0; i + j <= p; i++)
572 {
573 int k = p - i - j;
574 du(o,0) = ((dshape_x(i)* shape_l(k)) -
575 ( shape_x(i)*dshape_l(k)))*shape_y(j);
576 du(o,1) = ((dshape_y(j)* shape_l(k)) -
577 ( shape_y(j)*dshape_l(k)))*shape_x(i);
578 o++;
579 }
580
581 Ti.Mult(du, dshape);
582}
583
584void H1_TriangleElement::CalcHessian(const IntegrationPoint &ip,
585 DenseMatrix &ddshape) const
586{
587 const int p = order;
588#ifdef MFEM_THREAD_SAFE
589 Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
590 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
591 Vector ddshape_x(p + 1), ddshape_y(p + 1), ddshape_l(p + 1);
592 DenseMatrix ddu(dof, dim);
593#endif
594
595 poly1d.CalcBasis(p, ip.x, shape_x, dshape_x, ddshape_x);
596 poly1d.CalcBasis(p, ip.y, shape_y, dshape_y, ddshape_y);
597 poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l, dshape_l, ddshape_l);
598
599 for (int o = 0, j = 0; j <= p; j++)
600 for (int i = 0; i + j <= p; i++)
601 {
602 int k = p - i - j;
603 // u_xx, u_xy, u_yy
604 ddu(o,0) = ((ddshape_x(i) * shape_l(k)) - 2. * (dshape_x(i) * dshape_l(k)) +
605 (shape_x(i) * ddshape_l(k))) * shape_y(j);
606 ddu(o,1) = (((shape_x(i) * ddshape_l(k)) - dshape_x(i) * dshape_l(k)) * shape_y(
607 j)) + (((dshape_x(i) * shape_l(k)) - (shape_x(i) * dshape_l(k))) * dshape_y(j));
608 ddu(o,2) = ((ddshape_y(j) * shape_l(k)) - 2. * (dshape_y(j) * dshape_l(k)) +
609 (shape_y(j) * ddshape_l(k))) * shape_x(i);
610 o++;
611 }
612
613 Ti.Mult(ddu, ddshape);
614}
615
616
617H1_TetrahedronElement::H1_TetrahedronElement(const int p, const int btype)
618 : NodalFiniteElement(3, Geometry::TETRAHEDRON, ((p + 1)*(p + 2)*(p + 3))/6,
619 p, FunctionSpace::Pk)
620{
621 const real_t *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
622
623#ifndef MFEM_THREAD_SAFE
624 shape_x.SetSize(p + 1);
625 shape_y.SetSize(p + 1);
626 shape_z.SetSize(p + 1);
627 shape_l.SetSize(p + 1);
628 dshape_x.SetSize(p + 1);
629 dshape_y.SetSize(p + 1);
630 dshape_z.SetSize(p + 1);
631 dshape_l.SetSize(p + 1);
632 ddshape_x.SetSize(p + 1);
633 ddshape_y.SetSize(p + 1);
634 ddshape_z.SetSize(p + 1);
635 ddshape_l.SetSize(p + 1);
636 u.SetSize(dof);
637 du.SetSize(dof, dim);
638 ddu.SetSize(dof, (dim * (dim + 1)) / 2);
639#else
640 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
641#endif
642
643 auto tri = [](int k) { return (k*(k + 1))/2; };
644 auto tet = [](int k) { return (k*(k + 1)*(k + 2))/6; };
645 int ndof = tet(p+1);
646 auto idx = [tri, tet, p, ndof](int i, int j, int k)
647 {
648 return ndof - tet(p - k) - tri(p + 1 - k - j) + i;
649 };
650
651 lex_ordering.SetSize(dof);
652
653 // vertices
654 lex_ordering[idx(0,0,0)] = 0;
655 Nodes.IntPoint(0).Set3(cp[0], cp[0], cp[0]);
656 lex_ordering[idx(p,0,0)] = 1;
657 Nodes.IntPoint(1).Set3(cp[p], cp[0], cp[0]);
658 lex_ordering[idx(0,p,0)] = 2;
659 Nodes.IntPoint(2).Set3(cp[0], cp[p], cp[0]);
660 lex_ordering[idx(0,0,p)] = 3;
661 Nodes.IntPoint(3).Set3(cp[0], cp[0], cp[p]);
662
663 // edges (see Tetrahedron::edges in mesh/tetrahedron.cpp)
664 int o = 4;
665 for (int i = 1; i < p; i++) // (0,1)
666 {
667 lex_ordering[idx(i,0,0)] = o;
668 Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[0]);
669 }
670 for (int i = 1; i < p; i++) // (0,2)
671 {
672 lex_ordering[idx(0,i,0)] = o;
673 Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[0]);
674 }
675 for (int i = 1; i < p; i++) // (0,3)
676 {
677 lex_ordering[idx(0,0,i)] = o;
678 Nodes.IntPoint(o++).Set3(cp[0], cp[0], cp[i]);
679 }
680 for (int i = 1; i < p; i++) // (1,2)
681 {
682 lex_ordering[idx(p-i,i,0)] = o;
683 Nodes.IntPoint(o++).Set3(cp[p-i], cp[i], cp[0]);
684 }
685 for (int i = 1; i < p; i++) // (1,3)
686 {
687 lex_ordering[idx(p-i,0,i)] = o;
688 Nodes.IntPoint(o++).Set3(cp[p-i], cp[0], cp[i]);
689 }
690 for (int i = 1; i < p; i++) // (2,3)
691 {
692 lex_ordering[idx(0,p-i,i)] = o;
693 Nodes.IntPoint(o++).Set3(cp[0], cp[p-i], cp[i]);
694 }
695
696 // faces (see Mesh::GenerateFaces in mesh/mesh.cpp)
697 for (int j = 1; j < p; j++)
698 for (int i = 1; i + j < p; i++) // (1,2,3)
699 {
700 lex_ordering[idx(p-i-j,i,j)] = o;
701 real_t w = cp[i] + cp[j] + cp[p-i-j];
702 Nodes.IntPoint(o++).Set3(cp[p-i-j]/w, cp[i]/w, cp[j]/w);
703 }
704 for (int j = 1; j < p; j++)
705 for (int i = 1; i + j < p; i++) // (0,3,2)
706 {
707 lex_ordering[idx(0,j,i)] = o;
708 real_t w = cp[i] + cp[j] + cp[p-i-j];
709 Nodes.IntPoint(o++).Set3(cp[0], cp[j]/w, cp[i]/w);
710 }
711 for (int j = 1; j < p; j++)
712 for (int i = 1; i + j < p; i++) // (0,1,3)
713 {
714 lex_ordering[idx(i,0,j)] = o;
715 real_t w = cp[i] + cp[j] + cp[p-i-j];
716 Nodes.IntPoint(o++).Set3(cp[i]/w, cp[0], cp[j]/w);
717 }
718 for (int j = 1; j < p; j++)
719 for (int i = 1; i + j < p; i++) // (0,2,1)
720 {
721 lex_ordering[idx(j,i,0)] = o;
722 real_t w = cp[i] + cp[j] + cp[p-i-j];
723 Nodes.IntPoint(o++).Set3(cp[j]/w, cp[i]/w, cp[0]);
724 }
725
726 // interior
727 for (int k = 1; k < p; k++)
728 for (int j = 1; j + k < p; j++)
729 for (int i = 1; i + j + k < p; i++)
730 {
731 lex_ordering[idx(i,j,k)] = o;
732 real_t w = cp[i] + cp[j] + cp[k] + cp[p-i-j-k];
733 Nodes.IntPoint(o++).Set3(cp[i]/w, cp[j]/w, cp[k]/w);
734 }
735
736 DenseMatrix T(dof);
737 for (int m = 0; m < dof; m++)
738 {
739 IntegrationPoint &ip = Nodes.IntPoint(m);
740 poly1d.CalcBasis(p, ip.x, shape_x);
741 poly1d.CalcBasis(p, ip.y, shape_y);
742 poly1d.CalcBasis(p, ip.z, shape_z);
743 poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
744
745 o = 0;
746 for (int k = 0; k <= p; k++)
747 for (int j = 0; j + k <= p; j++)
748 for (int i = 0; i + j + k <= p; i++)
749 {
750 T(o++, m) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
751 }
752 }
753
754 Ti.Factor(T);
755 // mfem::out << "H1_TetrahedronElement(" << p << ") : "; Ti.TestInversion();
756}
757
758void H1_TetrahedronElement::CalcShape(const IntegrationPoint &ip,
759 Vector &shape) const
760{
761 const int p = order;
762
763#ifdef MFEM_THREAD_SAFE
764 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
765 Vector u(dof);
766#endif
767
768 poly1d.CalcBasis(p, ip.x, shape_x);
769 poly1d.CalcBasis(p, ip.y, shape_y);
770 poly1d.CalcBasis(p, ip.z, shape_z);
771 poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
772
773 for (int o = 0, k = 0; k <= p; k++)
774 for (int j = 0; j + k <= p; j++)
775 for (int i = 0; i + j + k <= p; i++)
776 {
777 u(o++) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
778 }
779
780 Ti.Mult(u, shape);
781}
782
783void H1_TetrahedronElement::CalcDShape(const IntegrationPoint &ip,
784 DenseMatrix &dshape) const
785{
786 const int p = order;
787
788#ifdef MFEM_THREAD_SAFE
789 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
790 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
791 DenseMatrix du(dof, dim);
792#endif
793
794 poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
795 poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
796 poly1d.CalcBasis(p, ip.z, shape_z, dshape_z);
797 poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l, dshape_l);
798
799 for (int o = 0, k = 0; k <= p; k++)
800 for (int j = 0; j + k <= p; j++)
801 for (int i = 0; i + j + k <= p; i++)
802 {
803 int l = p - i - j - k;
804 du(o,0) = ((dshape_x(i)* shape_l(l)) -
805 ( shape_x(i)*dshape_l(l)))*shape_y(j)*shape_z(k);
806 du(o,1) = ((dshape_y(j)* shape_l(l)) -
807 ( shape_y(j)*dshape_l(l)))*shape_x(i)*shape_z(k);
808 du(o,2) = ((dshape_z(k)* shape_l(l)) -
809 ( shape_z(k)*dshape_l(l)))*shape_x(i)*shape_y(j);
810 o++;
811 }
812
813 Ti.Mult(du, dshape);
814}
815
816void H1_TetrahedronElement::CalcHessian(const IntegrationPoint &ip,
817 DenseMatrix &ddshape) const
818{
819 const int p = order;
820
821#ifdef MFEM_THREAD_SAFE
822 Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
823 Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
824 Vector ddshape_x(p + 1), ddshape_y(p + 1), ddshape_z(p + 1), ddshape_l(p + 1);
825 DenseMatrix ddu(dof, ((dim + 1) * dim) / 2);
826#endif
827
828 poly1d.CalcBasis(p, ip.x, shape_x, dshape_x, ddshape_x);
829 poly1d.CalcBasis(p, ip.y, shape_y, dshape_y, ddshape_y);
830 poly1d.CalcBasis(p, ip.z, shape_z, dshape_z, ddshape_z);
831 poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l, dshape_l, ddshape_l);
832
833 for (int o = 0, k = 0; k <= p; k++)
834 for (int j = 0; j + k <= p; j++)
835 for (int i = 0; i + j + k <= p; i++)
836 {
837 // u_xx, u_xy, u_xz, u_yy, u_yz, u_zz
838 int l = p - i - j - k;
839 ddu(o,0) = ((ddshape_x(i) * shape_l(l)) - 2. * (dshape_x(i) * dshape_l(l)) +
840 (shape_x(i) * ddshape_l(l))) * shape_y(j) * shape_z(k);
841 ddu(o,1) = ((dshape_y(j) * ((dshape_x(i) * shape_l(l)) -
842 (shape_x(i) * dshape_l(l)))) +
843 (shape_y(j) * ((ddshape_l(l) * shape_x(i)) -
844 (dshape_x(i) * dshape_l(l)))))* shape_z(k);
845 ddu(o,2) = ((dshape_z(k) * ((dshape_x(i) * shape_l(l)) -
846 (shape_x(i) * dshape_l(l)))) +
847 (shape_z(k) * ((ddshape_l(l) * shape_x(i)) -
848 (dshape_x(i) * dshape_l(l)))))* shape_y(j);
849 ddu(o,3) = ((ddshape_y(j) * shape_l(l)) - 2. * (dshape_y(j) * dshape_l(l)) +
850 (shape_y(j) * ddshape_l(l))) * shape_x(i) * shape_z(k);
851 ddu(o,4) = ((dshape_z(k) * ((dshape_y(j) * shape_l(l)) -
852 (shape_y(j)*dshape_l(l))) ) +
853 (shape_z(k)* ((ddshape_l(l)*shape_y(j)) -
854 (dshape_y(j) * dshape_l(l)) ) ) )* shape_x(i);
855 ddu(o,5) = ((ddshape_z(k) * shape_l(l)) - 2. * (dshape_z(k) * dshape_l(l)) +
856 (shape_z(k) * ddshape_l(l))) * shape_y(j) * shape_x(i);
857 o++;
858 }
859 Ti.Mult(ddu, ddshape);
860}
861
862// TODO: use a FunctionSpace specific to wedges instead of Qk.
863H1_WedgeElement::H1_WedgeElement(const int p,
864 const int btype)
865 : NodalFiniteElement(3, Geometry::PRISM, ((p + 1)*(p + 1)*(p + 2))/2,
866 p, FunctionSpace::Qk),
867 TriangleFE(p, btype),
868 SegmentFE(p, btype)
869{
870#ifndef MFEM_THREAD_SAFE
871 t_shape.SetSize(TriangleFE.GetDof());
872 s_shape.SetSize(SegmentFE.GetDof());
873 t_dshape.SetSize(TriangleFE.GetDof(), 2);
874 s_dshape.SetSize(SegmentFE.GetDof(), 1);
875#endif
876
877 t_dof.SetSize(dof);
878 s_dof.SetSize(dof);
879
880 int p2p3 = 2*p + 3, ntri = ((p + 1)*(p + 2))/2;
881 auto idx = [p2p3,ntri](int i, int j, int k)
882 {
883 return k*ntri + ((p2p3-j)*j)/2+i;
884 };
885
886 lex_ordering.SetSize(dof);
887 int o = 0;
888
889 // Nodal DoFs
890 lex_ordering[idx(0,0,0)] = o++;
891 lex_ordering[idx(p,0,0)] = o++;
892 lex_ordering[idx(0,p,0)] = o++;
893 lex_ordering[idx(0,0,p)] = o++;
894 lex_ordering[idx(p,0,p)] = o++;
895 lex_ordering[idx(0,p,p)] = o++;
896 t_dof[0] = 0; s_dof[0] = 0;
897 t_dof[1] = 1; s_dof[1] = 0;
898 t_dof[2] = 2; s_dof[2] = 0;
899 t_dof[3] = 0; s_dof[3] = 1;
900 t_dof[4] = 1; s_dof[4] = 1;
901 t_dof[5] = 2; s_dof[5] = 1;
902
903 // Edge DoFs
904 int k = 0;
905 int ne = p-1;
906 for (int i=1; i<p; i++)
907 {
908 lex_ordering[idx(i,0,0)] = o + 0*ne + k;
909 lex_ordering[idx(p-i,i,0)] = o + 1*ne + k;
910 lex_ordering[idx(0,p-i,0)] = o + 2*ne + k;
911 lex_ordering[idx(i,0,p)] = o + 3*ne + k;
912 lex_ordering[idx(p-i,i,p)] = o + 4*ne + k;
913 lex_ordering[idx(0,p-i,p)] = o + 5*ne + k;
914 lex_ordering[idx(0,0,i)] = o + 6*ne + k;
915 lex_ordering[idx(p,0,i)] = o + 7*ne + k;
916 lex_ordering[idx(0,p,i)] = o + 8*ne + k;
917 t_dof[5 + 0 * ne + i] = 2 + 0 * ne + i; s_dof[5 + 0 * ne + i] = 0;
918 t_dof[5 + 1 * ne + i] = 2 + 1 * ne + i; s_dof[5 + 1 * ne + i] = 0;
919 t_dof[5 + 2 * ne + i] = 2 + 2 * ne + i; s_dof[5 + 2 * ne + i] = 0;
920 t_dof[5 + 3 * ne + i] = 2 + 0 * ne + i; s_dof[5 + 3 * ne + i] = 1;
921 t_dof[5 + 4 * ne + i] = 2 + 1 * ne + i; s_dof[5 + 4 * ne + i] = 1;
922 t_dof[5 + 5 * ne + i] = 2 + 2 * ne + i; s_dof[5 + 5 * ne + i] = 1;
923 t_dof[5 + 6 * ne + i] = 0; s_dof[5 + 6 * ne + i] = i + 1;
924 t_dof[5 + 7 * ne + i] = 1; s_dof[5 + 7 * ne + i] = i + 1;
925 t_dof[5 + 8 * ne + i] = 2; s_dof[5 + 8 * ne + i] = i + 1;
926 ++k;
927 }
928 o += 9*ne;
929
930 // Triangular Face DoFs
931 k=0;
932 int nt = (p-1)*(p-2)/2;
933 for (int j=1; j<p; j++)
934 {
935 for (int i=1; i<p-j; i++)
936 {
937 int l = j - p + (((2 * p - 1) - i) * i) / 2;
938 lex_ordering[idx(i,j,0)] = o+l;
939 lex_ordering[idx(i,j,p)] = o+nt+k;
940 t_dof[6 + 9 * ne + k] = 3 * p + l; s_dof[6 + 9 * ne + k] = 0;
941 t_dof[6 + 9 * ne + nt + k] = 3 * p + k; s_dof[6 + 9 * ne + nt + k] = 1;
942 k++;
943 }
944 }
945 o += 2*nt;
946
947 // Quadrilateral Face DoFs
948 k=0;
949 int nq = (p-1)*(p-1);
950 for (int j=1; j<p; j++)
951 {
952 for (int i=1; i<p; i++)
953 {
954 lex_ordering[idx(i,0,j)] = o+k;
955 lex_ordering[idx(p-i,i,j)] = o+nq+k;
956 lex_ordering[idx(0,p-i,j)] = o+2*nq+k;
957
958 t_dof[6 + 9 * ne + 2 * nt + 0 * nq + k] = 2 + 0 * ne + i;
959 t_dof[6 + 9 * ne + 2 * nt + 1 * nq + k] = 2 + 1 * ne + i;
960 t_dof[6 + 9 * ne + 2 * nt + 2 * nq + k] = 2 + 2 * ne + i;
961
962 s_dof[6 + 9 * ne + 2 * nt + 0 * nq + k] = 1 + j;
963 s_dof[6 + 9 * ne + 2 * nt + 1 * nq + k] = 1 + j;
964 s_dof[6 + 9 * ne + 2 * nt + 2 * nq + k] = 1 + j;
965
966 k++;
967 }
968 }
969 o += 3*nq;
970
971 // Interior DoFs
972 int m=0;
973 for (k=1; k<p; k++)
974 {
975 int l=0;
976 for (int j=1; j<p; j++)
977 {
978 for (int i=1; i+j<p; i++)
979 {
980 lex_ordering[idx(i,j,k)] = o++;
981 t_dof[6 + 9 * ne + 2 * nt + 3 * nq + m] = 3 * p + l;
982 s_dof[6 + 9 * ne + 2 * nt + 3 * nq + m] = 1 + k;
983 l++; m++;
984 }
985 }
986 }
987
988 // Define Nodes
989 const IntegrationRule & t_Nodes = TriangleFE.GetNodes();
990 const IntegrationRule & s_Nodes = SegmentFE.GetNodes();
991 for (int i=0; i<dof; i++)
992 {
993 Nodes.IntPoint(i).x = t_Nodes.IntPoint(t_dof[i]).x;
994 Nodes.IntPoint(i).y = t_Nodes.IntPoint(t_dof[i]).y;
995 Nodes.IntPoint(i).z = s_Nodes.IntPoint(s_dof[i]).x;
996 }
997}
998
999void H1_WedgeElement::CalcShape(const IntegrationPoint &ip,
1000 Vector &shape) const
1001{
1002#ifdef MFEM_THREAD_SAFE
1003 Vector t_shape(TriangleFE.GetDof());
1004 Vector s_shape(SegmentFE.GetDof());
1005#endif
1006
1007 IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
1008
1009 TriangleFE.CalcShape(ip, t_shape);
1010 SegmentFE.CalcShape(ipz, s_shape);
1011
1012 for (int i=0; i<dof; i++)
1013 {
1014 shape[i] = t_shape[t_dof[i]] * s_shape[s_dof[i]];
1015 }
1016}
1017
1018void H1_WedgeElement::CalcDShape(const IntegrationPoint &ip,
1019 DenseMatrix &dshape) const
1020{
1021#ifdef MFEM_THREAD_SAFE
1022 Vector t_shape(TriangleFE.GetDof());
1023 DenseMatrix t_dshape(TriangleFE.GetDof(), 2);
1024 Vector s_shape(SegmentFE.GetDof());
1025 DenseMatrix s_dshape(SegmentFE.GetDof(), 1);
1026#endif
1027
1028 IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
1029
1030 TriangleFE.CalcShape(ip, t_shape);
1031 TriangleFE.CalcDShape(ip, t_dshape);
1032 SegmentFE.CalcShape(ipz, s_shape);
1033 SegmentFE.CalcDShape(ipz, s_dshape);
1034
1035 for (int i=0; i<dof; i++)
1036 {
1037 dshape(i, 0) = t_dshape(t_dof[i],0) * s_shape[s_dof[i]];
1038 dshape(i, 1) = t_dshape(t_dof[i],1) * s_shape[s_dof[i]];
1039 dshape(i, 2) = t_shape[t_dof[i]] * s_dshape(s_dof[i],0);
1040 }
1041}
1042
1043}
mfem::Array::SetSize
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
mfem::DenseMatrixInverse::Factor
void Factor()
Factor the current DenseMatrix, *a.
Definition densemat.cpp:4185
mfem::DenseMatrixInverse::Mult
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition densemat.cpp:4231
mfem::DenseMatrix
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
mfem::DenseMatrix::SetSize
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:105
mfem::FiniteElement::VerifyNodal
static int VerifyNodal(int b_type)
Ensure that the BasisType of b_type nodal (satisfies the interpolation property).
Definition fe_base.hpp:647
mfem::FiniteElement::dof
int dof
Number of degrees of freedom.
Definition fe_base.hpp:248
mfem::FiniteElement::Nodes
IntegrationRule Nodes
Definition fe_base.hpp:251
mfem::FiniteElement::VerifyClosed
static int VerifyClosed(int b_type)
Ensure that the BasisType of b_type is closed (has Quadrature1D points on the boundary).
Definition fe_base.hpp:630
mfem::FiniteElement::GetNodes
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition fe_base.hpp:395
mfem::FiniteElement::GetDof
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition fe_base.hpp:329
mfem::FiniteElement::order
int order
Order/degree of the shape functions.
Definition fe_base.hpp:249
mfem::FiniteElement::dim
int dim
Dimension of reference space.
Definition fe_base.hpp:241
mfem::FunctionSpace
Describes the function space on each element.
Definition fe_base.hpp:222
mfem::H1_HexahedronElement::H1_HexahedronElement
H1_HexahedronElement(const int p, const int btype=BasisType::GaussLobatto)
Construct the H1_HexahedronElement of order p and BasisType btype.
Definition fe_h1.cpp:265
mfem::H1_HexahedronElement::CalcHessian
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &Hessian) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition fe_h1.cpp:338
mfem::H1_HexahedronElement::CalcShape
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_h1.cpp:293
mfem::H1_HexahedronElement::ProjectDelta
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition fe_h1.cpp:367
mfem::H1_HexahedronElement::CalcDShape
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_h1.cpp:314
mfem::H1_QuadrilateralElement::CalcHessian
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &Hessian) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition fe_h1.cpp:192
mfem::H1_QuadrilateralElement::CalcShape
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_h1.cpp:151
mfem::H1_QuadrilateralElement::ProjectDelta
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition fe_h1.cpp:216
mfem::H1_QuadrilateralElement::CalcDShape
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_h1.cpp:170
mfem::H1_QuadrilateralElement::H1_QuadrilateralElement
H1_QuadrilateralElement(const int p, const int btype=BasisType::GaussLobatto)
Construct the H1_QuadrilateralElement of order p and BasisType btype.
Definition fe_h1.cpp:125
mfem::H1_SegmentElement::CalcDShape
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_h1.cpp:59
mfem::H1_SegmentElement::H1_SegmentElement
H1_SegmentElement(const int p, const int btype=BasisType::GaussLobatto)
Construct the H1_SegmentElement of order p and BasisType btype.
Definition fe_h1.cpp:21
mfem::H1_SegmentElement::CalcHessian
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &Hessian) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition fe_h1.cpp:78
mfem::H1_SegmentElement::CalcShape
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_h1.cpp:40
mfem::H1_SegmentElement::ProjectDelta
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition fe_h1.cpp:97
mfem::H1_TetrahedronElement::CalcShape
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_h1.cpp:758
mfem::H1_TetrahedronElement::H1_TetrahedronElement
H1_TetrahedronElement(const int p, const int btype=BasisType::GaussLobatto)
Construct the H1_TetrahedronElement of order p and BasisType btype.
Definition fe_h1.cpp:617
mfem::H1_TetrahedronElement::CalcHessian
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &ddshape) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition fe_h1.cpp:816
mfem::H1_TetrahedronElement::CalcDShape
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_h1.cpp:783
mfem::H1_TriangleElement::CalcShape
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_h1.cpp:533
mfem::H1_TriangleElement::CalcHessian
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &ddshape) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition fe_h1.cpp:584
mfem::H1_TriangleElement::H1_TriangleElement
H1_TriangleElement(const int p, const int btype=BasisType::GaussLobatto)
Construct the H1_TriangleElement of order p and BasisType btype.
Definition fe_h1.cpp:451
mfem::H1_TriangleElement::CalcDShape
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_h1.cpp:555
mfem::H1_WedgeElement::CalcShape
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_h1.cpp:999
mfem::H1_WedgeElement::H1_WedgeElement
H1_WedgeElement(const int p, const int btype=BasisType::GaussLobatto)
Construct the H1_WedgeElement of order p and BasisType btype.
Definition fe_h1.cpp:863
mfem::H1_WedgeElement::CalcDShape
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_h1.cpp:1018
mfem::IntegrationPoint
Class for integration point with weight.
Definition intrules.hpp:35
mfem::IntegrationPoint::Set2
void Set2(const real_t x1, const real_t x2)
Definition intrules.hpp:89
mfem::IntegrationPoint::Set3
void Set3(const real_t x1, const real_t x2, const real_t x3)
Definition intrules.hpp:79
mfem::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
mfem::IntegrationRule::IntPoint
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition intrules.hpp:259
mfem::NodalFiniteElement
Class for standard nodal finite elements.
Definition fe_base.hpp:715
mfem::NodalFiniteElement::lex_ordering
Array< int > lex_ordering
Definition fe_base.hpp:720
mfem::Poly_1D::Basis::Eval
void Eval(const real_t x, Vector &u) const
Evaluate the basis functions at point x in [0,1].
Definition fe_base.cpp:1761
mfem::Poly_1D::CalcDelta
static real_t CalcDelta(const int p, const real_t x)
Evaluate a representation of a Delta function at point x.
Definition fe_base.hpp:1139
mfem::Poly_1D::ClosedPoints
const real_t * ClosedPoints(const int p, const int btype=BasisType::GaussLobatto)
Get coordinates of a closed (GaussLegendre) set of points if degree p.
Definition fe_base.hpp:1083
mfem::Poly_1D::CalcBasis
static void CalcBasis(const int p, const real_t x, real_t *u)
Evaluate the values of a hierarchical 1D basis at point x hierarchical = k-th basis function is degre...
Definition fe_base.hpp:1099
mfem::TensorBasisElement::basis1d
Poly_1D::Basis & basis1d
Definition fe_base.hpp:1200
mfem::Vector
Vector data type.
Definition vector.hpp:80
mfem::Vector::SetSize
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
fe_h1.hpp
mfem::SegmentFE
Linear1DFiniteElement SegmentFE
Definition segment.cpp:52
mfem::u
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
mfem::real_t
float real_t
Definition config.hpp:43
mfem::poly1d
Poly_1D poly1d
Definition fe.cpp:28
mfem::TriangleFE
MFEM_EXPORT Linear2DFiniteElement TriangleFE
Definition fe.cpp:32
p
real_t p(const Vector &x, real_t t)
Definition navier_mms.cpp:53
doxygen 1.11.0