46 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
49 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
53 "order in 'a' direction is higher than order "
56 "order in 'a' direction is higher than order "
59 "order in 'b' direction is higher than order "
88 int Q0 =
m_base[0]->GetNumPoints();
89 int Q1 =
m_base[1]->GetNumPoints();
90 int Q2 =
m_base[2]->GetNumPoints();
91 int Qtot = Q0 * Q1 * Q2;
98 bool Do_2 = (out_dxi2.size() > 0) ?
true :
false;
99 bool Do_1 = (out_dxi1.size() > 0) ?
true :
false;
116 eta_0 =
m_base[0]->GetZ();
117 eta_1 =
m_base[1]->GetZ();
118 eta_2 =
m_base[2]->GetZ();
124 for (
int k = 0; k < Q2; ++k)
126 for (
int j = 0; j < Q1; ++j, dEta0 += Q0)
128 Vmath::Smul(Q0, 2.0 / (1.0 - eta_1[j]), dEta0, 1, dEta0, 1);
130 fac = 1.0 / (1.0 - eta_2[k]);
131 Vmath::Smul(Q0 * Q1, fac, &out_dEta0[0] + k * Q0 * Q1, 1,
132 &out_dEta0[0] + k * Q0 * Q1, 1);
135 if (out_dxi0.size() > 0)
147 for (
int k = 0; k < Q1 * Q2; ++k)
149 Vmath::Vmul(Q0, &Fac0[0], 1, &out_dEta0[0] + k * Q0, 1,
150 &out_dEta0[0] + k * Q0, 1);
153 for (
int k = 0; k < Q2; ++k)
156 &out_dEta1[0] + k * Q0 * Q1, 1,
157 &out_dEta1[0] + k * Q0 * Q1, 1);
164 Vmath::Vadd(Qtot, out_dEta0, 1, out_dEta1, 1, out_dxi1, 1);
171 for (
int k = 0; k < Q2; ++k)
173 for (
int j = 0; j < Q1; ++j, dEta1 += Q0)
175 Vmath::Smul(Q0, (1.0 + eta_1[j]) / 2.0, dEta1, 1, dEta1, 1);
182 Vmath::Vadd(Qtot, out_dEta0, 1, out_dEta1, 1, out_dxi2, 1);
183 Vmath::Vadd(Qtot, out_dEta2, 1, out_dxi2, 1, out_dxi2, 1);
220 ASSERTL1(
false,
"input dir is out of range");
271 "Basis[1] is not a general tensor type");
275 "Basis[2] is not a general tensor type");
277 if (
m_base[0]->Collocation() &&
m_base[1]->Collocation() &&
282 inarray, 1, outarray, 1);
296 int nquad1 =
m_base[1]->GetNumPoints();
297 int nquad2 =
m_base[2]->GetNumPoints();
298 int order0 =
m_base[0]->GetNumModes();
299 int order1 =
m_base[1]->GetNumModes();
302 nquad2 * nquad1 * order0);
305 m_base[2]->GetBdata(), inarray, outarray, wsp,
true,
328 [[maybe_unused]]
bool doCheckCollDir0,
329 [[maybe_unused]]
bool doCheckCollDir1,
330 [[maybe_unused]]
bool doCheckCollDir2)
332 int nquad0 =
m_base[0]->GetNumPoints();
333 int nquad1 =
m_base[1]->GetNumPoints();
334 int nquad2 =
m_base[2]->GetNumPoints();
336 int order0 =
m_base[0]->GetNumModes();
337 int order1 =
m_base[1]->GetNumModes();
338 int order2 =
m_base[2]->GetNumModes();
342 tmp + nquad2 * order0 * (2 * order1 - order0 + 1) / 2;
344 int i, j, mode, mode1, cnt;
347 mode = mode1 = cnt = 0;
348 for (i = 0; i < order0; ++i)
350 for (j = 0; j < order1 - i; ++j, ++cnt)
353 base2.data() + mode * nquad2, nquad2,
354 inarray.data() + mode1, 1, 0.0,
355 tmp.data() + cnt * nquad2, 1);
356 mode += order2 - i - j;
357 mode1 += order2 - i - j;
360 for (j = order1 - i; j < order2 - i; ++j)
362 mode += order2 - i - j;
372 Blas::Daxpy(nquad2, inarray[1], base2.data() + nquad2, 1,
373 &tmp[0] + nquad2, 1);
376 Blas::Daxpy(nquad2, inarray[1], base2.data() + nquad2, 1,
377 &tmp[0] + order1 * nquad2, 1);
382 for (i = 0; i < order0; ++i)
384 Blas::Dgemm(
'N',
'T', nquad1, nquad2, order1 - i, 1.0,
385 base1.data() + mode * nquad1, nquad1,
386 tmp.data() + mode * nquad2, nquad2, 0.0,
387 tmp1.data() + i * nquad1 * nquad2, nquad1);
398 for (i = 0; i < nquad2; ++i)
400 Blas::Daxpy(nquad1, tmp[nquad2 + i], base1.data() + nquad1, 1,
401 &tmp1[nquad1 * nquad2] + i * nquad1, 1);
406 Blas::Dgemm(
'N',
'T', nquad0, nquad1 * nquad2, order0, 1.0, base0.data(),
407 nquad0, tmp1.data(), nquad1 * nquad2, 0.0, outarray.data(),
429 out = (*matsys) * in;
471 "Basis[1] is not a general tensor type");
475 "Basis[2] is not a general tensor type");
477 if (
m_base[0]->Collocation() &&
m_base[1]->Collocation())
497 int nquad0 =
m_base[0]->GetNumPoints();
498 int nquad1 =
m_base[1]->GetNumPoints();
499 int nquad2 =
m_base[2]->GetNumPoints();
500 int order0 =
m_base[0]->GetNumModes();
501 int order1 =
m_base[1]->GetNumModes();
504 nquad2 * order0 * (2 * order1 - order0 + 1) / 2);
506 if (multiplybyweights)
513 tmp, outarray, wsp,
true,
true,
true);
519 inarray, outarray, wsp,
true,
true,
true);
529 [[maybe_unused]]
bool doCheckCollDir0,
530 [[maybe_unused]]
bool doCheckCollDir1,
531 [[maybe_unused]]
bool doCheckCollDir2)
533 int nquad0 =
m_base[0]->GetNumPoints();
534 int nquad1 =
m_base[1]->GetNumPoints();
535 int nquad2 =
m_base[2]->GetNumPoints();
537 int order0 =
m_base[0]->GetNumModes();
538 int order1 =
m_base[1]->GetNumModes();
539 int order2 =
m_base[2]->GetNumModes();
544 int i, j, mode, mode1, cnt;
547 Blas::Dgemm(
'T',
'N', nquad1 * nquad2, order0, nquad0, 1.0, inarray.data(),
548 nquad0, base0.data(), nquad0, 0.0, tmp1.data(),
552 for (mode = i = 0; i < order0; ++i)
554 Blas::Dgemm(
'T',
'N', nquad2, order1 - i, nquad1, 1.0,
555 tmp1.data() + i * nquad1 * nquad2, nquad1,
556 base1.data() + mode * nquad1, nquad1, 0.0,
557 tmp2.data() + mode * nquad2, nquad2);
566 Blas::Dgemv(
'T', nquad1, nquad2, 1.0, tmp1.data() + nquad1 * nquad2,
567 nquad1, base1.data() + nquad1, 1, 1.0, tmp2.data() + nquad2,
572 mode = mode1 = cnt = 0;
573 for (i = 0; i < order0; ++i)
575 for (j = 0; j < order1 - i; ++j, ++cnt)
578 base2.data() + mode * nquad2, nquad2,
579 tmp2.data() + cnt * nquad2, 1, 0.0,
580 outarray.data() + mode1, 1);
581 mode += order2 - i - j;
582 mode1 += order2 - i - j;
585 for (j = order1 - i; j < order2 - i; ++j)
587 mode += order2 - i - j;
597 Blas::Ddot(nquad2, base2.data() + nquad2, 1, &tmp2[nquad2], 1);
600 outarray[1] +=
Blas::Ddot(nquad2, base2.data() + nquad2, 1,
601 &tmp2[nquad2 * order1], 1);
623 int nquad0 =
m_base[0]->GetNumPoints();
624 int nquad1 =
m_base[1]->GetNumPoints();
625 int nquad2 =
m_base[2]->GetNumPoints();
626 int nqtot = nquad0 * nquad1 * nquad2;
627 int nmodes0 =
m_base[0]->GetNumModes();
628 int nmodes1 =
m_base[1]->GetNumModes();
629 int wspsize = nquad0 + nquad1 + nquad2 + max(nqtot,
m_ncoeffs) +
630 nquad1 * nquad2 * nmodes0 +
631 nquad2 * nmodes0 * (2 * nmodes1 - nmodes0 + 1) / 2;
644 for (i = 0; i < nquad0; ++i)
646 gfac0[i] = 0.5 * (1 + z0[i]);
650 for (i = 0; i < nquad1; ++i)
652 gfac1[i] = 2.0 / (1 - z1[i]);
656 for (i = 0; i < nquad2; ++i)
658 gfac2[i] = 2.0 / (1 - z2[i]);
662 for (i = 0; i < nquad1 * nquad2; ++i)
664 Vmath::Smul(nquad0, gfac1[i % nquad1], &inarray[0] + i * nquad0, 1,
665 &tmp0[0] + i * nquad0, 1);
667 for (i = 0; i < nquad2; ++i)
669 Vmath::Smul(nquad0 * nquad1, gfac2[i], &tmp0[0] + i * nquad0 * nquad1,
670 1, &tmp0[0] + i * nquad0 * nquad1, 1);
681 m_base[2]->GetBdata(), tmp0, outarray, wsp,
false,
true,
true);
688 for (i = 0; i < nquad1 * nquad2; ++i)
690 Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp0[0] + i * nquad0, 1,
691 &tmp0[0] + i * nquad0, 1);
696 m_base[2]->GetBdata(), tmp0, tmp3, wsp,
false,
true,
true);
698 for (i = 0; i < nquad2; ++i)
701 &inarray[0] + i * nquad0 * nquad1, 1,
702 &tmp0[0] + i * nquad0 * nquad1, 1);
707 m_base[2]->GetBdata(), tmp0, outarray, wsp,
true,
false,
true);
716 for (i = 0; i < nquad1; ++i)
718 gfac1[i] = (1 + z1[i]) / 2;
721 for (i = 0; i < nquad1 * nquad2; ++i)
723 Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp0[0] + i * nquad0, 1,
724 &tmp0[0] + i * nquad0, 1);
728 m_base[2]->GetBdata(), tmp0, tmp3, wsp,
false,
true,
true);
730 for (i = 0; i < nquad2; ++i)
733 &inarray[0] + i * nquad0 * nquad1, 1,
734 &tmp0[0] + i * nquad0 * nquad1, 1);
736 for (i = 0; i < nquad1 * nquad2; ++i)
738 Vmath::Smul(nquad0, gfac1[i % nquad1], &tmp0[0] + i * nquad0, 1,
739 &tmp0[0] + i * nquad0, 1);
744 m_base[2]->GetBdata(), tmp0, tmp4, wsp,
true,
false,
true);
749 m_base[2]->GetDbdata(), tmp0, outarray, wsp,
true,
true,
false);
759 ASSERTL1(
false,
"input dir is out of range");
796 eta[0] = 2.0 * (1.0 + xi[0]) / d12 - 1.0;
797 eta[1] = 2.0 * (1.0 + xi[1]) / d2 - 1.0;
805 xi[1] = (1.0 + eta[1]) * (1.0 - xi[2]) * 0.5 - 1.0;
806 xi[0] = (1.0 + eta[0]) * (-xi[1] - xi[2]) * 0.5 - 1.0;
822 const int nm1 =
m_base[1]->GetNumModes();
823 const int nm2 =
m_base[2]->GetNumModes();
825 const int b = 2 * nm2 + 1;
826 const int mode0 = floor(0.5 * (b -
sqrt(b * b - 8.0 * mode / nm1)));
828 mode - nm1 * (mode0 * (nm2 - 1) + 1 - (mode0 - 2) * (mode0 - 1) / 2);
829 const int mode1 = tmp / (nm2 - mode0);
830 const int mode2 = tmp % (nm2 - mode0);
839 return StdExpansion::BaryEvaluateBasis<2>(coll[2], 1);
841 else if (mode0 == 0 && mode2 == 1)
843 return StdExpansion::BaryEvaluateBasis<1>(coll[1], 0) *
844 StdExpansion::BaryEvaluateBasis<2>(coll[2], 1);
846 else if (mode0 == 1 && mode1 == 1 && mode2 == 0)
848 return StdExpansion::BaryEvaluateBasis<0>(coll[0], 0) *
849 StdExpansion::BaryEvaluateBasis<1>(coll[1], 1);
853 return StdExpansion::BaryEvaluateBasis<0>(coll[0], mode0) *
854 StdExpansion::BaryEvaluateBasis<1>(coll[1], mode1) *
855 StdExpansion::BaryEvaluateBasis<2>(coll[2], mode2);
861 std::array<NekDouble, 3> &firstOrderDerivs)
868 if ((1 - coll[1]) < 1e-5 || (1 - coll[2]) < 1e-5)
872 EphysDeriv2(totPoints);
873 PhysDeriv(inarray, EphysDeriv0, EphysDeriv1, EphysDeriv2);
876 I[0] =
GetBase()[0]->GetI(coll);
877 I[1] =
GetBase()[1]->GetI(coll + 1);
878 I[2] =
GetBase()[2]->GetI(coll + 2);
886 std::array<NekDouble, 3> interDeriv;
890 NekDouble temp = 2.0 / ((1 - coll[1]) * (1 - coll[2]));
891 interDeriv[0] *= temp;
894 firstOrderDerivs[0] = 2 * interDeriv[0];
901 interDeriv[0] *= fac0;
904 fac0 = 2 / (1 - coll[2]);
905 interDeriv[1] *= fac0;
909 firstOrderDerivs[1] = interDeriv[0] + interDeriv[1];
912 fac0 = (1 + coll[1]) / 2;
913 interDeriv[1] *= fac0;
918 firstOrderDerivs[2] = interDeriv[0] + interDeriv[1] + interDeriv[2];
927 int nummodes[3] = {
m_base[0]->GetNumModes(),
m_base[1]->GetNumModes(),
928 m_base[2]->GetNumModes()};
933 numModes0 = nummodes[0];
934 numModes1 = nummodes[1];
939 numModes0 = nummodes[0];
940 numModes1 = nummodes[2];
946 numModes0 = nummodes[1];
947 numModes1 = nummodes[2];
981 "BasisType is not a boundary interior form");
984 "BasisType is not a boundary interior form");
987 "BasisType is not a boundary interior form");
989 int P =
m_base[0]->GetNumModes();
990 int Q =
m_base[1]->GetNumModes();
991 int R =
m_base[2]->GetNumModes();
1000 "BasisType is not a boundary interior form");
1003 "BasisType is not a boundary interior form");
1006 "BasisType is not a boundary interior form");
1008 int P =
m_base[0]->GetNumModes() - 1;
1009 int Q =
m_base[1]->GetNumModes() - 1;
1010 int R =
m_base[2]->GetNumModes() - 1;
1012 return (Q + 1) +
P * (1 + 2 * Q -
P) / 2
1013 + (R + 1) +
P * (1 + 2 * R -
P) / 2
1014 + 2 * (R + 1) + Q * (1 + 2 * R - Q);
1019 ASSERTL2((i >= 0) && (i <= 3),
"face id is out of range");
1020 int nFaceCoeffs = 0;
1021 int nummodesA, nummodesB,
P, Q;
1027 else if ((i == 1) || (i == 2))
1039 nFaceCoeffs = Q + 1 + (
P * (1 + 2 * Q -
P)) / 2;
1045 ASSERTL2((i >= 0) && (i <= 3),
"face id is out of range");
1046 int Pi =
m_base[0]->GetNumModes() - 2;
1047 int Qi =
m_base[1]->GetNumModes() - 2;
1048 int Ri =
m_base[2]->GetNumModes() - 2;
1052 return Pi * (2 * Qi - Pi - 1) / 2;
1056 return Pi * (2 * Ri - Pi - 1) / 2;
1060 return Qi * (2 * Ri - Qi - 1) / 2;
1066 ASSERTL2(i >= 0 && i <= 3,
"face id is out of range");
1070 return m_base[0]->GetNumPoints() *
m_base[1]->GetNumPoints();
1074 return m_base[0]->GetNumPoints() *
m_base[2]->GetNumPoints();
1078 return m_base[1]->GetNumPoints() *
m_base[2]->GetNumPoints();
1084 ASSERTL2((i >= 0) && (i <= 5),
"edge id is out of range");
1085 int P =
m_base[0]->GetNumModes();
1086 int Q =
m_base[1]->GetNumModes();
1087 int R =
m_base[2]->GetNumModes();
1093 else if (i == 1 || i == 2)
1106 ASSERTL2(i >= 0 && i <= 3,
"face id is out of range");
1107 ASSERTL2(j == 0 || j == 1,
"face direction is out of range");
1111 return m_base[j]->GetPointsKey();
1115 return m_base[2 * j]->GetPointsKey();
1119 return m_base[j + 1]->GetPointsKey();
1124 const std::vector<unsigned int> &nummodes,
int &modes_offset)
1127 nummodes[modes_offset], nummodes[modes_offset + 1],
1128 nummodes[modes_offset + 2]);
1138 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
1139 ASSERTL2(k == 0 || k == 1,
"face direction out of range");
1158 m_base[dir]->GetNumModes(), UseGLL);
1174 for (
int k = 0; k < Qz; ++k)
1176 for (
int j = 0; j < Qy; ++j)
1178 for (
int i = 0; i < Qx; ++i)
1180 int s = i + Qx * (j + Qy * k);
1182 (eta_x[i] + 1.0) * (1.0 - eta_y[j]) * (1.0 - eta_z[k]) / 4 -
1184 xi_y[s] = (eta_y[j] + 1.0) * (1.0 - eta_z[k]) / 2 - 1.0;
1206 "Mapping not defined for this type of basis");
1209 if (useCoeffPacking ==
true)
1211 switch (localVertexId)
1235 ASSERTL0(
false,
"Vertex ID must be between 0 and 3");
1242 switch (localVertexId)
1266 ASSERTL0(
false,
"Vertex ID must be between 0 and 3");
1286 "BasisType is not a boundary interior form");
1289 "BasisType is not a boundary interior form");
1292 "BasisType is not a boundary interior form");
1294 int P =
m_base[0]->GetNumModes();
1295 int Q =
m_base[1]->GetNumModes();
1296 int R =
m_base[2]->GetNumModes();
1300 if (outarray.size() != nIntCoeffs)
1306 for (
int i = 2; i <
P; ++i)
1308 for (
int j = 1; j < Q - i - 1; ++j)
1310 for (
int k = 1; k < R - i - j; ++k)
1312 outarray[idx++] =
GetMode(i, j, k);
1325 "BasisType is not a boundary interior form");
1328 "BasisType is not a boundary interior form");
1331 "BasisType is not a boundary interior form");
1333 int P =
m_base[0]->GetNumModes();
1334 int Q =
m_base[1]->GetNumModes();
1335 int R =
m_base[2]->GetNumModes();
1342 if (outarray.size() != nBnd)
1347 for (i = 0; i <
P; ++i)
1352 for (j = 0; j < Q - i; j++)
1354 for (k = 0; k < R - i - j; ++k)
1356 outarray[idx++] =
GetMode(i, j, k);
1364 for (k = 0; k < R - i; ++k)
1366 outarray[idx++] =
GetMode(i, 0, k);
1368 for (j = 1; j < Q - i; ++j)
1370 outarray[idx++] =
GetMode(i, j, 0);
1380 int P = 0, Q = 0, idx = 0;
1381 int nFaceCoeffs = 0;
1387 Q =
m_base[1]->GetNumModes();
1391 Q =
m_base[2]->GetNumModes();
1396 Q =
m_base[2]->GetNumModes();
1399 ASSERTL0(
false,
"fid must be between 0 and 3");
1402 nFaceCoeffs =
P * (2 * Q -
P + 1) / 2;
1404 if (maparray.size() != nFaceCoeffs)
1413 for (i = 0; i <
P; ++i)
1415 for (j = 0; j < Q - i; ++j)
1417 maparray[idx++] =
GetMode(i, j, 0);
1423 for (i = 0; i <
P; ++i)
1425 for (k = 0; k < Q - i; ++k)
1427 maparray[idx++] =
GetMode(i, 0, k);
1433 for (j = 0; j <
P - 1; ++j)
1435 for (k = 0; k < Q - 1 - j; ++k)
1437 maparray[idx++] =
GetMode(1, j, k);
1439 if (j == 0 && k == 0)
1441 maparray[idx++] =
GetMode(0, 0, 1);
1443 if (j == 0 && k == Q - 2)
1445 for (
int r = 0; r < Q - 1; ++r)
1447 maparray[idx++] =
GetMode(0, 1, r);
1455 for (j = 0; j <
P; ++j)
1457 for (k = 0; k < Q - j; ++k)
1459 maparray[idx++] =
GetMode(0, j, k);
1464 ASSERTL0(
false,
"Element map not available.");
1473 int nummodesA = 0, nummodesB = 0, i, j, k, idx;
1476 "Method only implemented for Modified_A BasisType (x "
1477 "direction), Modified_B BasisType (y direction), and "
1478 "Modified_C BasisType(z direction)");
1480 int nFaceCoeffs = 0;
1485 nummodesA =
m_base[0]->GetNumModes();
1486 nummodesB =
m_base[1]->GetNumModes();
1489 nummodesA =
m_base[0]->GetNumModes();
1490 nummodesB =
m_base[2]->GetNumModes();
1494 nummodesA =
m_base[1]->GetNumModes();
1495 nummodesB =
m_base[2]->GetNumModes();
1498 ASSERTL0(
false,
"fid must be between 0 and 3");
1507 nFaceCoeffs =
P * (2 * Q -
P + 1) / 2;
1510 if (maparray.size() != nFaceCoeffs)
1515 if (signarray.size() != nFaceCoeffs)
1521 fill(signarray.data(), signarray.data() + nFaceCoeffs, 1);
1528 int minPA = min(nummodesA,
P);
1529 int minQB = min(nummodesB, Q);
1531 for (j = 0; j < minPA; ++j)
1534 for (k = 0; k < minQB - j; ++k, ++cnt)
1536 maparray[idx++] = cnt;
1539 cnt += nummodesB - minQB;
1541 for (k = nummodesB - j; k < Q - j; ++k)
1543 signarray[idx] = 0.0;
1544 maparray[idx++] = maparray[0];
1548 for (j = nummodesA; j <
P; ++j)
1550 for (k = 0; k < Q - j; ++k)
1552 signarray[idx] = 0.0;
1553 maparray[idx++] = maparray[0];
1560 for (i = 0; i <
P; ++i)
1562 for (j = 0; j < Q - i; ++j, idx++)
1566 signarray[idx] = (i % 2 ? -1 : 1);
1571 swap(maparray[0], maparray[Q]);
1573 for (i = 1; i < Q - 1; ++i)
1575 swap(maparray[i + 1], maparray[Q + i]);
1588 const int P =
m_base[0]->GetNumModes();
1589 const int Q =
m_base[1]->GetNumModes();
1590 const int R =
m_base[2]->GetNumModes();
1594 if (maparray.size() != nEdgeIntCoeffs)
1600 fill(maparray.data(), maparray.data() + nEdgeIntCoeffs, 0);
1603 if (signarray.size() != nEdgeIntCoeffs)
1609 fill(signarray.data(), signarray.data() + nEdgeIntCoeffs, 1);
1615 for (i = 0; i <
P - 2; ++i)
1617 maparray[i] =
GetMode(i + 2, 0, 0);
1621 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1628 for (i = 0; i < Q - 2; ++i)
1630 maparray[i] =
GetMode(1, i + 1, 0);
1634 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1641 for (i = 0; i < Q - 2; ++i)
1643 maparray[i] =
GetMode(0, i + 2, 0);
1647 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1654 for (i = 0; i < R - 2; ++i)
1656 maparray[i] =
GetMode(0, 0, i + 2);
1660 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1667 for (i = 0; i < R - 2; ++i)
1669 maparray[i] =
GetMode(1, 0, i + 1);
1673 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1680 for (i = 0; i < R - 2; ++i)
1682 maparray[i] =
GetMode(0, 1, i + 1);
1686 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1693 ASSERTL0(
false,
"Edge not defined.");
1703 const int P =
m_base[0]->GetNumModes();
1704 const int Q =
m_base[1]->GetNumModes();
1705 const int R =
m_base[2]->GetNumModes();
1709 if (maparray.size() != nFaceIntCoeffs)
1714 if (signarray.size() != nFaceIntCoeffs)
1720 fill(signarray.data(), signarray.data() + nFaceIntCoeffs, 1);
1727 for (i = 2; i <
P; ++i)
1729 for (j = 1; j < Q - i; ++j)
1731 if ((
int)faceOrient == 7)
1733 signarray[idx] = (i % 2 ? -1 : 1);
1735 maparray[idx++] =
GetMode(i, j, 0);
1741 for (i = 2; i <
P; ++i)
1743 for (k = 1; k < R - i; ++k)
1745 if ((
int)faceOrient == 7)
1747 signarray[idx] = (i % 2 ? -1 : 1);
1749 maparray[idx++] =
GetMode(i, 0, k);
1755 for (j = 1; j < Q - 1; ++j)
1757 for (k = 1; k < R - 1 - j; ++k)
1759 if ((
int)faceOrient == 7)
1761 signarray[idx] = ((j + 1) % 2 ? -1 : 1);
1763 maparray[idx++] =
GetMode(1, j, k);
1769 for (j = 2; j < Q; ++j)
1771 for (k = 1; k < R - j; ++k)
1773 if ((
int)faceOrient == 7)
1775 signarray[idx] = (j % 2 ? -1 : 1);
1777 maparray[idx++] =
GetMode(0, j, k);
1782 ASSERTL0(
false,
"Face interior map not available.");
1800 int nq0 =
m_base[0]->GetNumPoints();
1801 int nq1 =
m_base[1]->GetNumPoints();
1802 int nq2 =
m_base[2]->GetNumPoints();
1812 nq = max(nq0, max(nq1, nq2));
1826 for (
int i = 0; i < nq; ++i)
1828 for (
int j = 0; j < nq - i; ++j)
1830 for (
int k = 0; k < nq - i - j; ++k, ++cnt)
1833 coords[cnt][0] = -1.0 + 2 * k / (
NekDouble)(nq - 1);
1834 coords[cnt][1] = -1.0 + 2 * j / (
NekDouble)(nq - 1);
1835 coords[cnt][2] = -1.0 + 2 * i / (
NekDouble)(nq - 1);
1840 for (
int i = 0; i < neq; ++i)
1844 I[0] =
m_base[0]->GetI(coll);
1845 I[1] =
m_base[1]->GetI(coll + 1);
1846 I[2] =
m_base[2]->GetI(coll + 2);
1850 for (
int k = 0; k < nq2; ++k)
1852 for (
int j = 0; j < nq1; ++j)
1855 fac = (I[1]->GetPtr())[j] * (I[2]->GetPtr())[k];
1859 Mat->GetRawPtr() + k * nq0 * nq1 * neq +
1923 const int Q =
m_base[1]->GetNumModes();
1924 const int R =
m_base[2]->GetNumModes();
1926 int i, j, q_hat, k_hat;
1930 for (i = 0; i < I; ++i)
1936 cnt += q_hat * (q_hat + 1) / 2 + k_hat * (Q - i);
1941 for (j = 0; j < J; ++j)
1959 int nquad0 =
m_base[0]->GetNumPoints();
1960 int nquad1 =
m_base[1]->GetNumPoints();
1961 int nquad2 =
m_base[2]->GetNumPoints();
1971 for (i = 0; i < nquad1 * nquad2; ++i)
1974 1, &outarray[0] + i * nquad0, 1);
1980 case LibUtilities::eGaussRadauMAlpha1Beta0:
1981 for (j = 0; j < nquad2; ++j)
1983 for (i = 0; i < nquad1; ++i)
1986 &outarray[0] + i * nquad0 + j * nquad0 * nquad1,
1993 for (j = 0; j < nquad2; ++j)
1995 for (i = 0; i < nquad1; ++i)
1998 &outarray[0] + i * nquad0 + j * nquad0 * nquad1,
2008 case LibUtilities::eGaussRadauMAlpha2Beta0:
2009 for (i = 0; i < nquad2; ++i)
2012 &outarray[0] + i * nquad0 * nquad1, 1);
2016 case LibUtilities::eGaussRadauMAlpha1Beta0:
2017 for (i = 0; i < nquad2; ++i)
2019 Blas::Dscal(nquad0 * nquad1, 0.25 * (1 - z2[i]) * w2[i],
2020 &outarray[0] + i * nquad0 * nquad1, 1);
2024 for (i = 0; i < nquad2; ++i)
2027 0.25 * (1 - z2[i]) * (1 - z2[i]) * w2[i],
2028 &outarray[0] + i * nquad0 * nquad1, 1);
2045 int qa =
m_base[0]->GetNumPoints();
2046 int qb =
m_base[1]->GetNumPoints();
2047 int qc =
m_base[2]->GetNumPoints();
2048 int nmodes_a =
m_base[0]->GetNumModes();
2049 int nmodes_b =
m_base[1]->GetNumModes();
2050 int nmodes_c =
m_base[2]->GetNumModes();
2064 int i, j, k, cnt = 0;
2067 OrthoExp.
FwdTrans(array, orthocoeffs);
2077 for (i = 0; i < nmodes_a; ++i)
2079 for (j = 0; j < nmodes_b - j; ++j)
2082 pow((1.0 * i) / (nmodes_a - 1), cutoff * nmodes_a),
2083 pow((1.0 * j) / (nmodes_b - 1), cutoff * nmodes_b));
2085 for (k = 0; k < nmodes_c - i - j; ++k)
2088 std::max(fac1, pow((1.0 * k) / (nmodes_c - 1),
2089 cutoff * nmodes_c));
2091 orthocoeffs[cnt] *= SvvDiffCoeff * fac;
2107 max_abc = max(max_abc, 0);
2110 for (i = 0; i < nmodes_a; ++i)
2112 for (j = 0; j < nmodes_b - j; ++j)
2114 int maxij = max(i, j);
2116 for (k = 0; k < nmodes_c - i - j; ++k)
2118 int maxijk = max(maxij, k);
2141 int cutoff_a = (int)(SVVCutOff * nmodes_a);
2142 int cutoff_b = (int)(SVVCutOff * nmodes_b);
2143 int cutoff_c = (int)(SVVCutOff * nmodes_c);
2144 int nmodes = min(min(nmodes_a, nmodes_b), nmodes_c);
2145 NekDouble cutoff = min(min(cutoff_a, cutoff_b), cutoff_c);
2149 for (i = 0; i < nmodes_a; ++i)
2151 for (j = 0; j < nmodes_b - i; ++j)
2153 for (k = 0; k < nmodes_c - i - j; ++k)
2155 if (i + j + k >= cutoff)
2157 orthocoeffs[cnt] *= ((SvvDiffCoeff)*exp(
2158 -(i + j + k - nmodes) * (i + j + k - nmodes) /
2159 ((
NekDouble)((i + j + k - cutoff + epsilon) *
2160 (i + j + k - cutoff + epsilon)))));
2164 orthocoeffs[cnt] *= 0.0;
2173 OrthoExp.
BwdTrans(orthocoeffs, array);
2180 int nquad0 =
m_base[0]->GetNumPoints();
2181 int nquad1 =
m_base[1]->GetNumPoints();
2182 int nquad2 =
m_base[2]->GetNumPoints();
2183 int nqtot = nquad0 * nquad1 * nquad2;
2184 int nmodes0 =
m_base[0]->GetNumModes();
2185 int nmodes1 =
m_base[1]->GetNumModes();
2186 int nmodes2 =
m_base[2]->GetNumModes();
2187 int numMax = nmodes0;
2209 bortho0, bortho1, bortho2);
2212 OrthoTetExp->FwdTrans(phys_tmp, coeff);
2218 for (
int u = 0; u < numMin; ++u)
2220 for (
int i = 0; i < numMin - u; ++i)
2223 tmp2 = coeff_tmp1 + cnt, 1);
2224 cnt += numMax - u - i;
2226 for (
int i = numMin; i < numMax - u; ++i)
2228 cnt += numMax - u - i;
2232 OrthoTetExp->BwdTrans(coeff_tmp1, phys_tmp);
2239 int np0 =
m_base[0]->GetNumPoints();
2240 int np1 =
m_base[1]->GetNumPoints();
2241 int np2 =
m_base[2]->GetNumPoints();
2242 int np = max(np0, max(np1, np2));
2253 for (
int i = 0; i < np - 1; ++i)
2255 planep1 += (np - i) * (np - i + 1) / 2;
2260 for (
int j = 0; j < np - i - 1; ++j)
2262 rowp1 += np - i - j;
2263 row1p1 += np - i - j - 1;
2264 for (
int k = 0; k < np - i - j - 2; ++k)
2266 conn[cnt++] = plane + row + k + 1;
2267 conn[cnt++] = plane + row + k;
2268 conn[cnt++] = plane + rowp1 + k;
2269 conn[cnt++] = planep1 + row1 + k;
2271 conn[cnt++] = plane + row + k + 1;
2272 conn[cnt++] = plane + rowp1 + k + 1;
2273 conn[cnt++] = planep1 + row1 + k + 1;
2274 conn[cnt++] = planep1 + row1 + k;
2276 conn[cnt++] = plane + rowp1 + k + 1;
2277 conn[cnt++] = plane + row + k + 1;
2278 conn[cnt++] = plane + rowp1 + k;
2279 conn[cnt++] = planep1 + row1 + k;
2281 conn[cnt++] = planep1 + row1 + k;
2282 conn[cnt++] = planep1 + row1p1 + k;
2283 conn[cnt++] = plane + rowp1 + k;
2284 conn[cnt++] = plane + rowp1 + k + 1;
2286 conn[cnt++] = planep1 + row1 + k;
2287 conn[cnt++] = planep1 + row1p1 + k;
2288 conn[cnt++] = planep1 + row1 + k + 1;
2289 conn[cnt++] = plane + rowp1 + k + 1;
2291 if (k < np - i - j - 3)
2293 conn[cnt++] = plane + rowp1 + k + 1;
2294 conn[cnt++] = planep1 + row1p1 + k + 1;
2295 conn[cnt++] = planep1 + row1 + k + 1;
2296 conn[cnt++] = planep1 + row1p1 + k;
2300 conn[cnt++] = plane + row + np - i - j - 1;
2301 conn[cnt++] = plane + row + np - i - j - 2;
2302 conn[cnt++] = plane + rowp1 + np - i - j - 2;
2303 conn[cnt++] = planep1 + row1 + np - i - j - 2;
2306 row1 += np - i - j - 1;
2308 plane += (np - i) * (np - i + 1) / 2;
#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.
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space.
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
int v_GetNtraces() const override
int v_GetTraceNcoeffs(const int i) const override
NekDouble v_PhysEvalFirstDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) override
void v_GetTraceNumModes(const int fid, int &numModes0, int &numModes1, Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2) override
int v_GetTraceIntNcoeffs(const int i) const override
void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray) override
void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_GetInteriorMap(Array< OneD, unsigned int > &outarray) override
int v_GetNedges() const override
int GetMode(const int i, const int j, const int k)
Compute the mode number in the expansion for a particular tensorial combination.
int v_GetEdgeNcoeffs(const int i) const override
bool v_IsBoundaryInteriorExpansion() const override
LibUtilities::PointsKey v_GetTracePointsKey(const int i, const int j) const override
NekDouble v_PhysEvaluateBasis(const Array< OneD, const NekDouble > &coords, int mode) final
void v_GetEdgeInteriorToElementMap(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
StdTetExp(const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const LibUtilities::BasisKey &Bc)
int v_GetTraceNumPoints(const int i) const override
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) 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
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey) override
DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey) override
void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta) override
void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray) override
void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true) override
void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
int v_NumBndryCoeffs() const override
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true) override
void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) 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
void v_GetCoords(Array< OneD, NekDouble > &coords_x, Array< OneD, NekDouble > &coords_y, Array< OneD, NekDouble > &coords_z) 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_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_GetElmtTraceToTraceMap(const unsigned int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient=eForwards, int P=-1, int Q=-1) override
int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false) override
int v_GetNverts() const override
void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_dx, Array< OneD, NekDouble > &out_dy, Array< OneD, NekDouble > &out_dz) override
Calculate the derivative of the physical points.
const LibUtilities::BasisKey v_GetTraceBasisKey(const int i, const int k, bool UseGLL=false) const override
LibUtilities::ShapeType v_DetShapeType() const override
void v_LocCollapsedToLocCoord(const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi) override
LibUtilities::ShapeType DetShapeType() const
int v_NumDGBndryCoeffs() 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
int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset) override
DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey) 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 getNumberOfBndCoefficients(int Na, int Nb, int Nc)
int getNumberOfCoefficients(int Na, int Nb, int Nc)
@ 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
@ eFactorSVVDGKerDiffCoeff
@ eFactorSVVPowerKerDiffCoeff
std::shared_ptr< StdTetExp > StdTetExpSharedPtr
LibUtilities::BasisKey EvaluateTriFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes, bool UseGLL)
const int kSVVDGFiltermodesmin
const int kSVVDGFiltermodesmax
const NekDouble kSVVDGFilter[9][11]
@ ePhysInterpToEquiSpaced
@ eDir1BwdDir1_Dir2FwdDir2
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 Zero(int n, T *x, const int incx)
Zero vector.
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
scalarT< T > sqrt(scalarT< T > in)