MFEM v4.8.0
Finite element discretization library
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
fe_pyramid.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// Finite Element classes on Pyramid shaped elements
13
14#include "fe_pyramid.hpp"
15
16namespace mfem
17{
18
19using namespace std;
20
21Vector FuentesPyramid::grad_lam1(real_t x, real_t y, real_t z)
22{
23 return Vector({CheckZ(z) ? - (1.0 - y - z) / (1.0 - z) : -0.5,
24 CheckZ(z) ? - (1.0 - x - z) / (1.0 - z) : -0.5,
25 CheckZ(z) ? x * y / ((1.0 - z) * (1.0 - z)) - 1.0 : -0.75});
26}
27
28Vector FuentesPyramid::grad_lam2(real_t x, real_t y, real_t z)
29{
30 return Vector({CheckZ(z) ? (1.0 - y - z) / (1.0 - z) : 0.5,
31 CheckZ(z) ? - x / (1.0 - z) : -0.5,
32 CheckZ(z) ? - x * y / ((1.0 - z) * (1.0 - z)) : -0.25});
33}
34
35Vector FuentesPyramid::grad_lam3(real_t x, real_t y, real_t z)
36{
37 return Vector({CheckZ(z) ? y / (1.0 - z) : 0.5,
38 CheckZ(z) ? x / (1.0 - z) : 0.5,
39 CheckZ(z) ? x * y / ((1.0 - z) * (1.0 - z)) : 0.25});
40}
41
42Vector FuentesPyramid::grad_lam4(real_t x, real_t y, real_t z)
43{
44 return Vector({CheckZ(z) ? - y / (1.0 - z) : -0.5,
45 CheckZ(z) ? (1.0 - x - z) / (1.0 - z) : 0.5,
46 CheckZ(z) ? - x * y / ((1.0 - z) * (1.0 - z)) : -0.25});
47}
48
49Vector FuentesPyramid::grad_lam5(real_t x, real_t y, real_t z)
50{
51 return Vector({0.0, 0.0, 1.0});
52}
53
54DenseMatrix FuentesPyramid::grad_lam15(real_t x, real_t y, real_t z)
55{
56 DenseMatrix dlam(2, 3);
57 dlam.SetRow(0, grad_lam1(x, y, z));
58 dlam.SetRow(1, grad_lam5(x, y, z));
59 return dlam;
60}
61
62DenseMatrix FuentesPyramid::grad_lam25(real_t x, real_t y, real_t z)
63{
64 DenseMatrix dlam(2, 3);
65 dlam.SetRow(0, grad_lam2(x, y, z));
66 dlam.SetRow(1, grad_lam5(x, y, z));
67 return dlam;
68}
69
70DenseMatrix FuentesPyramid::grad_lam35(real_t x, real_t y, real_t z)
71{
72 DenseMatrix dlam(2, 3);
73 dlam.SetRow(0, grad_lam3(x, y, z));
74 dlam.SetRow(1, grad_lam5(x, y, z));
75 return dlam;
76}
77
78DenseMatrix FuentesPyramid::grad_lam45(real_t x, real_t y, real_t z)
79{
80 DenseMatrix dlam(2, 3);
81 dlam.SetRow(0, grad_lam4(x, y, z));
82 dlam.SetRow(1, grad_lam5(x, y, z));
83 return dlam;
84}
85
86Vector FuentesPyramid::lam15_grad_lam15(real_t x, real_t y, real_t z)
87{
88 Vector lam = lam15(x, y, z);
89 Vector lamdlam(3);
90 add(lam(0), grad_lam5(x, y, z), -lam(1), grad_lam1(x, y, z), lamdlam);
91 return lamdlam;
92}
93
94Vector FuentesPyramid::lam25_grad_lam25(real_t x, real_t y, real_t z)
95{
96 Vector lam = lam25(x, y, z);
97 Vector lamdlam(3);
98 add(lam(0), grad_lam5(x, y, z), -lam(1), grad_lam2(x, y, z), lamdlam);
99 return lamdlam;
100}
101
102Vector FuentesPyramid::lam35_grad_lam35(real_t x, real_t y, real_t z)
103{
104 Vector lam = lam35(x, y, z);
105 Vector lamdlam(3);
106 add(lam(0), grad_lam5(x, y, z), -lam(1), grad_lam3(x, y, z), lamdlam);
107 return lamdlam;
108}
109
110Vector FuentesPyramid::lam45_grad_lam45(real_t x, real_t y, real_t z)
111{
112 Vector lam = lam45(x, y, z);
113 Vector lamdlam(3);
114 add(lam(0), grad_lam5(x, y, z), -lam(1), grad_lam4(x, y, z), lamdlam);
115 return lamdlam;
116}
117
118Vector FuentesPyramid::lam125_grad_lam125(real_t x, real_t y, real_t z)
119{
120 Vector lgl({-x * z / (one - z), y - one, z});
121 lgl *= (one - y - z) / (one - z);
122 return lgl;
123}
124
125Vector FuentesPyramid::lam235_grad_lam235(real_t x, real_t y, real_t z)
126{
127 Vector lgl({x, -y * z / (one - z), z});
128 lgl *= x / (one - z);
129 return lgl;
130}
131
132Vector FuentesPyramid::lam345_grad_lam345(real_t x, real_t y, real_t z)
133{
134 Vector lgl({-x * z / (one - z), y, z});
135 lgl *= y / (one - z);
136 return lgl;
137}
138
139Vector FuentesPyramid::lam435_grad_lam435(real_t x, real_t y, real_t z)
140{
141 Vector lgl({x * z / (one - z), -y, -z});
142 lgl *= y / (one - z);
143 return lgl;
144}
145
146Vector FuentesPyramid::lam415_grad_lam415(real_t x, real_t y, real_t z)
147{
148 Vector lgl({x - one, -y * z / (one - z), z});
149 lgl *= (one - x - z) / (one - z);
150 return lgl;
151}
152
153Vector FuentesPyramid::lam145_grad_lam145(real_t x, real_t y, real_t z)
154{
155 Vector lgl({one - x, y * z / (one - z), -z});
156 lgl *= (one - x - z) / (one - z);
157 return lgl;
158}
159
160real_t FuentesPyramid::div_lam125_grad_lam125(real_t x, real_t y, real_t z)
161{ return (1.0 - z - y) / (1.0 - z); }
162
165
168
169real_t FuentesPyramid::div_lam435_grad_lam435(real_t x, real_t y, real_t z)
170{ return -y / (1.0 - z); }
171
172real_t FuentesPyramid::div_lam415_grad_lam415(real_t x, real_t y, real_t z)
173{ return (1.0 - z - x) / (1.0 - z); }
174
175real_t FuentesPyramid::div_lam145_grad_lam145(real_t x, real_t y, real_t z)
176{ return -(1.0 - z - x) / (1.0 - z); }
177
178DenseMatrix FuentesPyramid::grad_mu01(real_t z)
179{
180 DenseMatrix dmu(2, 3);
181 dmu.SetRow(0, grad_mu0(z));
182 dmu.SetRow(1, grad_mu1(z));
183 return dmu;
184}
185
186Vector FuentesPyramid::grad_mu0(real_t z, const Vector xy, unsigned int ab)
187{
188 Vector dmu({0.0, 0.0, - xy[ab-1] / pow(1.0 - z, 2)});
189 dmu[ab-1] = -1.0 / (1.0 - z);
190 return dmu;
191}
192
193Vector FuentesPyramid::grad_mu1(real_t z, const Vector xy, unsigned int ab)
194{
195 Vector dmu({0.0, 0.0, xy[ab-1] / pow(1.0 - z, 2)});
196 dmu[ab-1] = 1.0 / (1.0 - z);
197 return dmu;
198}
199
200DenseMatrix FuentesPyramid::grad_mu01(real_t z, Vector xy, unsigned int ab)
201{
202 DenseMatrix dmu(2, 3);
203 dmu.SetRow(0, grad_mu0(z, xy, ab));
204 dmu.SetRow(1, grad_mu1(z, xy, ab));
205 return dmu;
206}
207
208Vector FuentesPyramid::mu01_grad_mu01(real_t z, Vector xy, unsigned int ab)
209{
210 Vector mu = mu01(z, xy, ab);
211 Vector mudmu(3);
212 add(mu(0), grad_mu1(z, xy, ab), -mu(1), grad_mu0(z, xy, ab), mudmu);
213 return mudmu;
214}
215
216Vector FuentesPyramid::grad_nu0(real_t z, const Vector xy, unsigned int ab)
217{
218 Vector dnu({0.0, 0.0, -1.0}); dnu[ab-1] = -1.0;
219 return dnu;
220}
221
222Vector FuentesPyramid::grad_nu1(real_t z, const Vector xy, unsigned int ab)
223{
224 Vector dnu({0.0, 0.0, 0.0}); dnu[ab-1] = 1.0;
225 return dnu;
226}
227
228Vector FuentesPyramid::grad_nu2(real_t z, const Vector xy, unsigned int ab)
229{
230 return Vector({0.0, 0.0, 1.0});
231}
232
233DenseMatrix FuentesPyramid::grad_nu01(real_t z, Vector xy, unsigned int ab)
234{
235 DenseMatrix dnu(2, 3);
236 dnu.SetRow(0, grad_nu0(z, xy, ab));
237 dnu.SetRow(1, grad_nu1(z, xy, ab));
238 return dnu;
239}
240
241DenseMatrix FuentesPyramid::grad_nu012(real_t z, Vector xy, unsigned int ab)
242{
243 DenseMatrix dnu(3, 3);
244 dnu.SetRow(0, grad_nu0(z, xy, ab));
245 dnu.SetRow(1, grad_nu1(z, xy, ab));
246 dnu.SetRow(2, grad_nu2(z, xy, ab));
247 return dnu;
248}
249
250DenseMatrix FuentesPyramid::grad_nu120(real_t z, Vector xy, unsigned int ab)
251{
252 DenseMatrix dnu(3, 3);
253 dnu.SetRow(0, grad_nu1(z, xy, ab));
254 dnu.SetRow(1, grad_nu2(z, xy, ab));
255 dnu.SetRow(2, grad_nu0(z, xy, ab));
256 return dnu;
257}
258
259Vector FuentesPyramid::nu01_grad_nu01(real_t z, Vector xy, unsigned int ab)
260{
261 Vector nu = nu01(z, xy, ab);
262 Vector nudnu(3);
263 add(nu(0), grad_nu1(z, xy, ab), -nu(1), grad_nu0(z, xy, ab), nudnu);
264 return nudnu;
265}
266
267Vector FuentesPyramid::nu12_grad_nu12(real_t z, Vector xy, unsigned int ab)
268{
269 Vector nu = nu12(z, xy, ab);
270 Vector nudnu(3);
271 add(nu(0), grad_nu2(z, xy, ab), -nu(1), grad_nu1(z, xy, ab), nudnu);
272 return nudnu;
273}
274
275Vector FuentesPyramid::nu012_grad_nu012(real_t z, Vector xy, unsigned int ab)
276{
277 Vector nu(nu012(z, xy, ab));
278 Vector dnu0(grad_nu0(z, xy, ab));
279 Vector dnu1(grad_nu1(z, xy, ab));
280 Vector dnu2(grad_nu2(z, xy, ab));
281
282 Vector v01(3), v12(3), v20(3);
283 dnu0.cross3D(dnu1, v01);
284 dnu1.cross3D(dnu2, v12);
285 dnu2.cross3D(dnu0, v20);
286
287 Vector nudnu(3);
288 add(nu(0), v12, nu(1), v20, nudnu);
289 nudnu.Add(nu(2), v01);
290 return nudnu;
291}
292
293void FuentesPyramid::CalcScaledLegendre(int p, real_t x, real_t t,
294 real_t *u)
295{
296 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
297 if (t > 0.0)
298 {
299 Poly_1D::CalcLegendre(p, x / t, u);
300 for (int i = 1; i <= p; i++)
301 {
302 u[i] *= pow(t, i);
303 }
304 }
305 else
306 {
307 // This assumes x = 0 as well as t = 0 since x \in [0,t]
308 u[0] = 1.0;
309 for (int i = 1; i <= p; i++) { u[i] = 0.0; }
310 }
311}
312
313void FuentesPyramid::CalcScaledLegendre(int p, real_t x, real_t t,
314 real_t *u,
315 real_t *dudx, real_t *dudt)
316{
317 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
318 if (t > 0.0)
319 {
320 Poly_1D::CalcLegendre(p, x / t, u, dudx);
321 dudx[0] = 0.0;
322 dudt[0] = - dudx[0] * x / t;
323 for (int i = 1; i <= p; i++)
324 {
325 u[i] *= pow(t, i);
326 dudx[i] *= pow(t, i - 1);
327 dudt[i] = (u[i] * i - dudx[i] * x) / t;
328 }
329 }
330 else
331 {
332 // This assumes x = 0 as well as t = 0 since x \in [0,t]
333 u[0] = 1.0;
334 dudx[0] = 0.0;
335 dudt[0] = 0.0;
336 if (p >=1)
337 {
338 u[1] = 0.0;
339 dudx[1] = 2.0;
340 dudt[1] = -1.0;
341 }
342 for (int i = 2; i <= p; i++)
343 {
344 u[i] = 0.0;
345 dudx[i] = 0.0;
346 dudt[i] = 0.0;
347 }
348 }
349}
350
351void FuentesPyramid::CalcScaledLegendre(int p, real_t x, real_t t,
352 Vector &u)
353{
354 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
355 MFEM_ASSERT(u.Size() >= p+1, "Size of u is too small");
356 CalcScaledLegendre(p, x, t, u.GetData());
357}
358
359void FuentesPyramid::CalcScaledLegendre(int p, real_t x, real_t t,
360 Vector &u, Vector &dudx, Vector &dudt)
361{
362 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
363 MFEM_ASSERT(u.Size() >= p+1, "Size of u is too small");
364 MFEM_ASSERT(dudx.Size() >= p+1, "Size of dudx is too small");
365 MFEM_ASSERT(dudt.Size() >= p+1, "Size of dudt is too small");
366 CalcScaledLegendre(p, x, t, u.GetData(), dudx.GetData(), dudt.GetData());
367}
368
369void FuentesPyramid::CalcIntegratedLegendre(int p, real_t x, real_t t,
370 real_t *u)
371{
372 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
373 if (t > 0.0)
374 {
375 CalcScaledLegendre(p, x, t, u);
376 for (int i = p; i >= 2; i--)
377 {
378 u[i] = (u[i] - t * t * u[i-2]) / (4.0 * i - 2.0);
379 }
380 if (p >= 1)
381 {
382 u[1] = x;
383 }
384 u[0] = 0.0;
385 }
386 else
387 {
388 for (int i = 0; i <= p; i++)
389 {
390 u[i] = 0.0;
391 }
392 }
393}
394
395void FuentesPyramid::CalcIntegratedLegendre(int p, real_t x, real_t t,
396 real_t *u,
397 real_t *dudx, real_t *dudt)
398{
399 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
400 if (t > 0.0)
401 {
402 CalcScaledLegendre(p, x, t, u, dudx, dudt);
403 for (int i = p; i >= 2; i--)
404 {
405 u[i] = (u[i] - t * t * u[i-2]) / (4.0 * i - 2.0);
406 dudx[i] = (dudx[i] - t * t * dudx[i-2]) / (4.0 * i - 2.0);
407 dudt[i] = (dudt[i] - t * t * dudt[i-2] - 2.0 * t * u[i-2]) /
408 (4.0 * i - 2.0);
409 }
410 if (p >= 1)
411 {
412 u[1] = x; dudx[1] = 1.0; dudt[1] = 0.0;
413 }
414 u[0] = 0.0; dudx[0] = 0.0; dudt[0] = 0.0;
415 }
416 else
417 {
418 for (int i = 0; i <= p; i++)
419 {
420 u[i] = 0.0;
421 dudx[i] = 0.0;
422 dudt[i] = 0.0;
423 }
424 }
425}
426
427void FuentesPyramid::CalcIntegratedLegendre(int p, real_t x,
428 real_t t, Vector &u)
429{
430 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
431 MFEM_ASSERT(u.Size() >= p+1, "Size of u is too small");
432 CalcIntegratedLegendre(p, x, t, u.GetData());
433}
434
435void FuentesPyramid::CalcIntegratedLegendre(int p, real_t x,
436 real_t t, Vector &u,
437 Vector &dudx, Vector &dudt)
438{
439 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
440 MFEM_ASSERT(u.Size() >= p+1, "Size of u is too small");
441 MFEM_ASSERT(dudx.Size() >= p+1, "Size of dudx is too small");
442 MFEM_ASSERT(dudt.Size() >= p+1, "Size of dudt is too small");
443 CalcIntegratedLegendre(p, x, t, u.GetData(),
444 dudx.GetData(), dudt.GetData());
445}
446
447void FuentesPyramid::CalcHomogenizedScaLegendre(int p, real_t s0, real_t s1,
448 real_t *u)
449{
450 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
451 CalcScaledLegendre(p, s1, s0 + s1, u);
452}
453
454void FuentesPyramid::CalcHomogenizedScaLegendre(int p,
455 real_t s0, real_t s1,
456 real_t *u,
457 real_t *duds0, real_t *duds1)
458{
459 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
460 CalcScaledLegendre(p, s1, s0+s1, u, duds1, duds0);
461 for (int i = 0; i <= p; i++) { duds1[i] += duds0[i]; }
462}
463
464void FuentesPyramid::CalcHomogenizedScaLegendre(int p, real_t s0, real_t s1,
465 Vector &u)
466{
467 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
468 MFEM_ASSERT(u.Size() >= p+1, "Size of u is too small");
469 CalcHomogenizedScaLegendre(p, s0, s1, u.GetData());
470}
471
472void FuentesPyramid::CalcHomogenizedScaLegendre(int p,
473 real_t s0, real_t s1,
474 Vector &u,
475 Vector &duds0, Vector &duds1)
476{
477 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
478 MFEM_ASSERT(u.Size() >= p+1, "Size of u is too small");
479 MFEM_ASSERT(duds0.Size() >= p+1, "Size of duds0 is too small");
480 MFEM_ASSERT(duds1.Size() >= p+1, "Size of duds1 is too small");
481 CalcHomogenizedScaLegendre(p, s0, s1, u.GetData(),
482 duds0.GetData(), duds1.GetData());
483}
484
485void FuentesPyramid::CalcHomogenizedIntLegendre(int p,
486 real_t t0, real_t t1,
487 real_t *u)
488{
489 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
490 CalcIntegratedLegendre(p, t1, t0 + t1, u);
491}
492
493void FuentesPyramid::CalcHomogenizedIntLegendre(int p,
494 real_t t0, real_t t1,
495 real_t *u,
496 real_t *dudt0, real_t *dudt1)
497{
498 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
499 CalcIntegratedLegendre(p, t1, t0+t1, u, dudt1, dudt0);
500 for (int i = 0; i <= p; i++) { dudt1[i] += dudt0[i]; }
501}
502
503void FuentesPyramid::CalcHomogenizedIntLegendre(int p,
504 real_t t0, real_t t1,
505 Vector &u)
506{
507 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
508 MFEM_ASSERT(u.Size() >= p+1, "Size of u is too small");
509 CalcHomogenizedIntLegendre(p, t0, t1, u.GetData());
510}
511
512void FuentesPyramid::CalcHomogenizedIntLegendre(int p,
513 real_t t0, real_t t1,
514 Vector &u,
515 Vector &dudt0, Vector &dudt1)
516{
517 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
518 MFEM_ASSERT(u.Size() >= p+1, "Size of u is too small");
519 MFEM_ASSERT(dudt0.Size() >= p+1, "Size of dudt0 is too small");
520 MFEM_ASSERT(dudt1.Size() >= p+1, "Size of dudt1 is too small");
521 CalcHomogenizedIntLegendre(p, t0, t1, u.GetData(),
522 dudt0.GetData(), dudt1.GetData());
523}
524
525void FuentesPyramid::CalcIntegratedJacobi(int p, real_t alpha,
526 real_t x, real_t t,
527 real_t *u)
528{
529 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
530 if (t > 0.0)
531 {
532 CalcScaledJacobi(p, alpha, x, t, u);
533 for (int i = p; i >= 2; i--)
534 {
535 real_t d0 = 2.0 * i + alpha;
536 real_t d1 = d0 - 1.0;
537 real_t d2 = d0 - 2.0;
538 real_t a = (alpha + i) / (d0 * d1);
539 real_t b = alpha / (d0 * d2);
540 real_t c = (real_t)(i - 1) / (d1 * d2);
541 u[i] = a * u[i] + b * t * u[i - 1] - c * t * t * u[i - 2];
542 }
543 if (p >= 1)
544 {
545 u[1] = x;
546 }
547 u[0] = 0.0;
548 }
549 else
550 {
551 u[0] = 1.0;
552 for (int i = 1; i <= p; i++)
553 {
554 u[i] = 0.0;
555 }
556 }
557}
558
559void FuentesPyramid::CalcIntegratedJacobi(int p, real_t alpha,
560 real_t x, real_t t,
561 real_t *u,
562 real_t *dudx,
563 real_t *dudt)
564{
565 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
566 CalcScaledJacobi(p, alpha, x, t, u, dudx, dudt);
567 for (int i = p; i >= 2; i--)
568 {
569 real_t d0 = 2.0 * i + alpha;
570 real_t d1 = d0 - 1.0;
571 real_t d2 = d0 - 2.0;
572 real_t a = (alpha + i) / (d0 * d1);
573 real_t b = alpha / (d0 * d2);
574 real_t c = (real_t)(i - 1) / (d1 * d2);
575 u[i] = a * u[i] + b * t * u[i - 1] - c * t * t * u[i - 2];
576 dudx[i] = a * dudx[i] + b * t * dudx[i - 1] - c * t * t * dudx[i - 2];
577 dudt[i] = a * dudt[i] + b * t * dudt[i - 1] + b * u[i - 1]
578 - c * t * t * dudt[i - 2] - 2.0 * c * t * u[i - 2];
579 }
580 if (p >= 1)
581 {
582 u[1] = x;
583 dudx[1] = 1.0;
584 dudt[1] = 0.0;
585 }
586 u[0] = 0.0;
587 dudx[0] = 0.0;
588 dudt[0] = 0.0;
589}
590
591void FuentesPyramid::CalcScaledJacobi(int p, real_t alpha,
592 real_t x, real_t t,
593 real_t *u)
594{
595 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
596
597 u[0] = 1.0;
598 if (p >= 1)
599 {
600 u[1] = (2.0 + alpha) * x - t;
601 }
602 for (int i = 2; i <= p; i++)
603 {
604 real_t a = 2.0 * i * (alpha + i) * (2.0 * i + alpha - 2.0);
605 real_t b = 2.0 * i + alpha - 1.0;
606 real_t c = (2.0 * i + alpha) * (2.0 * i + alpha - 2.0);
607 real_t d = 2.0 * (alpha + i - 1.0) * (i - 1) * (2.0 * i + alpha);
608 u[i] = (b * (c * (2.0 * x - t) + alpha * alpha * t) * u[i - 1]
609 - d * t * t * u[i - 2]) / a;
610 }
611}
612
613void FuentesPyramid::CalcScaledJacobi(int p, real_t alpha,
614 real_t x, real_t t,
615 real_t *u, real_t *dudx, real_t *dudt)
616{
617 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
618
619 u[0] = 1.0;
620 dudx[0] = 0.0;
621 dudt[0] = 0.0;
622 if (p >= 1)
623 {
624 u[1] = (2.0 + alpha) * x - t;
625 dudx[1] = 2.0 + alpha;
626 dudt[1] = -1.0;
627 }
628 for (int i = 2; i <= p; i++)
629 {
630 real_t a = 2.0 * i * (alpha + i) * (2.0 * i + alpha - 2.0);
631 real_t b = 2.0 * i + alpha - 1.0;
632 real_t c = (2.0 * i + alpha) * (2.0 * i + alpha - 2.0);
633 real_t d = 2.0 * (alpha + i - 1.0) * (i - 1) * (2.0 * i + alpha);
634 u[i] = (b * (c * (2.0 * x - t) + alpha * alpha * t) * u[i - 1]
635 - d * t * t * u[i - 2]) / a;
636 dudx[i] = (b * ((c * (2.0 * x - t) + alpha * alpha * t) * dudx[i - 1] +
637 2.0 * c * u[i - 1])
638 - d * t * t * dudx[i - 2]) / a;
639 dudt[i] = (b * ((c * (2.0 * x - t) + alpha * alpha * t) * dudt[i - 1] +
640 (alpha * alpha - c) * u[i - 1])
641 - d * t * t * dudt[i - 2] - 2.0 * d * t * u[i - 2]) / a;
642 }
643}
644
645void FuentesPyramid::CalcScaledJacobi(int p, real_t alpha,
646 real_t x, real_t t,
647 Vector &u)
648{
649 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
650 MFEM_ASSERT(u.Size() >= p+1, "Size of u is too small");
651 CalcScaledJacobi(p, alpha, x, t, u.GetData());
652}
653
654void FuentesPyramid::CalcScaledJacobi(int p, real_t alpha,
655 real_t x, real_t t,
656 Vector &u, Vector &dudx, Vector &dudt)
657{
658 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
659 MFEM_ASSERT(u.Size() >= p+1, "Size of u is too small");
660 MFEM_ASSERT(dudx.Size() >= p+1, "Size of dudx is too small");
661 MFEM_ASSERT(dudt.Size() >= p+1, "Size of dudt is too small");
662 CalcScaledJacobi(p, alpha, x, t, u.GetData(),
663 dudx.GetData(), dudt.GetData());
664}
665
666void FuentesPyramid::CalcHomogenizedScaJacobi(int p, real_t alpha,
667 real_t t0, real_t t1,
668 real_t *u,
669 real_t *dudt0, real_t *dudt1)
670{
671 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
672 CalcScaledJacobi(p, alpha, t1, t0+t1, u, dudt1, dudt0);
673 for (int i = 0; i <= p; i++) { dudt1[i] += dudt0[i]; }
674}
675
676void FuentesPyramid::CalcHomogenizedScaJacobi(int p, real_t alpha,
677 real_t t0, real_t t1,
678 Vector &u)
679{
680 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
681 MFEM_ASSERT(u.Size() >= p+1, "Size of u is too small");
682 CalcHomogenizedScaJacobi(p, alpha, t0, t1, u.GetData());
683}
684
685void FuentesPyramid::CalcHomogenizedScaJacobi(int p, real_t alpha,
686 real_t t0, real_t t1,
687 Vector &u,
688 Vector &dudt0, Vector &dudt1)
689{
690 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
691 MFEM_ASSERT(u.Size() >= p+1, "Size of u is too small");
692 MFEM_ASSERT(dudt0.Size() >= p+1, "Size of dudt0 is too small");
693 MFEM_ASSERT(dudt1.Size() >= p+1, "Size of dudt1 is too small");
694 CalcHomogenizedScaJacobi(p, alpha, t0, t1, u.GetData(),
695 dudt0.GetData(), dudt1.GetData());
696}
697
698void FuentesPyramid::CalcHomogenizedIntJacobi(int p, real_t alpha,
699 real_t t0, real_t t1,
700 real_t *u,
701 real_t *dudt0, real_t *dudt1)
702{
703 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
704 CalcIntegratedJacobi(p, alpha, t1, t0+t1, u, dudt1, dudt0);
705 for (int i = 0; i <= p; i++) { dudt1[i] += dudt0[i]; }
706}
707
708void FuentesPyramid::CalcHomogenizedIntJacobi(int p, real_t alpha,
709 real_t t0, real_t t1,
710 Vector &u)
711{
712 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
713 MFEM_ASSERT(u.Size() >= p+1, "Size of u is too small");
714 CalcHomogenizedIntJacobi(p, alpha, t0, t1, u.GetData());
715}
716
717void FuentesPyramid::CalcHomogenizedIntJacobi(int p, real_t alpha,
718 real_t t0, real_t t1,
719 Vector &u,
720 Vector &dudt0, Vector &dudt1)
721{
722 MFEM_ASSERT(p >= 0, "Polynomial order must be zero or larger");
723 MFEM_ASSERT(u.Size() >= p+1, "Size of u is too small");
724 MFEM_ASSERT(dudt0.Size() >= p+1, "Size of dudt0 is too small");
725 MFEM_ASSERT(dudt1.Size() >= p+1, "Size of dudt1 is too small");
726 CalcHomogenizedIntJacobi(p, alpha, t0, t1, u.GetData(),
727 dudt0.GetData(), dudt1.GetData());
728}
729
730void FuentesPyramid::phi_E(int p, real_t s0, real_t s1,
731 real_t *u)
732{
733 MFEM_ASSERT(p >= 2, "Polynomial order must be two or larger");
734 CalcHomogenizedIntLegendre(p, s0, s1, u);
735}
736
737void FuentesPyramid::phi_E(int p, real_t s0, real_t s1,
738 real_t *u, real_t *duds0, real_t *duds1)
739{
740 MFEM_ASSERT(p >= 2, "Polynomial order must be two or larger");
741 CalcHomogenizedIntLegendre(p, s0, s1, u, duds0, duds1);
742}
743
744void FuentesPyramid::phi_E(int p, Vector s, Vector &u)
745{
746 MFEM_ASSERT(p >= 2, "Polynomial order must be two or larger");
747 MFEM_ASSERT(u.Size() >= p+1, "Size of u is too small");
748 phi_E(p, s[0], s[1], u.GetData());
749}
750
751void FuentesPyramid::phi_E(int p, Vector s, Vector &u, DenseMatrix &duds)
752{
753 MFEM_ASSERT(p >= 2, "Polynomial order must be two or larger");
754 MFEM_ASSERT(s.Size() >= 2, "Size of s must be 2 or larger");
755 MFEM_ASSERT(u.Size() >= p+1, "Size of u is too small");
756 MFEM_ASSERT(duds.Height() >= p+1, "First dimension of duds is too small");
757 MFEM_ASSERT(duds.Width() >= 2,
758 "Second dimension of duds must be 2 or larger");
759 phi_E(p, s[0], s[1], u.GetData(), duds.GetColumn(0), duds.GetColumn(1));
760}
761
762void FuentesPyramid::phi_E(int p, Vector s, const DenseMatrix &grad_s,
763 Vector &u, DenseMatrix &grad_u) const
764{
765 MFEM_ASSERT(p >= 2, "Polynomial order must be two or larger");
766 MFEM_ASSERT(s.Size() >= 2, "Size of s must be 2 or larger");
767 MFEM_ASSERT(grad_s.Height() >= 2,
768 "First dimension of grad_s must be 2");
769 MFEM_ASSERT(grad_s.Width() >= 3,
770 "Second dimension of grad_s must be 3");
771 MFEM_ASSERT(u.Size() >= p+1, "Size of u is too small");
772 MFEM_ASSERT(grad_u.Height() >= p+1,
773 "First dimension of grad_u is too small");
774 MFEM_ASSERT(grad_u.Width() == grad_s.Width(),
775 "Second dimension of grad_u must match that of grad_s");
776#ifdef MFEM_THREAD_SAFE
777 DenseMatrix phi_E_mtmp;
778#endif
779 DenseMatrix &duds = phi_E_mtmp;
780 duds.SetSize(p + 1, grad_s.Height());
781
782 phi_E(p, s[0], s[1], u.GetData(), duds.GetColumn(0), duds.GetColumn(1));
783 Mult(duds, grad_s, grad_u);
784}
785
786void FuentesPyramid::phi_Q(int p, Vector s, Vector t, DenseMatrix &u) const
787{
788 MFEM_ASSERT(p >= 2, "Polynomial order must be two or larger");
789 MFEM_ASSERT(s.Size() >= 2, "Size of s must be 2 or larger");
790 MFEM_ASSERT(t.Size() >= 2, "Size of t must be 2 or larger");
791 MFEM_ASSERT(u.Height() >= p+1, "First dimension of u is too small");
792 MFEM_ASSERT(u.Width() >= p+1, "Second dimension of u is too small");
793
794#ifdef MFEM_THREAD_SAFE
795 Vector phi_Q_vtmp1;
796 Vector phi_Q_vtmp2;
797#endif
798 Vector &phi_E_i = phi_Q_vtmp1;
799 Vector &phi_E_j = phi_Q_vtmp2;
800
801 phi_E_i.SetSize(p+1);
802 phi_E(p, s, phi_E_i);
803
804 phi_E_j.SetSize(p+1);
805 phi_E(p, t, phi_E_j);
806
807 for (int j=0; j<=p; j++)
808 for (int i=0; i<=p; i++)
809 {
810 u(i,j) = phi_E_i[i] * phi_E_j[j];
811 }
812}
813
814void FuentesPyramid::phi_Q(int p, Vector s, const DenseMatrix &grad_s,
815 Vector t, const DenseMatrix &grad_t,
816 DenseMatrix &u, DenseTensor &grad_u) const
817{
818 MFEM_ASSERT(p >= 2, "Polynomial order must be two or larger");
819 MFEM_ASSERT(s.Size() >= 2, "Size of s must be 2 or larger");
820 MFEM_ASSERT(grad_s.Height() >= 2,
821 "First dimension of grad_s must be 2 or larger");
822 MFEM_ASSERT(grad_s.Width() >= 3,
823 "Second dimension of grad_s must be 3 or larger");
824 MFEM_ASSERT(t.Size() >= 2, "Size of t must be 2 or larger");
825 MFEM_ASSERT(grad_t.Height() >= 2,
826 "First dimension of grad_t must be 2 or larger");
827 MFEM_ASSERT(grad_t.Width() >= 3,
828 "Second dimension of grad_t must be 3 or larger");
829 MFEM_ASSERT(u.Height() >= p+1, "First dimension of u is too small");
830 MFEM_ASSERT(u.Width() >= p+1, "First dimension of u is too small");
831 MFEM_ASSERT(grad_u.SizeI() >= p+1,
832 "First dimension of grad_u is too small");
833 MFEM_ASSERT(grad_u.SizeJ() >= p+1,
834 "Second dimension of grad_u is too small");
835 MFEM_ASSERT(grad_u.SizeK() >= 3,
836 "Third dimension of grad_u must be 3 or larger");
837
838#ifdef MFEM_THREAD_SAFE
839 Vector phi_Q_vtmp1;
840 Vector phi_Q_vtmp2;
841 DenseMatrix phi_Q_mtmp1;
842 DenseMatrix phi_Q_mtmp2;
843#endif
844 Vector &phi_E_i = phi_Q_vtmp1;
845 Vector &phi_E_j = phi_Q_vtmp2;
846 DenseMatrix &dphi_E_i = phi_Q_mtmp1;
847 DenseMatrix &dphi_E_j = phi_Q_mtmp2;
848
849 phi_E_i.SetSize(p+1);
850 dphi_E_i.SetSize(p+1, grad_s.Width());
851 phi_E(p, s, grad_s, phi_E_i, dphi_E_i);
852
853 phi_E_j.SetSize(p+1);
854 dphi_E_j.SetSize(p+1, grad_t.Width());
855 phi_E(p, t, grad_t, phi_E_j, dphi_E_j);
856
857 for (int j=0; j<=p; j++)
858 for (int i=0; i<=p; i++)
859 {
860 u(i,j) = phi_E_i[i] * phi_E_j[j];
861
862 for (int k=0; k<3; k++)
863 grad_u(i,j,k) =
864 phi_E_i(i) * dphi_E_j(j,k) + dphi_E_i(i,k) * phi_E_j(j);
865 }
866}
867
868void FuentesPyramid::phi_T(int p, Vector s, DenseMatrix &u) const
869{
870 MFEM_ASSERT(p >= 3, "Polynomial order must be three or larger");
871 MFEM_ASSERT(s.Size() >= 3, "Size of s must be 3 or larger");
872 MFEM_ASSERT(u.Height() >= p, "First dimension of u is too small");
873 MFEM_ASSERT(u.Width() >= p-1, "Second dimension of u is too small");
874
875#ifdef MFEM_THREAD_SAFE
876 Vector phi_T_vtmp1;
877 Vector phi_T_vtmp2;
878#endif
879 Vector &phi_E_i = phi_T_vtmp1;
880 Vector &L_j = phi_T_vtmp2;
881
882 phi_E_i.SetSize(p);
883 phi_E(p-1, s, phi_E_i);
884
885 L_j.SetSize(p-1);
886
887 u = 0.0;
888 for (int i = 2; i < p; i++)
889 {
890 const real_t alpha = 2.0 * i;
891 CalcHomogenizedIntJacobi(p-2, alpha, s[0] + s[1], s[2], L_j);
892
893 for (int j = 1; i + j <= p; j++)
894 {
895 u(i,j) = phi_E_i[i] * L_j[j];
896 }
897 }
898}
899
900void FuentesPyramid::phi_T(int p, Vector s, const DenseMatrix &grad_s,
901 DenseMatrix &u, DenseTensor &grad_u) const
902{
903 MFEM_ASSERT(p >= 3, "Polynomial order must be three or larger");
904 MFEM_ASSERT(s.Size() >= 3, "Size of s must be 3 or larger");
905 MFEM_ASSERT(grad_s.Height() >= 3,
906 "First dimension of grad_s must be 2 or larger");
907 MFEM_ASSERT(grad_s.Width() >= 3,
908 "Second dimension of grad_s must be 3 or larger");
909 MFEM_ASSERT(u.Height() >= p, "First dimension of u is too small");
910 MFEM_ASSERT(u.Width() >= p-1, "Second dimension of u is too small");
911 MFEM_ASSERT(grad_u.SizeI() >= p,
912 "First dimension of grad_u is too small");
913 MFEM_ASSERT(grad_u.SizeJ() >= p-1,
914 "Second dimension of grad_u is too small");
915 MFEM_ASSERT(grad_u.SizeK() >= 3,
916 "Third dimension of grad_u must be 3 or larger");
917
918#ifdef MFEM_THREAD_SAFE
919 Vector phi_T_vtmp1;
920 Vector phi_T_vtmp2;
921 Vector phi_T_vtmp3;
922 Vector phi_T_vtmp4;
923 DenseMatrix phi_T_mtmp1;
924#endif
925 Vector &phi_E_i = phi_T_vtmp1;
926 DenseMatrix &dphi_E_i = phi_T_mtmp1;
927 Vector &L_j = phi_T_vtmp2;
928 Vector &dL_j_dx = phi_T_vtmp3;
929 Vector &dL_j_dt = phi_T_vtmp4;
930
931 phi_E_i.SetSize(p);
932 dphi_E_i.SetSize(p, 3);
933 phi_E(p-1, s, grad_s, phi_E_i, dphi_E_i);
934
935 L_j.SetSize(p-1);
936 dL_j_dx.SetSize(p-1);
937 dL_j_dt.SetSize(p-1);
938
939 u = 0.0;
940 grad_u = 0.0;
941 for (int i = 2; i < p; i++)
942 {
943 const real_t alpha = 2.0 * i;
944 CalcHomogenizedIntJacobi(p-2, alpha, s[0] + s[1], s[2], L_j,
945 dL_j_dx, dL_j_dt);
946
947 for (int j = 1; i + j <= p; j++)
948 {
949 u(i,j) = phi_E_i[i] * L_j[j];
950
951 for (int d=0; d<3; d++)
952 grad_u(i, j, d) = dphi_E_i(i, d) * L_j[j] +
953 phi_E_i[i] * (dL_j_dx[j] * (grad_s(0, d) +
954 grad_s(1, d)) +
955 dL_j_dt[j] * grad_s(2, d));
956 }
957 }
958}
959
960void FuentesPyramid::E_E(int p, Vector s, Vector sds, DenseMatrix &u) const
961{
962 MFEM_ASSERT(p >= 1, "Polynomial order must be one or larger");
963 MFEM_ASSERT(s.Size() >= 2, "Size of s must be 2 or larger");
964 MFEM_ASSERT(sds.Size() >= 3, "Size of sds must be 3 or larger");
965 MFEM_ASSERT(u.Height() >= p, "First dimension of u is too small");
966 MFEM_ASSERT(u.Width() >= 3, "Second dimension of u must be 3 or larger");
967
968#ifdef MFEM_THREAD_SAFE
969 Vector E_E_vtmp;
970#endif
971 Vector &P_i = E_E_vtmp;
972
973 P_i.SetSize(p);
974 CalcHomogenizedScaLegendre(p - 1, s[0], s[1], P_i);
975 for (int i=0; i<p; i++)
976 {
977 u(i,0) = P_i(i) * sds(0);
978 u(i,1) = P_i(i) * sds(1);
979 u(i,2) = P_i(i) * sds(2);
980 }
981}
982
983void FuentesPyramid::E_E(int p, Vector s, const DenseMatrix &grad_s,
984 DenseMatrix &u, DenseMatrix &curl_u) const
985{
986 MFEM_ASSERT(p >= 1, "Polynomial order must be one or larger");
987 MFEM_ASSERT(s.Size() >= 2, "Size of s must be 2 or larger");
988 MFEM_ASSERT(grad_s.Height() >= 2,
989 "First dimension of grad_s must be 2 or larger");
990 MFEM_ASSERT(grad_s.Width() >= 3,
991 "Second dimension of grad_s must be 3 or larger");
992 MFEM_ASSERT(u.Height() >= p, "First dimension of u is too small");
993 MFEM_ASSERT(u.Width() >= 3, "Second dimension of u must be 3 or larger");
994 MFEM_ASSERT(curl_u.Height() >= p,
995 "First dimension of curl_u is too small");
996 MFEM_ASSERT(curl_u.Width() >= 3,
997 "Second dimension of curl_u must be 3 or larger");
998
999#ifdef MFEM_THREAD_SAFE
1000 Vector E_E_vtmp;
1001#endif
1002 Vector &P_i = E_E_vtmp;
1003
1004 P_i.SetSize(p);
1005 CalcHomogenizedScaLegendre(p - 1, s[0], s[1], P_i);
1006
1007 Vector grad_s0({grad_s(0,0), grad_s(0,1), grad_s(0,2)});
1008 Vector grad_s1({grad_s(1,0), grad_s(1,1), grad_s(1,2)});
1009 Vector sds(3);
1010 add(s(0), grad_s1, -s(1), grad_s0, sds);
1011
1012 Vector dsxds(3);
1013 grad_s0.cross3D(grad_s1, dsxds);
1014
1015 for (int i=0; i<p; i++)
1016 {
1017 u(i,0) = P_i(i) * sds(0);
1018 u(i,1) = P_i(i) * sds(1);
1019 u(i,2) = P_i(i) * sds(2);
1020
1021 curl_u(i, 0) = (i + 2) * P_i(i) * dsxds(0);
1022 curl_u(i, 1) = (i + 2) * P_i(i) * dsxds(1);
1023 curl_u(i, 2) = (i + 2) * P_i(i) * dsxds(2);
1024 }
1025}
1026
1027void FuentesPyramid::E_Q(int p, Vector s, Vector sds, Vector t,
1028 DenseTensor &u) const
1029{
1030 MFEM_ASSERT(p >= 2, "Polynomial order must be two or larger");
1031 MFEM_ASSERT(s.Size() >= 2, "Size of s must be 2 or larger");
1032 MFEM_ASSERT(sds.Size() >= 3, "Size of sds must be 3 or larger");
1033 MFEM_ASSERT(t.Size() >= 2, "Size of t must be 2 or larger");
1034 MFEM_ASSERT(u.SizeI() >= p, "First dimension of u is too small");
1035 MFEM_ASSERT(u.SizeJ() >= p+1, "Second dimension of u is too small");
1036 MFEM_ASSERT(u.SizeK() >= 3, "Third dimension of u must be 3 or larger");
1037
1038#ifdef MFEM_THREAD_SAFE
1039 Vector E_Q_vtmp;
1040 DenseMatrix E_Q_mtmp1;
1041#endif
1042
1043 DenseMatrix &E_E_i = E_Q_mtmp1;
1044 Vector &phi_E_j = E_Q_vtmp;
1045
1046 E_E_i.SetSize(p, 3);
1047 E_E(p, s, sds, E_E_i);
1048
1049 phi_E_j.SetSize(p + 1);
1050 phi_E(p, t, phi_E_j);
1051
1052 for (int k=0; k<3; k++)
1053 {
1054 u(k).SetCol(0, 0.0);
1055 u(k).SetCol(1, 0.0);
1056 }
1057
1058 for (int j=2; j<=p; j++)
1059 for (int i=0; i<p; i++)
1060 for (int k=0; k<3; k++)
1061 {
1062 u(i, j, k) = phi_E_j(j) * E_E_i(i, k);
1063 }
1064}
1065
1066void FuentesPyramid::E_Q(int p, Vector s, const DenseMatrix &grad_s,
1067 Vector t, const DenseMatrix &grad_t,
1068 DenseTensor &u, DenseTensor &curl_u) const
1069{
1070 MFEM_ASSERT(p >= 2, "Polynomial order must be two or larger");
1071 MFEM_ASSERT(s.Size() >= 2, "Size of s must be 2 or larger");
1072 MFEM_ASSERT(grad_s.Height() >= 2,
1073 "First dimension of grad_s must be 2");
1074 MFEM_ASSERT(grad_s.Width() >= 3,
1075 "Second dimension of grad_s must be 3");
1076 MFEM_ASSERT(t.Size() >= 2, "Size of t must be 2 or larger");
1077 MFEM_ASSERT(grad_t.Height() >= 2,
1078 "First dimension of grad_t must be 2");
1079 MFEM_ASSERT(grad_t.Width() >= 3,
1080 "Second dimension of grad_t must be 3");
1081 MFEM_ASSERT(u.SizeI() >= p, "First dimension of u is too small");
1082 MFEM_ASSERT(u.SizeJ() >= p+1, "Second dimension of u is too small");
1083 MFEM_ASSERT(u.SizeK() >= 3, "Third dimension of u must be 3 or larger");
1084 MFEM_ASSERT(curl_u.SizeI() >= p, "First dimension of curl_u is too small");
1085 MFEM_ASSERT(curl_u.SizeJ() >= p+1,
1086 "Second dimension of curl_u is too small");
1087 MFEM_ASSERT(curl_u.SizeK() >= 3,
1088 "Third dimension of curl_u must be 3 or larger");
1089
1090#ifdef MFEM_THREAD_SAFE
1091 Vector E_Q_vtmp;
1092 DenseMatrix E_Q_mtmp1;
1093 DenseMatrix E_Q_mtmp2;
1094 DenseMatrix E_Q_mtmp3;
1095#endif
1096 Vector &phi_E_j = E_Q_vtmp;
1097 DenseMatrix &dphi_E_j = E_Q_mtmp1;
1098 DenseMatrix &E_E_i = E_Q_mtmp2;
1099 DenseMatrix &dE_E_i = E_Q_mtmp3;
1100
1101 phi_E_j.SetSize(p + 1);
1102 dphi_E_j.SetSize(p + 1, grad_t.Width());
1103 phi_E(p, t, grad_t, phi_E_j, dphi_E_j);
1104
1105 E_E_i.SetSize(p, 3);
1106 dE_E_i.SetSize(p, 3);
1107 E_E(p, s, grad_s, E_E_i, dE_E_i);
1108
1109 for (int k=0; k<3; k++)
1110 {
1111 u(k).SetCol(0, 0.0);
1112 u(k).SetCol(1, 0.0);
1113 curl_u(k).SetCol(0, 0.0);
1114 curl_u(k).SetCol(1, 0.0);
1115 }
1116
1117 for (int j=2; j<=p; j++)
1118 for (int i=0; i<p; i++)
1119 {
1120 for (int k=0; k<3; k++)
1121 {
1122 u(i, j, k) = phi_E_j(j) * E_E_i(i, k);
1123 }
1124
1125 curl_u(i, j, 0) = phi_E_j(j) * dE_E_i(i, 0)
1126 + dphi_E_j(j, 1) * E_E_i(i, 2)
1127 - dphi_E_j(j, 2) * E_E_i(i, 1);
1128 curl_u(i, j, 1) = phi_E_j(j) * dE_E_i(i, 1)
1129 + dphi_E_j(j, 2) * E_E_i(i, 0)
1130 - dphi_E_j(j, 0) * E_E_i(i, 2);
1131 curl_u(i, j, 2) = phi_E_j(j) * dE_E_i(i, 2)
1132 + dphi_E_j(j, 0) * E_E_i(i, 1)
1133 - dphi_E_j(j, 1) * E_E_i(i, 0);
1134 }
1135}
1136
1137void FuentesPyramid::E_T(int p, Vector s, Vector sds, DenseTensor &u) const
1138{
1139 MFEM_ASSERT(p >= 2, "Polynomial order must be two or larger");
1140 MFEM_ASSERT(s.Size() >= 3, "Size of s must be 3 or larger");
1141 MFEM_ASSERT(sds.Size() >= 3, "Size of sds must be 3 or larger");
1142 MFEM_ASSERT(u.SizeI() >= p - 1, "First dimension of u is too small");
1143 MFEM_ASSERT(u.SizeJ() >= p, "Second dimension of u is too small");
1144 MFEM_ASSERT(u.SizeK() >= 3, "Third dimension of u must be 3 or larger");
1145
1146#ifdef MFEM_THREAD_SAFE
1147 Vector E_T_vtmp1;
1148 DenseMatrix E_T_mtmp1;
1149#endif
1150 Vector &L_j = E_T_vtmp1;
1151 DenseMatrix &E_E_i = E_T_mtmp1;
1152
1153 E_E_i.SetSize(p - 1, 3);
1154 E_E(p - 1, s, sds, E_E_i);
1155
1156 L_j.SetSize(p);
1157 for (int i=0; i<p-1; i++)
1158 {
1159 const real_t alpha = 2.0 * i + 1.0;
1160 CalcHomogenizedIntJacobi(p - 1, alpha, s[0] + s[1], s[2], L_j);
1161
1162 u(i, 0, 0) = 0.0; u(i, 0, 1) = 0.0; u(i, 0, 2) = 0.0;
1163 for (int j=1; i+j<p; j++)
1164 for (int k=0; k<3; k++)
1165 {
1166 u(i, j, k) = L_j(j) * E_E_i(i, k);
1167 }
1168 }
1169}
1170
1171void FuentesPyramid::E_T(int p, Vector s, const DenseMatrix & grad_s,
1172 DenseTensor &u, DenseTensor &curl_u) const
1173{
1174 MFEM_ASSERT(p >= 2, "Polynomial order must be two or larger");
1175 MFEM_ASSERT(s.Size() >= 3, "Size of s must be 3 or larger");
1176 MFEM_ASSERT(grad_s.Height() >= 3,
1177 "First dimension of grad_s must be 3");
1178 MFEM_ASSERT(grad_s.Width() >= 3,
1179 "Second dimension of grad_s must be 3");
1180 MFEM_ASSERT(u.SizeI() >= p - 1, "First dimension of u is too small");
1181 MFEM_ASSERT(u.SizeJ() >= p, "Second dimension of u is too small");
1182 MFEM_ASSERT(u.SizeK() >= 3, "Third dimension of u must be 3 or larger");
1183 MFEM_ASSERT(curl_u.SizeI() >= p - 1,
1184 "First dimension of curl_u is too small");
1185 MFEM_ASSERT(curl_u.SizeJ() >= p,
1186 "Second dimension of curl_u is too small");
1187 MFEM_ASSERT(curl_u.SizeK() >= 3,
1188 "Third dimension of curl_u must be 3 or larger");
1189
1190#ifdef MFEM_THREAD_SAFE
1191 Vector E_T_vtmp1;
1192 Vector E_T_vtmp2;
1193 Vector E_T_vtmp3;
1194 DenseMatrix E_T_mtmp1;
1195 DenseMatrix E_T_mtmp2;
1196#endif
1197 Vector &L_j = E_T_vtmp1;
1198 Vector &dL_j_dx = E_T_vtmp2;
1199 Vector &dL_j_dt = E_T_vtmp3;
1200 DenseMatrix & E_E_i = E_T_mtmp1;
1201 DenseMatrix &dE_E_i = E_T_mtmp2;
1202
1203 Vector dL(3), grad_L(3);
1204
1205 E_E_i.SetSize(p - 1, 3);
1206 dE_E_i.SetSize(p - 1, 3);
1207 E_E(p - 1, s, grad_s, E_E_i, dE_E_i);
1208
1209 L_j.SetSize(p);
1210 dL_j_dx.SetSize(p);
1211 dL_j_dt.SetSize(p);
1212 for (int i=0; i<p-1; i++)
1213 {
1214 const real_t alpha = 2.0 * i + 1.0;
1215 CalcHomogenizedIntJacobi(p - 1, alpha, s[0] + s[1], s[2], L_j,
1216 dL_j_dx, dL_j_dt);
1217
1218 u(i, 0, 0) = 0.0; u(i, 0, 1) = 0.0; u(i, 0, 2) = 0.0;
1219 curl_u(i, 0, 0) = 0.0; curl_u(i, 0, 1) = 0.0; curl_u(i, 0, 2) = 0.0;
1220 for (int j=1; i+j<p; j++)
1221 {
1222 dL(0) = dL_j_dx(j); dL(1) = dL_j_dx(j); dL(2) = dL_j_dt(j);
1223
1224 grad_s.MultTranspose(dL, grad_L);
1225
1226 for (int k=0; k<3; k++)
1227 {
1228 u(i, j, k) = L_j(j) * E_E_i(i, k);
1229 curl_u(i, j, k) = L_j(j) * dE_E_i(i, k);
1230 }
1231 curl_u(i, j, 0) += grad_L(1) * E_E_i(i, 2) - grad_L(2) * E_E_i(i, 1);
1232 curl_u(i, j, 1) += grad_L(2) * E_E_i(i, 0) - grad_L(0) * E_E_i(i, 2);
1233 curl_u(i, j, 2) += grad_L(0) * E_E_i(i, 1) - grad_L(1) * E_E_i(i, 0);
1234 }
1235 }
1236}
1237
1238void FuentesPyramid::V_Q(int p, Vector s, Vector sds,
1239 Vector t, Vector tdt, DenseTensor &u) const
1240{
1241 MFEM_ASSERT(p >= 1, "Polynomial order must be one or larger");
1242 MFEM_ASSERT(s.Size() >= 2, "Size of s must be 2 or larger");
1243 MFEM_ASSERT(sds.Size() >= 3, "Size of sds must be 3 or larger");
1244 MFEM_ASSERT(t.Size() >= 2, "Size of t must be 2 or larger");
1245 MFEM_ASSERT(tdt.Size() >= 3, "Size of tdt must be 3 or larger");
1246 MFEM_ASSERT(u.SizeI() >= p, "First dimension of u is too small");
1247 MFEM_ASSERT(u.SizeJ() >= p, "Second dimension of u is too small");
1248 MFEM_ASSERT(u.SizeK() >= 3, "Third dimension of u must be 3 or larger");
1249
1250#ifdef MFEM_THREAD_SAFE
1251 DenseMatrix V_Q_mtmp1;
1252 DenseMatrix V_Q_mtmp2;
1253#endif
1254 DenseMatrix &E_E_i = V_Q_mtmp1;
1255 DenseMatrix &E_E_j = V_Q_mtmp2;
1256
1257 E_E_i.SetSize(p, 3);
1258 E_E(p, s, sds, E_E_i);
1259
1260 E_E_j.SetSize(p, 3);
1261 E_E(p, t, tdt, E_E_j);
1262
1263 for (int j=0; j<p; j++)
1264 for (int i=0; i<p; i++)
1265 {
1266 u(i, j, 0) = E_E_i(i, 1) * E_E_j(j, 2) - E_E_i(i, 2) * E_E_j(j, 1);
1267 u(i, j, 1) = E_E_i(i, 2) * E_E_j(j, 0) - E_E_i(i, 0) * E_E_j(j, 2);
1268 u(i, j, 2) = E_E_i(i, 0) * E_E_j(j, 1) - E_E_i(i, 1) * E_E_j(j, 0);
1269 }
1270}
1271
1272void FuentesPyramid::V_T(int p, Vector s, Vector sdsxds, DenseTensor &u) const
1273{
1274 MFEM_ASSERT(p >= 1, "Polynomial order must be one or larger");
1275 MFEM_ASSERT(s.Size() >= 2, "Size of s must be 2 or larger");
1276 MFEM_ASSERT(sdsxds.Size() >= 3, "Size of sdsxds must be 3 or larger");
1277 MFEM_ASSERT(u.SizeI() >= p, "First dimension of u is too small");
1278 MFEM_ASSERT(u.SizeJ() >= p, "Second dimension of u is too small");
1279 MFEM_ASSERT(u.SizeK() >= 3, "Third dimension of u must be 3 or larger");
1280
1281#ifdef MFEM_THREAD_SAFE
1282 Vector V_T_vtmp1;
1283 Vector V_T_vtmp2;
1284#endif
1285 Vector &P_i = V_T_vtmp1;
1286 Vector &P_j = V_T_vtmp2;
1287
1288 P_i.SetSize(p);
1289 CalcHomogenizedScaLegendre(p-1, s[0], s[1], P_i);
1290
1291 P_j.SetSize(p);
1292 for (int i=0; i<p; i++)
1293 {
1294 const real_t alpha = 2.0 * i + 1.0;
1295 CalcHomogenizedScaJacobi(p-1, alpha, s[0] + s[1], s[2], P_j);
1296 for (int j=0; i + j < p; j++)
1297 {
1298 const real_t vij = P_i(i) * P_j(j);
1299 u(i,j,0) = vij * sdsxds(0);
1300 u(i,j,1) = vij * sdsxds(1);
1301 u(i,j,2) = vij * sdsxds(2);
1302 }
1303 }
1304}
1305
1306void FuentesPyramid::V_T(int p, Vector s, Vector sdsxds, real_t dsdsxds,
1307 DenseTensor &u, DenseMatrix &du) const
1308{
1309 MFEM_ASSERT(p >= 1, "Polynomial order must be one or larger");
1310 MFEM_ASSERT(s.Size() >= 2, "Size of s must be 2 or larger");
1311 MFEM_ASSERT(sdsxds.Size() >= 3, "Size of sdsxds must be 3 or larger");
1312 MFEM_ASSERT(u.SizeI() >= p, "First dimension of u is too small");
1313 MFEM_ASSERT(u.SizeJ() >= p, "Second dimension of u is too small");
1314 MFEM_ASSERT(u.SizeK() >= 3, "Third dimension of u must be 3 or larger");
1315 MFEM_ASSERT(du.Height() >= p, "First dimension of du is too small");
1316 MFEM_ASSERT(du.Width() >= p, "Second dimension of du is too small");
1317
1318#ifdef MFEM_THREAD_SAFE
1319 Vector V_T_vtmp1;
1320 Vector V_T_vtmp2;
1321#endif
1322 Vector &P_i = V_T_vtmp1;
1323 Vector &P_j = V_T_vtmp2;
1324
1325 P_i.SetSize(p);
1326 CalcHomogenizedScaLegendre(p-1, s[0], s[1], P_i);
1327
1328 P_j.SetSize(p);
1329 for (int i=0; i<p; i++)
1330 {
1331 const real_t alpha = 2.0 * i + 1.0;
1332 CalcHomogenizedScaJacobi(p-1, alpha, s[0] + s[1], s[2], P_j);
1333 for (int j=0; i + j < p; j++)
1334 {
1335 const real_t vij = P_i(i) * P_j(j);
1336 u(i,j,0) = vij * sdsxds(0);
1337 u(i,j,1) = vij * sdsxds(1);
1338 u(i,j,2) = vij * sdsxds(2);
1339
1340 du(i,j) = (i+j+3) * vij * dsdsxds;
1341 }
1342 }
1343}
1344
1345void FuentesPyramid::VT_T(int p, Vector s, Vector sds, Vector sdsxds,
1346 real_t mu, Vector grad_mu, DenseTensor &u) const
1347{
1348 MFEM_ASSERT(p >= 1, "Polynomial order must be one or larger");
1349 MFEM_ASSERT(s.Size() >= 2, "Size of s must be 2 or larger");
1350 MFEM_ASSERT(sds.Size() >= 3, "Size of sds must be 3 or larger");
1351 MFEM_ASSERT(sdsxds.Size() >= 3, "Size of sdsxds must be 3 or larger");
1352 MFEM_ASSERT(grad_mu.Size() >= 3, "Size of grad_mu must be 3 or larger");
1353 MFEM_ASSERT(u.SizeI() >= p, "First dimension of u is too small");
1354 MFEM_ASSERT(u.SizeJ() >= p, "Second dimension of u is too small");
1355 MFEM_ASSERT(u.SizeK() >= 3, "Third dimension of u must be 3 or larger");
1356
1357#ifdef MFEM_THREAD_SAFE
1358 Vector VT_T_vtmp1;
1359 Vector VT_T_vtmp2;
1360 DenseMatrix VT_T_mtmp1;
1361 DenseTensor VT_T_ttmp1;
1362#endif
1363
1364 Vector ms({mu * s(0), mu * s(1), s(2)});
1365 Vector s2(s.GetData(), 2);
1366
1367 Vector &P_i = VT_T_vtmp1;
1368 P_i.SetSize(p);
1369 CalcHomogenizedScaLegendre(p-1, ms[0], ms[1], P_i);
1370
1371 DenseMatrix &EE0 = VT_T_mtmp1;
1372 EE0.SetSize(1,3);
1373 Vector EE(EE0.GetData(), 3);
1374 E_E(1, s2, sds, EE0);
1375
1376 Vector dmuxEE(3);
1377 grad_mu.cross3D(EE, dmuxEE);
1378
1379 DenseTensor &VT00 = VT_T_ttmp1;
1380 VT00.SetSize(1,1,3);
1381 V_T(1, s, sdsxds, VT00);
1382
1383 Vector &J_j = VT_T_vtmp2;
1384 J_j.SetSize(p);
1385
1386 u = 0.0;
1387
1388 for (int i=0; i<p; i++)
1389 {
1390 CalcHomogenizedScaJacobi(p-i-1, 2*i+1, ms[0] + ms[1], ms[2], J_j);
1391 for (int j=0; i+j<p; j++)
1392 for (int k=0; k<3; k++)
1393 u(i, j, k) = P_i(i) * J_j(j) *
1394 (mu * VT00(0,0,k) + s(2) * dmuxEE(k));
1395 }
1396}
1397
1398void FuentesPyramid::VT_T(int p, Vector s, Vector sds, Vector sdsxds,
1399 Vector grad_s2, real_t mu, Vector grad_mu,
1400 DenseTensor &u, DenseMatrix &du) const
1401{
1402 MFEM_ASSERT(p >= 1, "Polynomial order must be one or larger");
1403 MFEM_ASSERT(s.Size() >= 2, "Size of s must be 2 or larger");
1404 MFEM_ASSERT(sds.Size() >= 3, "Size of sds must be 3 or larger");
1405 MFEM_ASSERT(sdsxds.Size() >= 3, "Size of sdsxds must be 3 or larger");
1406 MFEM_ASSERT(grad_s2.Size() >= 3, "Size of grad_s2 must be 3 or larger");
1407 MFEM_ASSERT(grad_mu.Size() >= 3, "Size of grad_mu must be 3 or larger");
1408 MFEM_ASSERT(u.SizeI() >= p, "First dimension of u is too small");
1409 MFEM_ASSERT(u.SizeJ() >= p, "Second dimension of u is too small");
1410 MFEM_ASSERT(u.SizeK() >= 3, "Third dimension of u must be 3 or larger");
1411 MFEM_ASSERT(du.Height() >= p, "First dimension of du is too small");
1412 MFEM_ASSERT(du.Width() >= p, "Second dimension of du is too small");
1413
1414#ifdef MFEM_THREAD_SAFE
1415 Vector VT_T_vtmp1;
1416 Vector VT_T_vtmp2;
1417 DenseMatrix VT_T_mtmp1;
1418 DenseTensor VT_T_ttmp1;
1419#endif
1420
1421 Vector ms({mu * s(0), mu * s(1), s(2)});
1422 Vector s2(s.GetData(), 2);
1423
1424 Vector &P_i = VT_T_vtmp1;
1425 P_i.SetSize(p);
1426 CalcHomogenizedScaLegendre(p-1, ms[0], ms[1], P_i);
1427
1428 DenseMatrix &EE0 = VT_T_mtmp1;
1429 EE0.SetSize(1,3);
1430 Vector EE(EE0.GetData(), 3);
1431 E_E(1, s2, sds, EE0);
1432
1433 Vector dmuxEE(3);
1434 grad_mu.cross3D(EE, dmuxEE);
1435
1436 Vector EExds2(3);
1437 EE.cross3D(grad_s2, EExds2);
1438
1439 DenseTensor &VT00 = VT_T_ttmp1;
1440 VT00.SetSize(1,1,3);
1441 V_T(1, s, sdsxds, VT00);
1442
1443 Vector &J_j = VT_T_vtmp2;
1444 J_j.SetSize(p);
1445
1446 Vector EV(3);
1447
1448 u = 0.0;
1449 du = 0.0;
1450
1451 for (int i=0; i<p; i++)
1452 {
1453 CalcHomogenizedScaJacobi(p-i-1, 2*i+1, ms[0] + ms[1], ms[2], J_j);
1454 for (int j=0; i+j<p; j++)
1455 {
1456 for (int k=0; k<3; k++)
1457 {
1458 u(i, j, k) = P_i(i) * J_j(j) *
1459 (mu * VT00(0, 0, k) + s(2) * dmuxEE(k));
1460
1461 EV(k) = (i+j+3) * EExds2(k) - VT00(0, 0, k);
1462 }
1463
1464 du(i, j) = P_i(i) * J_j(j) * (grad_mu * EV);
1465 }
1466 }
1467}
1468
1469void FuentesPyramid::V_L(int p, Vector sx, const DenseMatrix &grad_sx,
1470 Vector sy, const DenseMatrix &grad_sy,
1471 real_t t, Vector grad_t,
1472 DenseTensor &u) const
1473{
1474 MFEM_ASSERT(p >= 2, "Polynomial order must be two or larger");
1475 MFEM_ASSERT(sx.Size() >= 2, "Size of sx must be 2 or larger");
1476 MFEM_ASSERT(grad_sx.Height() >= 2,
1477 "First dimension of grad_sx must be 2 or larger");
1478 MFEM_ASSERT(grad_sx.Width() >= 3,
1479 "Second dimension of grad_sx must be 3 or larger");
1480 MFEM_ASSERT(sy.Size() >= 2, "Size of sy must be 2 or larger");
1481 MFEM_ASSERT(grad_sy.Height() >= 2,
1482 "First dimension of grad_sy must be 2 or larger");
1483 MFEM_ASSERT(grad_sy.Width() >= 3,
1484 "Second dimension of grad_sy must be 3 or larger");
1485 MFEM_ASSERT(grad_t.Size() >= 3, "Size of grad_t must be 3 or larger");
1486 MFEM_ASSERT(u.SizeI() >= p+1, "First dimension of u is too small");
1487 MFEM_ASSERT(u.SizeJ() >= p+1, "Second dimension of u is too small");
1488 MFEM_ASSERT(u.SizeK() >= 3, "Third dimension of u must be 3 or larger");
1489
1490#ifdef MFEM_THREAD_SAFE
1491 Vector V_L_vtmp1;
1492 Vector V_L_vtmp2;
1493 DenseMatrix V_L_mtmp1;
1494 DenseMatrix V_L_mtmp2;
1495#endif
1496 Vector &phi_E_i = V_L_vtmp1;
1497 Vector &phi_E_j = V_L_vtmp2;
1498 DenseMatrix &dphi_E_i = V_L_mtmp1;
1499 DenseMatrix &dphi_E_j = V_L_mtmp2;
1500
1501 Vector grad_t3(grad_t.GetData(), 3);
1502
1503 phi_E_i.SetSize(p+1);
1504 dphi_E_i.SetSize(p+1, grad_sx.Width());
1505 phi_E(p, sx, grad_sx, phi_E_i, dphi_E_i);
1506
1507 phi_E_j.SetSize(p+1);
1508 dphi_E_j.SetSize(p+1, grad_sy.Width());
1509 phi_E(p, sy, grad_sy, phi_E_j, dphi_E_j);
1510
1511 Vector dphii(3);
1512 Vector dphij(3);
1513 Vector dphidphi(3);
1514 Vector phidphi(3);
1515 Vector dtphidphi(3);
1516
1517 for (int j=2; j<=p; j++)
1518 {
1519 for (int l=0; l<3; l++) { dphij[l] = dphi_E_j(j, l); }
1520 for (int i=2; i<=p; i++)
1521 {
1522 for (int l=0; l<3; l++) { dphii[l] = dphi_E_i(i, l); }
1523
1524 dphii.cross3D(dphij, dphidphi);
1525
1526 add(phi_E_i(i), dphij, -phi_E_j(j), dphii, phidphi);
1527
1528 grad_t3.cross3D(phidphi, dtphidphi);
1529
1530 for (int l=0; l<3; l++)
1531 {
1532 u(i, j, l) = t * (t * dphidphi(l) + dtphidphi(l));
1533 }
1534 }
1535 }
1536}
1537
1538void FuentesPyramid::V_R(int p, Vector s, const DenseMatrix &grad_s,
1539 real_t mu, Vector dmu, real_t t, Vector dt,
1540 DenseMatrix &u) const
1541{
1542 MFEM_ASSERT(p >= 2, "Polynomial order must be two or larger");
1543 MFEM_ASSERT(s.Size() >= 2, "Size of s must be 2 or larger");
1544 MFEM_ASSERT(grad_s.Height() >= 2,
1545 "First dimension of grad_s must be 2");
1546 MFEM_ASSERT(grad_s.Width() >= 3,
1547 "Second dimension of grad_s must be 3");
1548 MFEM_ASSERT(dmu.Size() >= 3, "Size of dmu must be 3 or larger");
1549 MFEM_ASSERT(dt.Size() >= 3, "Size of dt must be 3 or larger");
1550 MFEM_ASSERT(u.Height() >= p+1, "First dimension of u is too small");
1551 MFEM_ASSERT(u.Width() >= 3, "Second dimension of u is too small");
1552
1553#ifdef MFEM_THREAD_SAFE
1554 Vector V_R_vtmp;
1555 DenseMatrix V_R_mtmp;
1556#endif
1557 Vector &phi_E_i = V_R_vtmp;
1558 DenseMatrix &dphi_E_i = V_R_mtmp;
1559
1560 phi_E_i.SetSize(p+1);
1561 dphi_E_i.SetSize(p+1, grad_s.Width());
1562 phi_E(p, s, grad_s, phi_E_i, dphi_E_i);
1563
1564 u.SetRow(0, 0.0);
1565 u.SetRow(1, 0.0);
1566
1567 Vector dmu3(dmu.GetData(), 3);
1568 Vector dt3(dt.GetData(), 3);
1569 Vector dphit2(3);
1570 Vector dphixdmu(3);
1571 Vector dphi(3);
1572
1573 for (int i=2; i<=p; i++)
1574 {
1575 // dphi_E_i.GetRow(i, dphi);
1576 for (int l=0; l<3; l++) { dphi[l] = dphi_E_i(i, l); }
1577 add(t * t, dphi, 2.0 * t * phi_E_i(i), dt3, dphit2);
1578 dphit2.cross3D(dmu3, dphixdmu);
1579 // u.SetRow(i, dphixdmu);
1580 for (int l=0; l<3; l++) { u(i, l) = dphixdmu(l); }
1581 }
1582}
1583
1584} // namespace mfem
mfem::DenseMatrix
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
mfem::DenseMatrix::MultTranspose
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
Definition densemat.cpp:172
mfem::DenseMatrix::SetRow
void SetRow(int r, const real_t *row)
Definition densemat.cpp:2208
mfem::DenseMatrix::GetData
real_t * GetData() const
Returns the matrix data array.
Definition densemat.hpp:118
mfem::DenseMatrix::SetSize
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:108
mfem::DenseMatrix::GetColumn
void GetColumn(int c, Vector &col) const
Definition densemat.cpp:1361
mfem::DenseTensor
Rank 3 tensor (array of matrices)
Definition densemat.hpp:1112
mfem::DenseTensor::SetSize
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
Definition densemat.hpp:1163
mfem::DenseTensor::SizeJ
int SizeJ() const
Definition densemat.hpp:1158
mfem::DenseTensor::SizeI
int SizeI() const
Definition densemat.hpp:1157
mfem::DenseTensor::SizeK
int SizeK() const
Definition densemat.hpp:1159
mfem::FuentesPyramid::grad_mu01
static DenseMatrix grad_mu01(real_t z)
Definition fe_pyramid.cpp:178
mfem::FuentesPyramid::lam45_grad_lam45
static Vector lam45_grad_lam45(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:110
mfem::FuentesPyramid::phi_Q
void phi_Q(int p, Vector s, Vector t, DenseMatrix &u) const
Definition fe_pyramid.cpp:786
mfem::FuentesPyramid::CalcHomogenizedIntJacobi
static void CalcHomogenizedIntJacobi(int p, real_t alpha, real_t t0, real_t t1, real_t *u)
Definition fe_pyramid.hpp:390
mfem::FuentesPyramid::nu12
static Vector nu12(real_t z, Vector xy, unsigned int ab)
Definition fe_pyramid.hpp:207
mfem::FuentesPyramid::V_R
void V_R(int p, Vector s, const DenseMatrix &grad_s, real_t mu, Vector grad_mu, real_t t, Vector grad_t, DenseMatrix &u) const
Definition fe_pyramid.cpp:1538
mfem::FuentesPyramid::lam125_grad_lam125
static Vector lam125_grad_lam125(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:118
mfem::FuentesPyramid::lam145_grad_lam145
static Vector lam145_grad_lam145(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:153
mfem::FuentesPyramid::lam45
static Vector lam45(real_t x, real_t y, real_t z)
Definition fe_pyramid.hpp:115
mfem::FuentesPyramid::phi_E
static void phi_E(int p, real_t s0, real_t s1, real_t *u)
Definition fe_pyramid.cpp:730
mfem::FuentesPyramid::nu012
static Vector nu012(real_t z, Vector xy, unsigned int ab)
Definition fe_pyramid.hpp:209
mfem::FuentesPyramid::grad_nu120
static DenseMatrix grad_nu120(real_t z, Vector xy, unsigned int ab)
Definition fe_pyramid.cpp:250
mfem::FuentesPyramid::grad_lam25
static DenseMatrix grad_lam25(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:62
mfem::FuentesPyramid::E_E
void E_E(int p, Vector s, Vector sds, DenseMatrix &u) const
Definition fe_pyramid.cpp:960
mfem::FuentesPyramid::CalcIntegratedLegendre
static void CalcIntegratedLegendre(int p, real_t x, real_t t, real_t *u)
Integrated Legendre Polynomials.
Definition fe_pyramid.cpp:369
mfem::FuentesPyramid::nu01_grad_nu01
static Vector nu01_grad_nu01(real_t z, Vector xy, unsigned int ab)
Definition fe_pyramid.cpp:259
mfem::FuentesPyramid::lam35
static Vector lam35(real_t x, real_t y, real_t z)
Definition fe_pyramid.hpp:113
mfem::FuentesPyramid::grad_lam5
static Vector grad_lam5(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:49
mfem::FuentesPyramid::div_lam345_grad_lam345
static real_t div_lam345_grad_lam345(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:166
mfem::FuentesPyramid::lam25
static Vector lam25(real_t x, real_t y, real_t z)
Definition fe_pyramid.hpp:111
mfem::FuentesPyramid::nu01
static Vector nu01(real_t z, Vector xy, unsigned int ab)
Definition fe_pyramid.hpp:205
mfem::FuentesPyramid::VT_T
void VT_T(int p, Vector s, Vector sds, Vector sdsxds, real_t mu, Vector grad_mu, DenseTensor &u) const
Definition fe_pyramid.cpp:1345
mfem::FuentesPyramid::lam25_grad_lam25
static Vector lam25_grad_lam25(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:94
mfem::FuentesPyramid::lam415_grad_lam415
static Vector lam415_grad_lam415(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:146
mfem::FuentesPyramid::lam235_grad_lam235
static Vector lam235_grad_lam235(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:125
mfem::FuentesPyramid::mu01
static Vector mu01(real_t z)
Definition fe_pyramid.hpp:177
mfem::FuentesPyramid::grad_nu012
static DenseMatrix grad_nu012(real_t z, Vector xy, unsigned int ab)
Definition fe_pyramid.cpp:241
mfem::FuentesPyramid::CalcIntegratedJacobi
static void CalcIntegratedJacobi(int p, real_t alpha, real_t x, real_t t, real_t *u)
Integrated Jacobi Polynomials.
Definition fe_pyramid.cpp:525
mfem::FuentesPyramid::div_lam145_grad_lam145
static real_t div_lam145_grad_lam145(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:175
mfem::FuentesPyramid::V_L
void V_L(int p, Vector sx, const DenseMatrix &grad_sx, Vector sy, const DenseMatrix &grad_sy, real_t t, Vector grad_t, DenseTensor &u) const
Definition fe_pyramid.cpp:1469
mfem::FuentesPyramid::grad_lam4
static Vector grad_lam4(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:42
mfem::FuentesPyramid::grad_nu2
static Vector grad_nu2(real_t z, const Vector xy, unsigned int ab)
Definition fe_pyramid.cpp:228
mfem::FuentesPyramid::grad_mu1
static Vector grad_mu1(real_t z)
Definition fe_pyramid.hpp:174
mfem::FuentesPyramid::V_Q
void V_Q(int p, Vector s, Vector ds, Vector t, Vector dt, DenseTensor &u) const
Definition fe_pyramid.cpp:1238
mfem::FuentesPyramid::CalcScaledLegendre
static void CalcScaledLegendre(int p, real_t x, real_t t, real_t *u)
Shifted and Scaled Legendre Polynomials.
Definition fe_pyramid.cpp:293
mfem::FuentesPyramid::E_T
void E_T(int p, Vector s, Vector sds, DenseTensor &u) const
Definition fe_pyramid.cpp:1137
mfem::FuentesPyramid::CalcHomogenizedScaJacobi
static void CalcHomogenizedScaJacobi(int p, real_t alpha, real_t t0, real_t t1, real_t *u)
Definition fe_pyramid.hpp:373
mfem::FuentesPyramid::lam15_grad_lam15
static Vector lam15_grad_lam15(real_t x, real_t y, real_t z)
Computes .
Definition fe_pyramid.cpp:86
mfem::FuentesPyramid::lam35_grad_lam35
static Vector lam35_grad_lam35(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:102
mfem::FuentesPyramid::div_lam415_grad_lam415
static real_t div_lam415_grad_lam415(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:172
mfem::FuentesPyramid::V_T
void V_T(int p, Vector s, Vector sdsxds, DenseTensor &u) const
Definition fe_pyramid.cpp:1272
mfem::FuentesPyramid::lam435_grad_lam435
static Vector lam435_grad_lam435(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:139
mfem::FuentesPyramid::div_lam235_grad_lam235
static real_t div_lam235_grad_lam235(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:163
mfem::FuentesPyramid::grad_mu0
static Vector grad_mu0(real_t z)
Definition fe_pyramid.hpp:172
mfem::FuentesPyramid::one
static constexpr real_t one
Definition fe_pyramid.hpp:80
mfem::FuentesPyramid::grad_lam35
static DenseMatrix grad_lam35(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:70
mfem::FuentesPyramid::grad_lam15
static DenseMatrix grad_lam15(real_t x, real_t y, real_t z)
Gradients of the above two component vectors.
Definition fe_pyramid.cpp:54
mfem::FuentesPyramid::grad_nu1
static Vector grad_nu1(real_t z, const Vector xy, unsigned int ab)
Definition fe_pyramid.cpp:222
mfem::FuentesPyramid::E_Q
void E_Q(int p, Vector s, Vector ds, Vector t, DenseTensor &u) const
Definition fe_pyramid.cpp:1027
mfem::FuentesPyramid::grad_lam1
static Vector grad_lam1(real_t x, real_t y, real_t z)
Gradients of the "Affine" Coordinates.
Definition fe_pyramid.cpp:21
mfem::FuentesPyramid::lam15
static Vector lam15(real_t x, real_t y, real_t z)
Two component vectors associated with edges touching the apex.
Definition fe_pyramid.hpp:109
mfem::FuentesPyramid::mu01_grad_mu01
static Vector mu01_grad_mu01(real_t z, Vector xy, unsigned int ab)
Definition fe_pyramid.cpp:208
mfem::FuentesPyramid::nu012_grad_nu012
static Vector nu012_grad_nu012(real_t z, Vector xy, unsigned int ab)
Definition fe_pyramid.cpp:275
mfem::FuentesPyramid::grad_lam2
static Vector grad_lam2(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:28
mfem::FuentesPyramid::lam345_grad_lam345
static Vector lam345_grad_lam345(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:132
mfem::FuentesPyramid::CalcHomogenizedIntLegendre
static void CalcHomogenizedIntLegendre(int p, real_t t0, real_t t1, real_t *u)
Definition fe_pyramid.cpp:485
mfem::FuentesPyramid::div_lam125_grad_lam125
static real_t div_lam125_grad_lam125(real_t x, real_t y, real_t z)
Divergences of the above "normal" vector functions divided by 3.
Definition fe_pyramid.cpp:160
mfem::FuentesPyramid::grad_lam45
static DenseMatrix grad_lam45(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:78
mfem::FuentesPyramid::grad_nu01
static DenseMatrix grad_nu01(real_t z, Vector xy, unsigned int ab)
Definition fe_pyramid.cpp:233
mfem::FuentesPyramid::CheckZ
static bool CheckZ(real_t z)
Definition fe_pyramid.hpp:87
mfem::FuentesPyramid::CalcScaledJacobi
static void CalcScaledJacobi(int p, real_t alpha, real_t x, real_t t, real_t *u)
Shifted and Scaled Jacobi Polynomials.
Definition fe_pyramid.cpp:591
mfem::FuentesPyramid::CalcHomogenizedScaLegendre
static void CalcHomogenizedScaLegendre(int p, real_t s0, real_t s1, real_t *u)
Definition fe_pyramid.cpp:447
mfem::FuentesPyramid::grad_nu0
static Vector grad_nu0(real_t z, const Vector xy, unsigned int ab)
Definition fe_pyramid.cpp:216
mfem::FuentesPyramid::phi_T
void phi_T(int p, Vector nu, DenseMatrix &u) const
Definition fe_pyramid.cpp:868
mfem::FuentesPyramid::nu12_grad_nu12
static Vector nu12_grad_nu12(real_t z, Vector xy, unsigned int ab)
Definition fe_pyramid.cpp:267
mfem::FuentesPyramid::div_lam435_grad_lam435
static real_t div_lam435_grad_lam435(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:169
mfem::FuentesPyramid::grad_lam3
static Vector grad_lam3(real_t x, real_t y, real_t z)
Definition fe_pyramid.cpp:35
mfem::Operator::Height
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
mfem::Operator::Width
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
mfem::Poly_1D::CalcLegendre
static void CalcLegendre(const int p, const real_t x, real_t *u)
Definition fe_base.cpp:2245
mfem::Vector
Vector data type.
Definition vector.hpp:82
mfem::Vector::Size
int Size() const
Returns the size of the vector.
Definition vector.hpp:226
mfem::Vector::SetSize
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:558
mfem::Vector::GetData
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition vector.hpp:235
mfem::Vector::Add
Vector & Add(const real_t a, const Vector &Va)
(*this) += a * Va
Definition vector.cpp:322
mfem::Vector::cross3D
void cross3D(const Vector &vin, Vector &vout) const
Definition vector.cpp:616
alpha
const real_t alpha
Definition ex15.cpp:369
mu
real_t mu
Definition ex25.cpp:140
fe_pyramid.hpp
b
real_t b
Definition lissajous.cpp:42
a
real_t a
Definition lissajous.cpp:41
mfem::u
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
mfem::Mult
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition table.cpp:548
mfem::add
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition vector.cpp:391
mfem::real_t
float real_t
Definition config.hpp:43
p
real_t p(const Vector &x, real_t t)
Definition navier_mms.cpp:53
Generated on Wed Apr 9 2025 16:53:03 for MFEM by doxygen 1.11.0