48 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
51 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
56 std::string(
"PyrExpMatrix")),
57 m_staticCondMatrixManager(
std::bind(&
Expansion::CreateStaticCondMatrix,
58 this,
std::placeholders::_1),
59 std::string(
"PyrExpStaticCondMatrix"))
64 : StdExpansion(T), StdExpansion3D(T), StdPyrExp(T),
Expansion(T),
66 m_staticCondMatrixManager(T.m_staticCondMatrixManager)
96 int nquad0 =
m_base[0]->GetNumPoints();
97 int nquad1 =
m_base[1]->GetNumPoints();
98 int nquad2 =
m_base[2]->GetNumPoints();
106 (
NekDouble *)&inarray[0], 1, &tmp[0], 1);
111 (
NekDouble *)&inarray[0], 1, &tmp[0], 1);
115 return StdPyrExp::v_Integral(tmp);
127 int nquad0 =
m_base[0]->GetNumPoints();
128 int nquad1 =
m_base[1]->GetNumPoints();
129 int nquad2 =
m_base[2]->GetNumPoints();
136 StdPyrExp::v_PhysDeriv(inarray, diff0, diff1, diff2);
142 Vmath::Vmul(nquad0 * nquad1 * nquad2, &gmat[0][0], 1, &diff0[0], 1,
144 Vmath::Vvtvp(nquad0 * nquad1 * nquad2, &gmat[1][0], 1, &diff1[0], 1,
145 &out_d0[0], 1, &out_d0[0], 1);
146 Vmath::Vvtvp(nquad0 * nquad1 * nquad2, &gmat[2][0], 1, &diff2[0], 1,
147 &out_d0[0], 1, &out_d0[0], 1);
152 Vmath::Vmul(nquad0 * nquad1 * nquad2, &gmat[3][0], 1, &diff0[0], 1,
154 Vmath::Vvtvp(nquad0 * nquad1 * nquad2, &gmat[4][0], 1, &diff1[0], 1,
155 &out_d1[0], 1, &out_d1[0], 1);
156 Vmath::Vvtvp(nquad0 * nquad1 * nquad2, &gmat[5][0], 1, &diff2[0], 1,
157 &out_d1[0], 1, &out_d1[0], 1);
162 Vmath::Vmul(nquad0 * nquad1 * nquad2, &gmat[6][0], 1, &diff0[0], 1,
164 Vmath::Vvtvp(nquad0 * nquad1 * nquad2, &gmat[7][0], 1, &diff1[0], 1,
165 &out_d2[0], 1, &out_d2[0], 1);
166 Vmath::Vvtvp(nquad0 * nquad1 * nquad2, &gmat[8][0], 1, &diff2[0], 1,
167 &out_d2[0], 1, &out_d2[0], 1);
174 Vmath::Smul(nquad0 * nquad1 * nquad2, gmat[0][0], &diff0[0], 1,
176 Blas::Daxpy(nquad0 * nquad1 * nquad2, gmat[1][0], &diff1[0], 1,
178 Blas::Daxpy(nquad0 * nquad1 * nquad2, gmat[2][0], &diff2[0], 1,
184 Vmath::Smul(nquad0 * nquad1 * nquad2, gmat[3][0], &diff0[0], 1,
186 Blas::Daxpy(nquad0 * nquad1 * nquad2, gmat[4][0], &diff1[0], 1,
188 Blas::Daxpy(nquad0 * nquad1 * nquad2, gmat[5][0], &diff2[0], 1,
194 Vmath::Smul(nquad0 * nquad1 * nquad2, gmat[6][0], &diff0[0], 1,
196 Blas::Daxpy(nquad0 * nquad1 * nquad2, gmat[7][0], &diff1[0], 1,
198 Blas::Daxpy(nquad0 * nquad1 * nquad2, gmat[8][0], &diff2[0], 1,
224 if (
m_base[0]->Collocation() &&
m_base[1]->Collocation() &&
241 out = (*matsys) * in;
286 const int nquad0 =
m_base[0]->GetNumPoints();
287 const int nquad1 =
m_base[1]->GetNumPoints();
288 const int nquad2 =
m_base[2]->GetNumPoints();
289 const int order0 =
m_base[0]->GetNumModes();
290 const int order1 =
m_base[1]->GetNumModes();
294 if (multiplybyweights)
302 tmp, outarray, wsp,
true,
true,
true);
308 inarray, outarray, wsp,
true,
true,
true);
353 const int nquad0 =
m_base[0]->GetNumPoints();
354 const int nquad1 =
m_base[1]->GetNumPoints();
355 const int nquad2 =
m_base[2]->GetNumPoints();
356 const int order0 =
m_base[0]->GetNumModes();
357 const int order1 =
m_base[1]->GetNumModes();
358 const int nqtot = nquad0 * nquad1 * nquad2;
366 std::max(nqtot, order0 * nquad2 * (nquad1 + order1)));
378 m_base[2]->GetBdata(), tmp2, outarray, wsp,
382 m_base[2]->GetBdata(), tmp3, tmp6, wsp,
true,
388 m_base[2]->GetDbdata(), tmp4, tmp6, wsp,
true,
398 const int nquad0 =
m_base[0]->GetNumPoints();
399 const int nquad1 =
m_base[1]->GetNumPoints();
400 const int nquad2 =
m_base[2]->GetNumPoints();
401 const int order0 =
m_base[0]->GetNumModes();
402 const int order1 =
m_base[1]->GetNumModes();
403 const int nqtot = nquad0 * nquad1 * nquad2;
414 std::max(nqtot, order0 * nquad2 * (nquad1 + order1)));
428 Vmath::Vmul(nqtot, &df[3 * dir][0], 1, tmp1.data(), 1, tmp2.data(), 1);
429 Vmath::Vmul(nqtot, &df[3 * dir + 1][0], 1, tmp1.data(), 1, tmp3.data(),
431 Vmath::Vmul(nqtot, &df[3 * dir + 2][0], 1, tmp1.data(), 1, tmp4.data(),
436 Vmath::Smul(nqtot, df[3 * dir][0], tmp1.data(), 1, tmp2.data(), 1);
437 Vmath::Smul(nqtot, df[3 * dir + 1][0], tmp1.data(), 1, tmp3.data(), 1);
438 Vmath::Smul(nqtot, df[3 * dir + 2][0], tmp1.data(), 1, tmp4.data(), 1);
442 for (
int i = 0; i < nquad0; ++i)
444 gfac0[i] = 0.5 * (1 + z0[i]);
448 for (
int i = 0; i < nquad1; ++i)
450 gfac1[i] = 0.5 * (1 + z1[i]);
454 for (
int i = 0; i < nquad2; ++i)
456 gfac2[i] = 2.0 / (1 - z2[i]);
459 const int nq01 = nquad0 * nquad1;
461 for (
int i = 0; i < nquad2; ++i)
463 Vmath::Smul(nq01, gfac2[i], &tmp2[0] + i * nq01, 1, &tmp2[0] + i * nq01,
465 Vmath::Smul(nq01, gfac2[i], &tmp3[0] + i * nq01, 1, &tmp3[0] + i * nq01,
467 Vmath::Smul(nq01, gfac2[i], &tmp4[0] + i * nq01, 1, &tmp5[0] + i * nq01,
472 for (
int i = 0; i < nquad1 * nquad2; ++i)
474 Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp5[0] + i * nquad0, 1,
475 &wsp[0] + i * nquad0, 1);
478 Vmath::Vadd(nqtot, &tmp2[0], 1, &wsp[0], 1, &tmp2[0], 1);
481 for (
int i = 0; i < nquad1 * nquad2; ++i)
483 Vmath::Smul(nquad0, gfac1[i % nquad1], &tmp5[0] + i * nquad0, 1,
484 &tmp5[0] + i * nquad0, 1);
486 Vmath::Vadd(nqtot, &tmp3[0], 1, &tmp5[0], 1, &tmp3[0], 1);
497 m_base[2]->GetBasisKey());
503 m_base[0]->GetPointsKey());
505 m_base[1]->GetPointsKey());
507 m_base[2]->GetPointsKey());
522 ASSERTL1(Lcoords[0] <= -1.0 && Lcoords[0] >= 1.0 && Lcoords[1] <= -1.0 &&
523 Lcoords[1] >= 1.0 && Lcoords[2] <= -1.0 && Lcoords[2] >= 1.0,
524 "Local coordinates are not in region [-1,1]");
528 for (i = 0; i <
m_geom->GetCoordim(); ++i)
530 coords[i] =
m_geom->GetCoord(i, Lcoords);
542 const NekDouble *data,
const std::vector<unsigned int> &nummodes,
543 const int mode_offset,
NekDouble *coeffs,
544 std::vector<LibUtilities::BasisType> &fromType)
546 int data_order0 = nummodes[mode_offset];
547 int fillorder0 = min(
m_base[0]->GetNumModes(), data_order0);
548 int data_order1 = nummodes[mode_offset + 1];
549 int order1 =
m_base[1]->GetNumModes();
550 int fillorder1 = min(order1, data_order1);
551 int data_order2 = nummodes[mode_offset + 2];
552 int order2 =
m_base[2]->GetNumModes();
553 int fillorder2 = min(order2, data_order2);
560 data_order1 != fillorder1 || data_order2 != fillorder2)
567 m_base[0]->GetPointsKey()),
569 m_base[1]->GetPointsKey()),
571 m_base[2]->GetPointsKey()));
575 m_base[2]->GetBasisKey());
601 return StdExpansion3D::v_PhysEvaluate(Lcoord, physvals);
612 m_geom->GetLocCoords(coord, Lcoord);
614 return StdExpansion3D::v_PhysEvaluate(Lcoord, physvals);
620 std::array<NekDouble, 3> &firstOrderDerivs)
624 m_geom->GetLocCoords(coord, Lcoord);
625 return StdPyrExp::v_PhysEvalFirstDeriv(Lcoord, inarray, firstOrderDerivs);
634 int nquad0 =
m_base[0]->GetNumPoints();
635 int nquad1 =
m_base[1]->GetNumPoints();
636 int nquad2 =
m_base[2]->GetNumPoints();
646 if (outarray.size() != nq0 * nq1)
652 for (
int i = 0; i < nquad0 * nquad1; ++i)
661 if (outarray.size() != nq0 * nq1)
667 for (
int k = 0; k < nquad2; k++)
669 for (
int i = 0; i < nquad0; ++i)
671 outarray[k * nquad0 + i] = (nquad0 * nquad1 * k) + i;
679 if (outarray.size() != nq0 * nq1)
685 for (
int j = 0; j < nquad1 * nquad2; ++j)
687 outarray[j] = nquad0 - 1 + j * nquad0;
694 if (outarray.size() != nq0 * nq1)
700 for (
int k = 0; k < nquad2; k++)
702 for (
int i = 0; i < nquad0; ++i)
704 outarray[k * nquad0 + i] =
705 nquad0 * (nquad1 - 1) + (nquad0 * nquad1 * k) + i;
713 if (outarray.size() != nq0 * nq1)
719 for (
int j = 0; j < nquad1 * nquad2; ++j)
721 outarray[j] = j * nquad0;
725 ASSERTL0(
false,
"face value (> 4) is out of range");
736 for (
int i = 0; i < ptsKeys.size(); ++i)
748 geomFactors->GetDerivFactors(ptsKeys);
762 for (i = 0; i < vCoordDim; ++i)
767 size_t nqb = nq_face;
782 for (i = 0; i < vCoordDim; ++i)
784 normal[i][0] = -df[3 * i + 2][0];
790 for (i = 0; i < vCoordDim; ++i)
792 normal[i][0] = -df[3 * i + 1][0];
798 for (i = 0; i < vCoordDim; ++i)
800 normal[i][0] = df[3 * i][0] + df[3 * i + 2][0];
806 for (i = 0; i < vCoordDim; ++i)
808 normal[i][0] = df[3 * i + 1][0] + df[3 * i + 2][0];
814 for (i = 0; i < vCoordDim; ++i)
816 normal[i][0] = -df[3 * i][0];
821 ASSERTL0(
false,
"face is out of range (face < 4)");
826 for (i = 0; i < vCoordDim; ++i)
828 fac += normal[i][0] * normal[i][0];
830 fac = 1.0 /
sqrt(fac);
834 for (i = 0; i < vCoordDim; ++i)
836 Vmath::Fill(nq_face, fac * normal[i][0], normal[i], 1);
844 int nq0 = ptsKeys[0].GetNumPoints();
845 int nq1 = ptsKeys[1].GetNumPoints();
846 int nq2 = ptsKeys[2].GetNumPoints();
847 int nq01 = nq0 * nq1;
855 else if (face == 1 || face == 3)
877 for (j = 0; j < nq01; ++j)
879 normals[j] = -df[2][j] * jac[j];
880 normals[nqtot + j] = -df[5][j] * jac[j];
881 normals[2 * nqtot + j] = -df[8][j] * jac[j];
885 points0 = ptsKeys[0];
886 points1 = ptsKeys[1];
892 for (j = 0; j < nq0; ++j)
894 for (k = 0; k < nq2; ++k)
896 int tmp = j + nq01 * k;
897 normals[j + k * nq0] = -df[1][tmp] * jac[tmp];
898 normals[nqtot + j + k * nq0] = -df[4][tmp] * jac[tmp];
899 normals[2 * nqtot + j + k * nq0] =
900 -df[7][tmp] * jac[tmp];
901 faceJac[j + k * nq0] = jac[tmp];
905 points0 = ptsKeys[0];
906 points1 = ptsKeys[2];
912 for (j = 0; j < nq1; ++j)
914 for (k = 0; k < nq2; ++k)
916 int tmp = nq0 - 1 + nq0 * j + nq01 * k;
917 normals[j + k * nq1] =
918 (df[0][tmp] + df[2][tmp]) * jac[tmp];
919 normals[nqtot + j + k * nq1] =
920 (df[3][tmp] + df[5][tmp]) * jac[tmp];
921 normals[2 * nqtot + j + k * nq1] =
922 (df[6][tmp] + df[8][tmp]) * jac[tmp];
923 faceJac[j + k * nq1] = jac[tmp];
927 points0 = ptsKeys[1];
928 points1 = ptsKeys[2];
934 for (j = 0; j < nq0; ++j)
936 for (k = 0; k < nq2; ++k)
938 int tmp = nq0 * (nq1 - 1) + j + nq01 * k;
939 normals[j + k * nq0] =
940 (df[1][tmp] + df[2][tmp]) * jac[tmp];
941 normals[nqtot + j + k * nq0] =
942 (df[4][tmp] + df[5][tmp]) * jac[tmp];
943 normals[2 * nqtot + j + k * nq0] =
944 (df[7][tmp] + df[8][tmp]) * jac[tmp];
945 faceJac[j + k * nq0] = jac[tmp];
949 points0 = ptsKeys[0];
950 points1 = ptsKeys[2];
956 for (j = 0; j < nq1; ++j)
958 for (k = 0; k < nq2; ++k)
960 int tmp = j * nq0 + nq01 * k;
961 normals[j + k * nq1] = -df[0][tmp] * jac[tmp];
962 normals[nqtot + j + k * nq1] = -df[3][tmp] * jac[tmp];
963 normals[2 * nqtot + j + k * nq1] =
964 -df[6][tmp] * jac[tmp];
965 faceJac[j + k * nq1] = jac[tmp];
969 points0 = ptsKeys[1];
970 points1 = ptsKeys[2];
975 ASSERTL0(
false,
"face is out of range (face < 4)");
983 Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
986 for (i = 0; i < vCoordDim; ++i)
991 Vmath::Vmul(nq_face, work, 1, normal[i], 1, normal[i], 1);
998 Vmath::Vvtvp(nq_face, normal[i], 1, normal[i], 1, work, 1, work, 1);
1008 Vmath::Vmul(nq_face, normal[i], 1, work, 1, normal[i], 1);
1034 StdPyrExp::v_SVVLaplacianFilter(array, mkey);
1059 returnval = StdPyrExp::v_GenMatrix(mkey);
1073 return tmp->GetStdMatrix(mkey);
1105 const unsigned int dim = 3;
1111 for (
unsigned int i = 0; i < dim; ++i)
1113 for (
unsigned int j = i; j < dim; ++j)
1144 const unsigned int nquad0 =
m_base[0]->GetNumPoints();
1145 const unsigned int nquad1 =
m_base[1]->GetNumPoints();
1146 const unsigned int nquad2 =
m_base[2]->GetNumPoints();
1149 for (j = 0; j < nquad2; ++j)
1151 for (i = 0; i < nquad1; ++i)
1154 &h0[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1156 &h1[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1157 Vmath::Fill(nquad0, (1.0 + z1[i]) / (1.0 - z2[j]),
1158 &h2[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1161 for (i = 0; i < nquad0; i++)
1163 Blas::Dscal(nquad1 * nquad2, 1 + z0[i], &h1[0] + i, nquad0);
1172 Vmath::Vvtvvtp(nqtot, &df[0][0], 1, &h0[0], 1, &df[2][0], 1, &h1[0], 1,
1174 Vmath::Vvtvvtp(nqtot, &df[3][0], 1, &h0[0], 1, &df[5][0], 1, &h1[0], 1,
1176 Vmath::Vvtvvtp(nqtot, &df[6][0], 1, &h0[0], 1, &df[8][0], 1, &h1[0], 1,
1180 Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp1[0], 1, &wsp2[0], 1, &wsp2[0],
1182 Vmath::Vvtvp(nqtot, &wsp3[0], 1, &wsp3[0], 1, &g0[0], 1, &g0[0], 1);
1185 Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp1[0], 1, &df[5][0], 1, &wsp2[0],
1187 Vmath::Vvtvp(nqtot, &df[8][0], 1, &wsp3[0], 1, &g4[0], 1, &g4[0], 1);
1190 Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &h0[0], 1, &df[2][0], 1, &h2[0], 1,
1192 Vmath::Vvtvvtp(nqtot, &df[4][0], 1, &h0[0], 1, &df[5][0], 1, &h2[0], 1,
1194 Vmath::Vvtvvtp(nqtot, &df[7][0], 1, &h0[0], 1, &df[8][0], 1, &h2[0], 1,
1198 Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0],
1200 Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp6[0], 1, &g1[0], 1, &g1[0], 1);
1203 Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp4[0], 1, &wsp2[0], 1, &wsp5[0],
1205 Vmath::Vvtvp(nqtot, &wsp3[0], 1, &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1208 Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp4[0], 1, &df[5][0], 1, &wsp5[0],
1210 Vmath::Vvtvp(nqtot, &df[8][0], 1, &wsp6[0], 1, &g5[0], 1, &g5[0], 1);
1214 &df[5][0], 1, &g2[0], 1);
1215 Vmath::Vvtvp(nqtot, &df[8][0], 1, &df[8][0], 1, &g2[0], 1, &g2[0], 1);
1228 Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp1[0], 1, &wsp2[0], 1, &wsp2[0],
1230 Vmath::Vvtvp(nqtot, &wsp3[0], 1, &wsp3[0], 1, &g0[0], 1, &g0[0], 1);
1233 Vmath::Svtsvtp(nqtot, df[2][0], &wsp1[0], 1, df[5][0], &wsp2[0], 1,
1235 Vmath::Svtvp(nqtot, df[8][0], &wsp3[0], 1, &g4[0], 1, &g4[0], 1);
1246 Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0],
1248 Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp6[0], 1, &g1[0], 1, &g1[0], 1);
1251 Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp4[0], 1, &wsp2[0], 1, &wsp5[0],
1253 Vmath::Vvtvp(nqtot, &wsp3[0], 1, &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1256 Vmath::Svtsvtp(nqtot, df[2][0], &wsp4[0], 1, df[5][0], &wsp5[0], 1,
1258 Vmath::Svtvp(nqtot, df[8][0], &wsp6[0], 1, &g5[0], 1, &g5[0], 1);
1262 df[2][0] * df[2][0] + df[5][0] * df[5][0] +
1263 df[8][0] * df[8][0],
1267 for (
unsigned int i = 0; i < dim; ++i)
1269 for (
unsigned int j = i; j < dim; ++j)
1287 int nquad0 =
m_base[0]->GetNumPoints();
1288 int nquad1 =
m_base[1]->GetNumPoints();
1289 int nq2 =
m_base[2]->GetNumPoints();
1290 int nqtot = nquad0 * nquad1 * nq2;
1292 ASSERTL1(wsp.size() >= 6 * nqtot,
"Insufficient workspace size.");
1293 ASSERTL1(m_ncoeffs <= nqtot, "Workspace not set up for ncoeffs > nqtot
");
1295 const Array<OneD, const NekDouble> &base0 = m_base[0]->GetBdata();
1296 const Array<OneD, const NekDouble> &base1 = m_base[1]->GetBdata();
1297 const Array<OneD, const NekDouble> &base2 = m_base[2]->GetBdata();
1298 const Array<OneD, const NekDouble> &dbase0 = m_base[0]->GetDbdata();
1299 const Array<OneD, const NekDouble> &dbase1 = m_base[1]->GetDbdata();
1300 const Array<OneD, const NekDouble> &dbase2 = m_base[2]->GetDbdata();
1301 const Array<OneD, const NekDouble> &metric00 =
1302 m_metrics[eMetricLaplacian00];
1303 const Array<OneD, const NekDouble> &metric01 =
1304 m_metrics[eMetricLaplacian01];
1305 const Array<OneD, const NekDouble> &metric02 =
1306 m_metrics[eMetricLaplacian02];
1307 const Array<OneD, const NekDouble> &metric11 =
1308 m_metrics[eMetricLaplacian11];
1309 const Array<OneD, const NekDouble> &metric12 =
1310 m_metrics[eMetricLaplacian12];
1311 const Array<OneD, const NekDouble> &metric22 =
1312 m_metrics[eMetricLaplacian22];
1314 // Allocate temporary storage
1315 Array<OneD, NekDouble> wsp0(2 * nqtot, wsp);
1316 Array<OneD, NekDouble> wsp1(nqtot, wsp + 1 * nqtot);
1317 Array<OneD, NekDouble> wsp2(nqtot, wsp + 2 * nqtot);
1318 Array<OneD, NekDouble> wsp3(nqtot, wsp + 3 * nqtot);
1319 Array<OneD, NekDouble> wsp4(nqtot, wsp + 4 * nqtot);
1320 Array<OneD, NekDouble> wsp5(nqtot, wsp + 5 * nqtot);
1322 // LAPLACIAN MATRIX OPERATION
1323 // wsp1 = du_dxi1 = D_xi1 * inarray = D_xi1 * u
1324 // wsp2 = du_dxi2 = D_xi2 * inarray = D_xi2 * u
1325 // wsp2 = du_dxi3 = D_xi3 * inarray = D_xi3 * u
1326 StdExpansion3D::PhysTensorDeriv(inarray, wsp0, wsp1, wsp2);
1328 // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1329 // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1330 // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
1331 // especially for this purpose
1332 Vmath::Vvtvvtp(nqtot, &metric00[0], 1, &wsp0[0], 1, &metric01[0], 1,
1333 &wsp1[0], 1, &wsp3[0], 1);
1334 Vmath::Vvtvp(nqtot, &metric02[0], 1, &wsp2[0], 1, &wsp3[0], 1, &wsp3[0], 1);
1335 Vmath::Vvtvvtp(nqtot, &metric01[0], 1, &wsp0[0], 1, &metric11[0], 1,
1336 &wsp1[0], 1, &wsp4[0], 1);
1337 Vmath::Vvtvp(nqtot, &metric12[0], 1, &wsp2[0], 1, &wsp4[0], 1, &wsp4[0], 1);
1338 Vmath::Vvtvvtp(nqtot, &metric02[0], 1, &wsp0[0], 1, &metric12[0], 1,
1339 &wsp1[0], 1, &wsp5[0], 1);
1340 Vmath::Vvtvp(nqtot, &metric22[0], 1, &wsp2[0], 1, &wsp5[0], 1, &wsp5[0], 1);
1342 // outarray = m = (D_xi1 * B)^T * k
1343 // wsp1 = n = (D_xi2 * B)^T * l
1344 IProductWRTBase_SumFacKernel(dbase0, base1, base2, wsp3, outarray, wsp0,
1346 IProductWRTBase_SumFacKernel(base0, dbase1, base2, wsp4, wsp2, wsp0, true,
1348 Vmath::Vadd(m_ncoeffs, wsp2.data(), 1, outarray.data(), 1, outarray.data(),
1350 IProductWRTBase_SumFacKernel(base0, base1, dbase2, wsp5, wsp2, wsp0, true,
1352 Vmath::Vadd(m_ncoeffs, wsp2.data(), 1, outarray.data(), 1, outarray.data(),
1361void PyrExp::v_NormalTraceDerivFactors(
1362 Array<OneD, Array<OneD, NekDouble>> &d0factors,
1363 Array<OneD, Array<OneD, NekDouble>> &d1factors,
1364 Array<OneD, Array<OneD, NekDouble>> &d2factors)
1366 int nquad0 = GetNumPoints(0);
1367 int nquad1 = GetNumPoints(1);
1368 int nquad2 = GetNumPoints(2);
1370 const Array<TwoD, const NekDouble> &df =
1371 m_metricinfo->GetDerivFactors(GetPointsKeys());
1373 if (d0factors.size() != 5)
1375 d0factors = Array<OneD, Array<OneD, NekDouble>>(5);
1376 d1factors = Array<OneD, Array<OneD, NekDouble>>(5);
1377 d2factors = Array<OneD, Array<OneD, NekDouble>>(5);
1380 if (d0factors[0].size() != nquad0 * nquad1)
1382 d0factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1383 d1factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1384 d2factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1387 if (d0factors[1].size() != nquad0 * nquad2)
1389 d0factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1390 d0factors[3] = Array<OneD, NekDouble>(nquad0 * nquad2);
1391 d1factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1392 d1factors[3] = Array<OneD, NekDouble>(nquad0 * nquad2);
1393 d2factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1394 d2factors[3] = Array<OneD, NekDouble>(nquad0 * nquad2);
1397 if (d0factors[2].size() != nquad1 * nquad2)
1399 d0factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1400 d0factors[4] = Array<OneD, NekDouble>(nquad1 * nquad2);
1401 d1factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1402 d1factors[4] = Array<OneD, NekDouble>(nquad1 * nquad2);
1403 d2factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1404 d2factors[4] = Array<OneD, NekDouble>(nquad1 * nquad2);
1408 const Array<OneD, const Array<OneD, NekDouble>> &normal_0 =
1410 const Array<OneD, const Array<OneD, NekDouble>> &normal_1 =
1412 const Array<OneD, const Array<OneD, NekDouble>> &normal_2 =
1414 const Array<OneD, const Array<OneD, NekDouble>> &normal_3 =
1416 const Array<OneD, const Array<OneD, NekDouble>> &normal_4 =
1419 int ncoords = normal_0.size();
1421 // first gather together standard cartesian inner products
1422 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1425 for (int i = 0; i < nquad0 * nquad1; ++i)
1427 d0factors[0][i] = df[0][i] * normal_0[0][i];
1428 d1factors[0][i] = df[1][i] * normal_0[0][i];
1429 d2factors[0][i] = df[2][i] * normal_0[0][i];
1432 for (int n = 1; n < ncoords; ++n)
1434 for (int i = 0; i < nquad0 * nquad1; ++i)
1436 d0factors[0][i] += df[3 * n][i] * normal_0[n][i];
1437 d1factors[0][i] += df[3 * n + 1][i] * normal_0[n][i];
1438 d2factors[0][i] += df[3 * n + 2][i] * normal_0[n][i];
1443 for (int j = 0; j < nquad2; ++j)
1445 for (int i = 0; i < nquad0; ++i)
1447 d0factors[1][j * nquad0 + i] = df[0][j * nquad0 * nquad1 + i] *
1448 normal_1[0][j * nquad0 + i];
1449 d1factors[1][j * nquad0 + i] = df[1][j * nquad0 * nquad1 + i] *
1450 normal_1[0][j * nquad0 + i];
1451 d2factors[1][j * nquad0 + i] = df[2][j * nquad0 * nquad1 + i] *
1452 normal_1[0][j * nquad0 + i];
1454 d0factors[3][j * nquad0 + i] =
1455 df[0][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1456 normal_3[0][j * nquad0 + i];
1457 d1factors[3][j * nquad0 + i] =
1458 df[1][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1459 normal_3[0][j * nquad0 + i];
1460 d2factors[3][j * nquad0 + i] =
1461 df[2][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1462 normal_3[0][j * nquad0 + i];
1466 for (int n = 1; n < ncoords; ++n)
1468 for (int j = 0; j < nquad2; ++j)
1470 for (int i = 0; i < nquad0; ++i)
1472 d0factors[1][j * nquad0 + i] +=
1473 df[3 * n][j * nquad0 * nquad1 + i] *
1474 normal_1[0][j * nquad0 + i];
1475 d1factors[1][j * nquad0 + i] +=
1476 df[3 * n + 1][j * nquad0 * nquad1 + i] *
1477 normal_1[0][j * nquad0 + i];
1478 d2factors[1][j * nquad0 + i] +=
1479 df[3 * n + 2][j * nquad0 * nquad1 + i] *
1480 normal_1[0][j * nquad0 + i];
1482 d0factors[3][j * nquad0 + i] +=
1483 df[3 * n][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1484 normal_3[0][j * nquad0 + i];
1485 d1factors[3][j * nquad0 + i] +=
1486 df[3 * n + 1][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1487 normal_3[0][j * nquad0 + i];
1488 d2factors[3][j * nquad0 + i] +=
1489 df[3 * n + 2][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1490 normal_3[0][j * nquad0 + i];
1496 for (int j = 0; j < nquad2; ++j)
1498 for (int i = 0; i < nquad1; ++i)
1500 d0factors[2][j * nquad1 + i] =
1501 df[0][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1502 normal_2[0][j * nquad1 + i];
1503 d1factors[2][j * nquad1 + i] =
1504 df[1][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1505 normal_2[0][j * nquad1 + i];
1506 d2factors[2][j * nquad1 + i] =
1507 df[2][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1508 normal_2[0][j * nquad1 + i];
1510 d0factors[4][j * nquad1 + i] =
1511 df[0][j * nquad0 * nquad1 + i * nquad0] *
1512 normal_4[0][j * nquad1 + i];
1513 d1factors[4][j * nquad1 + i] =
1514 df[1][j * nquad0 * nquad1 + i * nquad0] *
1515 normal_4[0][j * nquad1 + i];
1516 d2factors[4][j * nquad1 + i] =
1517 df[2][j * nquad0 * nquad1 + i * nquad0] *
1518 normal_4[0][j * nquad1 + i];
1522 for (int n = 1; n < ncoords; ++n)
1524 for (int j = 0; j < nquad2; ++j)
1526 for (int i = 0; i < nquad1; ++i)
1528 d0factors[2][j * nquad1 + i] +=
1529 df[3 * n][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1530 normal_2[n][j * nquad1 + i];
1531 d1factors[2][j * nquad1 + i] +=
1533 [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1534 normal_2[n][j * nquad1 + i];
1535 d2factors[2][j * nquad1 + i] +=
1537 [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1538 normal_2[n][j * nquad1 + i];
1540 d0factors[4][j * nquad1 + i] +=
1541 df[3 * n][i * nquad0 + j * nquad0 * nquad1] *
1542 normal_4[n][j * nquad1 + i];
1543 d1factors[4][j * nquad1 + i] +=
1544 df[3 * n + 1][i * nquad0 + j * nquad0 * nquad1] *
1545 normal_4[n][j * nquad1 + i];
1546 d2factors[4][j * nquad1 + i] +=
1547 df[3 * n + 2][i * nquad0 + j * nquad0 * nquad1] *
1548 normal_4[n][j * nquad1 + i];
1556 for (int i = 0; i < nquad0 * nquad1; ++i)
1558 d0factors[0][i] = df[0][0] * normal_0[0][i];
1559 d1factors[0][i] = df[1][0] * normal_0[0][i];
1560 d2factors[0][i] = df[2][0] * normal_0[0][i];
1563 for (int n = 1; n < ncoords; ++n)
1565 for (int i = 0; i < nquad0 * nquad1; ++i)
1567 d0factors[0][i] += df[3 * n][0] * normal_0[n][i];
1568 d1factors[0][i] += df[3 * n + 1][0] * normal_0[n][i];
1569 d2factors[0][i] += df[3 * n + 2][0] * normal_0[n][i];
1574 for (int i = 0; i < nquad0 * nquad2; ++i)
1576 d0factors[1][i] = df[0][0] * normal_1[0][i];
1577 d0factors[3][i] = df[0][0] * normal_3[0][i];
1579 d1factors[1][i] = df[1][0] * normal_1[0][i];
1580 d1factors[3][i] = df[1][0] * normal_3[0][i];
1582 d2factors[1][i] = df[2][0] * normal_1[0][i];
1583 d2factors[3][i] = df[2][0] * normal_3[0][i];
1586 for (int n = 1; n < ncoords; ++n)
1588 for (int i = 0; i < nquad0 * nquad2; ++i)
1590 d0factors[1][i] += df[3 * n][0] * normal_1[n][i];
1591 d0factors[3][i] += df[3 * n][0] * normal_3[n][i];
1593 d1factors[1][i] += df[3 * n + 1][0] * normal_1[n][i];
1594 d1factors[3][i] += df[3 * n + 1][0] * normal_3[n][i];
1596 d2factors[1][i] += df[3 * n + 2][0] * normal_1[n][i];
1597 d2factors[3][i] += df[3 * n + 2][0] * normal_3[n][i];
1602 for (int i = 0; i < nquad1 * nquad2; ++i)
1604 d0factors[2][i] = df[0][0] * normal_2[0][i];
1605 d0factors[4][i] = df[0][0] * normal_4[0][i];
1607 d1factors[2][i] = df[1][0] * normal_2[0][i];
1608 d1factors[4][i] = df[1][0] * normal_4[0][i];
1610 d2factors[2][i] = df[2][0] * normal_2[0][i];
1611 d2factors[4][i] = df[2][0] * normal_4[0][i];
1614 for (int n = 1; n < ncoords; ++n)
1616 for (int i = 0; i < nquad1 * nquad2; ++i)
1618 d0factors[2][i] += df[3 * n][0] * normal_2[n][i];
1619 d0factors[4][i] += df[3 * n][0] * normal_4[n][i];
1621 d1factors[2][i] += df[3 * n + 1][0] * normal_2[n][i];
1622 d1factors[4][i] += df[3 * n + 1][0] * normal_4[n][i];
1624 d2factors[2][i] += df[3 * n + 2][0] * normal_2[n][i];
1625 d2factors[4][i] += df[3 * n + 2][0] * normal_4[n][i];
1631} // namespace Nektar::LocalRegions
#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 v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
void ComputeQuadratureMetric()
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
void v_ExtractDataToCoeffs(const NekDouble *data, const std::vector< unsigned int > &nummodes, const int mode_offset, NekDouble *coeffs, std::vector< LibUtilities::BasisType > &fromType) override
NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray) override
Integrate the physical point list inarray over pyramidic region and return the value.
void v_ComputeLaplacianMetric() override
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey) override
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_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey) override
NekDouble v_PhysEvalFirstDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) override
void v_ComputeTraceNormal(const int face) override
void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords) override
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Calculates the inner product .
StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const override
void v_DropLocMatrix(const MatrixKey &mkey) override
void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp) override
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
PyrExp(const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const LibUtilities::BasisKey &Bc, const SpatialDomains::PyrGeomSharedPtr &geom)
Constructor using BasisKey class for quadrature points and order definition.
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const override
DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey) override
void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray) override
void v_GetTracePhysMap(const int face, Array< OneD, int > &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
Calculate the derivative of the physical points.
NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coord, const Array< OneD, const NekDouble > &physvals) override
This function evaluates the expansion at a single (arbitrary) point of the domain.
void v_DropLocStaticCondMatrix(const MatrixKey &mkey) 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 basis B=base0*base1*base2 and put into out...
DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey) override
NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals) override
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true) 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)
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
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.
Array< OneD, LibUtilities::BasisSharedPtr > m_base
MatrixType GetMatrixType() const
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
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)
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
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< PyrGeom > PyrGeomSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
std::shared_ptr< StdPyrExp > StdPyrExpSharedPtr
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
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 Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
Svtsvtp (scalar times vector plus scalar times vector):
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 Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Svtvp (scalar times vector plus vector): z = alpha*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)