47 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
50 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
54 "order in 'a' direction is higher than order in 'c' direction");
80 int Qx =
m_base[0]->GetNumPoints();
81 int Qy =
m_base[1]->GetNumPoints();
82 int Qz =
m_base[2]->GetNumPoints();
83 int Qtot = Qx * Qy * Qz;
93 bool Do_1 = (out_dxi1.size() > 0) ?
true :
false;
94 bool Do_3 = (out_dxi3.size() > 0) ?
true :
false;
113 for (k = 0; k < Qz; ++k)
116 &dEta_bar1[0] + k * Qx * Qy, 1,
117 &out_dxi1[0] + k * Qx * Qy, 1);
124 for (k = 0; k < Qz; ++k)
127 &dEta_bar1[0] + k * Qx * Qy, 1,
128 &dEta_bar1[0] + k * Qx * Qy, 1);
132 for (i = 0; i < Qx; ++i)
134 Vmath::Svtvp(Qz * Qy, 1.0 + eta_x[i], &dEta_bar1[0] + i, Qx,
135 &out_dxi3[0] + i, Qx, &out_dxi3[0] + i, Qx);
169 ASSERTL1(
false,
"input dir is out of range");
225 "Basis[1] is not a general tensor type");
229 "Basis[2] is not a general tensor type");
231 if (
m_base[0]->Collocation() &&
m_base[1]->Collocation() &&
236 inarray, 1, outarray, 1);
247 int nquad1 =
m_base[1]->GetNumPoints();
248 int nquad2 =
m_base[2]->GetNumPoints();
249 int order0 =
m_base[0]->GetNumModes();
250 int order1 =
m_base[1]->GetNumModes();
253 nquad1 * nquad2 * order0);
256 m_base[2]->GetBdata(), inarray, outarray, wsp,
true,
266 [[maybe_unused]]
bool doCheckCollDir0,
267 [[maybe_unused]]
bool doCheckCollDir1,
268 [[maybe_unused]]
bool doCheckCollDir2)
271 int nquad0 =
m_base[0]->GetNumPoints();
272 int nquad1 =
m_base[1]->GetNumPoints();
273 int nquad2 =
m_base[2]->GetNumPoints();
274 int nummodes0 =
m_base[0]->GetNumModes();
275 int nummodes1 =
m_base[1]->GetNumModes();
276 int nummodes2 =
m_base[2]->GetNumModes();
280 for (i = mode = 0; i < nummodes0; ++i)
282 Blas::Dgemm(
'N',
'N', nquad2, nummodes1, nummodes2 - i, 1.0,
283 base2.data() + mode * nquad2, nquad2,
284 inarray.data() + mode * nummodes1, nummodes2 - i, 0.0,
285 tmp0.data() + i * nquad2 * nummodes1, nquad2);
286 mode += nummodes2 - i;
291 for (i = 0; i < nummodes1; i++)
294 base2.data() + nquad2, 1,
295 tmp0.data() + nquad2 * (nummodes1 + i), 1);
299 for (i = 0; i < nummodes0; i++)
301 Blas::Dgemm(
'N',
'T', nquad1, nquad2, nummodes1, 1.0, base1.data(),
302 nquad1, tmp0.data() + i * nquad2 * nummodes1, nquad2, 0.0,
303 tmp1.data() + i * nquad2 * nquad1, nquad1);
306 Blas::Dgemm(
'N',
'T', nquad0, nquad2 * nquad1, nummodes0, 1.0, base0.data(),
307 nquad0, tmp1.data(), nquad2 * nquad1, 0.0, outarray.data(),
335 out = (*matsys) * in;
373 "Basis[1] is not a general tensor type");
377 "Basis[2] is not a general tensor type");
379 if (
m_base[0]->Collocation() &&
m_base[1]->Collocation())
393 int nquad1 =
m_base[1]->GetNumPoints();
394 int nquad2 =
m_base[2]->GetNumPoints();
395 int order0 =
m_base[0]->GetNumModes();
396 int order1 =
m_base[1]->GetNumModes();
400 if (multiplybyweights)
407 tmp, outarray, wsp,
true,
true,
true);
413 inarray, outarray, wsp,
true,
true,
true);
423 [[maybe_unused]]
bool doCheckCollDir0,
424 [[maybe_unused]]
bool doCheckCollDir1,
425 [[maybe_unused]]
bool doCheckCollDir2)
429 const int nquad0 =
m_base[0]->GetNumPoints();
430 const int nquad1 =
m_base[1]->GetNumPoints();
431 const int nquad2 =
m_base[2]->GetNumPoints();
432 const int order0 =
m_base[0]->GetNumModes();
433 const int order1 =
m_base[1]->GetNumModes();
434 const int order2 =
m_base[2]->GetNumModes();
438 ASSERTL1(wsp.size() >= nquad1 * nquad2 * order0 + nquad2 * order0 * order1,
439 "Insufficient workspace size");
445 Blas::Dgemm(
'T',
'N', nquad1 * nquad2, order0, nquad0, 1.0, inarray.data(),
446 nquad0, base0.data(), nquad0, 0.0, tmp0.data(),
450 Blas::Dgemm(
'T',
'N', nquad2 * order0, order1, nquad1, 1.0, tmp0.data(),
451 nquad1, base1.data(), nquad1, 0.0, tmp1.data(),
455 for (mode = i = 0; i < order0; ++i)
457 Blas::Dgemm(
'T',
'N', order2 - i, order1, nquad2, 1.0,
458 base2.data() + mode * nquad2, nquad2,
459 tmp1.data() + i * nquad2, nquad2 * order0, 0.0,
460 outarray.data() + mode * order1, order2 - i);
468 for (i = 0; i < order1; ++i)
473 tmp1.data() + i * order0 * nquad2 + nquad2, 1);
493 ASSERTL0(dir >= 0 && dir <= 2,
"input dir is out of range");
496 int order0 =
m_base[0]->GetNumModes();
497 int order1 =
m_base[1]->GetNumModes();
498 int nquad0 =
m_base[0]->GetNumPoints();
499 int nquad1 =
m_base[1]->GetNumPoints();
500 int nquad2 =
m_base[2]->GetNumPoints();
510 for (i = 0; i < nquad0; ++i)
512 gfac0[i] = 0.5 * (1 + z0[i]);
516 for (i = 0; i < nquad2; ++i)
518 gfac2[i] = 2.0 / (1 - z2[i]);
524 for (i = 0; i < nquad2; ++i)
527 &inarray[0] + i * nquad0 * nquad1, 1,
528 &tmp0[0] + i * nquad0 * nquad1, 1);
539 m_base[2]->GetBdata(), tmp0, outarray, wsp,
true,
true,
true);
547 m_base[2]->GetBdata(), tmp0, outarray, wsp,
true,
true,
true);
556 for (i = 0; i < nquad1 * nquad2; ++i)
558 Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp0[0] + i * nquad0, 1,
559 &tmp0[0] + i * nquad0, 1);
564 m_base[2]->GetBdata(), tmp0, tmp1, wsp,
true,
true,
true);
569 m_base[2]->GetDbdata(), tmp0, outarray, wsp,
true,
true,
true);
599 eta[0] = 2.0 * (1.0 + xi[0]) / d2 - 1.0;
605 xi[0] = (1.0 + eta[0]) * (1.0 - eta[2]) * 0.5 - 1.0;
622 for (
int k = 0; k < Qz; ++k)
624 for (
int j = 0; j < Qy; ++j)
626 for (
int i = 0; i < Qx; ++i)
628 int s = i + Qx * (j + Qy * k);
629 xi_x[s] = (1.0 - eta_z[k]) * (1.0 + etaBar_x[i]) / 2.0 - 1.0;
643 const int nm1 =
m_base[1]->GetNumModes();
644 const int nm2 =
m_base[2]->GetNumModes();
645 const int b = 2 * nm2 + 1;
647 const int mode0 = floor(0.5 * (b -
sqrt(b * b - 8.0 * mode / nm1)));
649 mode - nm1 * (mode0 * (nm2 - 1) + 1 - (mode0 - 2) * (mode0 - 1) / 2);
650 const int mode1 = tmp / (nm2 - mode0);
651 const int mode2 = tmp % (nm2 - mode0);
653 if (mode0 == 0 && mode2 == 1 &&
657 return StdExpansion::BaryEvaluateBasis<1>(coll[1], mode1) *
658 StdExpansion::BaryEvaluateBasis<2>(coll[2], mode2);
662 return StdExpansion::BaryEvaluateBasis<0>(coll[0], mode0) *
663 StdExpansion::BaryEvaluateBasis<1>(coll[1], mode1) *
664 StdExpansion::BaryEvaluateBasis<2>(coll[2], mode2);
671 std::array<NekDouble, 3> &firstOrderDerivs)
680 if ((1 - coll[2]) < 1e-5)
684 EphysDeriv2(totPoints);
685 PhysDeriv(inarray, EphysDeriv0, EphysDeriv1, EphysDeriv2);
688 I[0] =
GetBase()[0]->GetI(coll);
689 I[1] =
GetBase()[1]->GetI(coll + 1);
690 I[2] =
GetBase()[2]->GetI(coll + 2);
700 NekDouble dEta_bar1 = firstOrderDerivs[0];
703 firstOrderDerivs[0] = fac * dEta_bar1;
706 fac = 1.0 / (1.0 - coll[2]);
707 dEta_bar1 = fac * dEta_bar1;
711 firstOrderDerivs[2] += fac * dEta_bar1;
726 int nummodes[3] = {
m_base[0]->GetNumModes(),
m_base[1]->GetNumModes(),
727 m_base[2]->GetNumModes()};
733 numModes0 = nummodes[0];
734 numModes1 = nummodes[1];
741 numModes0 = nummodes[1];
742 numModes1 = nummodes[2];
749 numModes0 = nummodes[0];
750 numModes1 = nummodes[2];
757 std::swap(numModes0, numModes1);
763 ASSERTL2(i >= 0 && i <= 8,
"edge id is out of range");
765 if (i == 0 || i == 2)
769 else if (i == 1 || i == 3 || i == 8)
811 "BasisType is not a boundary interior form");
814 "BasisType is not a boundary interior form");
817 "BasisType is not a boundary interior form");
819 int P =
m_base[0]->GetNumModes();
820 int Q =
m_base[1]->GetNumModes();
821 int R =
m_base[2]->GetNumModes();
830 "BasisType is not a boundary interior form");
833 "BasisType is not a boundary interior form");
836 "BasisType is not a boundary interior form");
838 int P =
m_base[0]->GetNumModes() - 1;
839 int Q =
m_base[1]->GetNumModes() - 1;
840 int R =
m_base[2]->GetNumModes() - 1;
842 return (
P + 1) * (Q + 1)
843 + 2 * (Q + 1) * (R + 1)
844 + 2 * (R + 1) +
P * (1 + 2 * R -
P);
849 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
854 else if (i == 1 || i == 3)
857 return Q + 1 + (
P * (1 + 2 * Q -
P)) / 2;
867 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
877 else if (i == 1 || i == 3)
879 return Pi * (2 * Ri - Pi - 1) / 2;
889 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
893 return m_base[0]->GetNumPoints() *
m_base[1]->GetNumPoints();
895 else if (i == 1 || i == 3)
897 return m_base[0]->GetNumPoints() *
m_base[2]->GetNumPoints();
901 return m_base[1]->GetNumPoints() *
m_base[2]->GetNumPoints();
908 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
909 ASSERTL2(j == 0 || j == 1,
"face direction is out of range");
913 return m_base[j]->GetPointsKey();
915 else if (i == 1 || i == 3)
917 return m_base[2 * j]->GetPointsKey();
921 return m_base[j + 1]->GetPointsKey();
929 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
930 ASSERTL2(k >= 0 && k <= 1,
"basis key id is out of range");
938 m_base[k]->GetNumModes());
945 m_base[k + 1]->GetNumModes());
952 m_base[2 * k]->GetNumModes(), UseGLL);
962 const std::vector<unsigned int> &nummodes,
int &modes_offset)
965 nummodes[modes_offset], nummodes[modes_offset + 1],
966 nummodes[modes_offset + 2]);
988 "Mapping not defined for this type of basis");
992 if (useCoeffPacking ==
true)
1015 ASSERTL0(
false,
"local vertex id must be between 0 and 5");
1041 ASSERTL0(
false,
"local vertex id must be between 0 and 5");
1052 "BasisType is not a boundary interior form");
1055 "BasisType is not a boundary interior form");
1058 "BasisType is not a boundary interior form");
1060 int P =
m_base[0]->GetNumModes() - 1,
p;
1061 int Q =
m_base[1]->GetNumModes() - 1,
q;
1062 int R =
m_base[2]->GetNumModes() - 1, r;
1066 if (outarray.size() != nIntCoeffs)
1074 for (
p = 2;
p <=
P; ++
p)
1076 for (
q = 2;
q <= Q; ++
q)
1078 for (r = 1; r <= R -
p; ++r)
1090 "BasisType is not a boundary interior form");
1093 "BasisType is not a boundary interior form");
1096 "BasisType is not a boundary interior form");
1098 int P =
m_base[0]->GetNumModes() - 1,
p;
1099 int Q =
m_base[1]->GetNumModes() - 1,
q;
1100 int R =
m_base[2]->GetNumModes() - 1, r;
1105 if (maparray.size() != nBnd)
1111 for (
p = 0;
p <=
P; ++
p)
1116 for (
q = 0;
q <= Q; ++
q)
1118 for (r = 0; r <= R -
p; ++r)
1128 for (
q = 0;
q <= Q; ++
q)
1132 for (r = 0; r <= R -
p; ++r)
1150 "Method only implemented if BasisType is identical"
1151 "in x and y directions");
1154 "Method only implemented for Modified_A BasisType"
1155 "(x and y direction) and Modified_B BasisType (z "
1157 int p,
q, r, idx = 0;
1164 Q =
m_base[1]->GetNumModes();
1169 Q =
m_base[2]->GetNumModes();
1174 Q =
m_base[2]->GetNumModes();
1177 ASSERTL0(
false,
"fid must be between 0 and 4");
1180 if (maparray.size() !=
P * Q)
1190 for (
q = 0;
q < Q; ++
q)
1192 for (
p = 0;
p <
P; ++
p)
1199 for (
p = 0;
p <
P; ++
p)
1201 for (r = 0; r < Q -
p; ++r)
1203 maparray[idx++] =
GetMode(
p, 0, r);
1208 for (
q = 0;
q <
P; ++
q)
1212 for (
q = 0;
q <
P; ++
q)
1216 for (r = 1; r < Q - 1; ++r)
1218 for (
q = 0;
q <
P; ++
q)
1220 maparray[(r + 1) *
P +
q] =
GetMode(1,
q, r);
1225 for (
p = 0;
p <
P; ++
p)
1227 for (r = 0; r < Q -
p; ++r)
1229 maparray[idx++] =
GetMode(
p, 1, r);
1234 for (r = 0; r < Q; ++r)
1236 for (
q = 0;
q <
P; ++
q)
1243 ASSERTL0(
false,
"Face to element map unavailable.");
1253 "Method only implemented if BasisType is identical"
1254 "in x and y directions");
1257 "Method only implemented for Modified_A BasisType"
1258 "(x and y direction) and Modified_B BasisType (z "
1261 int i, j, k,
p, r, nFaceCoeffs, idx = 0;
1262 int nummodesA = 0, nummodesB = 0;
1267 nummodesA =
m_base[0]->GetNumModes();
1268 nummodesB =
m_base[1]->GetNumModes();
1272 nummodesA =
m_base[0]->GetNumModes();
1273 nummodesB =
m_base[2]->GetNumModes();
1277 nummodesA =
m_base[1]->GetNumModes();
1278 nummodesB =
m_base[2]->GetNumModes();
1281 ASSERTL0(
false,
"fid must be between 0 and 4");
1290 else if (fid == 1 || fid == 3)
1292 nFaceCoeffs =
P * (2 * Q -
P + 1) / 2;
1296 nFaceCoeffs =
P * Q;
1300 if (maparray.size() != nFaceCoeffs)
1305 if (signarray.size() != nFaceCoeffs)
1311 fill(signarray.data(), signarray.data() + nFaceCoeffs, 1);
1314 int minPA = min(nummodesA,
P);
1315 int minQB = min(nummodesB, Q);
1317 if (fid == 1 || fid == 3)
1324 for (j = 0; j < minPA; ++j)
1327 for (k = 0; k < minQB - j; ++k, ++cnt)
1329 maparray[idx++] = cnt;
1332 cnt += nummodesB - minQB;
1335 for (k = nummodesB - j; k < Q - j; ++k)
1337 signarray[idx] = 0.0;
1338 maparray[idx++] = maparray[0];
1342 for (j = minPA; j < nummodesA; ++j)
1345 for (k = 0; k < minQB-j; ++k, ++cnt)
1347 maparray[idx++] = cnt;
1350 cnt += nummodesB-minQB;
1353 for (k = nummodesB-j; k < Q-j; ++k)
1355 signarray[idx] = 0.0;
1356 maparray[idx++] = maparray[0];
1360 for (j = nummodesA; j <
P; ++j)
1362 for (k = 0; k < Q - j; ++k)
1364 signarray[idx] = 0.0;
1365 maparray[idx++] = maparray[0];
1374 for (
p = 0;
p <
P; ++
p)
1376 for (r = 0; r < Q -
p; ++r, idx++)
1380 signarray[idx] =
p % 2 ? -1 : 1;
1385 swap(maparray[0], maparray[Q]);
1386 for (i = 1; i < Q - 1; ++i)
1388 swap(maparray[i + 1], maparray[Q + i]);
1398 for (i = 0; i < Q; i++)
1400 for (j = 0; j <
P; j++)
1404 arrayindx[i *
P + j] = i *
P + j;
1408 arrayindx[i *
P + j] = j * Q + i;
1415 for (j = 0; j <
P; ++j)
1418 for (k = 0; k < Q; k++)
1420 maparray[arrayindx[j + k *
P]] = j + k * nummodesA;
1423 for (k = nummodesB; k < Q; ++k)
1425 signarray[arrayindx[j + k *
P]] = 0.0;
1426 maparray[arrayindx[j + k *
P]] = maparray[0];
1430 for (j = nummodesA; j <
P; ++j)
1432 for (k = 0; k < Q; ++k)
1434 signarray[arrayindx[j + k *
P]] = 0.0;
1435 maparray[arrayindx[j + k *
P]] = maparray[0];
1450 for (i = 3; i < Q; i += 2)
1452 for (j = 0; j <
P; j++)
1454 signarray[arrayindx[i *
P + j]] *= -1;
1458 for (i = 0; i <
P; i++)
1460 swap(maparray[i], maparray[i +
P]);
1461 swap(signarray[i], signarray[i +
P]);
1466 for (i = 0; i < Q; i++)
1468 for (j = 3; j <
P; j += 2)
1470 signarray[arrayindx[i *
P + j]] *= -1;
1474 for (i = 0; i < Q; i++)
1476 swap(maparray[i], maparray[i + Q]);
1477 swap(signarray[i], signarray[i + Q]);
1489 for (i = 0; i < Q; i++)
1491 for (j = 3; j <
P; j += 2)
1493 signarray[arrayindx[i *
P + j]] *= -1;
1497 for (i = 0; i < Q; i++)
1499 swap(maparray[i *
P], maparray[i *
P + 1]);
1500 swap(signarray[i *
P], signarray[i *
P + 1]);
1505 for (i = 3; i < Q; i += 2)
1507 for (j = 0; j <
P; j++)
1509 signarray[arrayindx[i *
P + j]] *= -1;
1513 for (i = 0; i <
P; i++)
1515 swap(maparray[i * Q], maparray[i * Q + 1]);
1516 swap(signarray[i * Q], signarray[i * Q + 1]);
1529 const int P =
m_base[0]->GetNumModes() - 1;
1530 const int Q =
m_base[1]->GetNumModes() - 1;
1531 const int R =
m_base[2]->GetNumModes() - 1;
1534 if (maparray.size() != nEdgeIntCoeffs)
1539 if (signarray.size() != nEdgeIntCoeffs)
1545 fill(signarray.data(), signarray.data() + nEdgeIntCoeffs, 1);
1555 for (i = 2; i <=
P; ++i)
1557 maparray[i - 2] =
GetMode(i, 0, 0);
1562 for (i = 2; i <= Q; ++i)
1564 maparray[i - 2] =
GetMode(1, i, 0);
1571 for (i = 2; i <=
P; ++i)
1573 maparray[i - 2] =
GetMode(i, 1, 0);
1580 for (i = 2; i <= Q; ++i)
1582 maparray[i - 2] =
GetMode(0, i, 0);
1587 for (i = 2; i <= R; ++i)
1589 maparray[i - 2] =
GetMode(0, 0, i);
1594 for (i = 1; i <= R - 1; ++i)
1596 maparray[i - 1] =
GetMode(1, 0, i);
1601 for (i = 1; i <= R - 1; ++i)
1603 maparray[i - 1] =
GetMode(1, 1, i);
1608 for (i = 2; i <= R; ++i)
1610 maparray[i - 2] =
GetMode(0, 1, i);
1615 for (i = 2; i <= Q; ++i)
1617 maparray[i - 2] =
GetMode(0, i, 1);
1622 ASSERTL0(
false,
"Edge not defined.");
1628 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1639 const int P =
m_base[0]->GetNumModes() - 1;
1640 const int Q =
m_base[1]->GetNumModes() - 1;
1641 const int R =
m_base[2]->GetNumModes() - 1;
1643 int p,
q, r, idx = 0;
1649 if (maparray.size() != nFaceIntCoeffs)
1654 if (signarray.size() != nFaceIntCoeffs)
1660 fill(signarray.data(), signarray.data() + nFaceIntCoeffs, 1);
1666 if (fid != 1 && fid != 3)
1679 for (i = 0; i < nummodesB; i++)
1681 for (j = 0; j < nummodesA; j++)
1685 arrayindx[i * nummodesA + j] = i * nummodesA + j;
1689 arrayindx[i * nummodesA + j] = j * nummodesB + i;
1698 for (
q = 2;
q <= Q; ++
q)
1700 for (
p = 2;
p <=
P; ++
p)
1702 maparray[arrayindx[(
q - 2) * nummodesA + (
p - 2)]] =
1709 for (
p = 2;
p <=
P; ++
p)
1711 for (r = 1; r <= R -
p; ++r)
1715 signarray[idx] =
p % 2 ? -1 : 1;
1717 maparray[idx++] =
GetMode(
p, 0, r);
1723 for (r = 1; r <= R - 1; ++r)
1725 for (
q = 2;
q <= Q; ++
q)
1727 maparray[arrayindx[(r - 1) * nummodesA + (
q - 2)]] =
1734 for (
p = 2;
p <=
P; ++
p)
1736 for (r = 1; r <= R -
p; ++r)
1740 signarray[idx] =
p % 2 ? -1 : 1;
1742 maparray[idx++] =
GetMode(
p, 1, r);
1748 for (r = 2; r <= R; ++r)
1750 for (
q = 2;
q <= Q; ++
q)
1752 maparray[arrayindx[(r - 2) * nummodesA + (
q - 2)]] =
1759 ASSERTL0(
false,
"Face interior map not available.");
1764 if (fid == 1 || fid == 3)
1776 for (i = 1; i < nummodesB; i += 2)
1778 for (j = 0; j < nummodesA; j++)
1780 signarray[arrayindx[i * nummodesA + j]] *= -1;
1786 for (i = 0; i < nummodesB; i++)
1788 for (j = 1; j < nummodesA; j += 2)
1790 signarray[arrayindx[i * nummodesA + j]] *= -1;
1803 for (i = 0; i < nummodesB; i++)
1805 for (j = 1; j < nummodesA; j += 2)
1807 signarray[arrayindx[i * nummodesA + j]] *= -1;
1813 for (i = 1; i < nummodesB; i += 2)
1815 for (j = 0; j < nummodesA; j++)
1817 signarray[arrayindx[i * nummodesA + j]] *= -1;
1839 int nq0 =
m_base[0]->GetNumPoints();
1840 int nq1 =
m_base[1]->GetNumPoints();
1841 int nq2 =
m_base[2]->GetNumPoints();
1851 nq = max(nq0, max(nq1, nq2));
1864 for (
int i = 0; i < nq; ++i)
1866 for (
int j = 0; j < nq; ++j)
1868 for (
int k = 0; k < nq - i; ++k, ++cnt)
1871 coords[cnt][0] = -1.0 + 2 * k / (
NekDouble)(nq - 1);
1872 coords[cnt][1] = -1.0 + 2 * j / (
NekDouble)(nq - 1);
1873 coords[cnt][2] = -1.0 + 2 * i / (
NekDouble)(nq - 1);
1878 for (
int i = 0; i < neq; ++i)
1882 I[0] =
m_base[0]->GetI(coll);
1883 I[1] =
m_base[1]->GetI(coll + 1);
1884 I[2] =
m_base[2]->GetI(coll + 2);
1888 for (
int k = 0; k < nq2; ++k)
1890 for (
int j = 0; j < nq1; ++j)
1893 fac = (I[1]->GetPtr())[j] * (I[2]->GetPtr())[k];
1897 Mat->GetRawPtr() + k * nq0 * nq1 * neq +
1940 int Q =
m_base[1]->GetNumModes() - 1;
1941 int R =
m_base[2]->GetNumModes() - 1;
1945 (Q + 1) * (
p * R + 1 -
1946 (
p - 2) * (
p - 1) / 2);
1954 int nquad0 =
m_base[0]->GetNumPoints();
1955 int nquad1 =
m_base[1]->GetNumPoints();
1956 int nquad2 =
m_base[2]->GetNumPoints();
1965 for (i = 0; i < nquad1 * nquad2; ++i)
1967 Vmath::Vmul(nquad0, inarray.data() + i * nquad0, 1, w0.data(), 1,
1968 outarray.data() + i * nquad0, 1);
1972 for (j = 0; j < nquad2; ++j)
1974 for (i = 0; i < nquad1; ++i)
1977 &outarray[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1987 case LibUtilities::eGaussRadauMAlpha1Beta0:
1988 for (i = 0; i < nquad2; ++i)
1991 &outarray[0] + i * nquad0 * nquad1, 1);
1996 for (i = 0; i < nquad2; ++i)
1998 Blas::Dscal(nquad0 * nquad1, 0.5 * (1 - z2[i]) * w2[i],
1999 &outarray[0] + i * nquad0 * nquad1, 1);
2009 int qa =
m_base[0]->GetNumPoints();
2010 int qb =
m_base[1]->GetNumPoints();
2011 int qc =
m_base[2]->GetNumPoints();
2012 int nmodes_a =
m_base[0]->GetNumModes();
2013 int nmodes_b =
m_base[1]->GetNumModes();
2014 int nmodes_c =
m_base[2]->GetNumModes();
2026 int i, j, k, cnt = 0;
2029 OrthoExp.
FwdTrans(array, orthocoeffs);
2039 for (i = 0; i < nmodes_a; ++i)
2041 for (j = 0; j < nmodes_b; ++j)
2044 pow((1.0 * i) / (nmodes_a - 1), cutoff * nmodes_a),
2045 pow((1.0 * j) / (nmodes_b - 1), cutoff * nmodes_b));
2047 for (k = 0; k < nmodes_c - i; ++k)
2050 std::max(fac1, pow((1.0 * k) / (nmodes_c - 1),
2051 cutoff * nmodes_c));
2053 orthocoeffs[cnt] *= SvvDiffCoeff * fac;
2069 max_abc = max(max_abc, 0);
2072 for (i = 0; i < nmodes_a; ++i)
2074 for (j = 0; j < nmodes_b; ++j)
2076 int maxij = max(i, j);
2078 for (k = 0; k < nmodes_c - i; ++k)
2080 int maxijk = max(maxij, k);
2102 int cutoff_a = (int)(SVVCutOff * nmodes_a);
2103 int cutoff_b = (int)(SVVCutOff * nmodes_b);
2104 int cutoff_c = (int)(SVVCutOff * nmodes_c);
2108 int nmodes = min(min(nmodes_a, nmodes_b), nmodes_c);
2109 NekDouble cutoff = min(min(cutoff_a, cutoff_b), cutoff_c);
2112 for (i = 0; i < nmodes_a; ++i)
2114 for (j = 0; j < nmodes_b; ++j)
2116 for (k = 0; k < nmodes_c - i; ++k)
2118 if (j >= cutoff || i + k >= cutoff)
2122 exp(-(i + k - nmodes) * (i + k - nmodes) /
2123 ((
NekDouble)((i + k - cutoff + epsilon) *
2124 (i + k - cutoff + epsilon)))) *
2125 exp(-(j - nmodes) * (j - nmodes) /
2127 (j - cutoff + epsilon)))));
2131 orthocoeffs[cnt] *= 0.0;
2140 OrthoExp.
BwdTrans(orthocoeffs, array);
2147 int nquad0 =
m_base[0]->GetNumPoints();
2148 int nquad1 =
m_base[1]->GetNumPoints();
2149 int nquad2 =
m_base[2]->GetNumPoints();
2150 int nqtot = nquad0 * nquad1 * nquad2;
2151 int nmodes0 =
m_base[0]->GetNumModes();
2152 int nmodes1 =
m_base[1]->GetNumModes();
2153 int nmodes2 =
m_base[2]->GetNumModes();
2154 int numMax = nmodes0;
2175 bortho0, bortho1, bortho2);
2178 OrthoPrismExp->FwdTrans(phys_tmp, coeff);
2181 for (u = 0; u < numMin; ++u)
2183 for (i = 0; i < numMin; ++i)
2186 tmp2 = coeff_tmp1 + cnt, 1);
2190 for (i = numMin; i < numMax; ++i)
2196 OrthoPrismExp->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.
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 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)
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.
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
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
MatrixType GetMatrixType() const
NekDouble GetConstFactor(const ConstFactorType &factor) const
bool ConstFactorExists(const ConstFactorType &factor) const
Class representing a prismatic element in reference space.
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey) override
void v_GetElmtTraceToTraceMap(const unsigned int fid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation faceOrient, int P, int Q) override
void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_GetEdgeInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2) override
LibUtilities::PointsKey v_GetTracePointsKey(const int i, const int j) const override
const LibUtilities::BasisKey v_GetTraceBasisKey(const int i, const int k, bool UseGLL=false) const override
NekDouble v_PhysEvalFirstDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) override
int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false) override
int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset) override
int v_GetTraceIntNcoeffs(const int i) const override
StdPrismExp(const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const LibUtilities::BasisKey &Bc)
int v_NumBndryCoeffs() const override
void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray) override
NekDouble v_PhysEvaluateBasis(const Array< OneD, const NekDouble > &coords, int mode) final
void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true) override
LibUtilities::ShapeType v_DetShapeType() const override
Return Shape of region, using ShapeType enum list; i.e. prism.
void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Calculate the inner product of inarray with respect to the basis B=base0*base1*base2 and put into out...
int v_GetNverts() const override
int v_NumDGBndryCoeffs() 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_GetCoords(Array< OneD, NekDouble > &xi_x, Array< OneD, NekDouble > &xi_y, Array< OneD, NekDouble > &xi_z) override
void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray) override
DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey) override
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
int v_GetNtraces() const 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
int GetMode(int I, int J, int K)
Compute the local mode number in the expansion for a particular tensorial combination.
int v_GetTraceNcoeffs(const int i) const 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_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_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_GetEdgeNcoeffs(const int i) const override
void v_LocCollapsedToLocCoord(const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi) override
void v_GetInteriorMap(Array< OneD, unsigned int > &outarray) override
void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_GetTraceCoeffMap(const unsigned int fid, Array< OneD, unsigned int > &maparray) override
void v_GetTraceInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2) override
void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey) override
int v_GetTraceNumPoints(const int i) const override
bool v_IsBoundaryInteriorExpansion() const override
int v_GetNedges() const override
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Inner product of inarray over region with respect to the object's default expansion basis; output in ...
void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta) override
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.
@ eModified_B
Principle Modified Functions .
@ P
Monomial polynomials .
@ eOrtho_A
Principle Orthogonal Functions .
@ eModified_C
Principle Modified Functions .
@ eGLL_Lagrange
Lagrange for SEM basis .
@ eOrtho_C
Principle Orthogonal Functions .
@ eOrtho_B
Principle Orthogonal Functions .
@ eModified_A
Principle Modified Functions .
static const NekDouble kNekZeroTol
std::shared_ptr< StdPrismExp > StdPrismExpSharedPtr
@ eFactorSVVDGKerDiffCoeff
@ eFactorSVVPowerKerDiffCoeff
LibUtilities::BasisKey EvaluateTriFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes, bool UseGLL)
const int kSVVDGFiltermodesmin
const int kSVVDGFiltermodesmax
LibUtilities::BasisKey EvaluateQuadFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
const NekDouble kSVVDGFilter[9][11]
@ ePhysInterpToEquiSpaced
@ 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 Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Svtvp (scalar times vector plus vector): z = alpha*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)
scalarT< T > sqrt(scalarT< T > in)