48 :
StdExpansion(Ba.GetNumModes() * Bb.GetNumModes() * Bc.GetNumModes(), 3,
50 StdExpansion3D(Ba.GetNumModes() * Bb.GetNumModes() * Bc.GetNumModes(), Ba,
112 ASSERTL1(
false,
"input dir is out of range");
158 "Basis[1] is not a general tensor type");
162 "Basis[2] is not a general tensor type");
164 if (
m_base[0]->Collocation() &&
m_base[1]->Collocation() &&
169 inarray, 1, outarray, 1);
188 m_base[2]->GetBdata(), inarray, outarray, wsp,
true,
211 bool doCheckCollDir0,
bool doCheckCollDir1,
bool doCheckCollDir2)
213 int nquad0 =
m_base[0]->GetNumPoints();
214 int nquad1 =
m_base[1]->GetNumPoints();
215 int nquad2 =
m_base[2]->GetNumPoints();
216 int nmodes0 =
m_base[0]->GetNumModes();
217 int nmodes1 =
m_base[1]->GetNumModes();
218 int nmodes2 =
m_base[2]->GetNumModes();
221 bool colldir0 = doCheckCollDir0 ? (
m_base[0]->Collocation()) :
false;
222 bool colldir1 = doCheckCollDir1 ? (
m_base[1]->Collocation()) :
false;
223 bool colldir2 = doCheckCollDir2 ? (
m_base[2]->Collocation()) :
false;
227 if (colldir0 && colldir1 && colldir2)
234 ASSERTL1(wsp.size() >= nquad0 * nmodes2 * (nmodes1 + nquad1),
235 "Workspace size is not sufficient");
241 Blas::Dgemm(
'T',
'T', nmodes1 * nmodes2, nquad0, nmodes0, 1.0,
242 &inarray[0], nmodes0, base0.data(), nquad0, 0.0, &wsp[0],
244 Blas::Dgemm(
'T',
'T', nquad0 * nmodes2, nquad1, nmodes1, 1.0, &wsp[0],
245 nmodes1, base1.data(), nquad1, 0.0, &wsp2[0],
247 Blas::Dgemm(
'T',
'T', nquad0 * nquad1, nquad2, nmodes2, 1.0, &wsp2[0],
248 nmodes2, base2.data(), nquad2, 0.0, &outarray[0],
267 if ((
m_base[0]->Collocation()) && (
m_base[1]->Collocation()) &&
268 (
m_base[2]->Collocation()))
286 out = (*matsys) * in;
323 if (
m_base[0]->Collocation() &&
m_base[1]->Collocation() &&
341 int nquad0 =
m_base[0]->GetNumPoints();
342 int nquad1 =
m_base[1]->GetNumPoints();
343 int nquad2 =
m_base[2]->GetNumPoints();
344 int order0 =
m_base[0]->GetNumModes();
345 int order1 =
m_base[1]->GetNumModes();
348 order0 * order1 * nquad2);
350 if (multiplybyweights)
357 tmp, outarray, wsp,
true,
true,
true);
363 inarray, outarray, wsp,
true,
true,
true);
377 bool doCheckCollDir0,
bool doCheckCollDir1,
bool doCheckCollDir2)
379 int nquad0 =
m_base[0]->GetNumPoints();
380 int nquad1 =
m_base[1]->GetNumPoints();
381 int nquad2 =
m_base[2]->GetNumPoints();
382 int nmodes0 =
m_base[0]->GetNumModes();
383 int nmodes1 =
m_base[1]->GetNumModes();
384 int nmodes2 =
m_base[2]->GetNumModes();
386 bool colldir0 = doCheckCollDir0 ? (
m_base[0]->Collocation()) :
false;
387 bool colldir1 = doCheckCollDir1 ? (
m_base[1]->Collocation()) :
false;
388 bool colldir2 = doCheckCollDir2 ? (
m_base[2]->Collocation()) :
false;
390 if (colldir0 && colldir1 && colldir2)
396 ASSERTL1(wsp.size() >= nmodes0 * nquad2 * (nquad1 + nmodes1),
397 "Insufficient workspace size");
405 for (
int n = 0; n < nmodes0; ++n)
407 Vmath::Vcopy(nquad1 * nquad2, inarray.data() + n, nquad0,
408 tmp0.data() + nquad1 * nquad2 * n, 1);
413 Blas::Dgemm(
'T',
'N', nquad1 * nquad2, nmodes0, nquad0, 1.0,
414 inarray.data(), nquad0, base0.data(), nquad0, 0.0,
415 tmp0.data(), nquad1 * nquad2);
421 for (
int n = 0; n < nmodes1; ++n)
424 tmp1.data() + nquad2 * nmodes0 * n, 1);
429 Blas::Dgemm(
'T',
'N', nquad2 * nmodes0, nmodes1, nquad1, 1.0,
430 tmp0.data(), nquad1, base1.data(), nquad1, 0.0,
431 tmp1.data(), nquad2 * nmodes0);
437 for (
int n = 0; n < nmodes2; ++n)
439 Vmath::Vcopy(nmodes0 * nmodes1, tmp1.data() + n, nquad2,
440 outarray.data() + nmodes0 * nmodes1 * n, 1);
445 Blas::Dgemm(
'T',
'N', nmodes0 * nmodes1, nmodes2, nquad2, 1.0,
446 tmp1.data(), nquad2, base2.data(), nquad2, 0.0,
447 outarray.data(), nmodes0 * nmodes1);
463 ASSERTL0((dir == 0) || (dir == 1) || (dir == 2),
464 "input dir is out of range");
466 int nquad1 =
m_base[1]->GetNumPoints();
467 int nquad2 =
m_base[2]->GetNumPoints();
468 int order0 =
m_base[0]->GetNumModes();
469 int order1 =
m_base[1]->GetNumModes();
473 if (outarray.size() < inarray.size())
490 m_base[2]->GetBdata(), tmp, outarray, wsp,
false,
true,
true);
495 m_base[2]->GetBdata(), tmp, outarray, wsp,
true,
false,
true);
500 m_base[2]->GetDbdata(), tmp, outarray, wsp,
true,
true,
false);
526 int nquad0 =
m_base[0]->GetNumPoints();
527 int nquad1 =
m_base[1]->GetNumPoints();
528 int nquad2 =
m_base[2]->GetNumPoints();
534 int btmp0 =
m_base[0]->GetNumModes();
535 int btmp1 =
m_base[1]->GetNumModes();
536 int mode2 = mode / (btmp0 * btmp1);
537 int mode1 = (mode - mode2 * btmp0 * btmp1) / btmp0;
538 int mode0 = (mode - mode2 * btmp0 * btmp1) % btmp0;
540 ASSERTL2(mode == mode2 * btmp0 * btmp1 + mode1 * btmp0 + mode0,
541 "Mode lookup failed.");
543 "Calling argument mode is larger than total expansion "
546 for (
int i = 0; i < nquad1 * nquad2; ++i)
549 &outarray[0] + i * nquad0, 1);
552 for (
int j = 0; j < nquad2; ++j)
554 for (
int i = 0; i < nquad0; ++i)
557 &outarray[0] + i + j * nquad0 * nquad1, nquad0,
558 &outarray[0] + i + j * nquad0 * nquad1, nquad0);
562 for (
int i = 0; i < nquad2; i++)
564 Blas::Dscal(nquad0 * nquad1, base2[mode2 * nquad2 + i],
565 &outarray[0] + i * nquad0 * nquad1, 1);
579 const int nm0 =
m_base[0]->GetNumModes();
580 const int nm1 =
m_base[1]->GetNumModes();
581 const int mode2 = mode / (nm0 * nm1);
582 const int mode1 = (mode - mode2 * nm0 * nm1) / nm0;
583 const int mode0 = (mode - mode2 * nm0 * nm1) % nm0;
585 return StdExpansion::BaryEvaluateBasis<0>(coords[0], mode0) *
586 StdExpansion::BaryEvaluateBasis<1>(coords[1], mode1) *
587 StdExpansion::BaryEvaluateBasis<2>(coords[2], mode2);
614 "BasisType is not a boundary interior form");
617 "BasisType is not a boundary interior form");
620 "BasisType is not a boundary interior form");
622 int nmodes0 =
m_base[0]->GetNumModes();
623 int nmodes1 =
m_base[1]->GetNumModes();
624 int nmodes2 =
m_base[2]->GetNumModes();
626 return (2 * (nmodes0 * nmodes1 + nmodes0 * nmodes2 + nmodes1 * nmodes2) -
627 4 * (nmodes0 + nmodes1 + nmodes2) + 8);
634 "BasisType is not a boundary interior form");
637 "BasisType is not a boundary interior form");
640 "BasisType is not a boundary interior form");
642 int nmodes0 =
m_base[0]->GetNumModes();
643 int nmodes1 =
m_base[1]->GetNumModes();
644 int nmodes2 =
m_base[2]->GetNumModes();
646 return 2 * (nmodes0 * nmodes1 + nmodes0 * nmodes2 + nmodes1 * nmodes2);
651 ASSERTL2((i >= 0) && (i <= 5),
"face id is out of range");
652 if ((i == 0) || (i == 5))
656 else if ((i == 1) || (i == 3))
668 ASSERTL2((i >= 0) && (i <= 5),
"face id is out of range");
669 if ((i == 0) || (i == 5))
673 else if ((i == 1) || (i == 3))
685 ASSERTL2(i >= 0 && i <= 5,
"face id is out of range");
687 if (i == 0 || i == 5)
689 return m_base[0]->GetNumPoints() *
m_base[1]->GetNumPoints();
691 else if (i == 1 || i == 3)
693 return m_base[0]->GetNumPoints() *
m_base[2]->GetNumPoints();
697 return m_base[1]->GetNumPoints() *
m_base[2]->GetNumPoints();
704 ASSERTL2(i >= 0 && i <= 5,
"face id is out of range");
705 ASSERTL2(j == 0 || j == 1,
"face direction is out of range");
707 if (i == 0 || i == 5)
709 return m_base[j]->GetPointsKey();
711 else if (i == 1 || i == 3)
713 return m_base[2 * j]->GetPointsKey();
717 return m_base[j + 1]->GetPointsKey();
722 const std::vector<unsigned int> &nummodes,
int &modes_offset)
724 int nmodes = nummodes[modes_offset] * nummodes[modes_offset + 1] *
725 nummodes[modes_offset + 2];
732 const int i,
const int k, [[maybe_unused]]
bool UseGLL)
const
734 ASSERTL2(i >= 0 && i <= 5,
"face id is out of range");
735 ASSERTL2(k >= 0 && k <= 1,
"basis key id is out of range");
756 m_base[dir]->GetNumModes());
772 for (
int k = 0; k < Qz; ++k)
774 for (
int j = 0; j < Qy; ++j)
776 for (
int i = 0; i < Qx; ++i)
778 int s = i + Qx * (j + Qy * k);
790 int nummodes[3] = {
m_base[0]->GetNumModes(),
m_base[1]->GetNumModes(),
791 m_base[2]->GetNumModes()};
797 numModes0 = nummodes[0];
798 numModes1 = nummodes[1];
804 numModes0 = nummodes[0];
805 numModes1 = nummodes[2];
811 numModes0 = nummodes[1];
812 numModes1 = nummodes[2];
817 ASSERTL0(
false,
"fid out of range");
824 std::swap(numModes0, numModes1);
839 "BasisType is not a boundary interior form");
842 "BasisType is not a boundary interior form");
845 "BasisType is not a boundary interior form");
847 ASSERTL1((localVertexId >= 0) && (localVertexId < 8),
848 "local vertex id must be between 0 and 7");
855 int nummodes[3] = {
m_base[0]->GetNumModes(),
m_base[1]->GetNumModes(),
856 m_base[2]->GetNumModes()};
858 if (useCoeffPacking ==
true)
860 if (localVertexId > 3)
872 switch (localVertexId % 4)
919 if ((localVertexId % 4) % 3 > 0)
931 if (localVertexId % 4 > 1)
944 if (localVertexId > 3)
957 return r * nummodes[0] * nummodes[1] +
q * nummodes[0] +
p;
967 "BasisType is not a boundary interior form");
970 "BasisType is not a boundary interior form");
973 "BasisType is not a boundary interior form");
976 int nummodes[3] = {
m_base[0]->GetNumModes(),
m_base[1]->GetNumModes(),
977 m_base[2]->GetNumModes()};
981 if (outarray.size() != nIntCoeffs)
994 for (i = 0; i < 3; i++)
999 IntIdx[i][1] = nummodes[i];
1004 IntIdx[i][1] = nummodes[i] - 1;
1008 for (r = IntIdx[2][0]; r < IntIdx[2][1]; r++)
1010 for (
q = IntIdx[1][0];
q < IntIdx[1][1];
q++)
1012 for (
p = IntIdx[0][0];
p < IntIdx[0][1];
p++)
1015 r * nummodes[0] * nummodes[1] +
q * nummodes[0] +
p;
1028 "BasisType is not a boundary interior form");
1031 "BasisType is not a boundary interior form");
1034 "BasisType is not a boundary interior form");
1037 int nummodes[3] = {
m_base[0]->GetNumModes(),
m_base[1]->GetNumModes(),
1038 m_base[2]->GetNumModes()};
1042 if (outarray.size() != nBndCoeffs)
1056 for (i = 0; i < 3; i++)
1064 IntIdx[i][1] = nummodes[i];
1068 BndIdx[i][1] = nummodes[i] - 1;
1070 IntIdx[i][1] = nummodes[i] - 1;
1074 for (i = 0; i < 2; i++)
1077 for (
q = 0;
q < nummodes[1];
q++)
1079 for (
p = 0;
p < nummodes[0];
p++)
1082 r * nummodes[0] * nummodes[1] +
q * nummodes[0] +
p;
1087 for (r = IntIdx[2][0]; r < IntIdx[2][1]; r++)
1089 for (i = 0; i < 2; i++)
1092 for (
p = 0;
p < nummodes[0];
p++)
1095 r * nummodes[0] * nummodes[1] +
q * nummodes[0] +
p;
1099 for (
q = IntIdx[1][0];
q < IntIdx[1][1];
q++)
1101 for (i = 0; i < 2; i++)
1105 r * nummodes[0] * nummodes[1] +
q * nummodes[0] +
p;
1110 sort(outarray.data(), outarray.data() + nBndCoeffs);
1116 std::array<NekDouble, 3> &firstOrderDerivs)
1128 int nummodesA = 0, nummodesB = 0;
1132 "Method only implemented if BasisType is indentical in "
1136 "Method only implemented for Modified_A or "
1137 "GLL_Lagrange BasisType");
1139 const int nummodes0 =
m_base[0]->GetNumModes();
1140 const int nummodes1 =
m_base[1]->GetNumModes();
1141 const int nummodes2 =
m_base[2]->GetNumModes();
1147 nummodesA = nummodes0;
1148 nummodesB = nummodes1;
1152 nummodesA = nummodes0;
1153 nummodesB = nummodes2;
1157 nummodesA = nummodes1;
1158 nummodesB = nummodes2;
1161 ASSERTL0(
false,
"fid must be between 0 and 5");
1164 int nFaceCoeffs = nummodesA * nummodesB;
1166 if (maparray.size() != nFaceCoeffs)
1183 offset = nummodes0 * nummodes1;
1187 offset = (nummodes2 - 1) * nummodes0 * nummodes1;
1205 offset = nummodes0 * (nummodes1 - 1);
1206 jump1 = nummodes0 * nummodes1;
1212 jump1 = nummodes0 * nummodes1;
1223 offset = nummodes0 - 1;
1224 jump1 = nummodes0 * nummodes1;
1231 jump1 = nummodes0 * nummodes1;
1236 ASSERTL0(
false,
"fid must be between 0 and 5");
1239 for (i = 0; i < nummodesB; i++)
1241 for (j = 0; j < nummodesA; j++)
1243 maparray[i * nummodesA + j] = i * jump1 + j * jump2 + offset;
1254 int nummodesA = 0, nummodesB = 0;
1258 "Method only implemented if BasisType is indentical in "
1262 "Method only implemented for Modified_A or "
1263 "GLL_Lagrange BasisType");
1265 const int nummodes0 =
m_base[0]->GetNumModes();
1266 const int nummodes1 =
m_base[1]->GetNumModes();
1267 const int nummodes2 =
m_base[2]->GetNumModes();
1273 nummodesA = nummodes0;
1274 nummodesB = nummodes1;
1278 nummodesA = nummodes0;
1279 nummodesB = nummodes2;
1283 nummodesA = nummodes1;
1284 nummodesB = nummodes2;
1287 ASSERTL0(
false,
"fid must be between 0 and 5");
1299 if (modified ==
false)
1301 ASSERTL1((
P == nummodesA) || (Q == nummodesB),
1302 "Different trace space face dimention "
1303 "and element face dimention not possible for "
1304 "GLL-Lagrange bases");
1307 int nFaceCoeffs =
P * Q;
1309 if (maparray.size() != nFaceCoeffs)
1315 for (i = 0; i < nFaceCoeffs; ++i)
1320 if (signarray.size() != nFaceCoeffs)
1326 fill(signarray.data(), signarray.data() + nFaceCoeffs, 1);
1331 for (i = 0; i < Q; i++)
1333 for (j = 0; j <
P; j++)
1337 arrayindx[i *
P + j] = i *
P + j;
1341 arrayindx[i *
P + j] = j * Q + i;
1348 for (i = 0; i < nummodesB; i++)
1350 for (j = nummodesA; j <
P; j++)
1352 signarray[arrayindx[i *
P + j]] = 0.0;
1353 maparray[arrayindx[i *
P + j]] = maparray[0];
1357 for (i = nummodesB; i < Q; i++)
1359 for (j = 0; j <
P; j++)
1361 signarray[arrayindx[i *
P + j]] = 0.0;
1362 maparray[arrayindx[i *
P + j]] = maparray[0];
1368 for (i = 0; i < Q; i++)
1372 for (j = 0; j <
P; j++)
1374 maparray[arrayindx[i *
P + j]] = i * nummodesA + j;
1378 for (j = nummodesA; j <
P; j++)
1380 signarray[arrayindx[i *
P + j]] = 0.0;
1381 maparray[arrayindx[i *
P + j]] = maparray[0];
1386 for (i = nummodesB; i < Q; i++)
1388 for (j = 0; j <
P; j++)
1390 signarray[arrayindx[i *
P + j]] = 0.0;
1391 maparray[arrayindx[i *
P + j]] = maparray[0];
1405 for (i = 3; i < Q; i += 2)
1407 for (j = 0; j <
P; j++)
1409 signarray[arrayindx[i *
P + j]] *= -1;
1413 for (i = 0; i <
P; i++)
1415 swap(maparray[i], maparray[i +
P]);
1416 swap(signarray[i], signarray[i +
P]);
1421 for (i = 0; i <
P; i++)
1423 for (j = 0; j < Q / 2; j++)
1425 swap(maparray[i + j *
P],
1426 maparray[i +
P * Q -
P - j *
P]);
1427 swap(signarray[i + j *
P],
1428 signarray[i +
P * Q -
P - j *
P]);
1437 for (i = 0; i < Q; i++)
1439 for (j = 3; j <
P; j += 2)
1441 signarray[arrayindx[i *
P + j]] *= -1;
1445 for (i = 0; i < Q; i++)
1447 swap(maparray[i], maparray[i + Q]);
1448 swap(signarray[i], signarray[i + Q]);
1453 for (i = 0; i <
P; i++)
1455 for (j = 0; j < Q / 2; j++)
1457 swap(maparray[i * Q + j], maparray[i * Q + Q - 1 - j]);
1458 swap(signarray[i * Q + j],
1459 signarray[i * Q + Q - 1 - j]);
1475 for (i = 0; i < Q; i++)
1477 for (j = 3; j <
P; j += 2)
1479 signarray[arrayindx[i *
P + j]] *= -1;
1483 for (i = 0; i < Q; i++)
1485 swap(maparray[i *
P], maparray[i *
P + 1]);
1486 swap(signarray[i *
P], signarray[i *
P + 1]);
1491 for (i = 0; i < Q; i++)
1493 for (j = 0; j <
P / 2; j++)
1495 swap(maparray[i *
P + j], maparray[i *
P +
P - 1 - j]);
1496 swap(signarray[i *
P + j],
1497 signarray[i *
P +
P - 1 - j]);
1506 for (i = 3; i < Q; i += 2)
1508 for (j = 0; j <
P; j++)
1510 signarray[arrayindx[i *
P + j]] *= -1;
1514 for (i = 0; i <
P; i++)
1516 swap(maparray[i * Q], maparray[i * Q + 1]);
1517 swap(signarray[i * Q], signarray[i * Q + 1]);
1522 for (i = 0; i < Q; i++)
1524 for (j = 0; j <
P / 2; j++)
1526 swap(maparray[i + j * Q],
1527 maparray[i +
P * Q - Q - j * Q]);
1528 swap(signarray[i + j * Q],
1529 signarray[i +
P * Q - Q - j * Q]);
1549 "BasisType is not a boundary interior form");
1552 "BasisType is not a boundary interior form");
1555 "BasisType is not a boundary interior form");
1558 "local edge id must be between 0 and 11");
1562 if (maparray.size() != nEdgeIntCoeffs)
1567 if (signarray.size() != nEdgeIntCoeffs)
1573 fill(signarray.data(), signarray.data() + nEdgeIntCoeffs, 1);
1576 int nummodes[3] = {
m_base[0]->GetNumModes(),
m_base[1]->GetNumModes(),
1577 m_base[2]->GetNumModes()};
1582 bool reverseOrdering =
false;
1583 bool signChange =
false;
1585 int IdxRange[3][2] = {{0, 0}, {0, 0}, {0, 0}};
1605 IdxRange[2][0] = nummodes[2] - 1;
1606 IdxRange[2][1] = nummodes[2];
1623 IdxRange[2][1] = nummodes[2] - 1;
1627 reverseOrdering =
true;
1633 IdxRange[2][1] = nummodes[2];
1662 IdxRange[1][0] = nummodes[1] - 1;
1663 IdxRange[1][1] = nummodes[1];
1678 IdxRange[1][1] = nummodes[1] - 1;
1682 reverseOrdering =
true;
1688 IdxRange[1][1] = nummodes[1];
1703 IdxRange[1][1] = nummodes[1] - 1;
1707 reverseOrdering =
true;
1713 IdxRange[1][1] = nummodes[1];
1742 IdxRange[0][0] = nummodes[0] - 1;
1743 IdxRange[0][1] = nummodes[0];
1758 IdxRange[0][1] = nummodes[0] - 1;
1762 reverseOrdering =
true;
1768 IdxRange[0][1] = nummodes[0];
1783 IdxRange[0][1] = nummodes[0] - 1;
1787 reverseOrdering =
true;
1793 IdxRange[0][1] = nummodes[0];
1806 for (
int r = IdxRange[2][0]; r < IdxRange[2][1]; r++)
1808 for (
int q = IdxRange[1][0];
q < IdxRange[1][1];
q++)
1810 for (
int p = IdxRange[0][0];
p < IdxRange[0][1];
p++)
1813 r * nummodes[0] * nummodes[1] +
q * nummodes[0] +
p;
1818 if (reverseOrdering)
1820 reverse(maparray.data(), maparray.data() + nEdgeIntCoeffs);
1825 for (
int p = 1;
p < nEdgeIntCoeffs;
p += 2)
1842 "BasisType is not a boundary interior form");
1845 "BasisType is not a boundary interior form");
1848 "BasisType is not a boundary interior form");
1850 ASSERTL1((fid >= 0) && (fid < 6),
"local face id must be between 0 and 5");
1854 if (maparray.size() != nFaceIntCoeffs)
1859 if (signarray.size() != nFaceIntCoeffs)
1865 fill(signarray.data(), signarray.data() + nFaceIntCoeffs, 1);
1868 int nummodes[3] = {
m_base[0]->GetNumModes(),
m_base[1]->GetNumModes(),
1869 m_base[2]->GetNumModes()};
1884 nummodesA = nummodes[0];
1885 nummodesB = nummodes[1];
1891 nummodesA = nummodes[0];
1892 nummodesB = nummodes[2];
1898 nummodesA = nummodes[1];
1899 nummodesB = nummodes[2];
1907 for (
int i = 0; i < (nummodesB - 2); i++)
1909 for (
int j = 0; j < (nummodesA - 2); j++)
1913 arrayindx[i * (nummodesA - 2) + j] = i * (nummodesA - 2) + j;
1917 arrayindx[i * (nummodesA - 2) + j] = j * (nummodesB - 2) + i;
1944 IdxRange[2][0] = nummodes[2] - 1;
1945 IdxRange[2][1] = nummodes[2];
1962 IdxRange[2][0] = nummodes[2] - 2;
1969 IdxRange[2][1] = nummodes[2] - 1;
1976 IdxRange[2][1] = nummodes[2];
1981 for (
int i = 3; i < nummodes[2]; i += 2)
2005 IdxRange[1][0] = nummodes[1] - 1;
2006 IdxRange[1][1] = nummodes[1];
2024 IdxRange[1][0] = nummodes[1] - 2;
2031 IdxRange[1][1] = nummodes[1] - 1;
2038 IdxRange[1][1] = nummodes[1];
2043 for (
int i = 3; i < nummodes[1]; i += 2)
2057 IdxRange[1][0] = nummodes[1] - 2;
2064 IdxRange[1][1] = nummodes[1] - 1;
2071 IdxRange[1][1] = nummodes[1];
2076 for (
int i = 3; i < nummodes[1]; i += 2)
2098 IdxRange[0][0] = nummodes[0] - 1;
2099 IdxRange[0][1] = nummodes[0];
2116 IdxRange[0][0] = nummodes[0] - 2;
2123 IdxRange[0][1] = nummodes[0] - 1;
2130 IdxRange[0][1] = nummodes[0];
2135 for (
int i = 3; i < nummodes[0]; i += 2)
2146 for (
int r = IdxRange[2][0]; r != IdxRange[2][1]; r += Incr[2])
2148 for (
int q = IdxRange[1][0];
q != IdxRange[1][1];
q += Incr[1])
2150 for (
int p = IdxRange[0][0];
p != IdxRange[0][1];
p += Incr[0])
2152 maparray[arrayindx[cnt]] =
2153 r * nummodes[0] * nummodes[1] +
q * nummodes[0] +
p;
2154 signarray[arrayindx[cnt++]] = sign0[
p] * sign1[
q] * sign2[r];
2162 ASSERTL2((i >= 0) && (i <= 11),
"edge id is out of range");
2164 if ((i == 0) || (i == 2) || (i == 8) || (i == 10))
2168 else if ((i == 1) || (i == 3) || (i == 9) || (i == 11))
2189 int nq0 =
m_base[0]->GetNumPoints();
2190 int nq1 =
m_base[1]->GetNumPoints();
2191 int nq2 =
m_base[2]->GetNumPoints();
2201 nq = max(nq0, max(nq1, nq2));
2215 for (
int i = 0; i < nq; ++i)
2217 for (
int j = 0; j < nq; ++j)
2219 for (
int k = 0; k < nq; ++k, ++cnt)
2222 coords[cnt][0] = -1.0 + 2 * k / (
NekDouble)(nq - 1);
2223 coords[cnt][1] = -1.0 + 2 * j / (
NekDouble)(nq - 1);
2224 coords[cnt][2] = -1.0 + 2 * i / (
NekDouble)(nq - 1);
2229 for (
int i = 0; i < neq; ++i)
2233 I[0] =
m_base[0]->GetI(coll);
2234 I[1] =
m_base[1]->GetI(coll + 1);
2235 I[2] =
m_base[2]->GetI(coll + 2);
2239 for (
int k = 0; k < nq2; ++k)
2241 for (
int j = 0; j < nq1; ++j)
2244 fac = (I[1]->GetPtr())[j] * (I[2]->GetPtr())[k];
2248 Mat->GetRawPtr() + k * nq0 * nq1 * neq +
2312 int nquad0 =
m_base[0]->GetNumPoints();
2313 int nquad1 =
m_base[1]->GetNumPoints();
2314 int nquad2 =
m_base[2]->GetNumPoints();
2315 int nq01 = nquad0 * nquad1;
2316 int nq12 = nquad1 * nquad2;
2322 for (
int i = 0; i < nq12; ++i)
2324 Vmath::Vmul(nquad0, inarray.data() + i * nquad0, 1, w0.data(), 1,
2325 outarray.data() + i * nquad0, 1);
2328 for (
int i = 0; i < nq12; ++i)
2330 Vmath::Smul(nquad0, w1[i % nquad1], outarray.data() + i * nquad0, 1,
2331 outarray.data() + i * nquad0, 1);
2334 for (
int i = 0; i < nquad2; ++i)
2336 Vmath::Smul(nq01, w2[i], outarray.data() + i * nq01, 1,
2337 outarray.data() + i * nq01, 1);
2345 int qa =
m_base[0]->GetNumPoints();
2346 int qb =
m_base[1]->GetNumPoints();
2347 int qc =
m_base[2]->GetNumPoints();
2348 int nmodes_a =
m_base[0]->GetNumModes();
2349 int nmodes_b =
m_base[1]->GetNumModes();
2350 int nmodes_c =
m_base[2]->GetNumModes();
2365 OrthoExp.
FwdTrans(array, orthocoeffs);
2375 for (
int i = 0; i < nmodes_a; ++i)
2377 for (
int j = 0; j < nmodes_b; ++j)
2380 pow((1.0 * i) / (nmodes_a - 1), cutoff * nmodes_a),
2381 pow((1.0 * j) / (nmodes_b - 1), cutoff * nmodes_b));
2383 for (
int k = 0; k < nmodes_c; ++k)
2386 std::max(fac1, pow((1.0 * k) / (nmodes_c - 1),
2387 cutoff * nmodes_c));
2389 orthocoeffs[cnt] *= SvvDiffCoeff * fac;
2405 max_abc = max(max_abc, 0);
2408 for (
int i = 0; i < nmodes_a; ++i)
2410 for (
int j = 0; j < nmodes_b; ++j)
2412 int maxij = max(i, j);
2414 for (
int k = 0; k < nmodes_c; ++k)
2416 int maxijk = max(maxij, k);
2430 min(nmodes_a, nmodes_b));
2433 int nmodes = max(nmodes_a, nmodes_b);
2434 nmodes = max(nmodes, nmodes_c);
2437 for (
int j = cutoff; j < nmodes; ++j)
2439 fac[j] = fabs((j - nmodes) / ((
NekDouble)(j - cutoff + 1.0)));
2443 for (
int i = 0; i < nmodes_a; ++i)
2445 for (
int j = 0; j < nmodes_b; ++j)
2447 for (
int k = 0; k < nmodes_c; ++k)
2449 if ((i >= cutoff) || (j >= cutoff) || (k >= cutoff))
2451 orthocoeffs[i * nmodes_a * nmodes_b + j * nmodes_c +
2453 (SvvDiffCoeff * exp(-(fac[i] + fac[j] + fac[k])));
2457 orthocoeffs[i * nmodes_a * nmodes_b + j * nmodes_c +
2466 OrthoExp.
BwdTrans(orthocoeffs, array);
2475 int qa =
m_base[0]->GetNumPoints();
2476 int qb =
m_base[1]->GetNumPoints();
2477 int qc =
m_base[2]->GetNumPoints();
2478 int nmodesA =
m_base[0]->GetNumModes();
2479 int nmodesB =
m_base[1]->GetNumModes();
2480 int nmodesC =
m_base[2]->GetNumModes();
2481 int P = nmodesA - 1;
2482 int Q = nmodesB - 1;
2483 int R = nmodesC - 1;
2496 int Pcut = cutoff *
P;
2497 int Qcut = cutoff * Q;
2498 int Rcut = cutoff * R;
2502 OrthoExp.
FwdTrans(array, orthocoeffs);
2507 for (
int i = 0; i < nmodesA; ++i)
2509 for (
int j = 0; j < nmodesB; ++j)
2511 for (
int k = 0; k < nmodesC; ++k, ++index)
2514 if (i > Pcut || j > Qcut || k > Rcut)
2519 fac = max(max(fac1, fac2), fac3);
2520 fac = pow(fac, exponent);
2521 orthocoeffs[index] *= exp(-alpha * fac);
2528 OrthoExp.
BwdTrans(orthocoeffs, array);
2534 int np0 =
m_base[0]->GetNumPoints();
2535 int np1 =
m_base[1]->GetNumPoints();
2536 int np2 =
m_base[2]->GetNumPoints();
2537 int np = max(np0, max(np1, np2));
2545 for (
int i = 0; i < np - 1; ++i)
2547 for (
int j = 0; j < np - 1; ++j)
2550 for (
int k = 0; k < np - 1; ++k)
2552 conn[cnt++] = plane + row + k;
2553 conn[cnt++] = plane + row + k + 1;
2554 conn[cnt++] = plane + rowp1 + k;
2556 conn[cnt++] = plane + rowp1 + k + 1;
2557 conn[cnt++] = plane + rowp1 + k;
2558 conn[cnt++] = plane + row + k + 1;
#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.
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)
int GetEdgeNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th edge.
void v_HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
NekDouble BaryTensorDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs)
void v_LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
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.
void WeakDerivMatrixOp_MatFree(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
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.
void BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
this function calculates the inner product of a given function f with the different modes of the expa...
void LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space.
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
void MassMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Class representing a hexehedral element in reference space.
int v_NumDGBndryCoeffs() const override
void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey) override
void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2) override
Differentiation Methods.
void v_ExponentialFilter(Array< OneD, NekDouble > &array, const NekDouble alpha, const NekDouble exponent, const NekDouble cutoff) override
int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false) override
void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta) override
int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset) 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_GetTraceInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2) override
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey) override
int v_GetEdgeNcoeffs(const int i) const override
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multbyweights=true) override
StdHexExp(const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const LibUtilities::BasisKey &Bc)
int v_GetTraceNumPoints(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_GetTraceCoeffMap(const unsigned int fid, Array< OneD, unsigned int > &maparray) override
LibUtilities::ShapeType v_DetShapeType() const override
void v_GetTraceNumModes(const int fid, int &numModes0, int &numModes1, Orientation faceOrient=eDir1FwdDir1_Dir2FwdDir2) override
void v_LocCollapsedToLocCoord(const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi) override
void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey) override
int v_NumBndryCoeffs() const override
int v_GetNedges() const override
void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
int v_GetNtraces() const override
void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey) override
int v_GetNverts() const override
bool v_IsBoundaryInteriorExpansion() const override
void v_GetInteriorMap(Array< OneD, unsigned int > &outarray) override
void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray) override
void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray) override
void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true) override
void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey) override
int v_GetTraceIntNcoeffs(const int i) const override
void v_StdPhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2) override
void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_GetElmtTraceToTraceMap(const unsigned int fid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation faceOrient, int P, int Q) override
void v_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_GetEdgeInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2) override
void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey) override
const LibUtilities::BasisKey v_GetTraceBasisKey(const int i, const int k, bool useGLL=false) const override
void v_GetCoords(Array< OneD, NekDouble > &coords_x, Array< OneD, NekDouble > &coords_y, Array< OneD, NekDouble > &coords_z) override
DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey) override
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) 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_GetTraceNcoeffs(const int i) const override
void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
MatrixType GetMatrixType() const
NekDouble GetConstFactor(const ConstFactorType &factor) const
bool ConstFactorExists(const ConstFactorType &factor) const
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
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...
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
const int kSVVDGFiltermodesmin
const int kSVVDGFiltermodesmax
LibUtilities::BasisKey EvaluateQuadFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
const NekDouble kSVVDGFilter[9][11]
@ ePhysInterpToEquiSpaced
@ eDir1BwdDir2_Dir2BwdDir1
@ eDir1FwdDir1_Dir2FwdDir2
@ eDir1BwdDir1_Dir2BwdDir2
@ eDir1BwdDir2_Dir2FwdDir1
@ eDir1FwdDir1_Dir2BwdDir2
@ eDir1BwdDir1_Dir2FwdDir2
@ eDir1FwdDir2_Dir2FwdDir1
@ eDir1FwdDir2_Dir2BwdDir1
std::vector< double > q(NPUPPER *NPUPPER)
static Array< OneD, NekDouble > NullNekDouble1DArray
std::shared_ptr< DNekMat > DNekMatSharedPtr
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)