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.get() + mode * nquad2, nquad2,
284 inarray.get() + mode * nummodes1, nummodes2 - i, 0.0,
285 tmp0.get() + i * nquad2 * nummodes1, nquad2);
286 mode += nummodes2 - i;
291 for (i = 0; i < nummodes1; i++)
294 base2.get() + nquad2, 1,
295 tmp0.get() + nquad2 * (nummodes1 + i), 1);
299 for (i = 0; i < nummodes0; i++)
301 Blas::Dgemm(
'N',
'T', nquad1, nquad2, nummodes1, 1.0, base1.get(),
302 nquad1, tmp0.get() + i * nquad2 * nummodes1, nquad2, 0.0,
303 tmp1.get() + i * nquad2 * nquad1, nquad1);
306 Blas::Dgemm(
'N',
'T', nquad0, nquad2 * nquad1, nummodes0, 1.0, base0.get(),
307 nquad0, tmp1.get(), nquad2 * nquad1, 0.0, outarray.get(),
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.get(),
446 nquad0, base0.get(), nquad0, 0.0, tmp0.get(), nquad1 * nquad2);
449 Blas::Dgemm(
'T',
'N', nquad2 * order0, order1, nquad1, 1.0, tmp0.get(),
450 nquad1, base1.get(), nquad1, 0.0, tmp1.get(), nquad2 * order0);
453 for (mode = i = 0; i < order0; ++i)
455 Blas::Dgemm(
'T',
'N', order2 - i, order1, nquad2, 1.0,
456 base2.get() + mode * nquad2, nquad2,
457 tmp1.get() + i * nquad2, nquad2 * order0, 0.0,
458 outarray.get() + mode * order1, order2 - i);
466 for (i = 0; i < order1; ++i)
471 tmp1.get() + i * order0 * nquad2 + nquad2, 1);
491 ASSERTL0(dir >= 0 && dir <= 2,
"input dir is out of range");
494 int order0 =
m_base[0]->GetNumModes();
495 int order1 =
m_base[1]->GetNumModes();
496 int nquad0 =
m_base[0]->GetNumPoints();
497 int nquad1 =
m_base[1]->GetNumPoints();
498 int nquad2 =
m_base[2]->GetNumPoints();
508 for (i = 0; i < nquad0; ++i)
510 gfac0[i] = 0.5 * (1 + z0[i]);
514 for (i = 0; i < nquad2; ++i)
516 gfac2[i] = 2.0 / (1 - z2[i]);
522 for (i = 0; i < nquad2; ++i)
525 &inarray[0] + i * nquad0 * nquad1, 1,
526 &tmp0[0] + i * nquad0 * nquad1, 1);
537 m_base[2]->GetBdata(), tmp0, outarray, wsp,
true,
true,
true);
545 m_base[2]->GetBdata(), tmp0, outarray, wsp,
true,
true,
true);
554 for (i = 0; i < nquad1 * nquad2; ++i)
556 Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp0[0] + i * nquad0, 1,
557 &tmp0[0] + i * nquad0, 1);
562 m_base[2]->GetBdata(), tmp0, tmp1, wsp,
true,
true,
true);
567 m_base[2]->GetDbdata(), tmp0, outarray, wsp,
true,
true,
true);
597 eta[0] = 2.0 * (1.0 + xi[0]) / d2 - 1.0;
603 xi[0] = (1.0 + eta[0]) * (1.0 - eta[2]) * 0.5 - 1.0;
620 for (
int k = 0; k < Qz; ++k)
622 for (
int j = 0; j < Qy; ++j)
624 for (
int i = 0; i < Qx; ++i)
626 int s = i + Qx * (j + Qy * k);
627 xi_x[s] = (1.0 - eta_z[k]) * (1.0 + etaBar_x[i]) / 2.0 - 1.0;
638 std::array<NekDouble, 3> &firstOrderDerivs)
645 if ((1 - coll[2]) < 1e-5)
649 EphysDeriv2(totPoints);
650 PhysDeriv(inarray, EphysDeriv0, EphysDeriv1, EphysDeriv2);
653 I[0] =
GetBase()[0]->GetI(coll);
654 I[1] =
GetBase()[1]->GetI(coll + 1);
655 I[2] =
GetBase()[2]->GetI(coll + 2);
665 NekDouble dEta_bar1 = firstOrderDerivs[0];
668 firstOrderDerivs[0] = fac * dEta_bar1;
671 fac = 1.0 / (1.0 - coll[2]);
672 dEta_bar1 = fac * dEta_bar1;
676 firstOrderDerivs[2] += fac * dEta_bar1;
694 const int nm1 =
m_base[1]->GetNumModes();
695 const int nm2 =
m_base[2]->GetNumModes();
696 const int b = 2 * nm2 + 1;
698 const int mode0 = floor(0.5 * (b -
sqrt(b * b - 8.0 * mode / nm1)));
700 mode - nm1 * (mode0 * (nm2 - 1) + 1 - (mode0 - 2) * (mode0 - 1) / 2);
701 const int mode1 = tmp / (nm2 - mode0);
702 const int mode2 = tmp % (nm2 - mode0);
704 if (mode0 == 0 && mode2 == 1 &&
708 return StdExpansion::BaryEvaluateBasis<1>(coll[1], mode1) *
709 StdExpansion::BaryEvaluateBasis<2>(coll[2], mode2);
713 return StdExpansion::BaryEvaluateBasis<0>(coll[0], mode0) *
714 StdExpansion::BaryEvaluateBasis<1>(coll[1], mode1) *
715 StdExpansion::BaryEvaluateBasis<2>(coll[2], mode2);
722 int nummodes[3] = {
m_base[0]->GetNumModes(),
m_base[1]->GetNumModes(),
723 m_base[2]->GetNumModes()};
729 numModes0 = nummodes[0];
730 numModes1 = nummodes[1];
737 numModes0 = nummodes[1];
738 numModes1 = nummodes[2];
745 numModes0 = nummodes[0];
746 numModes1 = nummodes[2];
753 std::swap(numModes0, numModes1);
759 ASSERTL2(i >= 0 && i <= 8,
"edge id is out of range");
761 if (i == 0 || i == 2)
765 else if (i == 1 || i == 3 || i == 8)
807 "BasisType is not a boundary interior form");
810 "BasisType is not a boundary interior form");
813 "BasisType is not a boundary interior form");
815 int P =
m_base[0]->GetNumModes();
816 int Q =
m_base[1]->GetNumModes();
817 int R =
m_base[2]->GetNumModes();
826 "BasisType is not a boundary interior form");
829 "BasisType is not a boundary interior form");
832 "BasisType is not a boundary interior form");
834 int P =
m_base[0]->GetNumModes() - 1;
835 int Q =
m_base[1]->GetNumModes() - 1;
836 int R =
m_base[2]->GetNumModes() - 1;
838 return (
P + 1) * (Q + 1)
839 + 2 * (Q + 1) * (R + 1)
840 + 2 * (R + 1) +
P * (1 + 2 * R -
P);
845 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
850 else if (i == 1 || i == 3)
853 return Q + 1 + (
P * (1 + 2 * Q -
P)) / 2;
863 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
873 else if (i == 1 || i == 3)
875 return Pi * (2 * Ri - Pi - 1) / 2;
885 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
889 return m_base[0]->GetNumPoints() *
m_base[1]->GetNumPoints();
891 else if (i == 1 || i == 3)
893 return m_base[0]->GetNumPoints() *
m_base[2]->GetNumPoints();
897 return m_base[1]->GetNumPoints() *
m_base[2]->GetNumPoints();
904 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
905 ASSERTL2(j == 0 || j == 1,
"face direction is out of range");
909 return m_base[j]->GetPointsKey();
911 else if (i == 1 || i == 3)
913 return m_base[2 * j]->GetPointsKey();
917 return m_base[j + 1]->GetPointsKey();
924 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
925 ASSERTL2(k >= 0 && k <= 1,
"basis key id is out of range");
933 m_base[k]->GetNumModes());
940 m_base[k + 1]->GetNumModes());
947 m_base[2 * k]->GetNumModes());
957 const std::vector<unsigned int> &nummodes,
int &modes_offset)
960 nummodes[modes_offset], nummodes[modes_offset + 1],
961 nummodes[modes_offset + 2]);
983 "Mapping not defined for this type of basis");
987 if (useCoeffPacking ==
true)
1010 ASSERTL0(
false,
"local vertex id must be between 0 and 5");
1036 ASSERTL0(
false,
"local vertex id must be between 0 and 5");
1047 "BasisType is not a boundary interior form");
1050 "BasisType is not a boundary interior form");
1053 "BasisType is not a boundary interior form");
1055 int P =
m_base[0]->GetNumModes() - 1,
p;
1056 int Q =
m_base[1]->GetNumModes() - 1,
q;
1057 int R =
m_base[2]->GetNumModes() - 1, r;
1061 if (outarray.size() != nIntCoeffs)
1069 for (
p = 2;
p <=
P; ++
p)
1071 for (
q = 2;
q <= Q; ++
q)
1073 for (r = 1; r <= R -
p; ++r)
1085 "BasisType is not a boundary interior form");
1088 "BasisType is not a boundary interior form");
1091 "BasisType is not a boundary interior form");
1093 int P =
m_base[0]->GetNumModes() - 1,
p;
1094 int Q =
m_base[1]->GetNumModes() - 1,
q;
1095 int R =
m_base[2]->GetNumModes() - 1, r;
1100 if (maparray.size() != nBnd)
1106 for (
p = 0;
p <=
P; ++
p)
1111 for (
q = 0;
q <= Q; ++
q)
1113 for (r = 0; r <= R -
p; ++r)
1123 for (
q = 0;
q <= Q; ++
q)
1127 for (r = 0; r <= R -
p; ++r)
1145 "Method only implemented if BasisType is identical"
1146 "in x and y directions");
1149 "Method only implemented for Modified_A BasisType"
1150 "(x and y direction) and Modified_B BasisType (z "
1152 int p,
q, r, idx = 0;
1159 Q =
m_base[1]->GetNumModes();
1164 Q =
m_base[2]->GetNumModes();
1169 Q =
m_base[2]->GetNumModes();
1172 ASSERTL0(
false,
"fid must be between 0 and 4");
1175 if (maparray.size() !=
P * Q)
1185 for (
q = 0;
q < Q; ++
q)
1187 for (
p = 0;
p <
P; ++
p)
1194 for (
p = 0;
p <
P; ++
p)
1196 for (r = 0; r < Q -
p; ++r)
1198 maparray[idx++] =
GetMode(
p, 0, r);
1203 for (
q = 0;
q <
P; ++
q)
1207 for (
q = 0;
q <
P; ++
q)
1211 for (r = 1; r < Q - 1; ++r)
1213 for (
q = 0;
q <
P; ++
q)
1215 maparray[(r + 1) *
P +
q] =
GetMode(1,
q, r);
1220 for (
p = 0;
p <
P; ++
p)
1222 for (r = 0; r < Q -
p; ++r)
1224 maparray[idx++] =
GetMode(
p, 1, r);
1229 for (r = 0; r < Q; ++r)
1231 for (
q = 0;
q <
P; ++
q)
1238 ASSERTL0(
false,
"Face to element map unavailable.");
1248 "Method only implemented if BasisType is identical"
1249 "in x and y directions");
1252 "Method only implemented for Modified_A BasisType"
1253 "(x and y direction) and Modified_B BasisType (z "
1256 int i, j, k,
p, r, nFaceCoeffs, idx = 0;
1257 int nummodesA = 0, nummodesB = 0;
1262 nummodesA =
m_base[0]->GetNumModes();
1263 nummodesB =
m_base[1]->GetNumModes();
1267 nummodesA =
m_base[0]->GetNumModes();
1268 nummodesB =
m_base[2]->GetNumModes();
1272 nummodesA =
m_base[1]->GetNumModes();
1273 nummodesB =
m_base[2]->GetNumModes();
1276 ASSERTL0(
false,
"fid must be between 0 and 4");
1285 else if (fid == 1 || fid == 3)
1287 nFaceCoeffs =
P * (2 * Q -
P + 1) / 2;
1291 nFaceCoeffs =
P * Q;
1295 if (maparray.size() != nFaceCoeffs)
1300 if (signarray.size() != nFaceCoeffs)
1306 fill(signarray.get(), signarray.get() + nFaceCoeffs, 1);
1309 int minPA = min(nummodesA,
P);
1310 int minQB = min(nummodesB, Q);
1312 if (fid == 1 || fid == 3)
1319 for (j = 0; j < minPA; ++j)
1322 for (k = 0; k < minQB - j; ++k, ++cnt)
1324 maparray[idx++] = cnt;
1327 cnt += nummodesB - minQB;
1330 for (k = nummodesB - j; k < Q - j; ++k)
1332 signarray[idx] = 0.0;
1333 maparray[idx++] = maparray[0];
1337 for (j = minPA; j < nummodesA; ++j)
1340 for (k = 0; k < minQB-j; ++k, ++cnt)
1342 maparray[idx++] = cnt;
1345 cnt += nummodesB-minQB;
1348 for (k = nummodesB-j; k < Q-j; ++k)
1350 signarray[idx] = 0.0;
1351 maparray[idx++] = maparray[0];
1355 for (j = nummodesA; j <
P; ++j)
1357 for (k = 0; k < Q - j; ++k)
1359 signarray[idx] = 0.0;
1360 maparray[idx++] = maparray[0];
1369 for (
p = 0;
p <
P; ++
p)
1371 for (r = 0; r < Q -
p; ++r, idx++)
1375 signarray[idx] =
p % 2 ? -1 : 1;
1380 swap(maparray[0], maparray[Q]);
1381 for (i = 1; i < Q - 1; ++i)
1383 swap(maparray[i + 1], maparray[Q + i]);
1393 for (i = 0; i < Q; i++)
1395 for (j = 0; j <
P; j++)
1399 arrayindx[i *
P + j] = i *
P + j;
1403 arrayindx[i *
P + j] = j * Q + i;
1410 for (j = 0; j <
P; ++j)
1413 for (k = 0; k < Q; k++)
1415 maparray[arrayindx[j + k *
P]] = j + k * nummodesA;
1418 for (k = nummodesB; k < Q; ++k)
1420 signarray[arrayindx[j + k *
P]] = 0.0;
1421 maparray[arrayindx[j + k *
P]] = maparray[0];
1425 for (j = nummodesA; j <
P; ++j)
1427 for (k = 0; k < Q; ++k)
1429 signarray[arrayindx[j + k *
P]] = 0.0;
1430 maparray[arrayindx[j + k *
P]] = maparray[0];
1445 for (i = 3; i < Q; i += 2)
1447 for (j = 0; j <
P; j++)
1449 signarray[arrayindx[i *
P + j]] *= -1;
1453 for (i = 0; i <
P; i++)
1455 swap(maparray[i], maparray[i +
P]);
1456 swap(signarray[i], signarray[i +
P]);
1461 for (i = 0; i < Q; i++)
1463 for (j = 3; j <
P; j += 2)
1465 signarray[arrayindx[i *
P + j]] *= -1;
1469 for (i = 0; i < Q; i++)
1471 swap(maparray[i], maparray[i + Q]);
1472 swap(signarray[i], signarray[i + Q]);
1484 for (i = 0; i < Q; i++)
1486 for (j = 3; j <
P; j += 2)
1488 signarray[arrayindx[i *
P + j]] *= -1;
1492 for (i = 0; i < Q; i++)
1494 swap(maparray[i *
P], maparray[i *
P + 1]);
1495 swap(signarray[i *
P], signarray[i *
P + 1]);
1500 for (i = 3; i < Q; i += 2)
1502 for (j = 0; j <
P; j++)
1504 signarray[arrayindx[i *
P + j]] *= -1;
1508 for (i = 0; i <
P; i++)
1510 swap(maparray[i * Q], maparray[i * Q + 1]);
1511 swap(signarray[i * Q], signarray[i * Q + 1]);
1524 const int P =
m_base[0]->GetNumModes() - 1;
1525 const int Q =
m_base[1]->GetNumModes() - 1;
1526 const int R =
m_base[2]->GetNumModes() - 1;
1529 if (maparray.size() != nEdgeIntCoeffs)
1534 if (signarray.size() != nEdgeIntCoeffs)
1540 fill(signarray.get(), signarray.get() + nEdgeIntCoeffs, 1);
1550 for (i = 2; i <=
P; ++i)
1552 maparray[i - 2] =
GetMode(i, 0, 0);
1557 for (i = 2; i <= Q; ++i)
1559 maparray[i - 2] =
GetMode(1, i, 0);
1566 for (i = 2; i <=
P; ++i)
1568 maparray[i - 2] =
GetMode(i, 1, 0);
1575 for (i = 2; i <= Q; ++i)
1577 maparray[i - 2] =
GetMode(0, i, 0);
1582 for (i = 2; i <= R; ++i)
1584 maparray[i - 2] =
GetMode(0, 0, i);
1589 for (i = 1; i <= R - 1; ++i)
1591 maparray[i - 1] =
GetMode(1, 0, i);
1596 for (i = 1; i <= R - 1; ++i)
1598 maparray[i - 1] =
GetMode(1, 1, i);
1603 for (i = 2; i <= R; ++i)
1605 maparray[i - 2] =
GetMode(0, 1, i);
1610 for (i = 2; i <= Q; ++i)
1612 maparray[i - 2] =
GetMode(0, i, 1);
1617 ASSERTL0(
false,
"Edge not defined.");
1623 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1634 const int P =
m_base[0]->GetNumModes() - 1;
1635 const int Q =
m_base[1]->GetNumModes() - 1;
1636 const int R =
m_base[2]->GetNumModes() - 1;
1638 int p,
q, r, idx = 0;
1644 if (maparray.size() != nFaceIntCoeffs)
1649 if (signarray.size() != nFaceIntCoeffs)
1655 fill(signarray.get(), signarray.get() + nFaceIntCoeffs, 1);
1661 if (fid != 1 && fid != 3)
1674 for (i = 0; i < nummodesB; i++)
1676 for (j = 0; j < nummodesA; j++)
1680 arrayindx[i * nummodesA + j] = i * nummodesA + j;
1684 arrayindx[i * nummodesA + j] = j * nummodesB + i;
1693 for (
q = 2;
q <= Q; ++
q)
1695 for (
p = 2;
p <=
P; ++
p)
1697 maparray[arrayindx[(
q - 2) * nummodesA + (
p - 2)]] =
1704 for (
p = 2;
p <=
P; ++
p)
1706 for (r = 1; r <= R -
p; ++r)
1710 signarray[idx] =
p % 2 ? -1 : 1;
1712 maparray[idx++] =
GetMode(
p, 0, r);
1718 for (r = 1; r <= R - 1; ++r)
1720 for (
q = 2;
q <= Q; ++
q)
1722 maparray[arrayindx[(r - 1) * nummodesA + (
q - 2)]] =
1729 for (
p = 2;
p <=
P; ++
p)
1731 for (r = 1; r <= R -
p; ++r)
1735 signarray[idx] =
p % 2 ? -1 : 1;
1737 maparray[idx++] =
GetMode(
p, 1, r);
1743 for (r = 2; r <= R; ++r)
1745 for (
q = 2;
q <= Q; ++
q)
1747 maparray[arrayindx[(r - 2) * nummodesA + (
q - 2)]] =
1754 ASSERTL0(
false,
"Face interior map not available.");
1759 if (fid == 1 || fid == 3)
1771 for (i = 1; i < nummodesB; i += 2)
1773 for (j = 0; j < nummodesA; j++)
1775 signarray[arrayindx[i * nummodesA + j]] *= -1;
1781 for (i = 0; i < nummodesB; i++)
1783 for (j = 1; j < nummodesA; j += 2)
1785 signarray[arrayindx[i * nummodesA + j]] *= -1;
1798 for (i = 0; i < nummodesB; i++)
1800 for (j = 1; j < nummodesA; j += 2)
1802 signarray[arrayindx[i * nummodesA + j]] *= -1;
1808 for (i = 1; i < nummodesB; i += 2)
1810 for (j = 0; j < nummodesA; j++)
1812 signarray[arrayindx[i * nummodesA + j]] *= -1;
1834 int nq0 =
m_base[0]->GetNumPoints();
1835 int nq1 =
m_base[1]->GetNumPoints();
1836 int nq2 =
m_base[2]->GetNumPoints();
1846 nq = max(nq0, max(nq1, nq2));
1859 for (
int i = 0; i < nq; ++i)
1861 for (
int j = 0; j < nq; ++j)
1863 for (
int k = 0; k < nq - i; ++k, ++cnt)
1866 coords[cnt][0] = -1.0 + 2 * k / (
NekDouble)(nq - 1);
1867 coords[cnt][1] = -1.0 + 2 * j / (
NekDouble)(nq - 1);
1868 coords[cnt][2] = -1.0 + 2 * i / (
NekDouble)(nq - 1);
1873 for (
int i = 0; i < neq; ++i)
1877 I[0] =
m_base[0]->GetI(coll);
1878 I[1] =
m_base[1]->GetI(coll + 1);
1879 I[2] =
m_base[2]->GetI(coll + 2);
1883 for (
int k = 0; k < nq2; ++k)
1885 for (
int j = 0; j < nq1; ++j)
1888 fac = (I[1]->GetPtr())[j] * (I[2]->GetPtr())[k];
1892 Mat->GetRawPtr() + k * nq0 * nq1 * neq +
1935 int Q =
m_base[1]->GetNumModes() - 1;
1936 int R =
m_base[2]->GetNumModes() - 1;
1940 (Q + 1) * (
p * R + 1 -
1941 (
p - 2) * (
p - 1) / 2);
1949 int nquad0 =
m_base[0]->GetNumPoints();
1950 int nquad1 =
m_base[1]->GetNumPoints();
1951 int nquad2 =
m_base[2]->GetNumPoints();
1960 for (i = 0; i < nquad1 * nquad2; ++i)
1962 Vmath::Vmul(nquad0, inarray.get() + i * nquad0, 1, w0.get(), 1,
1963 outarray.get() + i * nquad0, 1);
1967 for (j = 0; j < nquad2; ++j)
1969 for (i = 0; i < nquad1; ++i)
1972 &outarray[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1982 case LibUtilities::eGaussRadauMAlpha1Beta0:
1983 for (i = 0; i < nquad2; ++i)
1986 &outarray[0] + i * nquad0 * nquad1, 1);
1991 for (i = 0; i < nquad2; ++i)
1993 Blas::Dscal(nquad0 * nquad1, 0.5 * (1 - z2[i]) * w2[i],
1994 &outarray[0] + i * nquad0 * nquad1, 1);
2004 int qa =
m_base[0]->GetNumPoints();
2005 int qb =
m_base[1]->GetNumPoints();
2006 int qc =
m_base[2]->GetNumPoints();
2007 int nmodes_a =
m_base[0]->GetNumModes();
2008 int nmodes_b =
m_base[1]->GetNumModes();
2009 int nmodes_c =
m_base[2]->GetNumModes();
2021 int i, j, k, cnt = 0;
2024 OrthoExp.
FwdTrans(array, orthocoeffs);
2034 for (
int i = 0; i < nmodes_a; ++i)
2036 for (
int j = 0; j < nmodes_b; ++j)
2039 pow((1.0 * i) / (nmodes_a - 1), cutoff * nmodes_a),
2040 pow((1.0 * j) / (nmodes_b - 1), cutoff * nmodes_b));
2042 for (
int k = 0; k < nmodes_c - i; ++k)
2045 std::max(fac1, pow((1.0 * k) / (nmodes_c - 1),
2046 cutoff * nmodes_c));
2048 orthocoeffs[cnt] *= SvvDiffCoeff * fac;
2064 max_abc = max(max_abc, 0);
2067 for (
int i = 0; i < nmodes_a; ++i)
2069 for (
int j = 0; j < nmodes_b; ++j)
2071 int maxij = max(i, j);
2073 for (
int k = 0; k < nmodes_c - i; ++k)
2075 int maxijk = max(maxij, k);
2097 int cutoff_a = (int)(SVVCutOff * nmodes_a);
2098 int cutoff_b = (int)(SVVCutOff * nmodes_b);
2099 int cutoff_c = (int)(SVVCutOff * nmodes_c);
2103 int nmodes = min(min(nmodes_a, nmodes_b), nmodes_c);
2104 NekDouble cutoff = min(min(cutoff_a, cutoff_b), cutoff_c);
2107 for (i = 0; i < nmodes_a; ++i)
2109 for (j = 0; j < nmodes_b; ++j)
2111 for (k = 0; k < nmodes_c - i; ++k)
2113 if (j >= cutoff || i + k >= cutoff)
2117 exp(-(i + k - nmodes) * (i + k - nmodes) /
2118 ((
NekDouble)((i + k - cutoff + epsilon) *
2119 (i + k - cutoff + epsilon)))) *
2120 exp(-(j - nmodes) * (j - nmodes) /
2122 (j - cutoff + epsilon)))));
2126 orthocoeffs[cnt] *= 0.0;
2135 OrthoExp.
BwdTrans(orthocoeffs, array);
2142 int nquad0 =
m_base[0]->GetNumPoints();
2143 int nquad1 =
m_base[1]->GetNumPoints();
2144 int nquad2 =
m_base[2]->GetNumPoints();
2145 int nqtot = nquad0 * nquad1 * nquad2;
2146 int nmodes0 =
m_base[0]->GetNumModes();
2147 int nmodes1 =
m_base[1]->GetNumModes();
2148 int nmodes2 =
m_base[2]->GetNumModes();
2149 int numMax = nmodes0;
2170 bortho0, bortho1, bortho2);
2173 OrthoPrismExp->FwdTrans(phys_tmp, coeff);
2176 for (u = 0; u < numMin; ++u)
2178 for (i = 0; i < numMin; ++i)
2181 tmp2 = coeff_tmp1 + cnt, 1);
2185 for (i = numMin; i < numMax; ++i)
2191 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
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
const LibUtilities::BasisKey v_GetTraceBasisKey(const int i, const int k) const 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...
NekDouble v_PhysEvaluate(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) 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_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
LibUtilities::BasisKey EvaluateTriFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
@ eFactorSVVDGKerDiffCoeff
@ eFactorSVVPowerKerDiffCoeff
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)