47 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
50 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
55 "order in 'a' direction is higher "
56 "than order in 'c' direction");
58 "order in 'b' direction is higher "
59 "than order in 'c' direction");
62 "Expected basis type in 'c' direction to be ModifiedPyr_C or "
89 int Qx =
m_base[0]->GetNumPoints();
90 int Qy =
m_base[1]->GetNumPoints();
91 int Qz =
m_base[2]->GetNumPoints();
100 eta_y =
m_base[1]->GetZ();
101 eta_z =
m_base[2]->GetZ();
105 if (out_dxi1.size() > 0)
107 for (k = 0, n = 0; k < Qz; ++k)
110 for (j = 0; j < Qy; ++j)
112 for (i = 0; i < Qx; ++i, ++n)
114 out_dxi1[n] = fac * dEta_bar1[n];
120 if (out_dxi2.size() > 0)
122 for (k = 0, n = 0; k < Qz; ++k)
125 for (j = 0; j < Qy; ++j)
127 for (i = 0; i < Qx; ++i, ++n)
129 out_dxi2[n] = fac * dXi2[n];
135 if (out_dxi3.size() > 0)
137 for (k = 0, n = 0; k < Qz; ++k)
140 for (j = 0; j < Qy; ++j)
143 for (i = 0; i < Qx; ++i, ++n)
145 out_dxi3[n] = (1.0 + eta_x[i]) * fac * dEta_bar1[n] +
146 fac1 * fac * dXi2[n] + dEta3[n];
182 ASSERTL1(
false,
"input dir is out of range");
231 if (
m_base[0]->Collocation() &&
m_base[1]->Collocation() &&
236 inarray, 1, outarray, 1);
250 int nquad0 =
m_base[0]->GetNumPoints();
251 int nquad1 =
m_base[1]->GetNumPoints();
252 int nquad2 =
m_base[2]->GetNumPoints();
253 int order0 =
m_base[0]->GetNumModes();
254 int order1 =
m_base[1]->GetNumModes();
257 nquad2 * nquad1 * nquad0);
260 m_base[2]->GetBdata(), inarray, outarray, wsp,
true,
270 [[maybe_unused]]
bool doCheckCollDir0,
271 [[maybe_unused]]
bool doCheckCollDir1,
272 [[maybe_unused]]
bool doCheckCollDir2)
274 int nquad0 =
m_base[0]->GetNumPoints();
275 int nquad1 =
m_base[1]->GetNumPoints();
276 int nquad2 =
m_base[2]->GetNumPoints();
278 int order0 =
m_base[0]->GetNumModes();
279 int order1 =
m_base[1]->GetNumModes();
280 int order2 =
m_base[2]->GetNumModes();
285 int i, j, mode, mode1, cnt;
288 mode = mode1 = cnt = 0;
289 for (i = 0; i < order0; ++i)
291 for (j = 0; j < order1; ++j, ++cnt)
293 int ijmax = max(i, j);
295 base2.data() + mode * nquad2, nquad2,
296 inarray.data() + mode1, 1, 0.0,
297 tmp.data() + cnt * nquad2, 1);
298 mode += order2 - ijmax;
299 mode1 += order2 - ijmax;
302 for (j = order1; j < order2; ++j)
304 int ijmax = max(i, j);
305 mode += order2 - ijmax;
317 Blas::Daxpy(nquad2, inarray[1], base2.data() + nquad2, 1,
318 &tmp[0] + nquad2, 1);
321 Blas::Daxpy(nquad2, inarray[1], base2.data() + nquad2, 1,
322 &tmp[0] + order1 * nquad2, 1);
325 Blas::Daxpy(nquad2, inarray[1], base2.data() + nquad2, 1,
326 &tmp[0] + order1 * nquad2 + nquad2, 1);
331 for (i = 0; i < order0; ++i)
333 Blas::Dgemm(
'N',
'T', nquad1, nquad2, order1, 1.0, base1.data(), nquad1,
334 tmp.data() + mode * nquad2, nquad2, 0.0,
335 tmp1.data() + i * nquad1 * nquad2, nquad1);
340 Blas::Dgemm(
'N',
'T', nquad0, nquad1 * nquad2, order0, 1.0, base0.data(),
341 nquad0, tmp1.data(), nquad1 * nquad2, 0.0, outarray.data(),
371 out = (*imatsys) * in;
396 if (
m_base[0]->Collocation() &&
m_base[1]->Collocation() &&
412 int nquad1 =
m_base[1]->GetNumPoints();
413 int nquad2 =
m_base[2]->GetNumPoints();
414 int order0 =
m_base[0]->GetNumModes();
415 int order1 =
m_base[1]->GetNumModes();
419 if (multiplybyweights)
427 tmp, outarray, wsp,
true,
true,
true);
433 inarray, outarray, wsp,
true,
true,
true);
443 [[maybe_unused]]
bool doCheckCollDir0,
444 [[maybe_unused]]
bool doCheckCollDir1,
445 [[maybe_unused]]
bool doCheckCollDir2)
447 int nquad0 =
m_base[0]->GetNumPoints();
448 int nquad1 =
m_base[1]->GetNumPoints();
449 int nquad2 =
m_base[2]->GetNumPoints();
451 int order0 =
m_base[0]->GetNumModes();
452 int order1 =
m_base[1]->GetNumModes();
453 int order2 =
m_base[2]->GetNumModes();
455 ASSERTL1(wsp.size() >= nquad1 * nquad2 * order0 + nquad2 * order0 * order1,
456 "Insufficient workspace size");
461 int i, j, mode, mode1, cnt;
464 Blas::Dgemm(
'T',
'N', nquad1 * nquad2, order0, nquad0, 1.0, inarray.data(),
465 nquad0, base0.data(), nquad0, 0.0, tmp1.data(),
469 for (mode = i = 0; i < order0; ++i)
472 tmp1.data() + i * nquad1 * nquad2, nquad1, base1.data(),
473 nquad1, 0.0, tmp2.data() + mode * nquad2, nquad2);
478 mode = mode1 = cnt = 0;
479 for (i = 0; i < order0; ++i)
481 for (j = 0; j < order1; ++j, ++cnt)
483 int ijmax = max(i, j);
486 base2.data() + mode * nquad2, nquad2,
487 tmp2.data() + cnt * nquad2, 1, 0.0,
488 outarray.data() + mode1, 1);
489 mode += order2 - ijmax;
490 mode1 += order2 - ijmax;
494 for (j = order1; j < order2; ++j)
496 int ijmax = max(i, j);
497 mode += order2 - ijmax;
507 Blas::Ddot(nquad2, base2.data() + nquad2, 1, &tmp2[nquad2], 1);
510 outarray[1] +=
Blas::Ddot(nquad2, base2.data() + nquad2, 1,
511 &tmp2[nquad2 * order1], 1);
514 outarray[1] +=
Blas::Ddot(nquad2, base2.data() + nquad2, 1,
515 &tmp2[nquad2 * order1 + nquad2], 1);
537 int nquad0 =
m_base[0]->GetNumPoints();
538 int nquad1 =
m_base[1]->GetNumPoints();
539 int nquad2 =
m_base[2]->GetNumPoints();
540 int nqtot = nquad0 * nquad1 * nquad2;
547 int order0 =
m_base[0]->GetNumModes();
548 int order1 =
m_base[1]->GetNumModes();
551 nquad2 * order0 * order1);
558 for (i = 0; i < nquad0; ++i)
560 gfac0[i] = 0.5 * (1 + z0[i]);
564 for (i = 0; i < nquad1; ++i)
566 gfac1[i] = 0.5 * (1 + z1[i]);
570 for (i = 0; i < nquad2; ++i)
572 gfac2[i] = 2.0 / (1 - z2[i]);
576 const int nq01 = nquad0 * nquad1;
577 for (i = 0; i < nquad2; ++i)
579 Vmath::Smul(nq01, gfac2[i], &inarray[0] + i * nq01, 1,
580 &tmp0[0] + i * nq01, 1);
591 m_base[2]->GetBdata(), tmp0, outarray, wsp,
false,
true,
true);
598 m_base[2]->GetBdata(), tmp0, outarray, wsp,
true,
false,
true);
607 for (i = 0; i < nquad1 * nquad2; ++i)
609 Vmath::Vmul(nquad0, tmp0.data() + i * nquad0, 1, gfac0.data(),
610 1, tmp0.data() + i * nquad0, 1);
614 m_base[2]->GetBdata(), tmp0, tmp3, wsp,
false,
true,
true);
617 for (i = 0; i < nquad2; ++i)
619 Vmath::Smul(nq01, gfac2[i], &inarray[0] + i * nq01, 1,
620 &tmp0[0] + i * nq01, 1);
622 for (i = 0; i < nquad1 * nquad2; ++i)
624 Vmath::Smul(nquad0, gfac1[i % nquad1], &tmp0[0] + i * nquad0, 1,
625 &tmp0[0] + i * nquad0, 1);
631 m_base[2]->GetBdata(), tmp0, tmp4, wsp,
true,
false,
true);
636 m_base[2]->GetDbdata(), tmp0, outarray, wsp,
true,
true,
false);
646 ASSERTL1(
false,
"input dir is out of range");
672 eta[1] = 2.0 * (1.0 + xi[1]) / d2 - 1.0;
673 eta[0] = 2.0 * (1.0 + xi[0]) / d2 - 1.0;
679 xi[0] = (1.0 + eta[0]) * (1.0 - eta[2]) * 0.5 - 1.0;
680 xi[1] = (1.0 + eta[1]) * (1.0 - eta[2]) * 0.5 - 1.0;
696 for (
int k = 0; k < Qz; ++k)
698 for (
int j = 0; j < Qy; ++j)
700 for (
int i = 0; i < Qx; ++i)
702 int s = i + Qx * (j + Qy * k);
705 xi_y[s] = (1.0 + eta_y[j]) * (1.0 - eta_z[k]) / 2.0 - 1.0;
706 xi_x[s] = (1.0 + etaBar_x[i]) * (1.0 - eta_z[k]) / 2.0 - 1.0;
718 const int nm0 =
m_base[0]->GetNumModes();
719 const int nm1 =
m_base[1]->GetNumModes();
720 const int nm2 =
m_base[2]->GetNumModes();
722 int mode0 = 0, mode1 = 0, mode2 = 0, cnt = 0;
725 for (mode0 = 0; mode0 < nm0; ++mode0)
727 for (mode1 = 0; mode1 < nm1; ++mode1)
729 int maxpq = max(mode0, mode1);
730 for (mode2 = 0; mode2 < nm2 - maxpq; ++mode2, ++cnt)
750 for (
int j = nm1; j < nm2; ++j)
752 int ijmax = max(mode0, j);
753 mode2 += nm2 - ijmax;
759 return StdExpansion::BaryEvaluateBasis<2>(coll[2], 1);
763 return StdExpansion::BaryEvaluateBasis<0>(coll[0], mode0) *
764 StdExpansion::BaryEvaluateBasis<1>(coll[1], mode1) *
765 StdExpansion::BaryEvaluateBasis<2>(coll[2], mode2);
772 std::array<NekDouble, 3> &firstOrderDerivs)
779 if ((1 - coll[2]) < 1e-5)
783 EphysDeriv2(totPoints);
784 PhysDeriv(inarray, EphysDeriv0, EphysDeriv1, EphysDeriv2);
787 I[0] =
GetBase()[0]->GetI(coll);
788 I[1] =
GetBase()[1]->GetI(coll + 1);
789 I[2] =
GetBase()[2]->GetI(coll + 2);
797 std::array<NekDouble, 3> interDeriv;
802 firstOrderDerivs[0] = fac * interDeriv[0];
803 firstOrderDerivs[1] = fac * interDeriv[1];
804 firstOrderDerivs[2] = ((1.0 + coll[0]) / (1.0 - coll[2])) * interDeriv[0] +
805 ((1.0 + coll[1]) / (1.0 - coll[2])) * interDeriv[1] +
821 int nummodes[3] = {
m_base[0]->GetNumModes(),
m_base[1]->GetNumModes(),
822 m_base[2]->GetNumModes()};
828 numModes0 = nummodes[0];
829 numModes1 = nummodes[1];
835 numModes0 = nummodes[0];
836 numModes1 = nummodes[2];
842 numModes0 = nummodes[1];
843 numModes1 = nummodes[2];
850 std::swap(numModes0, numModes1);
882 "BasisType is not a boundary interior form");
885 "BasisType is not a boundary interior form");
888 "BasisType is not a boundary interior form");
890 int P =
m_base[0]->GetNumModes();
891 int Q =
m_base[1]->GetNumModes();
892 int R =
m_base[2]->GetNumModes();
901 "BasisType is not a boundary interior form");
904 "BasisType is not a boundary interior form");
907 "BasisType is not a boundary interior form");
909 int P =
m_base[0]->GetNumModes() - 1;
910 int Q =
m_base[1]->GetNumModes() - 1;
911 int R =
m_base[2]->GetNumModes() - 1;
913 return (
P + 1) * (Q + 1)
914 + 2 * (R + 1) +
P * (1 + 2 * R -
P)
915 + 2 * (R + 1) + Q * (1 + 2 * R - Q);
920 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
926 else if (i == 1 || i == 3)
929 return Q + 1 + (
P * (1 + 2 * Q -
P)) / 2;
934 return Q + 1 + (
P * (1 + 2 * Q -
P)) / 2;
940 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
942 int P =
m_base[0]->GetNumModes() - 1;
943 int Q =
m_base[1]->GetNumModes() - 1;
944 int R =
m_base[2]->GetNumModes() - 1;
948 return (
P - 1) * (Q - 1);
950 else if (i == 1 || i == 3)
952 return (
P - 1) * (2 * (R - 1) - (
P - 1) - 1) / 2;
956 return (Q - 1) * (2 * (R - 1) - (Q - 1) - 1) / 2;
962 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
966 return m_base[0]->GetNumPoints() *
m_base[1]->GetNumPoints();
968 else if (i == 1 || i == 3)
970 return m_base[0]->GetNumPoints() *
m_base[2]->GetNumPoints();
974 return m_base[1]->GetNumPoints() *
m_base[2]->GetNumPoints();
980 ASSERTL2(i >= 0 && i <= 7,
"edge id is out of range");
982 if (i == 0 || i == 2)
986 else if (i == 1 || i == 3)
1000 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
1001 ASSERTL2(k >= 0 && k <= 1,
"basis key id is out of range");
1009 m_base[k]->GetNumModes());
1016 m_base[2 * k]->GetNumModes(), UseGLL);
1023 m_base[k + 1]->GetNumModes(), UseGLL);
1032 const std::vector<unsigned int> &nummodes,
int &modes_offset)
1035 nummodes[modes_offset], nummodes[modes_offset + 1],
1036 nummodes[modes_offset + 2]);
1047 "Mapping not defined for this type of basis");
1051 if (useCoeffPacking ==
true)
1071 ASSERTL0(
false,
"local vertex id must be between 0 and 4");
1094 ASSERTL0(
false,
"local vertex id must be between 0 and 4");
1105 "BasisType is not a boundary interior form");
1108 "BasisType is not a boundary interior form");
1111 "BasisType is not a boundary interior form");
1113 int P =
m_base[0]->GetNumModes() - 1,
p;
1114 int Q =
m_base[1]->GetNumModes() - 1,
q;
1115 int R =
m_base[2]->GetNumModes() - 1, r;
1119 if (outarray.size() != nIntCoeffs)
1127 for (
p = 2;
p <=
P; ++
p)
1129 for (
q = 2;
q <= Q; ++
q)
1131 int maxpq = max(
p,
q);
1132 for (r = 1; r <= R - maxpq; ++r)
1144 "BasisType is not a boundary interior form");
1147 "BasisType is not a boundary interior form");
1150 "BasisType is not a boundary interior form");
1152 int P =
m_base[0]->GetNumModes() - 1,
p;
1153 int Q =
m_base[1]->GetNumModes() - 1,
q;
1154 int R =
m_base[2]->GetNumModes() - 1, r;
1159 if (maparray.size() != nBnd)
1165 for (
p = 0;
p <=
P; ++
p)
1170 for (
q = 0;
q <= Q; ++
q)
1172 int maxpq = max(
p,
q);
1173 for (r = 0; r <= R - maxpq; ++r)
1183 for (
q = 0;
q <= Q; ++
q)
1187 for (r = 0; r <= R -
p; ++r)
1205 "Method only implemented if BasisType is identical"
1206 "in x and y directions");
1209 "Method only implemented for Modified_A BasisType"
1210 "(x and y direction) and ModifiedPyr_C BasisType (z "
1213 int p,
q, r,
P = 0, Q = 0, idx = 0;
1215 int order0 =
m_base[0]->GetNumModes();
1216 int order1 =
m_base[1]->GetNumModes();
1217 int order2 =
m_base[2]->GetNumModes();
1236 ASSERTL0(
false,
"fid must be between 0 and 4");
1239 if (maparray.size() !=
P * Q)
1250 for (
q = 0;
q < Q; ++
q)
1252 for (
p = 0;
p <
P; ++
p)
1260 for (
p = 0;
p <
P; ++
p)
1262 for (r = 0; r < Q -
p; ++r)
1264 maparray[idx++] =
GetMode(
p, 0, r);
1270 maparray[idx++] =
GetMode(1, 0, 0);
1271 maparray[idx++] =
GetMode(0, 0, 1);
1272 for (r = 1; r < Q - 1; ++r)
1274 maparray[idx++] =
GetMode(1, 0, r);
1277 for (
q = 1;
q <
P; ++
q)
1279 for (r = 0; r < Q -
q; ++r)
1281 maparray[idx++] =
GetMode(1,
q, r);
1287 maparray[idx++] =
GetMode(0, 1, 0);
1288 maparray[idx++] =
GetMode(0, 0, 1);
1289 for (r = 1; r < Q - 1; ++r)
1291 maparray[idx++] =
GetMode(0, 1, r);
1294 for (
p = 1;
p <
P; ++
p)
1296 for (r = 0; r < Q -
p; ++r)
1298 maparray[idx++] =
GetMode(
p, 1, r);
1304 for (
q = 0;
q <
P; ++
q)
1306 for (r = 0; r < Q -
q; ++r)
1308 maparray[idx++] =
GetMode(0,
q, r);
1314 ASSERTL0(
false,
"Face to element map unavailable.");
1324 "Method only implemented if BasisType is identical"
1325 "in x and y directions");
1328 "Method only implemented for Modified_A BasisType"
1329 "(x and y direction) and ModifiedPyr_C BasisType (z "
1332 int i, j, k,
p, r, nFaceCoeffs;
1333 int nummodesA = 0, nummodesB = 0;
1335 int order0 =
m_base[0]->GetNumModes();
1336 int order1 =
m_base[1]->GetNumModes();
1337 int order2 =
m_base[2]->GetNumModes();
1356 ASSERTL0(
false,
"fid must be between 0 and 4");
1367 nFaceCoeffs =
P * (2 * Q -
P + 1) / 2;
1371 nFaceCoeffs =
P * Q;
1375 if (maparray.size() != nFaceCoeffs)
1380 if (signarray.size() != nFaceCoeffs)
1386 fill(signarray.data(), signarray.data() + nFaceCoeffs, 1);
1396 int minPA = min(nummodesA,
P);
1397 int minQB = min(nummodesB, Q);
1399 for (j = 0; j < minPA; ++j)
1402 for (k = 0; k < minQB - j; ++k, ++cnt)
1404 maparray[idx++] = cnt;
1407 cnt += nummodesB - minQB;
1409 for (k = nummodesB - j; k < Q - j; ++k)
1411 signarray[idx] = 0.0;
1412 maparray[idx++] = maparray[0];
1417 for (j = minPA; j < nummodesA; ++j)
1420 for (k = 0; k < minQB-j; ++k, ++cnt)
1422 maparray[idx++] = cnt;
1425 cnt += nummodesB-minQB;
1427 for (k = nummodesB-j; k < Q-j; ++k)
1429 signarray[idx] = 0.0;
1430 maparray[idx++] = maparray[0];
1434 for (j = nummodesA; j <
P; ++j)
1436 for (k = 0; k < Q - j; ++k)
1438 signarray[idx] = 0.0;
1439 maparray[idx++] = maparray[0];
1447 swap(maparray[0], maparray[Q]);
1448 for (i = 1; i < Q - 1; ++i)
1450 swap(maparray[i + 1], maparray[Q + i]);
1454 for (
p = 0;
p <
P; ++
p)
1456 for (r = 0; r < Q -
p; ++r, idx++)
1460 signarray[idx] =
p % 2 ? -1 : 1;
1473 for (i = 0; i < Q; i++)
1475 for (j = 0; j <
P; j++)
1479 arrayindx[i *
P + j] = i *
P + j;
1483 arrayindx[i *
P + j] = j * Q + i;
1490 for (j = 0; j <
P; ++j)
1493 for (k = 0; k < Q; k++)
1495 maparray[arrayindx[j + k *
P]] = j + k * nummodesA;
1498 for (k = nummodesB; k < Q; ++k)
1500 signarray[arrayindx[j + k *
P]] = 0.0;
1501 maparray[arrayindx[j + k *
P]] = maparray[0];
1505 for (j = nummodesA; j <
P; ++j)
1507 for (k = 0; k < Q; ++k)
1509 signarray[arrayindx[j + k *
P]] = 0.0;
1510 maparray[arrayindx[j + k *
P]] = maparray[0];
1525 for (i = 3; i < Q; i += 2)
1527 for (j = 0; j <
P; j++)
1529 signarray[arrayindx[i *
P + j]] *= -1;
1533 for (i = 0; i <
P; i++)
1535 swap(maparray[i], maparray[i +
P]);
1536 swap(signarray[i], signarray[i +
P]);
1541 for (i = 0; i < Q; i++)
1543 for (j = 3; j <
P; j += 2)
1545 signarray[arrayindx[i *
P + j]] *= -1;
1549 for (i = 0; i < Q; i++)
1551 swap(maparray[i], maparray[i + Q]);
1552 swap(signarray[i], signarray[i + Q]);
1564 for (i = 0; i < Q; i++)
1566 for (j = 3; j <
P; j += 2)
1568 signarray[arrayindx[i *
P + j]] *= -1;
1572 for (i = 0; i < Q; i++)
1574 swap(maparray[i *
P], maparray[i *
P + 1]);
1575 swap(signarray[i *
P], signarray[i *
P + 1]);
1580 for (i = 3; i < Q; i += 2)
1582 for (j = 0; j <
P; j++)
1584 signarray[arrayindx[i *
P + j]] *= -1;
1588 for (i = 0; i <
P; i++)
1590 swap(maparray[i * Q], maparray[i * Q + 1]);
1591 swap(signarray[i * Q], signarray[i * Q + 1]);
1604 const int P =
m_base[0]->GetNumModes() - 1;
1605 const int Q =
m_base[1]->GetNumModes() - 1;
1606 const int R =
m_base[2]->GetNumModes() - 1;
1609 if (maparray.size() != nEdgeIntCoeffs)
1614 if (signarray.size() != nEdgeIntCoeffs)
1620 fill(signarray.data(), signarray.data() + nEdgeIntCoeffs, 1);
1630 for (i = 2; i <=
P; ++i)
1632 maparray[i - 2] =
GetMode(i, 0, 0);
1637 for (i = 2; i <= Q; ++i)
1639 maparray[i - 2] =
GetMode(1, i, 0);
1643 for (i = 2; i <=
P; ++i)
1645 maparray[i - 2] =
GetMode(i, 1, 0);
1650 for (i = 2; i <= Q; ++i)
1652 maparray[i - 2] =
GetMode(0, i, 0);
1656 for (i = 2; i <= R; ++i)
1658 maparray[i - 2] =
GetMode(0, 0, i);
1663 for (i = 1; i <= R - 1; ++i)
1665 maparray[i - 1] =
GetMode(1, 0, i);
1669 for (i = 1; i <= R - 1; ++i)
1671 maparray[i - 1] =
GetMode(1, 1, i);
1676 for (i = 1; i <= R - 1; ++i)
1678 maparray[i - 1] =
GetMode(0, 1, i);
1682 ASSERTL0(
false,
"Edge not defined.");
1688 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1699 const int P =
m_base[0]->GetNumModes() - 1;
1700 const int Q =
m_base[1]->GetNumModes() - 1;
1701 const int R =
m_base[2]->GetNumModes() - 1;
1703 int p,
q, r, idx = 0;
1708 if (maparray.size() != nFaceIntCoeffs)
1713 if (signarray.size() != nFaceIntCoeffs)
1719 fill(signarray.data(), signarray.data() + nFaceIntCoeffs, 1);
1730 for (i = 0; i < nummodesB; i++)
1732 for (j = 0; j < nummodesA; j++)
1736 arrayindx[i * nummodesA + j] = i * nummodesA + j;
1740 arrayindx[i * nummodesA + j] = j * nummodesB + i;
1749 for (
q = 2;
q <= Q; ++
q)
1751 for (
p = 2;
p <=
P; ++
p)
1753 maparray[arrayindx[(
q - 2) * nummodesA + (
p - 2)]] =
1759 for (
p = 2;
p <=
P; ++
p)
1761 for (r = 1; r <= R -
p; ++r)
1763 if ((
int)faceOrient == 7)
1765 signarray[idx] =
p % 2 ? -1 : 1;
1767 maparray[idx++] =
GetMode(
p, 0, r);
1772 for (
q = 2;
q <= Q; ++
q)
1774 for (r = 1; r <= R -
q; ++r)
1776 if ((
int)faceOrient == 7)
1778 signarray[idx] =
q % 2 ? -1 : 1;
1780 maparray[idx++] =
GetMode(1,
q, r);
1786 for (
p = 2;
p <=
P; ++
p)
1788 for (r = 1; r <= R -
p; ++r)
1790 if ((
int)faceOrient == 7)
1792 signarray[idx] =
p % 2 ? -1 : 1;
1794 maparray[idx++] =
GetMode(
p, 1, r);
1800 for (
q = 2;
q <= Q; ++
q)
1802 for (r = 1; r <= R -
q; ++r)
1804 if ((
int)faceOrient == 7)
1806 signarray[idx] =
q % 2 ? -1 : 1;
1808 maparray[idx++] =
GetMode(0,
q, r);
1813 ASSERTL0(
false,
"Face interior map not available.");
1823 if (faceOrient == 6 || faceOrient == 8 || faceOrient == 11 ||
1828 for (i = 1; i < nummodesB; i += 2)
1830 for (j = 0; j < nummodesA; j++)
1832 signarray[arrayindx[i * nummodesA + j]] *= -1;
1838 for (i = 0; i < nummodesB; i++)
1840 for (j = 1; j < nummodesA; j += 2)
1842 signarray[arrayindx[i * nummodesA + j]] *= -1;
1848 if (faceOrient == 7 || faceOrient == 8 || faceOrient == 10 ||
1853 for (i = 0; i < nummodesB; i++)
1855 for (j = 1; j < nummodesA; j += 2)
1857 signarray[arrayindx[i * nummodesA + j]] *= -1;
1863 for (i = 1; i < nummodesB; i += 2)
1865 for (j = 0; j < nummodesA; j++)
1867 signarray[arrayindx[i * nummodesA + j]] *= -1;
1914 const int Q =
m_base[1]->GetNumModes() - 1;
1915 const int R =
m_base[2]->GetNumModes() - 1;
1921 for (i = 0; i < I; ++i)
1927 cnt += (R + 1 - i) * (Q + 1) - l * (l + 1) / 2;
1931 l = max(0, J - 1 - I);
1932 cnt += (R + 1 - I) * J - l * (l + 1) / 2;
1946 int nquad0 =
m_base[0]->GetNumPoints();
1947 int nquad1 =
m_base[1]->GetNumPoints();
1948 int nquad2 =
m_base[2]->GetNumPoints();
1957 for (i = 0; i < nquad1 * nquad2; ++i)
1959 Vmath::Vmul(nquad0, inarray.data() + i * nquad0, 1, w0.data(), 1,
1960 outarray.data() + i * nquad0, 1);
1964 for (j = 0; j < nquad2; ++j)
1966 for (i = 0; i < nquad1; ++i)
1969 &outarray[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1979 case LibUtilities::eGaussRadauMAlpha2Beta0:
1980 for (i = 0; i < nquad2; ++i)
1983 &outarray[0] + i * nquad0 * nquad1, 1);
1988 for (i = 0; i < nquad2; ++i)
1991 0.25 * (1 - z2[i]) * (1 - z2[i]) * w2[i],
1992 &outarray[0] + i * nquad0 * nquad1, 1);
2002 int qa =
m_base[0]->GetNumPoints();
2003 int qb =
m_base[1]->GetNumPoints();
2004 int qc =
m_base[2]->GetNumPoints();
2005 int nmodes_a =
m_base[0]->GetNumModes();
2006 int nmodes_b =
m_base[1]->GetNumModes();
2007 int nmodes_c =
m_base[2]->GetNumModes();
2019 int i, j, k, cnt = 0;
2022 OrthoExp.
FwdTrans(array, orthocoeffs);
2032 for (i = 0; i < nmodes_a; ++i)
2034 for (j = 0; j < nmodes_b; ++j)
2036 int maxij = max(i, j);
2038 pow((1.0 * i) / (nmodes_a - 1), cutoff * nmodes_a),
2039 pow((1.0 * j) / (nmodes_b - 1), cutoff * nmodes_b));
2041 for (k = 0; k < nmodes_c - maxij; ++k)
2044 std::max(fac1, pow((1.0 * k) / (nmodes_c - 1),
2045 cutoff * nmodes_c));
2047 orthocoeffs[cnt] *= SvvDiffCoeff * fac;
2063 max_abc = max(max_abc, 0);
2066 for (i = 0; i < nmodes_a; ++i)
2068 for (j = 0; j < nmodes_b; ++j)
2070 int maxij = max(i, j);
2072 for (k = 0; k < nmodes_c - maxij; ++k)
2074 int maxijk = max(maxij, k);
2096 int cutoff_a = (int)(SVVCutOff * nmodes_a);
2097 int cutoff_b = (int)(SVVCutOff * nmodes_b);
2098 int cutoff_c = (int)(SVVCutOff * nmodes_c);
2102 int nmodes = min(min(nmodes_a, nmodes_b), nmodes_c);
2103 NekDouble cutoff = min(min(cutoff_a, cutoff_b), cutoff_c);
2105 for (i = 0; i < nmodes_a; ++i)
2107 for (j = 0; j < nmodes_b; ++j)
2109 int maxij = max(i, j);
2110 for (k = 0; k < nmodes_c - maxij; ++k)
2112 if (j + k >= cutoff || i + k >= cutoff)
2116 exp(-(i + k - nmodes) * (i + k - nmodes) /
2117 ((
NekDouble)((i + k - cutoff + epsilon) *
2118 (i + k - cutoff + epsilon)))) *
2119 exp(-(j - nmodes) * (j - nmodes) /
2121 (j - cutoff + epsilon)))));
2125 orthocoeffs[cnt] *= 0.0;
2134 OrthoExp.
BwdTrans(orthocoeffs, array);
2141 int nquad0 =
m_base[0]->GetNumPoints();
2142 int nquad1 =
m_base[1]->GetNumPoints();
2143 int nquad2 =
m_base[2]->GetNumPoints();
2144 int nqtot = nquad0 * nquad1 * nquad2;
2145 int nmodes0 =
m_base[0]->GetNumModes();
2146 int nmodes1 =
m_base[1]->GetNumModes();
2147 int nmodes2 =
m_base[2]->GetNumModes();
2148 int numMax = nmodes0;
2169 bortho0, bortho1, bortho2);
2172 OrthoPyrExp->FwdTrans(phys_tmp, coeff);
2175 for (u = 0; u < numMin; ++u)
2177 for (i = 0; i < numMin; ++i)
2180 int maxui = max(u, i);
2182 tmp2 = coeff_tmp1 + cnt, 1);
2183 cnt += nmodes2 - maxui;
2186 for (i = numMin; i < nmodes1; ++i)
2188 int maxui = max(u, i);
2189 cnt += numMax - maxui;
2193 OrthoPyrExp->BwdTrans(coeff_tmp1, phys_tmp);
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Describes the specification for a Basis.
BasisType GetBasisType() const
Return type of expansion basis.
int GetNumModes() const
Returns the order of the basis.
Defines a specification for a set of points.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
void IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
NekDouble BaryTensorDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs)
void PhysTensorDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray_d1, Array< OneD, NekDouble > &outarray_d2, Array< OneD, NekDouble > &outarray_d3)
Calculate the 3D derivative in the local tensor/collapsed coordinate at the physical points.
The base class for all shapes.
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
int NumBndryCoeffs(void) const
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
void LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Convert local cartesian coordinate xi into local collapsed coordinates eta.
const Array< OneD, const LibUtilities::BasisSharedPtr > & GetBase() const
This function gets the shared point to basis.
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
NekDouble PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &physvals)
This function evaluates the expansion at a single (arbitrary) point of the domain.
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space.
int GetTraceNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th trace.
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space.
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Array< OneD, LibUtilities::BasisSharedPtr > m_base
NekDouble GetConstFactor(const ConstFactorType &factor) const
bool ConstFactorExists(const ConstFactorType &factor) const
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey) override
int v_NumBndryCoeffs() const override
NekDouble v_PhysEvaluateBasis(const Array< OneD, const NekDouble > &coords, int mode) final
void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
int v_NumDGBndryCoeffs() const override
void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Forward transform from physical quadrature space stored in inarray and evaluate the expansion coeffic...
void v_GetElmtTraceToTraceMap(const unsigned int fid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation faceOrient, int P, int Q) override
int v_GetTraceNcoeffs(const int i) const override
void v_LocCollapsedToLocCoord(const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi) override
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey) override
int v_GetEdgeNcoeffs(const int i) const override
void v_BwdTrans_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2) override
int v_GetNtraces() const override
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray) override
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true) override
int v_GetNedges() const override
int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset) override
int GetMode(int I, int J, int K)
Compute the mode number in the expansion for a particular tensorial combination.
int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false) override
int v_GetTraceNumPoints(const int i) const override
void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Backward transformation is evaluated at the quadrature points.
void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray) override
void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Inner product of inarray over region with respect to the expansion basis m_base[0]->GetBdata(),...
void v_GetInteriorMap(Array< OneD, unsigned int > &outarray) override
void v_GetEdgeInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2) override
void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta) override
DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey) override
int v_GetNverts() const override
void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_GetTraceInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2) override
void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2) override
void v_StdPhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2) override
LibUtilities::ShapeType v_DetShapeType() const override
void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2) override
Calculate the derivative of the physical points.
void v_GetTraceNumModes(const int fid, int &numModes0, int &numModes1, Orientation faceOrient=eDir1FwdDir1_Dir2FwdDir2) override
void v_GetTraceCoeffMap(const unsigned int fid, Array< OneD, unsigned int > &maparray) override
NekDouble v_PhysEvalFirstDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) override
void v_GetCoords(Array< OneD, NekDouble > &xi_x, Array< OneD, NekDouble > &xi_y, Array< OneD, NekDouble > &xi_z) override
const LibUtilities::BasisKey v_GetTraceBasisKey(const int i, const int k, bool UseGLL=false) const override
int v_GetTraceIntNcoeffs(const int i) const override
static void Dgemv(const char &trans, const int &m, const int &n, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = alpha A x plus beta y where A[m x n].
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
static double Ddot(const int &n, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: output = .
static void Dgemm(const char &transa, const char &transb, const int &m, const int &n, const int &k, const double &alpha, const double *a, const int &lda, const double *b, const int &ldb, const double &beta, double *c, const int &ldc)
BLAS level 3: Matrix-matrix multiply C = A x B where op(A)[m x k], op(B)[k x n], C[m x n] DGEMM perfo...
static void Daxpy(const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: y = alpha x plus y.
int getNumberOfCoefficients(int Na, int Nb, int Nc)
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
int getNumberOfCoefficients(int Na)
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
@ P
Monomial polynomials .
@ eOrtho_A
Principle Orthogonal Functions .
@ eGLL_Lagrange
Lagrange for SEM basis .
@ eModifiedPyr_C
Principle Modified Functions.
@ eModified_A
Principle Modified Functions .
@ eOrthoPyr_C
Principle Orthogonal Functions .
static const NekDouble kNekZeroTol
@ eFactorSVVDGKerDiffCoeff
@ eFactorSVVPowerKerDiffCoeff
LibUtilities::BasisKey EvaluateTriFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes, bool UseGLL)
const int kSVVDGFiltermodesmin
std::shared_ptr< StdPyrExp > StdPyrExpSharedPtr
const int kSVVDGFiltermodesmax
LibUtilities::BasisKey EvaluateQuadFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
const NekDouble kSVVDGFilter[9][11]
@ eDir1BwdDir2_Dir2BwdDir1
@ eDir1BwdDir1_Dir2BwdDir2
@ eDir1BwdDir2_Dir2FwdDir1
@ eDir1FwdDir1_Dir2BwdDir2
@ eDir1BwdDir1_Dir2FwdDir2
@ eDir1FwdDir2_Dir2FwdDir1
@ eDir1FwdDir2_Dir2BwdDir1
std::vector< double > q(NPUPPER *NPUPPER)
static Array< OneD, NekDouble > NullNekDouble1DArray
std::shared_ptr< DNekMat > DNekMatSharedPtr
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)