36 #include <boost/core/ignore_unused.hpp>
46 StdTetExp::StdTetExp()
54 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
57 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
61 "order in 'a' direction is higher than order "
64 "order in 'a' direction is higher than order "
67 "order in 'b' direction is higher than order "
104 int Q0 =
m_base[0]->GetNumPoints();
105 int Q1 =
m_base[1]->GetNumPoints();
106 int Q2 =
m_base[2]->GetNumPoints();
107 int Qtot = Q0 * Q1 * Q2;
114 bool Do_2 = (out_dxi2.size() > 0) ?
true :
false;
115 bool Do_1 = (out_dxi1.size() > 0) ?
true :
false;
132 eta_0 =
m_base[0]->GetZ();
133 eta_1 =
m_base[1]->GetZ();
134 eta_2 =
m_base[2]->GetZ();
140 for (
int k = 0; k < Q2; ++k)
142 for (
int j = 0; j < Q1; ++j, dEta0 += Q0)
144 Vmath::Smul(Q0, 2.0 / (1.0 - eta_1[j]), dEta0, 1, dEta0, 1);
146 fac = 1.0 / (1.0 - eta_2[k]);
147 Vmath::Smul(Q0 * Q1, fac, &out_dEta0[0] + k * Q0 * Q1, 1,
148 &out_dEta0[0] + k * Q0 * Q1, 1);
151 if (out_dxi0.size() > 0)
163 for (
int k = 0; k < Q1 * Q2; ++k)
165 Vmath::Vmul(Q0, &Fac0[0], 1, &out_dEta0[0] + k * Q0, 1,
166 &out_dEta0[0] + k * Q0, 1);
169 for (
int k = 0; k < Q2; ++k)
172 &out_dEta1[0] + k * Q0 * Q1, 1,
173 &out_dEta1[0] + k * Q0 * Q1, 1);
180 Vmath::Vadd(Qtot, out_dEta0, 1, out_dEta1, 1, out_dxi1, 1);
187 for (
int k = 0; k < Q2; ++k)
189 for (
int j = 0; j < Q1; ++j, dEta1 += Q0)
191 Vmath::Smul(Q0, (1.0 + eta_1[j]) / 2.0, dEta1, 1, dEta1, 1);
198 Vmath::Vadd(Qtot, out_dEta0, 1, out_dEta1, 1, out_dxi2, 1);
199 Vmath::Vadd(Qtot, out_dEta2, 1, out_dxi2, 1, out_dxi2, 1);
236 ASSERTL1(
false,
"input dir is out of range");
287 "Basis[1] is not a general tensor type");
291 "Basis[2] is not a general tensor type");
293 if (
m_base[0]->Collocation() &&
m_base[1]->Collocation() &&
298 inarray, 1, outarray, 1);
312 int nquad1 =
m_base[1]->GetNumPoints();
313 int nquad2 =
m_base[2]->GetNumPoints();
314 int order0 =
m_base[0]->GetNumModes();
315 int order1 =
m_base[1]->GetNumModes();
318 nquad2 * nquad1 * order0);
321 m_base[2]->GetBdata(), inarray, outarray, wsp,
true,
344 bool doCheckCollDir0,
bool doCheckCollDir1,
bool doCheckCollDir2)
346 boost::ignore_unused(doCheckCollDir0, doCheckCollDir1, doCheckCollDir2);
348 int nquad0 =
m_base[0]->GetNumPoints();
349 int nquad1 =
m_base[1]->GetNumPoints();
350 int nquad2 =
m_base[2]->GetNumPoints();
352 int order0 =
m_base[0]->GetNumModes();
353 int order1 =
m_base[1]->GetNumModes();
354 int order2 =
m_base[2]->GetNumModes();
358 tmp + nquad2 * order0 * (2 * order1 - order0 + 1) / 2;
360 int i, j, mode, mode1, cnt;
363 mode = mode1 = cnt = 0;
364 for (i = 0; i < order0; ++i)
366 for (j = 0; j < order1 - i; ++j, ++cnt)
369 base2.get() + mode * nquad2, nquad2,
370 inarray.get() + mode1, 1, 0.0, tmp.get() + cnt * nquad2,
372 mode += order2 - i - j;
373 mode1 += order2 - i - j;
376 for (j = order1 - i; j < order2 - i; ++j)
378 mode += order2 - i - j;
388 Blas::Daxpy(nquad2, inarray[1], base2.get() + nquad2, 1,
389 &tmp[0] + nquad2, 1);
392 Blas::Daxpy(nquad2, inarray[1], base2.get() + nquad2, 1,
393 &tmp[0] + order1 * nquad2, 1);
398 for (i = 0; i < order0; ++i)
400 Blas::Dgemm(
'N',
'T', nquad1, nquad2, order1 - i, 1.0,
401 base1.get() + mode * nquad1, nquad1,
402 tmp.get() + mode * nquad2, nquad2, 0.0,
403 tmp1.get() + i * nquad1 * nquad2, nquad1);
414 for (i = 0; i < nquad2; ++i)
416 Blas::Daxpy(nquad1, tmp[nquad2 + i], base1.get() + nquad1, 1,
417 &tmp1[nquad1 * nquad2] + i * nquad1, 1);
422 Blas::Dgemm(
'N',
'T', nquad0, nquad1 * nquad2, order0, 1.0, base0.get(),
423 nquad0, tmp1.get(), nquad1 * nquad2, 0.0, outarray.get(),
445 out = (*matsys) * in;
487 "Basis[1] is not a general tensor type");
491 "Basis[2] is not a general tensor type");
493 if (
m_base[0]->Collocation() &&
m_base[1]->Collocation())
512 inarray.get(), 1, 0.0, outarray.get(), 1);
525 int nquad0 =
m_base[0]->GetNumPoints();
526 int nquad1 =
m_base[1]->GetNumPoints();
527 int nquad2 =
m_base[2]->GetNumPoints();
528 int order0 =
m_base[0]->GetNumModes();
529 int order1 =
m_base[1]->GetNumModes();
532 nquad2 * order0 * (2 * order1 - order0 + 1) / 2);
534 if (multiplybyweights)
541 tmp, outarray, wsp,
true,
true,
true);
547 inarray, outarray, wsp,
true,
true,
true);
557 bool doCheckCollDir0,
bool doCheckCollDir1,
bool doCheckCollDir2)
559 boost::ignore_unused(doCheckCollDir0, doCheckCollDir1, doCheckCollDir2);
561 int nquad0 =
m_base[0]->GetNumPoints();
562 int nquad1 =
m_base[1]->GetNumPoints();
563 int nquad2 =
m_base[2]->GetNumPoints();
565 int order0 =
m_base[0]->GetNumModes();
566 int order1 =
m_base[1]->GetNumModes();
567 int order2 =
m_base[2]->GetNumModes();
572 int i, j, mode, mode1, cnt;
575 Blas::Dgemm(
'T',
'N', nquad1 * nquad2, order0, nquad0, 1.0, inarray.get(),
576 nquad0, base0.get(), nquad0, 0.0, tmp1.get(), nquad1 * nquad2);
579 for (mode = i = 0; i < order0; ++i)
581 Blas::Dgemm(
'T',
'N', nquad2, order1 - i, nquad1, 1.0,
582 tmp1.get() + i * nquad1 * nquad2, nquad1,
583 base1.get() + mode * nquad1, nquad1, 0.0,
584 tmp2.get() + mode * nquad2, nquad2);
593 Blas::Dgemv(
'T', nquad1, nquad2, 1.0, tmp1.get() + nquad1 * nquad2,
594 nquad1, base1.get() + nquad1, 1, 1.0, tmp2.get() + nquad2,
599 mode = mode1 = cnt = 0;
600 for (i = 0; i < order0; ++i)
602 for (j = 0; j < order1 - i; ++j, ++cnt)
605 base2.get() + mode * nquad2, nquad2,
606 tmp2.get() + cnt * nquad2, 1, 0.0,
607 outarray.get() + mode1, 1);
608 mode += order2 - i - j;
609 mode1 += order2 - i - j;
612 for (j = order1 - i; j < order2 - i; ++j)
614 mode += order2 - i - j;
624 Blas::Ddot(nquad2, base2.get() + nquad2, 1, &tmp2[nquad2], 1);
627 outarray[1] +=
Blas::Ddot(nquad2, base2.get() + nquad2, 1,
628 &tmp2[nquad2 * order1], 1);
643 ASSERTL0((dir == 0) || (dir == 1) || (dir == 2),
644 "input dir is out of range");
666 inarray.get(), 1, 0.0, outarray.get(), 1);
680 int nquad0 =
m_base[0]->GetNumPoints();
681 int nquad1 =
m_base[1]->GetNumPoints();
682 int nquad2 =
m_base[2]->GetNumPoints();
683 int nqtot = nquad0 * nquad1 * nquad2;
684 int nmodes0 =
m_base[0]->GetNumModes();
685 int nmodes1 =
m_base[1]->GetNumModes();
686 int wspsize = nquad0 + nquad1 + nquad2 + max(nqtot,
m_ncoeffs) +
687 nquad1 * nquad2 * nmodes0 +
688 nquad2 * nmodes0 * (2 * nmodes1 - nmodes0 + 1) / 2;
701 for (i = 0; i < nquad0; ++i)
703 gfac0[i] = 0.5 * (1 + z0[i]);
707 for (i = 0; i < nquad1; ++i)
709 gfac1[i] = 2.0 / (1 - z1[i]);
713 for (i = 0; i < nquad2; ++i)
715 gfac2[i] = 2.0 / (1 - z2[i]);
719 for (i = 0; i < nquad1 * nquad2; ++i)
721 Vmath::Smul(nquad0, gfac1[i % nquad1], &inarray[0] + i * nquad0, 1,
722 &tmp0[0] + i * nquad0, 1);
724 for (i = 0; i < nquad2; ++i)
726 Vmath::Smul(nquad0 * nquad1, gfac2[i], &tmp0[0] + i * nquad0 * nquad1,
727 1, &tmp0[0] + i * nquad0 * nquad1, 1);
738 m_base[2]->GetBdata(), tmp0, outarray, wsp,
false,
true,
true);
745 for (i = 0; i < nquad1 * nquad2; ++i)
747 Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp0[0] + i * nquad0, 1,
748 &tmp0[0] + i * nquad0, 1);
753 m_base[2]->GetBdata(), tmp0, tmp3, wsp,
false,
true,
true);
755 for (i = 0; i < nquad2; ++i)
758 &inarray[0] + i * nquad0 * nquad1, 1,
759 &tmp0[0] + i * nquad0 * nquad1, 1);
764 m_base[2]->GetBdata(), tmp0, outarray, wsp,
true,
false,
true);
773 for (i = 0; i < nquad1; ++i)
775 gfac1[i] = (1 + z1[i]) / 2;
778 for (i = 0; i < nquad1 * nquad2; ++i)
780 Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp0[0] + i * nquad0, 1,
781 &tmp0[0] + i * nquad0, 1);
785 m_base[2]->GetBdata(), tmp0, tmp3, wsp,
false,
true,
true);
787 for (i = 0; i < nquad2; ++i)
790 &inarray[0] + i * nquad0 * nquad1, 1,
791 &tmp0[0] + i * nquad0 * nquad1, 1);
793 for (i = 0; i < nquad1 * nquad2; ++i)
795 Vmath::Smul(nquad0, gfac1[i % nquad1], &tmp0[0] + i * nquad0, 1,
796 &tmp0[0] + i * nquad0, 1);
801 m_base[2]->GetBdata(), tmp0, tmp4, wsp,
true,
false,
true);
806 m_base[2]->GetDbdata(), tmp0, outarray, wsp,
true,
true,
false);
816 ASSERTL1(
false,
"input dir is out of range");
853 eta[0] = 2.0 * (1.0 + xi[0]) / d12 - 1.0;
854 eta[1] = 2.0 * (1.0 + xi[1]) / d2 - 1.0;
861 xi[1] = (1.0 + eta[0]) * (1.0 - eta[2]) * 0.5 - 1.0;
862 xi[0] = (1.0 + eta[0]) * (-xi[1] - eta[2]) * 0.5 - 1.0;
879 const int nm1 =
m_base[1]->GetNumModes();
880 const int nm2 =
m_base[2]->GetNumModes();
882 const int b = 2 * nm2 + 1;
883 const int mode0 = floor(0.5 * (b -
sqrt(b * b - 8.0 * mode / nm1)));
885 mode - nm1 * (mode0 * (nm2 - 1) + 1 - (mode0 - 2) * (mode0 - 1) / 2);
886 const int mode1 = tmp / (nm2 - mode0);
887 const int mode2 = tmp % (nm2 - mode0);
896 return StdExpansion::BaryEvaluateBasis<2>(coll[2], 1);
898 else if (mode0 == 0 && mode2 == 1)
900 return StdExpansion::BaryEvaluateBasis<1>(coll[1], 0) *
901 StdExpansion::BaryEvaluateBasis<2>(coll[2], 1);
903 else if (mode0 == 1 && mode1 == 1 && mode2 == 0)
905 return StdExpansion::BaryEvaluateBasis<0>(coll[0], 0) *
906 StdExpansion::BaryEvaluateBasis<1>(coll[1], 1);
910 return StdExpansion::BaryEvaluateBasis<0>(coll[0], mode0) *
911 StdExpansion::BaryEvaluateBasis<1>(coll[1], mode1) *
912 StdExpansion::BaryEvaluateBasis<2>(coll[2], mode2);
917 int &numModes0,
int &numModes1,
920 boost::ignore_unused(faceOrient);
922 int nummodes[3] = {
m_base[0]->GetNumModes(),
m_base[1]->GetNumModes(),
923 m_base[2]->GetNumModes()};
928 numModes0 = nummodes[0];
929 numModes1 = nummodes[1];
934 numModes0 = nummodes[0];
935 numModes1 = nummodes[2];
941 numModes0 = nummodes[1];
942 numModes1 = nummodes[2];
976 "BasisType is not a boundary interior form");
979 "BasisType is not a boundary interior form");
982 "BasisType is not a boundary interior form");
984 int P =
m_base[0]->GetNumModes();
985 int Q =
m_base[1]->GetNumModes();
986 int R =
m_base[2]->GetNumModes();
995 "BasisType is not a boundary interior form");
998 "BasisType is not a boundary interior form");
1001 "BasisType is not a boundary interior form");
1003 int P =
m_base[0]->GetNumModes() - 1;
1004 int Q =
m_base[1]->GetNumModes() - 1;
1005 int R =
m_base[2]->GetNumModes() - 1;
1007 return (Q + 1) +
P * (1 + 2 * Q -
P) / 2
1008 + (R + 1) +
P * (1 + 2 * R -
P) / 2
1009 + 2 * (R + 1) + Q * (1 + 2 * R - Q);
1014 ASSERTL2((i >= 0) && (i <= 3),
"face id is out of range");
1015 int nFaceCoeffs = 0;
1016 int nummodesA, nummodesB,
P, Q;
1022 else if ((i == 1) || (i == 2))
1034 nFaceCoeffs = Q + 1 + (
P * (1 + 2 * Q -
P)) / 2;
1040 ASSERTL2((i >= 0) && (i <= 3),
"face id is out of range");
1041 int Pi =
m_base[0]->GetNumModes() - 2;
1042 int Qi =
m_base[1]->GetNumModes() - 2;
1043 int Ri =
m_base[2]->GetNumModes() - 2;
1047 return Pi * (2 * Qi - Pi - 1) / 2;
1051 return Pi * (2 * Ri - Pi - 1) / 2;
1055 return Qi * (2 * Ri - Qi - 1) / 2;
1061 int Pi =
m_base[0]->GetNumModes() - 2;
1062 int Qi =
m_base[1]->GetNumModes() - 2;
1063 int Ri =
m_base[2]->GetNumModes() - 2;
1065 return Pi * (2 * Qi - Pi - 1) / 2 + Pi * (2 * Ri - Pi - 1) / 2 +
1066 Qi * (2 * Ri - Qi - 1);
1071 ASSERTL2(i >= 0 && i <= 3,
"face id is out of range");
1075 return m_base[0]->GetNumPoints() *
m_base[1]->GetNumPoints();
1079 return m_base[0]->GetNumPoints() *
m_base[2]->GetNumPoints();
1083 return m_base[1]->GetNumPoints() *
m_base[2]->GetNumPoints();
1089 ASSERTL2((i >= 0) && (i <= 5),
"edge id is out of range");
1090 int P =
m_base[0]->GetNumModes();
1091 int Q =
m_base[1]->GetNumModes();
1092 int R =
m_base[2]->GetNumModes();
1098 else if (i == 1 || i == 2)
1111 ASSERTL2(i >= 0 && i <= 3,
"face id is out of range");
1112 ASSERTL2(j == 0 || j == 1,
"face direction is out of range");
1116 return m_base[j]->GetPointsKey();
1120 return m_base[2 * j]->GetPointsKey();
1124 return m_base[j + 1]->GetPointsKey();
1129 const std::vector<unsigned int> &nummodes,
int &modes_offset)
1132 nummodes[modes_offset], nummodes[modes_offset + 1],
1133 nummodes[modes_offset + 2]);
1142 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
1143 ASSERTL2(k == 0 || k == 1,
"face direction out of range");
1162 m_base[dir]->GetNumModes());
1181 for (
int k = 0; k < Qz; ++k)
1183 for (
int j = 0; j < Qy; ++j)
1185 for (
int i = 0; i < Qx; ++i)
1187 int s = i + Qx * (j + Qy * k);
1189 (eta_x[i] + 1.0) * (1.0 - eta_y[j]) * (1.0 - eta_z[k]) / 4 -
1191 xi_y[s] = (eta_y[j] + 1.0) * (1.0 - eta_z[k]) / 2 - 1.0;
1213 "Mapping not defined for this type of basis");
1216 if (useCoeffPacking ==
true)
1218 switch (localVertexId)
1242 ASSERTL0(
false,
"Vertex ID must be between 0 and 3");
1249 switch (localVertexId)
1273 ASSERTL0(
false,
"Vertex ID must be between 0 and 3");
1293 "BasisType is not a boundary interior form");
1296 "BasisType is not a boundary interior form");
1299 "BasisType is not a boundary interior form");
1301 int P =
m_base[0]->GetNumModes();
1302 int Q =
m_base[1]->GetNumModes();
1303 int R =
m_base[2]->GetNumModes();
1307 if (outarray.size() != nIntCoeffs)
1313 for (
int i = 2; i <
P - 2; ++i)
1315 for (
int j = 1; j < Q - i - 1; ++j)
1317 for (
int k = 1; k < R - i - j; ++k)
1319 outarray[idx++] =
GetMode(i, j, k);
1332 "BasisType is not a boundary interior form");
1335 "BasisType is not a boundary interior form");
1338 "BasisType is not a boundary interior form");
1340 int P =
m_base[0]->GetNumModes();
1341 int Q =
m_base[1]->GetNumModes();
1342 int R =
m_base[2]->GetNumModes();
1349 if (outarray.size() != nBnd)
1354 for (i = 0; i <
P; ++i)
1359 for (j = 0; j < Q - i; j++)
1361 for (k = 0; k < R - i - j; ++k)
1363 outarray[idx++] =
GetMode(i, j, k);
1371 for (k = 0; k < R - i; ++k)
1373 outarray[idx++] =
GetMode(i, 0, k);
1375 for (j = 1; j < Q - i; ++j)
1377 outarray[idx++] =
GetMode(i, j, 0);
1387 int P = 0, Q = 0, idx = 0;
1388 int nFaceCoeffs = 0;
1394 Q =
m_base[1]->GetNumModes();
1398 Q =
m_base[2]->GetNumModes();
1403 Q =
m_base[2]->GetNumModes();
1406 ASSERTL0(
false,
"fid must be between 0 and 3");
1409 nFaceCoeffs =
P * (2 * Q -
P + 1) / 2;
1411 if (maparray.size() != nFaceCoeffs)
1420 for (i = 0; i <
P; ++i)
1422 for (j = 0; j < Q - i; ++j)
1424 maparray[idx++] =
GetMode(i, j, 0);
1430 for (i = 0; i <
P; ++i)
1432 for (k = 0; k < Q - i; ++k)
1434 maparray[idx++] =
GetMode(i, 0, k);
1440 for (j = 0; j <
P - 1; ++j)
1442 for (k = 0; k < Q - 1 - j; ++k)
1444 maparray[idx++] =
GetMode(1, j, k);
1446 if (j == 0 && k == 0)
1448 maparray[idx++] =
GetMode(0, 0, 1);
1450 if (j == 0 && k == Q - 2)
1452 for (
int r = 0; r < Q - 1; ++r)
1454 maparray[idx++] =
GetMode(0, 1, r);
1462 for (j = 0; j <
P; ++j)
1464 for (k = 0; k < Q - j; ++k)
1466 maparray[idx++] =
GetMode(0, j, k);
1471 ASSERTL0(
false,
"Element map not available.");
1480 int nummodesA = 0, nummodesB = 0, i, j, k, idx;
1483 "Method only implemented for Modified_A BasisType (x "
1484 "direction), Modified_B BasisType (y direction), and "
1485 "Modified_C BasisType(z direction)");
1487 int nFaceCoeffs = 0;
1492 nummodesA =
m_base[0]->GetNumModes();
1493 nummodesB =
m_base[1]->GetNumModes();
1496 nummodesA =
m_base[0]->GetNumModes();
1497 nummodesB =
m_base[2]->GetNumModes();
1501 nummodesA =
m_base[1]->GetNumModes();
1502 nummodesB =
m_base[2]->GetNumModes();
1505 ASSERTL0(
false,
"fid must be between 0 and 3");
1514 nFaceCoeffs =
P * (2 * Q -
P + 1) / 2;
1517 if (maparray.size() != nFaceCoeffs)
1522 if (signarray.size() != nFaceCoeffs)
1528 fill(signarray.get(), signarray.get() + nFaceCoeffs, 1);
1535 int minPA = min(nummodesA,
P);
1536 int minQB = min(nummodesB, Q);
1538 for (j = 0; j < minPA; ++j)
1541 for (k = 0; k < minQB - j; ++k, ++cnt)
1543 maparray[idx++] = cnt;
1546 cnt += nummodesB - minQB;
1548 for (k = nummodesB - j; k < Q - j; ++k)
1550 signarray[idx] = 0.0;
1551 maparray[idx++] = maparray[0];
1555 for (j = nummodesA; j <
P; ++j)
1557 for (k = 0; k < Q - j; ++k)
1559 signarray[idx] = 0.0;
1560 maparray[idx++] = maparray[0];
1567 for (i = 0; i <
P; ++i)
1569 for (j = 0; j < Q - i; ++j, idx++)
1573 signarray[idx] = (i % 2 ? -1 : 1);
1578 swap(maparray[0], maparray[Q]);
1580 for (i = 1; i < Q - 1; ++i)
1582 swap(maparray[i + 1], maparray[Q + i]);
1595 const int P =
m_base[0]->GetNumModes();
1596 const int Q =
m_base[1]->GetNumModes();
1597 const int R =
m_base[2]->GetNumModes();
1601 if (maparray.size() != nEdgeIntCoeffs)
1607 fill(maparray.get(), maparray.get() + nEdgeIntCoeffs, 0);
1610 if (signarray.size() != nEdgeIntCoeffs)
1616 fill(signarray.get(), signarray.get() + nEdgeIntCoeffs, 1);
1622 for (i = 0; i <
P - 2; ++i)
1624 maparray[i] =
GetMode(i + 2, 0, 0);
1628 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1635 for (i = 0; i < Q - 2; ++i)
1637 maparray[i] =
GetMode(1, i + 1, 0);
1641 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1648 for (i = 0; i < Q - 2; ++i)
1650 maparray[i] =
GetMode(0, i + 2, 0);
1654 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1661 for (i = 0; i < R - 2; ++i)
1663 maparray[i] =
GetMode(0, 0, i + 2);
1667 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1674 for (i = 0; i < R - 2; ++i)
1676 maparray[i] =
GetMode(1, 0, i + 1);
1680 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1687 for (i = 0; i < R - 2; ++i)
1689 maparray[i] =
GetMode(0, 1, i + 1);
1693 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1700 ASSERTL0(
false,
"Edge not defined.");
1710 const int P =
m_base[0]->GetNumModes();
1711 const int Q =
m_base[1]->GetNumModes();
1712 const int R =
m_base[2]->GetNumModes();
1716 if (maparray.size() != nFaceIntCoeffs)
1721 if (signarray.size() != nFaceIntCoeffs)
1727 fill(signarray.get(), signarray.get() + nFaceIntCoeffs, 1);
1734 for (i = 2; i <
P - 1; ++i)
1736 for (j = 1; j < Q - i; ++j)
1738 if ((
int)faceOrient == 7)
1740 signarray[idx] = (i % 2 ? -1 : 1);
1742 maparray[idx++] =
GetMode(i, j, 0);
1748 for (i = 2; i <
P; ++i)
1750 for (k = 1; k < R - i; ++k)
1752 if ((
int)faceOrient == 7)
1754 signarray[idx] = (i % 2 ? -1 : 1);
1756 maparray[idx++] =
GetMode(i, 0, k);
1762 for (j = 1; j < Q - 2; ++j)
1764 for (k = 1; k < R - 1 - j; ++k)
1766 if ((
int)faceOrient == 7)
1768 signarray[idx] = ((j + 1) % 2 ? -1 : 1);
1770 maparray[idx++] =
GetMode(1, j, k);
1776 for (j = 2; j < Q - 1; ++j)
1778 for (k = 1; k < R - j; ++k)
1780 if ((
int)faceOrient == 7)
1782 signarray[idx] = (j % 2 ? -1 : 1);
1784 maparray[idx++] =
GetMode(0, j, k);
1789 ASSERTL0(
false,
"Face interior map not available.");
1807 int nq0 =
m_base[0]->GetNumPoints();
1808 int nq1 =
m_base[1]->GetNumPoints();
1809 int nq2 =
m_base[2]->GetNumPoints();
1819 nq = max(nq0, max(nq1, nq2));
1833 for (
int i = 0; i < nq; ++i)
1835 for (
int j = 0; j < nq - i; ++j)
1837 for (
int k = 0; k < nq - i - j; ++k, ++cnt)
1840 coords[cnt][0] = -1.0 + 2 * k / (
NekDouble)(nq - 1);
1841 coords[cnt][1] = -1.0 + 2 * j / (
NekDouble)(nq - 1);
1842 coords[cnt][2] = -1.0 + 2 * i / (
NekDouble)(nq - 1);
1847 for (
int i = 0; i < neq; ++i)
1851 I[0] =
m_base[0]->GetI(coll);
1852 I[1] =
m_base[1]->GetI(coll + 1);
1853 I[2] =
m_base[2]->GetI(coll + 2);
1857 for (
int k = 0; k < nq2; ++k)
1859 for (
int j = 0; j < nq1; ++j)
1862 fac = (I[1]->GetPtr())[j] * (I[2]->GetPtr())[k];
1866 Mat->GetRawPtr() + k * nq0 * nq1 * neq +
1916 const int Q =
m_base[1]->GetNumModes();
1917 const int R =
m_base[2]->GetNumModes();
1919 int i, j, q_hat, k_hat;
1923 for (i = 0; i < I; ++i)
1926 q_hat = min(Q, R - i);
1928 k_hat = max(R - Q - i, 0);
1929 cnt += q_hat * (q_hat + 1) / 2 + k_hat * Q;
1934 for (j = 0; j < J; ++j)
1952 int nquad0 =
m_base[0]->GetNumPoints();
1953 int nquad1 =
m_base[1]->GetNumPoints();
1954 int nquad2 =
m_base[2]->GetNumPoints();
1964 for (i = 0; i < nquad1 * nquad2; ++i)
1967 1, &outarray[0] + i * nquad0, 1);
1973 case LibUtilities::eGaussRadauMAlpha1Beta0:
1974 for (j = 0; j < nquad2; ++j)
1976 for (i = 0; i < nquad1; ++i)
1979 &outarray[0] + i * nquad0 + j * nquad0 * nquad1,
1986 for (j = 0; j < nquad2; ++j)
1988 for (i = 0; i < nquad1; ++i)
1991 &outarray[0] + i * nquad0 + j * nquad0 * nquad1,
2001 case LibUtilities::eGaussRadauMAlpha2Beta0:
2002 for (i = 0; i < nquad2; ++i)
2005 &outarray[0] + i * nquad0 * nquad1, 1);
2009 case LibUtilities::eGaussRadauMAlpha1Beta0:
2010 for (i = 0; i < nquad2; ++i)
2012 Blas::Dscal(nquad0 * nquad1, 0.25 * (1 - z2[i]) * w2[i],
2013 &outarray[0] + i * nquad0 * nquad1, 1);
2017 for (i = 0; i < nquad2; ++i)
2020 0.25 * (1 - z2[i]) * (1 - z2[i]) * w2[i],
2021 &outarray[0] + i * nquad0 * nquad1, 1);
2038 int qa =
m_base[0]->GetNumPoints();
2039 int qb =
m_base[1]->GetNumPoints();
2040 int qc =
m_base[2]->GetNumPoints();
2041 int nmodes_a =
m_base[0]->GetNumModes();
2042 int nmodes_b =
m_base[1]->GetNumModes();
2043 int nmodes_c =
m_base[2]->GetNumModes();
2057 int i, j, k, cnt = 0;
2060 OrthoExp.
FwdTrans(array, orthocoeffs);
2070 for (
int i = 0; i < nmodes_a; ++i)
2072 for (
int j = 0; j < nmodes_b - j; ++j)
2075 pow((1.0 * i) / (nmodes_a - 1), cutoff * nmodes_a),
2076 pow((1.0 * j) / (nmodes_b - 1), cutoff * nmodes_b));
2078 for (
int k = 0; k < nmodes_c - i - j; ++k)
2081 std::max(fac1, pow((1.0 * k) / (nmodes_c - 1),
2082 cutoff * nmodes_c));
2084 orthocoeffs[cnt] *= SvvDiffCoeff * fac;
2100 max_abc = max(max_abc, 0);
2103 for (
int i = 0; i < nmodes_a; ++i)
2105 for (
int j = 0; j < nmodes_b - j; ++j)
2107 int maxij = max(i, j);
2109 for (
int k = 0; k < nmodes_c - i - j; ++k)
2111 int maxijk = max(maxij, k);
2134 int cutoff_a = (int)(SVVCutOff * nmodes_a);
2135 int cutoff_b = (int)(SVVCutOff * nmodes_b);
2136 int cutoff_c = (int)(SVVCutOff * nmodes_c);
2137 int nmodes = min(min(nmodes_a, nmodes_b), nmodes_c);
2138 NekDouble cutoff = min(min(cutoff_a, cutoff_b), cutoff_c);
2142 for (i = 0; i < nmodes_a; ++i)
2144 for (j = 0; j < nmodes_b - i; ++j)
2146 for (k = 0; k < nmodes_c - i - j; ++k)
2148 if (i + j + k >= cutoff)
2150 orthocoeffs[cnt] *= ((SvvDiffCoeff)*exp(
2151 -(i + j + k - nmodes) * (i + j + k - nmodes) /
2152 ((
NekDouble)((i + j + k - cutoff + epsilon) *
2153 (i + j + k - cutoff + epsilon)))));
2157 orthocoeffs[cnt] *= 0.0;
2166 OrthoExp.
BwdTrans(orthocoeffs, array);
2173 int nquad0 =
m_base[0]->GetNumPoints();
2174 int nquad1 =
m_base[1]->GetNumPoints();
2175 int nquad2 =
m_base[2]->GetNumPoints();
2176 int nqtot = nquad0 * nquad1 * nquad2;
2177 int nmodes0 =
m_base[0]->GetNumModes();
2178 int nmodes1 =
m_base[1]->GetNumModes();
2179 int nmodes2 =
m_base[2]->GetNumModes();
2180 int numMax = nmodes0;
2202 bortho0, bortho1, bortho2);
2205 OrthoTetExp->FwdTrans(phys_tmp, coeff);
2211 for (
int u = 0; u < numMin; ++u)
2213 for (
int i = 0; i < numMin - u; ++i)
2216 tmp2 = coeff_tmp1 + cnt, 1);
2217 cnt += numMax - u - i;
2219 for (
int i = numMin; i < numMax - u; ++i)
2221 cnt += numMax - u - i;
2225 OrthoTetExp->BwdTrans(coeff_tmp1, phys_tmp);
2232 boost::ignore_unused(standard);
2234 int np0 =
m_base[0]->GetNumPoints();
2235 int np1 =
m_base[1]->GetNumPoints();
2236 int np2 =
m_base[2]->GetNumPoints();
2237 int np = max(np0, max(np1, np2));
2248 for (
int i = 0; i < np - 1; ++i)
2250 planep1 += (np - i) * (np - i + 1) / 2;
2255 for (
int j = 0; j < np - i - 1; ++j)
2257 rowp1 += np - i - j;
2258 row1p1 += np - i - j - 1;
2259 for (
int k = 0; k < np - i - j - 2; ++k)
2261 conn[cnt++] = plane + row + k + 1;
2262 conn[cnt++] = plane + row + k;
2263 conn[cnt++] = plane + rowp1 + k;
2264 conn[cnt++] = planep1 + row1 + k;
2266 conn[cnt++] = plane + row + k + 1;
2267 conn[cnt++] = plane + rowp1 + k + 1;
2268 conn[cnt++] = planep1 + row1 + k + 1;
2269 conn[cnt++] = planep1 + row1 + k;
2271 conn[cnt++] = plane + rowp1 + k + 1;
2272 conn[cnt++] = plane + row + k + 1;
2273 conn[cnt++] = plane + rowp1 + k;
2274 conn[cnt++] = planep1 + row1 + k;
2276 conn[cnt++] = planep1 + row1 + k;
2277 conn[cnt++] = planep1 + row1p1 + k;
2278 conn[cnt++] = plane + rowp1 + k;
2279 conn[cnt++] = plane + rowp1 + k + 1;
2281 conn[cnt++] = planep1 + row1 + k;
2282 conn[cnt++] = planep1 + row1p1 + k;
2283 conn[cnt++] = planep1 + row1 + k + 1;
2284 conn[cnt++] = plane + rowp1 + k + 1;
2286 if (k < np - i - j - 3)
2288 conn[cnt++] = plane + rowp1 + k + 1;
2289 conn[cnt++] = planep1 + row1p1 + k + 1;
2290 conn[cnt++] = planep1 + row1 + k + 1;
2291 conn[cnt++] = planep1 + row1p1 + k;
2295 conn[cnt++] = plane + row + np - i - j - 1;
2296 conn[cnt++] = plane + row + np - i - j - 2;
2297 conn[cnt++] = plane + rowp1 + np - i - j - 2;
2298 conn[cnt++] = planep1 + row1 + np - i - j - 2;
2301 row1 += np - i - j - 1;
2303 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)
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.
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
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.
Array< OneD, LibUtilities::BasisSharedPtr > m_base
MatrixType GetMatrixType() const
NekDouble GetConstFactor(const ConstFactorType &factor) const
bool ConstFactorExists(const ConstFactorType &factor) const
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual int v_NumBndryCoeffs() const
virtual int v_GetNedges() const
virtual void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true)
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_IProductWRTBase_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_StdPhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
virtual 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)
int GetMode(const int i, const int j, const int k)
Compute the mode number in the expansion for a particular tensorial combination.
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
virtual int v_GetNtraces() const
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
NekDouble v_PhysEvaluateBasis(const Array< OneD, const NekDouble > &coords, int mode) final
virtual int v_GetTraceNumPoints(const int i) const
virtual void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual int v_GetEdgeNcoeffs(const int i) const
virtual int v_GetTotalTraceIntNcoeffs() const
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_GetTraceNumModes(const int fid, int &numModes0, int &numModes1, Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2)
virtual int v_GetTraceNcoeffs(const int i) const
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_x, Array< OneD, NekDouble > &coords_y, Array< OneD, NekDouble > &coords_z)
virtual bool v_IsBoundaryInteriorExpansion()
virtual void v_GetTraceCoeffMap(const unsigned int fid, Array< OneD, unsigned int > &maparray)
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
virtual void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
virtual int v_GetTraceIntNcoeffs(const int i) const
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
virtual int v_GetNverts() const
virtual LibUtilities::PointsKey v_GetTracePointsKey(const int i, const int j) const
virtual void v_GetEdgeInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2)
virtual 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)
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_dx, Array< OneD, NekDouble > &out_dy, Array< OneD, NekDouble > &out_dz)
Calculate the derivative of the physical points.
virtual void v_GetTraceInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2)
virtual 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)
virtual void v_LocCollapsedToLocCoord(const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi)
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
LibUtilities::ShapeType DetShapeType() const
virtual const LibUtilities::BasisKey v_GetTraceBasisKey(const int i, const int k) const
virtual int v_NumDGBndryCoeffs() const
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
virtual LibUtilities::ShapeType v_DetShapeType() const
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 = A x 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 .
@ 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
The above copyright notice and this permission notice shall be included.
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)