61 : StdExpansion(Ba.GetNumModes() * Bb.GetNumModes() * Bc.GetNumModes(), 3,
63 StdExpansion3D(Ba.GetNumModes() * Bb.GetNumModes() * Bc.GetNumModes(), Ba,
68 std::string(
"HexExpMatrix")),
69 m_staticCondMatrixManager(
std::bind(&
Expansion::CreateStaticCondMatrix,
70 this,
std::placeholders::_1),
71 std::string(
"HexExpStaticCondMatrix"))
81 : StdExpansion(T), StdExpansion3D(T), StdHexExp(T),
Expansion(T),
83 m_staticCondMatrixManager(T.m_staticCondMatrixManager)
103 int nquad0 =
m_base[0]->GetNumPoints();
104 int nquad1 =
m_base[1]->GetNumPoints();
105 int nquad2 =
m_base[2]->GetNumPoints();
115 (
NekDouble *)&inarray[0], 1, &tmp[0], 1);
120 (
NekDouble *)&inarray[0], 1, &tmp[0], 1);
124 returnVal = StdHexExp::v_Integral(tmp);
147 int nquad0 =
m_base[0]->GetNumPoints();
148 int nquad1 =
m_base[1]->GetNumPoints();
149 int nquad2 =
m_base[2]->GetNumPoints();
150 int ntot = nquad0 * nquad1 * nquad2;
158 StdHexExp::v_PhysDeriv(inarray, Diff0, Diff1, Diff2);
164 Vmath::Vmul(ntot, &df[0][0], 1, &Diff0[0], 1, &out_d0[0], 1);
165 Vmath::Vvtvp(ntot, &df[1][0], 1, &Diff1[0], 1, &out_d0[0], 1,
167 Vmath::Vvtvp(ntot, &df[2][0], 1, &Diff2[0], 1, &out_d0[0], 1,
173 Vmath::Vmul(ntot, &df[3][0], 1, &Diff0[0], 1, &out_d1[0], 1);
174 Vmath::Vvtvp(ntot, &df[4][0], 1, &Diff1[0], 1, &out_d1[0], 1,
176 Vmath::Vvtvp(ntot, &df[5][0], 1, &Diff2[0], 1, &out_d1[0], 1,
182 Vmath::Vmul(ntot, &df[6][0], 1, &Diff0[0], 1, &out_d2[0], 1);
183 Vmath::Vvtvp(ntot, &df[7][0], 1, &Diff1[0], 1, &out_d2[0], 1,
185 Vmath::Vvtvp(ntot, &df[8][0], 1, &Diff2[0], 1, &out_d2[0], 1,
193 Vmath::Smul(ntot, df[0][0], &Diff0[0], 1, &out_d0[0], 1);
194 Blas::Daxpy(ntot, df[1][0], &Diff1[0], 1, &out_d0[0], 1);
195 Blas::Daxpy(ntot, df[2][0], &Diff2[0], 1, &out_d0[0], 1);
200 Vmath::Smul(ntot, df[3][0], &Diff0[0], 1, &out_d1[0], 1);
201 Blas::Daxpy(ntot, df[4][0], &Diff1[0], 1, &out_d1[0], 1);
202 Blas::Daxpy(ntot, df[5][0], &Diff2[0], 1, &out_d1[0], 1);
207 Vmath::Smul(ntot, df[6][0], &Diff0[0], 1, &out_d2[0], 1);
208 Blas::Daxpy(ntot, df[7][0], &Diff1[0], 1, &out_d2[0], 1);
209 Blas::Daxpy(ntot, df[8][0], &Diff2[0], 1, &out_d2[0], 1);
249 ASSERTL1(
false,
"input dir is out of range");
262 int nquad0 =
m_base[0]->GetNumPoints();
263 int nquad1 =
m_base[1]->GetNumPoints();
264 int nquad2 =
m_base[2]->GetNumPoints();
265 int ntot = nquad0 * nquad1 * nquad2;
273 StdHexExp::v_PhysDeriv(inarray, Diff0, Diff1, Diff2);
278 Vmath::Vmul(ntot, &dfdir[0][0], 1, &Diff0[0], 1, &outarray[0], 1);
279 Vmath::Vvtvp(ntot, &dfdir[1][0], 1, &Diff1[0], 1, &outarray[0], 1,
281 Vmath::Vvtvp(ntot, &dfdir[2][0], 1, &Diff2[0], 1, &outarray[0], 1,
300 if (
m_base[0]->Collocation() &&
m_base[1]->Collocation() &&
317 out = (*matsys) * in;
374 int nquad0 =
m_base[0]->GetNumPoints();
375 int nquad1 =
m_base[1]->GetNumPoints();
376 int nquad2 =
m_base[2]->GetNumPoints();
377 int order0 =
m_base[0]->GetNumModes();
378 int order1 =
m_base[1]->GetNumModes();
381 order0 * order1 * nquad2);
383 if (multiplybyweights)
390 tmp, outarray, wsp,
true,
true,
true);
396 inarray, outarray, wsp,
true,
true,
true);
431 ASSERTL1((dir == 0) || (dir == 1) || (dir == 2),
"Invalid direction.");
433 const int nq0 =
m_base[0]->GetNumPoints();
434 const int nq1 =
m_base[1]->GetNumPoints();
435 const int nq2 =
m_base[2]->GetNumPoints();
436 const int nq = nq0 * nq1 * nq2;
437 const int nm0 =
m_base[0]->GetNumModes();
438 const int nm1 =
m_base[1]->GetNumModes();
458 m_base[2]->GetBdata(), tmp2, outarray, wsp,
462 m_base[2]->GetBdata(), tmp3, tmp5, wsp,
true,
467 m_base[2]->GetDbdata(), tmp4, tmp5, wsp,
true,
476 ASSERTL1((dir == 0) || (dir == 1) || (dir == 2),
"Invalid direction.");
478 const int nq0 =
m_base[0]->GetNumPoints();
479 const int nq1 =
m_base[1]->GetNumPoints();
480 const int nq2 =
m_base[2]->GetNumPoints();
481 const int nq = nq0 * nq1 * nq2;
496 Vmath::Vmul(nq, &df[3 * dir][0], 1, tmp1.data(), 1, tmp2.data(), 1);
497 Vmath::Vmul(nq, &df[3 * dir + 1][0], 1, tmp1.data(), 1, tmp3.data(), 1);
498 Vmath::Vmul(nq, &df[3 * dir + 2][0], 1, tmp1.data(), 1, tmp4.data(), 1);
502 Vmath::Smul(nq, df[3 * dir][0], tmp1.data(), 1, tmp2.data(), 1);
503 Vmath::Smul(nq, df[3 * dir + 1][0], tmp1.data(), 1, tmp3.data(), 1);
504 Vmath::Smul(nq, df[3 * dir + 2][0], tmp1.data(), 1, tmp4.data(), 1);
520 const int nq0 =
m_base[0]->GetNumPoints();
521 const int nq1 =
m_base[1]->GetNumPoints();
522 const int nq2 =
m_base[2]->GetNumPoints();
523 const int nq = nq0 * nq1 * nq2;
524 const int nm0 =
m_base[0]->GetNumModes();
525 const int nm1 =
m_base[1]->GetNumModes();
543 Vmath::Vmul(nq, &dfdir[0][0], 1, tmp1.data(), 1, tmp2.data(), 1);
544 Vmath::Vmul(nq, &dfdir[1][0], 1, tmp1.data(), 1, tmp3.data(), 1);
545 Vmath::Vmul(nq, &dfdir[2][0], 1, tmp1.data(), 1, tmp4.data(), 1);
548 m_base[2]->GetBdata(), tmp2, outarray, wsp,
552 m_base[2]->GetBdata(), tmp3, tmp5, wsp,
true,
558 m_base[2]->GetDbdata(), tmp4, tmp5, wsp,
true,
578 return StdExpansion3D::v_PhysEvaluate(Lcoord, physvals);
587 m_geom->GetLocCoords(coord, Lcoord);
588 return StdExpansion3D::v_PhysEvaluate(Lcoord, physvals);
594 std::array<NekDouble, 3> &firstOrderDerivs)
598 m_geom->GetLocCoords(coord, Lcoord);
599 return StdHexExp::v_PhysEvalFirstDeriv(Lcoord, inarray, firstOrderDerivs);
606 m_base[2]->GetBasisKey());
612 m_base[0]->GetPointsKey());
614 m_base[1]->GetPointsKey());
616 m_base[2]->GetPointsKey());
634 ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[0] <= 1.0 && Lcoords[1] >= -1.0 &&
635 Lcoords[1] <= 1.0 && Lcoords[2] >= -1.0 && Lcoords[2] <= 1.0,
636 "Local coordinates are not in region [-1,1]");
640 for (i = 0; i <
m_geom->GetCoordim(); ++i)
642 coords[i] =
m_geom->GetCoord(i, Lcoords);
664 const NekDouble *data,
const std::vector<unsigned int> &nummodes,
665 const int mode_offset,
NekDouble *coeffs,
666 std::vector<LibUtilities::BasisType> &fromType)
668 int data_order0 = nummodes[mode_offset];
669 int fillorder0 = min(
m_base[0]->GetNumModes(), data_order0);
670 int data_order1 = nummodes[mode_offset + 1];
671 int order1 =
m_base[1]->GetNumModes();
672 int fillorder1 = min(order1, data_order1);
673 int data_order2 = nummodes[mode_offset + 2];
674 int order2 =
m_base[2]->GetNumModes();
675 int fillorder2 = min(order2, data_order2);
687 m_base[0]->GetPointsKey()),
689 m_base[1]->GetPointsKey()),
691 m_base[2]->GetPointsKey()));
694 m_base[2]->GetBasisKey());
716 "Extraction routine not set up for this basis");
718 "Extraction routine not set up for this basis");
721 for (j = 0; j < fillorder0; ++j)
723 for (i = 0; i < fillorder1; ++i)
725 Vmath::Vcopy(fillorder2, &data[cnt], 1, &coeffs[cnt1], 1);
731 for (i = fillorder1; i < data_order1; ++i)
736 for (i = fillorder1; i < order1; ++i)
761 ASSERTL0(
false,
"basis is either not set up or not "
768 int nquad0 =
m_base[0]->GetNumPoints();
769 int nquad1 =
m_base[1]->GetNumPoints();
770 int nquad2 =
m_base[2]->GetNumPoints();
782 if (outarray.size() != nq0 * nq1)
787 for (
int i = 0; i < nquad0 * nquad1; ++i)
797 if (outarray.size() != nq0 * nq1)
803 for (
int k = 0; k < nquad2; k++)
805 for (
int i = 0; i < nquad0; ++i)
807 outarray[k * nquad0 + i] = nquad0 * nquad1 * k + i;
816 if (outarray.size() != nq0 * nq1)
821 for (
int i = 0; i < nquad1 * nquad2; i++)
823 outarray[i] = nquad0 - 1 + i * nquad0;
831 if (outarray.size() != nq0 * nq1)
836 for (
int k = 0; k < nquad2; k++)
838 for (
int i = 0; i < nquad0; i++)
840 outarray[k * nquad0 + i] =
841 (nquad0 * (nquad1 - 1)) + (k * nquad0 * nquad1) + i;
850 if (outarray.size() != nq0 * nq1)
855 for (
int i = 0; i < nquad1 * nquad2; i++)
857 outarray[i] = i * nquad0;
864 if (outarray.size() != nq0 * nq1)
869 for (
int i = 0; i < nquad0 * nquad1; i++)
871 outarray[i] = nquad0 * nquad1 * (nquad2 - 1) + i;
876 ASSERTL0(
false,
"face value (> 5) is out of range");
889 for (i = 0; i < ptsKeys.size(); ++i)
900 geomFactors->GetDerivFactors(ptsKeys);
913 for (i = 0; i < vCoordDim; ++i)
918 size_t nqb = nq_face;
932 for (i = 0; i < vCoordDim; ++i)
934 normal[i][0] = -df[3 * i + 2][0];
938 for (i = 0; i < vCoordDim; ++i)
940 normal[i][0] = -df[3 * i + 1][0];
944 for (i = 0; i < vCoordDim; ++i)
946 normal[i][0] = df[3 * i][0];
950 for (i = 0; i < vCoordDim; ++i)
952 normal[i][0] = df[3 * i + 1][0];
956 for (i = 0; i < vCoordDim; ++i)
958 normal[i][0] = -df[3 * i][0];
962 for (i = 0; i < vCoordDim; ++i)
964 normal[i][0] = df[3 * i + 2][0];
968 ASSERTL0(
false,
"face is out of range (edge < 5)");
973 for (i = 0; i < vCoordDim; ++i)
975 fac += normal[i][0] * normal[i][0];
977 fac = 1.0 /
sqrt(fac);
980 for (i = 0; i < vCoordDim; ++i)
982 Vmath::Fill(nq_face, fac * normal[i][0], normal[i], 1);
989 int nqe0 = ptsKeys[0].GetNumPoints();
990 int nqe1 = ptsKeys[1].GetNumPoints();
991 int nqe2 = ptsKeys[2].GetNumPoints();
992 int nqe01 = nqe0 * nqe1;
993 int nqe02 = nqe0 * nqe2;
994 int nqe12 = nqe1 * nqe2;
997 if (face == 0 || face == 5)
1001 else if (face == 1 || face == 3)
1022 for (j = 0; j < nqe; ++j)
1024 normals[j] = -df[2][j] * jac[j];
1025 normals[nqe + j] = -df[5][j] * jac[j];
1026 normals[2 * nqe + j] = -df[8][j] * jac[j];
1027 faceJac[j] = jac[j];
1030 points0 = ptsKeys[0];
1031 points1 = ptsKeys[1];
1034 for (j = 0; j < nqe0; ++j)
1036 for (k = 0; k < nqe2; ++k)
1038 int idx = j + nqe01 * k;
1039 normals[j + k * nqe0] = -df[1][idx] * jac[idx];
1040 normals[nqe + j + k * nqe0] = -df[4][idx] * jac[idx];
1041 normals[2 * nqe + j + k * nqe0] =
1042 -df[7][idx] * jac[idx];
1043 faceJac[j + k * nqe0] = jac[idx];
1046 points0 = ptsKeys[0];
1047 points1 = ptsKeys[2];
1050 for (j = 0; j < nqe1; ++j)
1052 for (k = 0; k < nqe2; ++k)
1054 int idx = nqe0 - 1 + nqe0 * j + nqe01 * k;
1055 normals[j + k * nqe1] = df[0][idx] * jac[idx];
1056 normals[nqe + j + k * nqe1] = df[3][idx] * jac[idx];
1057 normals[2 * nqe + j + k * nqe1] = df[6][idx] * jac[idx];
1058 faceJac[j + k * nqe1] = jac[idx];
1061 points0 = ptsKeys[1];
1062 points1 = ptsKeys[2];
1065 for (j = 0; j < nqe0; ++j)
1067 for (k = 0; k < nqe2; ++k)
1069 int idx = nqe0 * (nqe1 - 1) + j + nqe01 * k;
1070 normals[j + k * nqe0] = df[1][idx] * jac[idx];
1071 normals[nqe + j + k * nqe0] = df[4][idx] * jac[idx];
1072 normals[2 * nqe + j + k * nqe0] = df[7][idx] * jac[idx];
1073 faceJac[j + k * nqe0] = jac[idx];
1076 points0 = ptsKeys[0];
1077 points1 = ptsKeys[2];
1080 for (j = 0; j < nqe1; ++j)
1082 for (k = 0; k < nqe2; ++k)
1084 int idx = j * nqe0 + nqe01 * k;
1085 normals[j + k * nqe1] = -df[0][idx] * jac[idx];
1086 normals[nqe + j + k * nqe1] = -df[3][idx] * jac[idx];
1087 normals[2 * nqe + j + k * nqe1] =
1088 -df[6][idx] * jac[idx];
1089 faceJac[j + k * nqe1] = jac[idx];
1092 points0 = ptsKeys[1];
1093 points1 = ptsKeys[2];
1096 for (j = 0; j < nqe01; ++j)
1098 int idx = j + nqe01 * (nqe2 - 1);
1099 normals[j] = df[2][idx] * jac[idx];
1100 normals[nqe + j] = df[5][idx] * jac[idx];
1101 normals[2 * nqe + j] = df[8][idx] * jac[idx];
1102 faceJac[j] = jac[idx];
1104 points0 = ptsKeys[0];
1105 points1 = ptsKeys[1];
1108 ASSERTL0(
false,
"face is out of range (face < 5)");
1117 Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
1125 Vmath::Vmul(nq_face, work, 1, normal[i], 1, normal[i], 1);
1132 Vmath::Vvtvp(nq_face, normal[i], 1, normal[i], 1, work, 1, work, 1);
1142 Vmath::Vmul(nq_face, normal[i], 1, work, 1, normal[i], 1);
1154 StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
1169 StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
1177 StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
1184 StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(inarray, outarray, mkey);
1191 StdExpansion::MassLevelCurvatureMatrixOp_MatFree(inarray, outarray, mkey);
1216 int n_coeffs = inarray.size();
1217 int nmodes0 =
m_base[0]->GetNumModes();
1218 int nmodes1 =
m_base[1]->GetNumModes();
1219 int nmodes2 =
m_base[2]->GetNumModes();
1220 int numMax = nmodes0;
1248 int cnt = 0, cnt2 = 0;
1250 for (
int u = 0; u < numMin + 1; ++u)
1252 for (
int i = 0; i < numMin; ++i)
1255 tmp2 = coeff_tmp1 + cnt, 1);
1261 tmp4 = coeff_tmp2 + cnt2, 1);
1263 cnt2 = u * nmodes0 * nmodes1;
1291 StdHexExp::v_SVVLaplacianFilter(array, mkey);
1316 returnval = StdHexExp::v_GenMatrix(mkey);
1331 return tmp->GetStdMatrix(mkey);
1365 int nquad0 =
m_base[0]->GetNumPoints();
1366 int nquad1 =
m_base[1]->GetNumPoints();
1367 int nquad2 =
m_base[2]->GetNumPoints();
1368 int nqtot = nquad0 * nquad1 * nquad2;
1370 ASSERTL1(wsp.size() >= 6 * nqtot,
"Insufficient workspace size.");
1399 StdExpansion3D::PhysTensorDeriv(inarray, wsp0, wsp1, wsp2);
1405 Vmath::Vvtvvtp(nqtot, &metric00[0], 1, &wsp0[0], 1, &metric01[0], 1,
1406 &wsp1[0], 1, &wsp3[0], 1);
1407 Vmath::Vvtvp(nqtot, &metric02[0], 1, &wsp2[0], 1, &wsp3[0], 1, &wsp3[0], 1);
1408 Vmath::Vvtvvtp(nqtot, &metric01[0], 1, &wsp0[0], 1, &metric11[0], 1,
1409 &wsp1[0], 1, &wsp4[0], 1);
1410 Vmath::Vvtvp(nqtot, &metric12[0], 1, &wsp2[0], 1, &wsp4[0], 1, &wsp4[0], 1);
1411 Vmath::Vvtvvtp(nqtot, &metric02[0], 1, &wsp0[0], 1, &metric12[0], 1,
1412 &wsp1[0], 1, &wsp5[0], 1);
1413 Vmath::Vvtvp(nqtot, &metric22[0], 1, &wsp2[0], 1, &wsp5[0], 1, &wsp5[0], 1);
1438 const unsigned int dim = 3;
1444 for (
unsigned int i = 0; i < dim; ++i)
1446 for (
unsigned int j = i; j < dim; ++j)
1483 if (d0factors.size() != 6)
1490 if (d0factors[0].size() != nquad0 * nquad1)
1500 if (d0factors[1].size() != nquad0 * nquad2)
1510 if (d0factors[2].size() != nquad1 * nquad2)
1534 int ncoords = normal_0.size();
1539 for (
int i = 0; i < nquad0 * nquad1; ++i)
1541 d0factors[0][i] = df[0][i] * normal_0[0][i];
1542 d1factors[0][i] = df[1][i] * normal_0[0][i];
1543 d2factors[0][i] = df[2][i] * normal_0[0][i];
1546 df[0][nquad0 * nquad1 * (nquad2 - 1) + i] * normal_5[0][i];
1548 df[1][nquad0 * nquad1 * (nquad2 - 1) + i] * normal_5[0][i];
1550 df[2][nquad0 * nquad1 * (nquad2 - 1) + i] * normal_5[0][i];
1553 for (
int n = 1; n < ncoords; ++n)
1555 for (
int i = 0; i < nquad0 * nquad1; ++i)
1557 d0factors[0][i] += df[3 * n][i] * normal_0[n][i];
1558 d1factors[0][i] += df[3 * n + 1][i] * normal_0[n][i];
1559 d2factors[0][i] += df[3 * n + 2][i] * normal_0[n][i];
1562 df[3 * n][nquad0 * nquad1 * (nquad2 - 1) + i] *
1565 df[3 * n + 1][nquad0 * nquad1 * (nquad2 - 1) + i] *
1568 df[3 * n + 2][nquad0 * nquad1 * (nquad2 - 1) + i] *
1574 for (
int j = 0; j < nquad2; ++j)
1576 for (
int i = 0; i < nquad0; ++i)
1578 d0factors[1][j * nquad0 + i] = df[0][j * nquad0 * nquad1 + i] *
1579 normal_1[0][j * nquad0 + i];
1580 d1factors[1][j * nquad0 + i] = df[1][j * nquad0 * nquad1 + i] *
1581 normal_1[0][j * nquad0 + i];
1582 d2factors[1][j * nquad0 + i] = df[2][j * nquad0 * nquad1 + i] *
1583 normal_1[0][j * nquad0 + i];
1585 d0factors[3][j * nquad0 + i] =
1586 df[0][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1587 normal_3[0][j * nquad0 + i];
1588 d1factors[3][j * nquad0 + i] =
1589 df[1][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1590 normal_3[0][j * nquad0 + i];
1591 d2factors[3][j * nquad0 + i] =
1592 df[2][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1593 normal_3[0][j * nquad0 + i];
1597 for (
int n = 1; n < ncoords; ++n)
1599 for (
int j = 0; j < nquad2; ++j)
1601 for (
int i = 0; i < nquad0; ++i)
1603 d0factors[1][j * nquad0 + i] +=
1604 df[3 * n][j * nquad0 * nquad1 + i] *
1605 normal_1[0][j * nquad0 + i];
1606 d1factors[1][j * nquad0 + i] +=
1607 df[3 * n + 1][j * nquad0 * nquad1 + i] *
1608 normal_1[0][j * nquad0 + i];
1609 d2factors[1][j * nquad0 + i] +=
1610 df[3 * n + 2][j * nquad0 * nquad1 + i] *
1611 normal_1[0][j * nquad0 + i];
1613 d0factors[3][j * nquad0 + i] +=
1614 df[3 * n][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1615 normal_3[0][j * nquad0 + i];
1616 d1factors[3][j * nquad0 + i] +=
1617 df[3 * n + 1][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1618 normal_3[0][j * nquad0 + i];
1619 d2factors[3][j * nquad0 + i] +=
1620 df[3 * n + 2][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1621 normal_3[0][j * nquad0 + i];
1627 for (
int j = 0; j < nquad2; ++j)
1629 for (
int i = 0; i < nquad1; ++i)
1631 d0factors[2][j * nquad1 + i] =
1632 df[0][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1633 normal_2[0][j * nquad1 + i];
1634 d1factors[2][j * nquad1 + i] =
1635 df[1][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1636 normal_2[0][j * nquad1 + i];
1637 d2factors[2][j * nquad1 + i] =
1638 df[2][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1639 normal_2[0][j * nquad1 + i];
1641 d0factors[4][j * nquad1 + i] =
1642 df[0][j * nquad0 * nquad1 + i * nquad0] *
1643 normal_4[0][j * nquad1 + i];
1644 d1factors[4][j * nquad1 + i] =
1645 df[1][j * nquad0 * nquad1 + i * nquad0] *
1646 normal_4[0][j * nquad1 + i];
1647 d2factors[4][j * nquad1 + i] =
1648 df[2][j * nquad0 * nquad1 + i * nquad0] *
1649 normal_4[0][j * nquad1 + i];
1653 for (
int n = 1; n < ncoords; ++n)
1655 for (
int j = 0; j < nquad2; ++j)
1657 for (
int i = 0; i < nquad1; ++i)
1659 d0factors[2][j * nquad1 + i] +=
1660 df[3 * n][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1661 normal_2[n][j * nquad1 + i];
1662 d1factors[2][j * nquad1 + i] +=
1664 [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1665 normal_2[n][j * nquad1 + i];
1666 d2factors[2][j * nquad1 + i] +=
1668 [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1669 normal_2[n][j * nquad1 + i];
1671 d0factors[4][j * nquad1 + i] +=
1672 df[3 * n][i * nquad0 + j * nquad0 * nquad1] *
1673 normal_4[n][j * nquad1 + i];
1674 d1factors[4][j * nquad1 + i] +=
1675 df[3 * n + 1][i * nquad0 + j * nquad0 * nquad1] *
1676 normal_4[n][j * nquad1 + i];
1677 d2factors[4][j * nquad1 + i] +=
1678 df[3 * n + 2][i * nquad0 + j * nquad0 * nquad1] *
1679 normal_4[n][j * nquad1 + i];
1687 for (
int i = 0; i < nquad0 * nquad1; ++i)
1689 d0factors[0][i] = df[0][0] * normal_0[0][i];
1690 d0factors[5][i] = df[0][0] * normal_5[0][i];
1692 d1factors[0][i] = df[1][0] * normal_0[0][i];
1693 d1factors[5][i] = df[1][0] * normal_5[0][i];
1695 d2factors[0][i] = df[2][0] * normal_0[0][i];
1696 d2factors[5][i] = df[2][0] * normal_5[0][i];
1699 for (
int n = 1; n < ncoords; ++n)
1701 for (
int i = 0; i < nquad0 * nquad1; ++i)
1703 d0factors[0][i] += df[3 * n][0] * normal_0[n][i];
1704 d0factors[5][i] += df[3 * n][0] * normal_5[n][i];
1706 d1factors[0][i] += df[3 * n + 1][0] * normal_0[n][i];
1707 d1factors[5][i] += df[3 * n + 1][0] * normal_5[n][i];
1709 d2factors[0][i] += df[3 * n + 2][0] * normal_0[n][i];
1710 d2factors[5][i] += df[3 * n + 2][0] * normal_5[n][i];
1715 for (
int i = 0; i < nquad0 * nquad2; ++i)
1717 d0factors[1][i] = df[0][0] * normal_1[0][i];
1718 d0factors[3][i] = df[0][0] * normal_3[0][i];
1720 d1factors[1][i] = df[1][0] * normal_1[0][i];
1721 d1factors[3][i] = df[1][0] * normal_3[0][i];
1723 d2factors[1][i] = df[2][0] * normal_1[0][i];
1724 d2factors[3][i] = df[2][0] * normal_3[0][i];
1727 for (
int n = 1; n < ncoords; ++n)
1729 for (
int i = 0; i < nquad0 * nquad2; ++i)
1731 d0factors[1][i] += df[3 * n][0] * normal_1[n][i];
1732 d0factors[3][i] += df[3 * n][0] * normal_3[n][i];
1734 d1factors[1][i] += df[3 * n + 1][0] * normal_1[n][i];
1735 d1factors[3][i] += df[3 * n + 1][0] * normal_3[n][i];
1737 d2factors[1][i] += df[3 * n + 2][0] * normal_1[n][i];
1738 d2factors[3][i] += df[3 * n + 2][0] * normal_3[n][i];
1743 for (
int i = 0; i < nquad1 * nquad2; ++i)
1745 d0factors[2][i] = df[0][0] * normal_2[0][i];
1746 d0factors[4][i] = df[0][0] * normal_4[0][i];
1748 d1factors[2][i] = df[1][0] * normal_2[0][i];
1749 d1factors[4][i] = df[1][0] * normal_4[0][i];
1751 d2factors[2][i] = df[2][0] * normal_2[0][i];
1752 d2factors[4][i] = df[2][0] * normal_4[0][i];
1755 for (
int n = 1; n < ncoords; ++n)
1757 for (
int i = 0; i < nquad1 * nquad2; ++i)
1759 d0factors[2][i] += df[3 * n][0] * normal_2[n][i];
1760 d0factors[4][i] += df[3 * n][0] * normal_4[n][i];
1762 d1factors[2][i] += df[3 * n + 1][0] * normal_2[n][i];
1763 d1factors[4][i] += df[3 * n + 1][0] * normal_4[n][i];
1765 d2factors[2][i] += df[3 * n + 2][0] * normal_2[n][i];
1766 d2factors[4][i] += df[3 * n + 2][0] * normal_4[n][i];
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Describes the specification for a Basis.
int GetNumPoints() const
Return points order at which basis is defined.
PointsKey GetPointsKey() const
Return distribution of points.
Defines a specification for a set of points.
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
std::map< int, NormalVector > m_traceNormals
std::map< int, Array< OneD, NekDouble > > m_elmtBndNormDirElmtLen
the element length in each element boundary(Vertex, edge or face) normal direction calculated based o...
SpatialDomains::GeometrySharedPtr GetGeom() const
SpatialDomains::GeometrySharedPtr m_geom
void ComputeLaplacianMetric()
void ComputeGmatcdotMF(const Array< TwoD, const NekDouble > &df, const Array< OneD, const NekDouble > &direction, Array< OneD, Array< OneD, NekDouble > > &dfdir)
void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
void ComputeQuadratureMetric()
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
const NormalVector & GetTraceNormal(const int id)
void v_WeakDirectionalDerivMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray) override
Integrate the physical point list inarray over region.
void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
NekDouble v_PhysEvalFirstDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) override
void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey) override
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
void v_DropLocStaticCondMatrix(const MatrixKey &mkey) override
DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey) override
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true) override
Calculate the inner product of inarray with respect to the given basis B = base0 * base1 * base2.
void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Forward transform from physical quadrature space stored in inarray and evaluate the expansion coeffic...
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Calculates the inner product .
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey) override
void v_ComputeTraceNormal(const int face) override
HexExp(const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const LibUtilities::BasisKey &Bc, const SpatialDomains::HexGeomSharedPtr &geom)
Constructor using BasisKey class for quadrature points and order definition.
void v_ComputeLaplacianMetric() override
void v_ExtractDataToCoeffs(const NekDouble *data, const std::vector< unsigned int > &nummodes, const int mode_offset, NekDouble *coeffs, std::vector< LibUtilities::BasisType > &fromType) override
void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp) override
void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2) override
Calculate the derivative of the physical points.
void v_DropLocMatrix(const MatrixKey &mkey) override
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords) override
Retrieves the physical coordinates of a given set of reference coordinates.
NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &physvals) override
This function evaluates the expansion at a single (arbitrary) point of the domain.
void v_IProductWRTDirectionalDerivBase_SumFac(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray) override
void v_PhysDirectionalDeriv(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &direction, Array< OneD, NekDouble > &out) override
Physical derivative along a direction vector.
NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals) override
DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey) override
void v_GetTracePhysMap(const int face, Array< OneD, int > &outarray) override
LibUtilities::ShapeType v_DetShapeType() const override
Return the region shape using the enum-list of ShapeType.
void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const override
StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const override
void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Calculate the inner product of inarray with respect to the elements basis.
void v_NormalTraceDerivFactors(Array< OneD, Array< OneD, NekDouble > > &factors, Array< OneD, Array< OneD, NekDouble > > &d0factors, Array< OneD, Array< OneD, NekDouble > > &d1factors) override
: This method gets all of the factors which are required as part of the Gradient Jump Penalty stabili...
void v_MassLevelCurvatureMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
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 v_HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
void v_LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
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.
const LibUtilities::PointsKeyVector GetPointsKeys() const
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...
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)
const LibUtilities::BasisKey GetTraceBasisKey(const int i, int k=-1, bool UseGLL=false) const
This function returns the basis key belonging to the i-th trace.
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
Class representing a hexehedral element in reference space.
MatrixType GetMatrixType() const
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.
void Interp3D(const BasisKey &fbasis0, const BasisKey &fbasis1, const BasisKey &fbasis2, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, const BasisKey &tbasis2, Array< OneD, NekDouble > &to)
this function interpolates a 3D function evaluated at the quadrature points of the 3D basis,...
void InterpCoeff3D(const BasisKey &fbasis0, const BasisKey &fbasis1, const BasisKey &fbasis2, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, const BasisKey &tbasis2, Array< OneD, NekDouble > &to)
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis,...
std::vector< PointsKey > PointsKeyVector
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
@ eOrtho_A
Principle Orthogonal Functions .
@ eGLL_Lagrange
Lagrange for SEM basis .
@ eModified_A
Principle Modified Functions .
std::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
GeomType
Indicates the type of element geometry.
@ eRegular
Geometry is straight-sided with constant geometric factors.
@ eMovingRegular
Currently unused.
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< HexGeom > HexGeomSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
@ eInvLaplacianWithUnityMean
std::shared_ptr< StdHexExp > StdHexExpSharedPtr
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
static Array< OneD, NekDouble > NullNekDouble1DArray
std::shared_ptr< DNekMat > DNekMatSharedPtr
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
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 Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*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 Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/x.
void Vdiv(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 Zero(int n, T *x, const int incx)
Zero vector.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
scalarT< T > sqrt(scalarT< T > in)