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.get() + mode * nquad2, nquad2,
354 inarray.get() + mode1, 1, 0.0, tmp.get() + cnt * nquad2,
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.get() + nquad2, 1,
373 &tmp[0] + nquad2, 1);
376 Blas::Daxpy(nquad2, inarray[1], base2.get() + 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.get() + mode * nquad1, nquad1,
386 tmp.get() + mode * nquad2, nquad2, 0.0,
387 tmp1.get() + i * nquad1 * nquad2, nquad1);
398 for (i = 0; i < nquad2; ++i)
400 Blas::Daxpy(nquad1, tmp[nquad2 + i], base1.get() + nquad1, 1,
401 &tmp1[nquad1 * nquad2] + i * nquad1, 1);
406 Blas::Dgemm(
'N',
'T', nquad0, nquad1 * nquad2, order0, 1.0, base0.get(),
407 nquad0, tmp1.get(), nquad1 * nquad2, 0.0, outarray.get(),
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.get(),
548 nquad0, base0.get(), nquad0, 0.0, tmp1.get(), nquad1 * nquad2);
551 for (mode = i = 0; i < order0; ++i)
553 Blas::Dgemm(
'T',
'N', nquad2, order1 - i, nquad1, 1.0,
554 tmp1.get() + i * nquad1 * nquad2, nquad1,
555 base1.get() + mode * nquad1, nquad1, 0.0,
556 tmp2.get() + mode * nquad2, nquad2);
565 Blas::Dgemv(
'T', nquad1, nquad2, 1.0, tmp1.get() + nquad1 * nquad2,
566 nquad1, base1.get() + nquad1, 1, 1.0, tmp2.get() + nquad2,
571 mode = mode1 = cnt = 0;
572 for (i = 0; i < order0; ++i)
574 for (j = 0; j < order1 - i; ++j, ++cnt)
577 base2.get() + mode * nquad2, nquad2,
578 tmp2.get() + cnt * nquad2, 1, 0.0,
579 outarray.get() + mode1, 1);
580 mode += order2 - i - j;
581 mode1 += order2 - i - j;
584 for (j = order1 - i; j < order2 - i; ++j)
586 mode += order2 - i - j;
596 Blas::Ddot(nquad2, base2.get() + nquad2, 1, &tmp2[nquad2], 1);
599 outarray[1] +=
Blas::Ddot(nquad2, base2.get() + nquad2, 1,
600 &tmp2[nquad2 * order1], 1);
622 int nquad0 =
m_base[0]->GetNumPoints();
623 int nquad1 =
m_base[1]->GetNumPoints();
624 int nquad2 =
m_base[2]->GetNumPoints();
625 int nqtot = nquad0 * nquad1 * nquad2;
626 int nmodes0 =
m_base[0]->GetNumModes();
627 int nmodes1 =
m_base[1]->GetNumModes();
628 int wspsize = nquad0 + nquad1 + nquad2 + max(nqtot,
m_ncoeffs) +
629 nquad1 * nquad2 * nmodes0 +
630 nquad2 * nmodes0 * (2 * nmodes1 - nmodes0 + 1) / 2;
643 for (i = 0; i < nquad0; ++i)
645 gfac0[i] = 0.5 * (1 + z0[i]);
649 for (i = 0; i < nquad1; ++i)
651 gfac1[i] = 2.0 / (1 - z1[i]);
655 for (i = 0; i < nquad2; ++i)
657 gfac2[i] = 2.0 / (1 - z2[i]);
661 for (i = 0; i < nquad1 * nquad2; ++i)
663 Vmath::Smul(nquad0, gfac1[i % nquad1], &inarray[0] + i * nquad0, 1,
664 &tmp0[0] + i * nquad0, 1);
666 for (i = 0; i < nquad2; ++i)
668 Vmath::Smul(nquad0 * nquad1, gfac2[i], &tmp0[0] + i * nquad0 * nquad1,
669 1, &tmp0[0] + i * nquad0 * nquad1, 1);
680 m_base[2]->GetBdata(), tmp0, outarray, wsp,
false,
true,
true);
687 for (i = 0; i < nquad1 * nquad2; ++i)
689 Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp0[0] + i * nquad0, 1,
690 &tmp0[0] + i * nquad0, 1);
695 m_base[2]->GetBdata(), tmp0, tmp3, wsp,
false,
true,
true);
697 for (i = 0; i < nquad2; ++i)
700 &inarray[0] + i * nquad0 * nquad1, 1,
701 &tmp0[0] + i * nquad0 * nquad1, 1);
706 m_base[2]->GetBdata(), tmp0, outarray, wsp,
true,
false,
true);
715 for (i = 0; i < nquad1; ++i)
717 gfac1[i] = (1 + z1[i]) / 2;
720 for (i = 0; i < nquad1 * nquad2; ++i)
722 Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp0[0] + i * nquad0, 1,
723 &tmp0[0] + i * nquad0, 1);
727 m_base[2]->GetBdata(), tmp0, tmp3, wsp,
false,
true,
true);
729 for (i = 0; i < nquad2; ++i)
732 &inarray[0] + i * nquad0 * nquad1, 1,
733 &tmp0[0] + i * nquad0 * nquad1, 1);
735 for (i = 0; i < nquad1 * nquad2; ++i)
737 Vmath::Smul(nquad0, gfac1[i % nquad1], &tmp0[0] + i * nquad0, 1,
738 &tmp0[0] + i * nquad0, 1);
743 m_base[2]->GetBdata(), tmp0, tmp4, wsp,
true,
false,
true);
748 m_base[2]->GetDbdata(), tmp0, outarray, wsp,
true,
true,
false);
758 ASSERTL1(
false,
"input dir is out of range");
795 eta[0] = 2.0 * (1.0 + xi[0]) / d12 - 1.0;
796 eta[1] = 2.0 * (1.0 + xi[1]) / d2 - 1.0;
803 xi[1] = (1.0 + eta[0]) * (1.0 - eta[2]) * 0.5 - 1.0;
804 xi[0] = (1.0 + eta[0]) * (-xi[1] - eta[2]) * 0.5 - 1.0;
821 const int nm1 =
m_base[1]->GetNumModes();
822 const int nm2 =
m_base[2]->GetNumModes();
824 const int b = 2 * nm2 + 1;
825 const int mode0 = floor(0.5 * (b -
sqrt(b * b - 8.0 * mode / nm1)));
827 mode - nm1 * (mode0 * (nm2 - 1) + 1 - (mode0 - 2) * (mode0 - 1) / 2);
828 const int mode1 = tmp / (nm2 - mode0);
829 const int mode2 = tmp % (nm2 - mode0);
838 return StdExpansion::BaryEvaluateBasis<2>(coll[2], 1);
840 else if (mode0 == 0 && mode2 == 1)
842 return StdExpansion::BaryEvaluateBasis<1>(coll[1], 0) *
843 StdExpansion::BaryEvaluateBasis<2>(coll[2], 1);
845 else if (mode0 == 1 && mode1 == 1 && mode2 == 0)
847 return StdExpansion::BaryEvaluateBasis<0>(coll[0], 0) *
848 StdExpansion::BaryEvaluateBasis<1>(coll[1], 1);
852 return StdExpansion::BaryEvaluateBasis<0>(coll[0], mode0) *
853 StdExpansion::BaryEvaluateBasis<1>(coll[1], mode1) *
854 StdExpansion::BaryEvaluateBasis<2>(coll[2], mode2);
859 std::array<NekDouble, 3> &firstOrderDerivs)
866 if ((1 - coll[1]) < 1e-5 || (1 - coll[2]) < 1e-5)
870 EphysDeriv2(totPoints);
871 PhysDeriv(inarray, EphysDeriv0, EphysDeriv1, EphysDeriv2);
874 I[0] =
GetBase()[0]->GetI(coll);
875 I[1] =
GetBase()[1]->GetI(coll + 1);
876 I[2] =
GetBase()[2]->GetI(coll + 2);
884 std::array<NekDouble, 3> interDeriv;
888 NekDouble temp = 2.0 / ((1 - coll[1]) * (1 - coll[2]));
889 interDeriv[0] *= temp;
892 firstOrderDerivs[0] = 2 * interDeriv[0];
899 interDeriv[0] *= fac0;
902 fac0 = 2 / (1 - coll[2]);
903 interDeriv[1] *= fac0;
907 firstOrderDerivs[1] = interDeriv[0] + interDeriv[1];
910 fac0 = (1 + coll[1]) / 2;
911 interDeriv[1] *= fac0;
916 firstOrderDerivs[2] = interDeriv[0] + interDeriv[1] + interDeriv[2];
925 int nummodes[3] = {
m_base[0]->GetNumModes(),
m_base[1]->GetNumModes(),
926 m_base[2]->GetNumModes()};
931 numModes0 = nummodes[0];
932 numModes1 = nummodes[1];
937 numModes0 = nummodes[0];
938 numModes1 = nummodes[2];
944 numModes0 = nummodes[1];
945 numModes1 = nummodes[2];
979 "BasisType is not a boundary interior form");
982 "BasisType is not a boundary interior form");
985 "BasisType is not a boundary interior form");
987 int P =
m_base[0]->GetNumModes();
988 int Q =
m_base[1]->GetNumModes();
989 int R =
m_base[2]->GetNumModes();
998 "BasisType is not a boundary interior form");
1001 "BasisType is not a boundary interior form");
1004 "BasisType is not a boundary interior form");
1006 int P =
m_base[0]->GetNumModes() - 1;
1007 int Q =
m_base[1]->GetNumModes() - 1;
1008 int R =
m_base[2]->GetNumModes() - 1;
1010 return (Q + 1) +
P * (1 + 2 * Q -
P) / 2
1011 + (R + 1) +
P * (1 + 2 * R -
P) / 2
1012 + 2 * (R + 1) + Q * (1 + 2 * R - Q);
1017 ASSERTL2((i >= 0) && (i <= 3),
"face id is out of range");
1018 int nFaceCoeffs = 0;
1019 int nummodesA, nummodesB,
P, Q;
1025 else if ((i == 1) || (i == 2))
1037 nFaceCoeffs = Q + 1 + (
P * (1 + 2 * Q -
P)) / 2;
1043 ASSERTL2((i >= 0) && (i <= 3),
"face id is out of range");
1044 int Pi =
m_base[0]->GetNumModes() - 2;
1045 int Qi =
m_base[1]->GetNumModes() - 2;
1046 int Ri =
m_base[2]->GetNumModes() - 2;
1050 return Pi * (2 * Qi - Pi - 1) / 2;
1054 return Pi * (2 * Ri - Pi - 1) / 2;
1058 return Qi * (2 * Ri - Qi - 1) / 2;
1064 ASSERTL2(i >= 0 && i <= 3,
"face id is out of range");
1068 return m_base[0]->GetNumPoints() *
m_base[1]->GetNumPoints();
1072 return m_base[0]->GetNumPoints() *
m_base[2]->GetNumPoints();
1076 return m_base[1]->GetNumPoints() *
m_base[2]->GetNumPoints();
1082 ASSERTL2((i >= 0) && (i <= 5),
"edge id is out of range");
1083 int P =
m_base[0]->GetNumModes();
1084 int Q =
m_base[1]->GetNumModes();
1085 int R =
m_base[2]->GetNumModes();
1091 else if (i == 1 || i == 2)
1104 ASSERTL2(i >= 0 && i <= 3,
"face id is out of range");
1105 ASSERTL2(j == 0 || j == 1,
"face direction is out of range");
1109 return m_base[j]->GetPointsKey();
1113 return m_base[2 * j]->GetPointsKey();
1117 return m_base[j + 1]->GetPointsKey();
1122 const std::vector<unsigned int> &nummodes,
int &modes_offset)
1125 nummodes[modes_offset], nummodes[modes_offset + 1],
1126 nummodes[modes_offset + 2]);
1135 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
1136 ASSERTL2(k == 0 || k == 1,
"face direction out of range");
1155 m_base[dir]->GetNumModes());
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 - 2; ++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.get(), signarray.get() + 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.get(), maparray.get() + nEdgeIntCoeffs, 0);
1603 if (signarray.size() != nEdgeIntCoeffs)
1609 fill(signarray.get(), signarray.get() + 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.get(), signarray.get() + nFaceIntCoeffs, 1);
1727 for (i = 2; i <
P - 1; ++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 - 2; ++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 - 1; ++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 +
1909 const int Q =
m_base[1]->GetNumModes();
1910 const int R =
m_base[2]->GetNumModes();
1912 int i, j, q_hat, k_hat;
1916 for (i = 0; i < I; ++i)
1919 q_hat = min(Q, R - i);
1921 k_hat = max(R - Q - i, 0);
1922 cnt += q_hat * (q_hat + 1) / 2 + k_hat * Q;
1927 for (j = 0; j < J; ++j)
1945 int nquad0 =
m_base[0]->GetNumPoints();
1946 int nquad1 =
m_base[1]->GetNumPoints();
1947 int nquad2 =
m_base[2]->GetNumPoints();
1957 for (i = 0; i < nquad1 * nquad2; ++i)
1960 1, &outarray[0] + i * nquad0, 1);
1966 case LibUtilities::eGaussRadauMAlpha1Beta0:
1967 for (j = 0; j < nquad2; ++j)
1969 for (i = 0; i < nquad1; ++i)
1972 &outarray[0] + i * nquad0 + j * nquad0 * nquad1,
1979 for (j = 0; j < nquad2; ++j)
1981 for (i = 0; i < nquad1; ++i)
1984 &outarray[0] + i * nquad0 + j * nquad0 * nquad1,
1994 case LibUtilities::eGaussRadauMAlpha2Beta0:
1995 for (i = 0; i < nquad2; ++i)
1998 &outarray[0] + i * nquad0 * nquad1, 1);
2002 case LibUtilities::eGaussRadauMAlpha1Beta0:
2003 for (i = 0; i < nquad2; ++i)
2005 Blas::Dscal(nquad0 * nquad1, 0.25 * (1 - z2[i]) * w2[i],
2006 &outarray[0] + i * nquad0 * nquad1, 1);
2010 for (i = 0; i < nquad2; ++i)
2013 0.25 * (1 - z2[i]) * (1 - z2[i]) * w2[i],
2014 &outarray[0] + i * nquad0 * nquad1, 1);
2031 int qa =
m_base[0]->GetNumPoints();
2032 int qb =
m_base[1]->GetNumPoints();
2033 int qc =
m_base[2]->GetNumPoints();
2034 int nmodes_a =
m_base[0]->GetNumModes();
2035 int nmodes_b =
m_base[1]->GetNumModes();
2036 int nmodes_c =
m_base[2]->GetNumModes();
2050 int i, j, k, cnt = 0;
2053 OrthoExp.
FwdTrans(array, orthocoeffs);
2063 for (
int i = 0; i < nmodes_a; ++i)
2065 for (
int j = 0; j < nmodes_b - j; ++j)
2068 pow((1.0 * i) / (nmodes_a - 1), cutoff * nmodes_a),
2069 pow((1.0 * j) / (nmodes_b - 1), cutoff * nmodes_b));
2071 for (
int k = 0; k < nmodes_c - i - j; ++k)
2074 std::max(fac1, pow((1.0 * k) / (nmodes_c - 1),
2075 cutoff * nmodes_c));
2077 orthocoeffs[cnt] *= SvvDiffCoeff * fac;
2093 max_abc = max(max_abc, 0);
2096 for (
int i = 0; i < nmodes_a; ++i)
2098 for (
int j = 0; j < nmodes_b - j; ++j)
2100 int maxij = max(i, j);
2102 for (
int k = 0; k < nmodes_c - i - j; ++k)
2104 int maxijk = max(maxij, k);
2127 int cutoff_a = (int)(SVVCutOff * nmodes_a);
2128 int cutoff_b = (int)(SVVCutOff * nmodes_b);
2129 int cutoff_c = (int)(SVVCutOff * nmodes_c);
2130 int nmodes = min(min(nmodes_a, nmodes_b), nmodes_c);
2131 NekDouble cutoff = min(min(cutoff_a, cutoff_b), cutoff_c);
2135 for (i = 0; i < nmodes_a; ++i)
2137 for (j = 0; j < nmodes_b - i; ++j)
2139 for (k = 0; k < nmodes_c - i - j; ++k)
2141 if (i + j + k >= cutoff)
2143 orthocoeffs[cnt] *= ((SvvDiffCoeff)*exp(
2144 -(i + j + k - nmodes) * (i + j + k - nmodes) /
2145 ((
NekDouble)((i + j + k - cutoff + epsilon) *
2146 (i + j + k - cutoff + epsilon)))));
2150 orthocoeffs[cnt] *= 0.0;
2159 OrthoExp.
BwdTrans(orthocoeffs, array);
2166 int nquad0 =
m_base[0]->GetNumPoints();
2167 int nquad1 =
m_base[1]->GetNumPoints();
2168 int nquad2 =
m_base[2]->GetNumPoints();
2169 int nqtot = nquad0 * nquad1 * nquad2;
2170 int nmodes0 =
m_base[0]->GetNumModes();
2171 int nmodes1 =
m_base[1]->GetNumModes();
2172 int nmodes2 =
m_base[2]->GetNumModes();
2173 int numMax = nmodes0;
2195 bortho0, bortho1, bortho2);
2198 OrthoTetExp->FwdTrans(phys_tmp, coeff);
2204 for (
int u = 0; u < numMin; ++u)
2206 for (
int i = 0; i < numMin - u; ++i)
2209 tmp2 = coeff_tmp1 + cnt, 1);
2210 cnt += numMax - u - i;
2212 for (
int i = numMin; i < numMax - u; ++i)
2214 cnt += numMax - u - i;
2218 OrthoTetExp->BwdTrans(coeff_tmp1, phys_tmp);
2225 int np0 =
m_base[0]->GetNumPoints();
2226 int np1 =
m_base[1]->GetNumPoints();
2227 int np2 =
m_base[2]->GetNumPoints();
2228 int np = max(np0, max(np1, np2));
2239 for (
int i = 0; i < np - 1; ++i)
2241 planep1 += (np - i) * (np - i + 1) / 2;
2246 for (
int j = 0; j < np - i - 1; ++j)
2248 rowp1 += np - i - j;
2249 row1p1 += np - i - j - 1;
2250 for (
int k = 0; k < np - i - j - 2; ++k)
2252 conn[cnt++] = plane + row + k + 1;
2253 conn[cnt++] = plane + row + k;
2254 conn[cnt++] = plane + rowp1 + k;
2255 conn[cnt++] = planep1 + row1 + k;
2257 conn[cnt++] = plane + row + k + 1;
2258 conn[cnt++] = plane + rowp1 + k + 1;
2259 conn[cnt++] = planep1 + row1 + k + 1;
2260 conn[cnt++] = planep1 + row1 + k;
2262 conn[cnt++] = plane + rowp1 + k + 1;
2263 conn[cnt++] = plane + row + k + 1;
2264 conn[cnt++] = plane + rowp1 + k;
2265 conn[cnt++] = planep1 + row1 + k;
2267 conn[cnt++] = planep1 + row1 + k;
2268 conn[cnt++] = planep1 + row1p1 + k;
2269 conn[cnt++] = plane + rowp1 + k;
2270 conn[cnt++] = plane + rowp1 + k + 1;
2272 conn[cnt++] = planep1 + row1 + k;
2273 conn[cnt++] = planep1 + row1p1 + k;
2274 conn[cnt++] = planep1 + row1 + k + 1;
2275 conn[cnt++] = plane + rowp1 + k + 1;
2277 if (k < np - i - j - 3)
2279 conn[cnt++] = plane + rowp1 + k + 1;
2280 conn[cnt++] = planep1 + row1p1 + k + 1;
2281 conn[cnt++] = planep1 + row1 + k + 1;
2282 conn[cnt++] = planep1 + row1p1 + k;
2286 conn[cnt++] = plane + row + np - i - j - 1;
2287 conn[cnt++] = plane + row + np - i - j - 2;
2288 conn[cnt++] = plane + rowp1 + np - i - j - 2;
2289 conn[cnt++] = planep1 + row1 + np - i - j - 2;
2292 row1 += np - i - j - 1;
2294 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
void v_GetTraceNumModes(const int fid, int &numModes0, int &numModes1, Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2) override
NekDouble v_PhysEvaluate(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) 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
const LibUtilities::BasisKey v_GetTraceBasisKey(const int i, const int k) 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.
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)
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
LibUtilities::BasisKey EvaluateTriFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
@ eFactorSVVDGKerDiffCoeff
@ eFactorSVVPowerKerDiffCoeff
std::shared_ptr< StdTetExp > StdTetExpSharedPtr
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)