50 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
53 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
58 std::string(
"PrismExpMatrix")),
59 m_staticCondMatrixManager(
std::bind(&
Expansion::CreateStaticCondMatrix,
60 this,
std::placeholders::_1),
61 std::string(
"PrismExpStaticCondMatrix"))
66 : StdExpansion(T), StdExpansion3D(T), StdPrismExp(T),
Expansion(T),
68 m_staticCondMatrixManager(T.m_staticCondMatrixManager)
98 int nquad0 =
m_base[0]->GetNumPoints();
99 int nquad1 =
m_base[1]->GetNumPoints();
100 int nquad2 =
m_base[2]->GetNumPoints();
108 (
NekDouble *)&inarray[0], 1, &tmp[0], 1);
113 (
NekDouble *)&inarray[0], 1, &tmp[0], 1);
117 return StdPrismExp::v_Integral(tmp);
136 StdPrismExp::v_PhysDeriv(inarray, diff0, diff1, diff2);
142 Vmath::Vmul(nqtot, &df[0][0], 1, &diff0[0], 1, &out_d0[0], 1);
143 Vmath::Vvtvp(nqtot, &df[1][0], 1, &diff1[0], 1, &out_d0[0], 1,
145 Vmath::Vvtvp(nqtot, &df[2][0], 1, &diff2[0], 1, &out_d0[0], 1,
151 Vmath::Vmul(nqtot, &df[3][0], 1, &diff0[0], 1, &out_d1[0], 1);
152 Vmath::Vvtvp(nqtot, &df[4][0], 1, &diff1[0], 1, &out_d1[0], 1,
154 Vmath::Vvtvp(nqtot, &df[5][0], 1, &diff2[0], 1, &out_d1[0], 1,
160 Vmath::Vmul(nqtot, &df[6][0], 1, &diff0[0], 1, &out_d2[0], 1);
161 Vmath::Vvtvp(nqtot, &df[7][0], 1, &diff1[0], 1, &out_d2[0], 1,
163 Vmath::Vvtvp(nqtot, &df[8][0], 1, &diff2[0], 1, &out_d2[0], 1,
171 Vmath::Smul(nqtot, df[0][0], &diff0[0], 1, &out_d0[0], 1);
172 Blas::Daxpy(nqtot, df[1][0], &diff1[0], 1, &out_d0[0], 1);
173 Blas::Daxpy(nqtot, df[2][0], &diff2[0], 1, &out_d0[0], 1);
178 Vmath::Smul(nqtot, df[3][0], &diff0[0], 1, &out_d1[0], 1);
179 Blas::Daxpy(nqtot, df[4][0], &diff1[0], 1, &out_d1[0], 1);
180 Blas::Daxpy(nqtot, df[5][0], &diff2[0], 1, &out_d1[0], 1);
185 Vmath::Smul(nqtot, df[6][0], &diff0[0], 1, &out_d2[0], 1);
186 Blas::Daxpy(nqtot, df[7][0], &diff1[0], 1, &out_d2[0], 1);
187 Blas::Daxpy(nqtot, df[8][0], &diff2[0], 1, &out_d2[0], 1);
212 if (
m_base[0]->Collocation() &&
m_base[1]->Collocation() &&
229 out = (*matsys) * in;
271 const int nquad0 =
m_base[0]->GetNumPoints();
272 const int nquad1 =
m_base[1]->GetNumPoints();
273 const int nquad2 =
m_base[2]->GetNumPoints();
274 const int order0 =
m_base[0]->GetNumModes();
275 const int order1 =
m_base[1]->GetNumModes();
279 if (multiplybyweights)
287 tmp, outarray, wsp,
true,
true,
true);
293 inarray, outarray, wsp,
true,
true,
true);
338 const int nquad0 =
m_base[0]->GetNumPoints();
339 const int nquad1 =
m_base[1]->GetNumPoints();
340 const int nquad2 =
m_base[2]->GetNumPoints();
341 const int order0 =
m_base[0]->GetNumModes();
342 const int order1 =
m_base[1]->GetNumModes();
343 const int nqtot = nquad0 * nquad1 * nquad2;
362 m_base[2]->GetBdata(), tmp2, outarray, wsp,
366 m_base[2]->GetBdata(), tmp3, tmp6, wsp,
true,
372 m_base[2]->GetDbdata(), tmp4, tmp6, wsp,
true,
382 const int nquad0 =
m_base[0]->GetNumPoints();
383 const int nquad1 =
m_base[1]->GetNumPoints();
384 const int nquad2 =
m_base[2]->GetNumPoints();
385 const int order0 =
m_base[0]->GetNumModes();
386 const int order1 =
m_base[1]->GetNumModes();
387 const int nqtot = nquad0 * nquad1 * nquad2;
410 Vmath::Vmul(nqtot, &df[3 * dir][0], 1, tmp1.data(), 1, tmp2.data(), 1);
411 Vmath::Vmul(nqtot, &df[3 * dir + 1][0], 1, tmp1.data(), 1, tmp3.data(),
413 Vmath::Vmul(nqtot, &df[3 * dir + 2][0], 1, tmp1.data(), 1, tmp4.data(),
418 Vmath::Smul(nqtot, df[3 * dir][0], tmp1.data(), 1, tmp2.data(), 1);
419 Vmath::Smul(nqtot, df[3 * dir + 1][0], tmp1.data(), 1, tmp3.data(), 1);
420 Vmath::Smul(nqtot, df[3 * dir + 2][0], tmp1.data(), 1, tmp4.data(), 1);
424 for (
int i = 0; i < nquad0; ++i)
426 gfac0[i] = 0.5 * (1 + z0[i]);
430 for (
int i = 0; i < nquad2; ++i)
432 gfac2[i] = 2.0 / (1 - z2[i]);
435 const int nq01 = nquad0 * nquad1;
437 for (
int i = 0; i < nquad2; ++i)
439 Vmath::Smul(nq01, gfac2[i], &tmp2[0] + i * nq01, 1, &tmp2[0] + i * nq01,
441 Vmath::Smul(nq01, gfac2[i], &tmp4[0] + i * nq01, 1, &tmp5[0] + i * nq01,
445 for (
int i = 0; i < nquad1 * nquad2; ++i)
447 Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp5[0] + i * nquad0, 1,
448 &tmp5[0] + i * nquad0, 1);
451 Vmath::Vadd(nqtot, &tmp2[0], 1, &tmp5[0], 1, &tmp2[0], 1);
462 m_base[2]->GetBasisKey());
468 m_base[0]->GetPointsKey());
470 m_base[1]->GetPointsKey());
472 m_base[2]->GetPointsKey());
475 bkey0, bkey1, bkey2);
487 ASSERTL1(Lcoords[0] <= -1.0 && Lcoords[0] >= 1.0 && Lcoords[1] <= -1.0 &&
488 Lcoords[1] >= 1.0 && Lcoords[2] <= -1.0 && Lcoords[2] >= 1.0,
489 "Local coordinates are not in region [-1,1]");
493 for (i = 0; i <
m_geom->GetCoordim(); ++i)
495 coords[i] =
m_geom->GetCoord(i, Lcoords);
516 return StdExpansion3D::v_PhysEvaluate(Lcoord, physvals);
526 m_geom->GetLocCoords(coord, Lcoord);
528 return StdExpansion3D::v_PhysEvaluate(Lcoord, physvals);
534 std::array<NekDouble, 3> &firstOrderDerivs)
538 m_geom->GetLocCoords(coord, Lcoord);
539 return StdPrismExp::v_PhysEvalFirstDeriv(Lcoord, inarray, firstOrderDerivs);
547 const NekDouble *data,
const std::vector<unsigned int> &nummodes,
548 const int mode_offset,
NekDouble *coeffs,
549 [[maybe_unused]] std::vector<LibUtilities::BasisType> &fromType)
551 int data_order0 = nummodes[mode_offset];
552 int fillorder0 = min(
m_base[0]->GetNumModes(), data_order0);
553 int data_order1 = nummodes[mode_offset + 1];
554 int order1 =
m_base[1]->GetNumModes();
555 int fillorder1 = min(order1, data_order1);
556 int data_order2 = nummodes[mode_offset + 2];
557 int order2 =
m_base[2]->GetNumModes();
558 int fillorder2 = min(order2, data_order2);
569 "Extraction routine not set up for this basis");
571 "Extraction routine not set up for this basis");
574 for (j = 0; j < fillorder0; ++j)
576 for (i = 0; i < fillorder1; ++i)
578 Vmath::Vcopy(fillorder2 - j, &data[cnt], 1, &coeffs[cnt1],
580 cnt += data_order2 - j;
585 for (i = fillorder1; i < data_order1; ++i)
587 cnt += data_order2 - j;
590 for (i = fillorder1; i < order1; ++i)
598 ASSERTL0(
false,
"basis is either not set up or not "
605 int nquad0 =
m_base[0]->GetNumPoints();
606 int nquad1 =
m_base[1]->GetNumPoints();
607 int nquad2 =
m_base[2]->GetNumPoints();
616 if (outarray.size() != nq0 * nq1)
622 for (
int i = 0; i < nquad0 * nquad1; ++i)
631 if (outarray.size() != nq0 * nq1)
637 for (
int k = 0; k < nquad2; k++)
639 for (
int i = 0; i < nquad0; ++i)
641 outarray[k * nquad0 + i] = (nquad0 * nquad1 * k) + i;
650 if (outarray.size() != nq0 * nq1)
656 for (
int j = 0; j < nquad1 * nquad2; ++j)
658 outarray[j] = nquad0 - 1 + j * nquad0;
664 if (outarray.size() != nq0 * nq1)
670 for (
int k = 0; k < nquad2; k++)
672 for (
int i = 0; i < nquad0; ++i)
674 outarray[k * nquad0 + i] =
675 nquad0 * (nquad1 - 1) + (nquad0 * nquad1 * k) + i;
683 if (outarray.size() != nq0 * nq1)
689 for (
int j = 0; j < nquad1 * nquad2; ++j)
691 outarray[j] = j * nquad0;
695 ASSERTL0(
false,
"face value (> 4) is out of range");
710 for (
int i = 0; i < ptsKeys.size(); ++i)
722 geomFactors->GetDerivFactors(ptsKeys);
725 int nq0 = ptsKeys[0].GetNumPoints();
726 int nq1 = ptsKeys[1].GetNumPoints();
727 int nq2 = ptsKeys[2].GetNumPoints();
728 int nq01 = nq0 * nq1;
742 for (i = 0; i < vCoordDim; ++i)
747 size_t nqb = nq_face;
762 for (i = 0; i < vCoordDim; ++i)
764 normal[i][0] = -df[3 * i + 2][0];
771 for (i = 0; i < vCoordDim; ++i)
773 normal[i][0] = -df[3 * i + 1][0];
779 for (i = 0; i < vCoordDim; ++i)
781 normal[i][0] = df[3 * i][0] + df[3 * i + 2][0];
787 for (i = 0; i < vCoordDim; ++i)
789 normal[i][0] = df[3 * i + 1][0];
795 for (i = 0; i < vCoordDim; ++i)
797 normal[i][0] = -df[3 * i][0];
802 ASSERTL0(
false,
"face is out of range (face < 4)");
807 for (i = 0; i < vCoordDim; ++i)
809 fac += normal[i][0] * normal[i][0];
811 fac = 1.0 /
sqrt(fac);
815 for (i = 0; i < vCoordDim; ++i)
817 Vmath::Fill(nq_face, fac * normal[i][0], normal[i], 1);
830 else if (face == 1 || face == 3)
852 for (j = 0; j < nq01; ++j)
854 normals[j] = -df[2][j] * jac[j];
855 normals[nqtot + j] = -df[5][j] * jac[j];
856 normals[2 * nqtot + j] = -df[8][j] * jac[j];
860 points0 = ptsKeys[0];
861 points1 = ptsKeys[1];
867 for (j = 0; j < nq0; ++j)
869 for (k = 0; k < nq2; ++k)
871 int tmp = j + nq01 * k;
872 normals[j + k * nq0] = -df[1][tmp] * jac[tmp];
873 normals[nqtot + j + k * nq0] = -df[4][tmp] * jac[tmp];
874 normals[2 * nqtot + j + k * nq0] =
875 -df[7][tmp] * jac[tmp];
876 faceJac[j + k * nq0] = jac[tmp];
880 points0 = ptsKeys[0];
881 points1 = ptsKeys[2];
887 for (j = 0; j < nq1; ++j)
889 for (k = 0; k < nq2; ++k)
891 int tmp = nq0 - 1 + nq0 * j + nq01 * k;
892 normals[j + k * nq1] =
893 (df[0][tmp] + df[2][tmp]) * jac[tmp];
894 normals[nqtot + j + k * nq1] =
895 (df[3][tmp] + df[5][tmp]) * jac[tmp];
896 normals[2 * nqtot + j + k * nq1] =
897 (df[6][tmp] + df[8][tmp]) * jac[tmp];
898 faceJac[j + k * nq1] = jac[tmp];
902 points0 = ptsKeys[1];
903 points1 = ptsKeys[2];
909 for (j = 0; j < nq0; ++j)
911 for (k = 0; k < nq2; ++k)
913 int tmp = nq0 * (nq1 - 1) + j + nq01 * k;
914 normals[j + k * nq0] = df[1][tmp] * jac[tmp];
915 normals[nqtot + j + k * nq0] = df[4][tmp] * jac[tmp];
916 normals[2 * nqtot + j + k * nq0] =
917 df[7][tmp] * jac[tmp];
918 faceJac[j + k * nq0] = jac[tmp];
922 points0 = ptsKeys[0];
923 points1 = ptsKeys[2];
929 for (j = 0; j < nq1; ++j)
931 for (k = 0; k < nq2; ++k)
933 int tmp = j * nq0 + nq01 * k;
934 normals[j + k * nq1] = -df[0][tmp] * jac[tmp];
935 normals[nqtot + j + k * nq1] = -df[3][tmp] * jac[tmp];
936 normals[2 * nqtot + j + k * nq1] =
937 -df[6][tmp] * jac[tmp];
938 faceJac[j + k * nq1] = jac[tmp];
942 points0 = ptsKeys[1];
943 points1 = ptsKeys[2];
948 ASSERTL0(
false,
"face is out of range (face < 4)");
956 Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
959 for (i = 0; i < vCoordDim; ++i)
964 Vmath::Vmul(nq_face, work, 1, normal[i], 1, normal[i], 1);
971 Vmath::Vvtvp(nq_face, normal[i], 1, normal[i], 1, work, 1, work, 1);
981 Vmath::Vmul(nq_face, normal[i], 1, work, 1, normal[i], 1);
990 StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
1005 StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
1036 StdPrismExp::v_SVVLaplacianFilter(array, mkey);
1062 returnval = StdPrismExp::v_GenMatrix(mkey);
1078 return tmp->GetStdMatrix(mkey);
1124 int nquad0 =
m_base[0]->GetNumPoints();
1125 int nquad1 =
m_base[1]->GetNumPoints();
1126 int nquad2 =
m_base[2]->GetNumPoints();
1127 int nqtot = nquad0 * nquad1 * nquad2;
1161 StdExpansion3D::PhysTensorDeriv(inarray, wsp1, wsp2, wsp3);
1170 for (i = 0; i < nquad2; ++i)
1173 &h0[0] + i * nquad0 * nquad1, 1);
1175 &h1[0] + i * nquad0 * nquad1, 1);
1177 for (i = 0; i < nquad0; i++)
1179 Blas::Dscal(nquad1 * nquad2, 0.5 * (1 + z0[i]), &h1[0] + i, nquad0);
1188 Vmath::Vvtvvtp(nqtot, &df[0][0], 1, &h0[0], 1, &df[2][0], 1, &h1[0], 1,
1191 Vmath::Vvtvvtp(nqtot, &df[3][0], 1, &h0[0], 1, &df[5][0], 1, &h1[0], 1,
1194 Vmath::Vvtvvtp(nqtot, &df[6][0], 1, &h0[0], 1, &df[8][0], 1, &h1[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, &g0[0], 1, &g0[0], 1);
1203 Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &wsp4[0], 1, &df[4][0], 1, &wsp5[0],
1205 Vmath::Vvtvp(nqtot, &df[7][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, &g4[0], 1, &g4[0], 1);
1215 &df[4][0], 1, &g1[0], 1);
1216 Vmath::Vvtvp(nqtot, &df[7][0], 1, &df[7][0], 1, &g1[0], 1, &g1[0], 1);
1220 &df[5][0], 1, &g2[0], 1);
1221 Vmath::Vvtvp(nqtot, &df[8][0], 1, &df[8][0], 1, &g2[0], 1, &g2[0], 1);
1225 &df[5][0], 1, &g5[0], 1);
1226 Vmath::Vvtvp(nqtot, &df[7][0], 1, &df[8][0], 1, &g5[0], 1, &g5[0], 1);
1241 Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0],
1243 Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1246 Vmath::Svtsvtp(nqtot, df[1][0], &wsp4[0], 1, df[4][0], &wsp5[0], 1,
1248 Vmath::Svtvp(nqtot, df[7][0], &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1251 Vmath::Svtsvtp(nqtot, df[2][0], &wsp4[0], 1, df[5][0], &wsp5[0], 1,
1253 Vmath::Svtvp(nqtot, df[8][0], &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1258 df[1][0] * df[1][0] + df[4][0] * df[4][0] +
1259 df[7][0] * df[7][0],
1264 df[2][0] * df[2][0] + df[5][0] * df[5][0] +
1265 df[8][0] * df[8][0],
1270 df[1][0] * df[2][0] + df[4][0] * df[5][0] +
1271 df[7][0] * df[8][0],
1276 Vmath::Vvtvvtp(nqtot, &g0[0], 1, &wsp1[0], 1, &g3[0], 1, &wsp2[0], 1,
1278 Vmath::Vvtvp(nqtot, &g4[0], 1, &wsp3[0], 1, &wsp7[0], 1, &wsp7[0], 1);
1279 Vmath::Vvtvvtp(nqtot, &g1[0], 1, &wsp2[0], 1, &g3[0], 1, &wsp1[0], 1,
1281 Vmath::Vvtvp(nqtot, &g5[0], 1, &wsp3[0], 1, &wsp8[0], 1, &wsp8[0], 1);
1282 Vmath::Vvtvvtp(nqtot, &g2[0], 1, &wsp3[0], 1, &g4[0], 1, &wsp1[0], 1,
1284 Vmath::Vvtvp(nqtot, &g5[0], 1, &wsp2[0], 1, &wsp9[0], 1, &wsp9[0], 1);
1311 int np0 =
m_base[0]->GetNumPoints();
1312 int np1 =
m_base[1]->GetNumPoints();
1313 int np2 =
m_base[2]->GetNumPoints();
1314 int np = max(np0, max(np1, np2));
1316 bool standard =
true;
1318 int vid0 =
m_geom->GetVid(0);
1319 int vid1 =
m_geom->GetVid(1);
1320 int vid2 =
m_geom->GetVid(4);
1324 if ((vid2 < vid1) && (vid2 < vid0))
1332 else if ((vid1 < vid2) && (vid1 < vid0))
1340 else if ((vid0 < vid2) && (vid0 < vid1))
1361 rot[0] = (0 + rotate) % 3;
1362 rot[1] = (1 + rotate) % 3;
1363 rot[2] = (2 + rotate) % 3;
1366 for (
int i = 0; i < np - 1; ++i)
1368 planep1 += (np - i) * np;
1373 if (standard ==
false)
1375 for (
int j = 0; j < np - 1; ++j)
1378 row1p1 += np - i - 1;
1379 for (
int k = 0; k < np - i - 2; ++k)
1382 prismpt[rot[0]] = plane + row + k;
1383 prismpt[rot[1]] = plane + row + k + 1;
1384 prismpt[rot[2]] = planep1 + row1 + k;
1386 prismpt[3 + rot[0]] = plane + rowp1 + k;
1387 prismpt[3 + rot[1]] = plane + rowp1 + k + 1;
1388 prismpt[3 + rot[2]] = planep1 + row1p1 + k;
1390 conn[cnt++] = prismpt[0];
1391 conn[cnt++] = prismpt[1];
1392 conn[cnt++] = prismpt[3];
1393 conn[cnt++] = prismpt[2];
1395 conn[cnt++] = prismpt[5];
1396 conn[cnt++] = prismpt[2];
1397 conn[cnt++] = prismpt[3];
1398 conn[cnt++] = prismpt[4];
1400 conn[cnt++] = prismpt[3];
1401 conn[cnt++] = prismpt[1];
1402 conn[cnt++] = prismpt[4];
1403 conn[cnt++] = prismpt[2];
1406 prismpt[rot[0]] = planep1 + row1 + k + 1;
1407 prismpt[rot[1]] = planep1 + row1 + k;
1408 prismpt[rot[2]] = plane + row + k + 1;
1410 prismpt[3 + rot[0]] = planep1 + row1p1 + k + 1;
1411 prismpt[3 + rot[1]] = planep1 + row1p1 + k;
1412 prismpt[3 + rot[2]] = plane + rowp1 + k + 1;
1414 conn[cnt++] = prismpt[0];
1415 conn[cnt++] = prismpt[1];
1416 conn[cnt++] = prismpt[2];
1417 conn[cnt++] = prismpt[5];
1419 conn[cnt++] = prismpt[5];
1420 conn[cnt++] = prismpt[0];
1421 conn[cnt++] = prismpt[4];
1422 conn[cnt++] = prismpt[1];
1424 conn[cnt++] = prismpt[3];
1425 conn[cnt++] = prismpt[4];
1426 conn[cnt++] = prismpt[0];
1427 conn[cnt++] = prismpt[5];
1431 prismpt[rot[0]] = plane + row + np - i - 2;
1432 prismpt[rot[1]] = plane + row + np - i - 1;
1433 prismpt[rot[2]] = planep1 + row1 + np - i - 2;
1435 prismpt[3 + rot[0]] = plane + rowp1 + np - i - 2;
1436 prismpt[3 + rot[1]] = plane + rowp1 + np - i - 1;
1437 prismpt[3 + rot[2]] = planep1 + row1p1 + np - i - 2;
1439 conn[cnt++] = prismpt[0];
1440 conn[cnt++] = prismpt[1];
1441 conn[cnt++] = prismpt[3];
1442 conn[cnt++] = prismpt[2];
1444 conn[cnt++] = prismpt[5];
1445 conn[cnt++] = prismpt[2];
1446 conn[cnt++] = prismpt[3];
1447 conn[cnt++] = prismpt[4];
1449 conn[cnt++] = prismpt[3];
1450 conn[cnt++] = prismpt[1];
1451 conn[cnt++] = prismpt[4];
1452 conn[cnt++] = prismpt[2];
1460 for (
int j = 0; j < np - 1; ++j)
1463 row1p1 += np - i - 1;
1464 for (
int k = 0; k < np - i - 2; ++k)
1467 prismpt[rot[0]] = plane + row + k;
1468 prismpt[rot[1]] = plane + row + k + 1;
1469 prismpt[rot[2]] = planep1 + row1 + k;
1471 prismpt[3 + rot[0]] = plane + rowp1 + k;
1472 prismpt[3 + rot[1]] = plane + rowp1 + k + 1;
1473 prismpt[3 + rot[2]] = planep1 + row1p1 + k;
1475 conn[cnt++] = prismpt[0];
1476 conn[cnt++] = prismpt[1];
1477 conn[cnt++] = prismpt[4];
1478 conn[cnt++] = prismpt[2];
1480 conn[cnt++] = prismpt[4];
1481 conn[cnt++] = prismpt[3];
1482 conn[cnt++] = prismpt[0];
1483 conn[cnt++] = prismpt[2];
1485 conn[cnt++] = prismpt[3];
1486 conn[cnt++] = prismpt[4];
1487 conn[cnt++] = prismpt[5];
1488 conn[cnt++] = prismpt[2];
1491 prismpt[rot[0]] = planep1 + row1 + k + 1;
1492 prismpt[rot[1]] = planep1 + row1 + k;
1493 prismpt[rot[2]] = plane + row + k + 1;
1495 prismpt[3 + rot[0]] = planep1 + row1p1 + k + 1;
1496 prismpt[3 + rot[1]] = planep1 + row1p1 + k;
1497 prismpt[3 + rot[2]] = plane + rowp1 + k + 1;
1499 conn[cnt++] = prismpt[0];
1500 conn[cnt++] = prismpt[2];
1501 conn[cnt++] = prismpt[1];
1502 conn[cnt++] = prismpt[5];
1504 conn[cnt++] = prismpt[3];
1505 conn[cnt++] = prismpt[5];
1506 conn[cnt++] = prismpt[0];
1507 conn[cnt++] = prismpt[1];
1509 conn[cnt++] = prismpt[5];
1510 conn[cnt++] = prismpt[3];
1511 conn[cnt++] = prismpt[4];
1512 conn[cnt++] = prismpt[1];
1516 prismpt[rot[0]] = plane + row + np - i - 2;
1517 prismpt[rot[1]] = plane + row + np - i - 1;
1518 prismpt[rot[2]] = planep1 + row1 + np - i - 2;
1520 prismpt[3 + rot[0]] = plane + rowp1 + np - i - 2;
1521 prismpt[3 + rot[1]] = plane + rowp1 + np - i - 1;
1522 prismpt[3 + rot[2]] = planep1 + row1p1 + np - i - 2;
1524 conn[cnt++] = prismpt[0];
1525 conn[cnt++] = prismpt[1];
1526 conn[cnt++] = prismpt[4];
1527 conn[cnt++] = prismpt[2];
1529 conn[cnt++] = prismpt[4];
1530 conn[cnt++] = prismpt[3];
1531 conn[cnt++] = prismpt[0];
1532 conn[cnt++] = prismpt[2];
1534 conn[cnt++] = prismpt[3];
1535 conn[cnt++] = prismpt[4];
1536 conn[cnt++] = prismpt[5];
1537 conn[cnt++] = prismpt[2];
1543 plane += (np - i) * np;
1564 if (d0factors.size() != 5)
1571 if (d0factors[0].size() != nquad0 * nquad1)
1578 if (d0factors[1].size() != nquad0 * nquad2)
1588 if (d0factors[2].size() != nquad1 * nquad2)
1610 int ncoords = normal_0.size();
1616 for (
int i = 0; i < nquad0 * nquad1; ++i)
1618 d0factors[0][i] = df[0][i] * normal_0[0][i];
1619 d1factors[0][i] = df[1][i] * normal_0[0][i];
1620 d2factors[0][i] = df[2][i] * normal_0[0][i];
1623 for (
int n = 1; n < ncoords; ++n)
1625 for (
int i = 0; i < nquad0 * nquad1; ++i)
1627 d0factors[0][i] += df[3 * n][i] * normal_0[n][i];
1628 d1factors[0][i] += df[3 * n + 1][i] * normal_0[n][i];
1629 d2factors[0][i] += df[3 * n + 2][i] * normal_0[n][i];
1634 for (
int j = 0; j < nquad2; ++j)
1636 for (
int i = 0; i < nquad0; ++i)
1638 d0factors[1][j * nquad0 + i] = df[0][j * nquad0 * nquad1 + i] *
1639 normal_1[0][j * nquad0 + i];
1640 d1factors[1][j * nquad0 + i] = df[1][j * nquad0 * nquad1 + i] *
1641 normal_1[0][j * nquad0 + i];
1642 d2factors[1][j * nquad0 + i] = df[2][j * nquad0 * nquad1 + i] *
1643 normal_1[0][j * nquad0 + i];
1645 d0factors[3][j * nquad0 + i] =
1646 df[0][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1647 normal_3[0][j * nquad0 + i];
1648 d1factors[3][j * nquad0 + i] =
1649 df[1][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1650 normal_3[0][j * nquad0 + i];
1651 d2factors[3][j * nquad0 + i] =
1652 df[2][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1653 normal_3[0][j * nquad0 + i];
1657 for (
int n = 1; n < ncoords; ++n)
1659 for (
int j = 0; j < nquad2; ++j)
1661 for (
int i = 0; i < nquad0; ++i)
1663 d0factors[1][j * nquad0 + i] +=
1664 df[3 * n][j * nquad0 * nquad1 + i] *
1665 normal_1[n][j * nquad0 + i];
1666 d1factors[1][j * nquad0 + i] +=
1667 df[3 * n + 1][j * nquad0 * nquad1 + i] *
1668 normal_1[n][j * nquad0 + i];
1669 d2factors[1][j * nquad0 + i] +=
1670 df[3 * n + 2][j * nquad0 * nquad1 + i] *
1671 normal_1[n][j * nquad0 + i];
1673 d0factors[3][j * nquad0 + i] +=
1674 df[3 * n][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1675 normal_3[n][j * nquad0 + i];
1676 d1factors[3][j * nquad0 + i] +=
1677 df[3 * n + 1][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1678 normal_3[n][j * nquad0 + i];
1679 d2factors[3][j * nquad0 + i] +=
1680 df[3 * n + 2][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1681 normal_3[n][j * nquad0 + i];
1687 for (
int j = 0; j < nquad2; ++j)
1689 for (
int i = 0; i < nquad1; ++i)
1691 d0factors[2][j * nquad1 + i] =
1692 df[0][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1693 normal_2[0][j * nquad1 + i];
1694 d1factors[2][j * nquad1 + i] =
1695 df[1][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1696 normal_2[0][j * nquad1 + i];
1697 d2factors[2][j * nquad1 + i] =
1698 df[2][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1699 normal_2[0][j * nquad1 + i];
1701 d0factors[4][j * nquad1 + i] =
1702 df[0][j * nquad0 * nquad1 + i * nquad0] *
1703 normal_4[0][j * nquad1 + i];
1704 d1factors[4][j * nquad1 + i] =
1705 df[1][j * nquad0 * nquad1 + i * nquad0] *
1706 normal_4[0][j * nquad1 + i];
1707 d2factors[4][j * nquad1 + i] =
1708 df[2][j * nquad0 * nquad1 + i * nquad0] *
1709 normal_4[0][j * nquad1 + i];
1713 for (
int n = 1; n < ncoords; ++n)
1715 for (
int j = 0; j < nquad2; ++j)
1717 for (
int i = 0; i < nquad1; ++i)
1719 d0factors[2][j * nquad1 + i] +=
1720 df[3 * n][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1721 normal_2[n][j * nquad1 + i];
1722 d1factors[2][j * nquad1 + i] +=
1724 [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1725 normal_2[n][j * nquad1 + i];
1726 d2factors[2][j * nquad1 + i] +=
1728 [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1729 normal_2[n][j * nquad1 + i];
1731 d0factors[4][j * nquad1 + i] +=
1732 df[3 * n][j * nquad0 * nquad1 + i * nquad0] *
1733 normal_4[n][j * nquad1 + i];
1734 d1factors[4][j * nquad1 + i] +=
1735 df[3 * n + 1][j * nquad0 * nquad1 + i * nquad0] *
1736 normal_4[n][j * nquad1 + i];
1737 d2factors[4][j * nquad1 + i] +=
1738 df[3 * n + 2][j * nquad0 * nquad1 + i * nquad0] *
1739 normal_4[n][j * nquad1 + i];
1747 for (
int i = 0; i < nquad0 * nquad1; ++i)
1749 d0factors[0][i] = df[0][0] * normal_0[0][i];
1750 d1factors[0][i] = df[1][0] * normal_0[0][i];
1751 d2factors[0][i] = df[2][0] * normal_0[0][i];
1754 for (
int n = 1; n < ncoords; ++n)
1756 for (
int i = 0; i < nquad0 * nquad1; ++i)
1758 d0factors[0][i] += df[3 * n][0] * normal_0[n][i];
1759 d1factors[0][i] += df[3 * n + 1][0] * normal_0[n][i];
1760 d2factors[0][i] += df[3 * n + 2][0] * normal_0[n][i];
1765 for (
int i = 0; i < nquad0 * nquad2; ++i)
1767 d0factors[1][i] = df[0][0] * normal_1[0][i];
1768 d0factors[3][i] = df[0][0] * normal_3[0][i];
1770 d1factors[1][i] = df[1][0] * normal_1[0][i];
1771 d1factors[3][i] = df[1][0] * normal_3[0][i];
1773 d2factors[1][i] = df[2][0] * normal_1[0][i];
1774 d2factors[3][i] = df[2][0] * normal_3[0][i];
1777 for (
int n = 1; n < ncoords; ++n)
1779 for (
int i = 0; i < nquad0 * nquad2; ++i)
1781 d0factors[1][i] += df[3 * n][0] * normal_1[n][i];
1782 d0factors[3][i] += df[3 * n][0] * normal_3[n][i];
1784 d1factors[1][i] += df[3 * n + 1][0] * normal_1[n][i];
1785 d1factors[3][i] += df[3 * n + 1][0] * normal_3[n][i];
1787 d2factors[1][i] += df[3 * n + 2][0] * normal_1[n][i];
1788 d2factors[3][i] += df[3 * n + 2][0] * normal_3[n][i];
1793 for (
int i = 0; i < nquad1 * nquad2; ++i)
1795 d0factors[2][i] = df[0][0] * normal_2[0][i];
1796 d0factors[4][i] = df[0][0] * normal_4[0][i];
1798 d1factors[2][i] = df[1][0] * normal_2[0][i];
1799 d1factors[4][i] = df[1][0] * normal_4[0][i];
1801 d2factors[2][i] = df[2][0] * normal_2[0][i];
1802 d2factors[4][i] = df[2][0] * normal_4[0][i];
1805 for (
int n = 1; n < ncoords; ++n)
1807 for (
int i = 0; i < nquad1 * nquad2; ++i)
1809 d0factors[2][i] += df[3 * n][0] * normal_2[n][i];
1810 d0factors[4][i] += df[3 * n][0] * normal_4[n][i];
1812 d1factors[2][i] += df[3 * n + 1][0] * normal_2[n][i];
1813 d1factors[4][i] += df[3 * n + 1][0] * normal_4[n][i];
1815 d2factors[2][i] += df[3 * n + 2][0] * normal_2[n][i];
1816 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 v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
const NormalVector & GetTraceNormal(const int id)
DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey) override
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const override
void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
void v_DropLocStaticCondMatrix(const MatrixKey &mkey) override
void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, 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_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true) 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
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Calculates the inner product .
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.
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true) 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_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords) override
Get the coordinates #coords at the local coordinates #Lcoords.
void v_NormalTraceDerivFactors(Array< OneD, Array< OneD, NekDouble > > &d0factors, Array< OneD, Array< OneD, NekDouble > > &d1factors, Array< OneD, Array< OneD, NekDouble > > &d2factors) override
: This method gets all of the factors which are required as part of the Gradient Jump Penalty stabili...
void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray) override
StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const override
void v_ComputeTraceNormal(const int face) override
Get the normals along specficied face Get the face normals interplated to a points0 x points 0 type d...
void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp) override
Calculate the Laplacian multiplication in a matrix-free manner.
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals) 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...
PrismExp(const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const LibUtilities::BasisKey &Bc, const SpatialDomains::PrismGeomSharedPtr &geom)
Constructor using BasisKey class for quadrature points and order definition.
void v_GetTracePhysMap(const int face, Array< OneD, int > &outarray) override
void v_DropLocMatrix(const MatrixKey &mkey) override
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_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey) override
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey) override
NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray) override
Integrate the physical point list inarray over prismatic region and return the value.
void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) 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
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 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.
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
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
@ eModified_B
Principle Modified Functions .
@ eModified_A
Principle Modified Functions .
std::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
std::shared_ptr< PrismGeom > PrismGeomSharedPtr
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< StdPrismExp > StdPrismExpSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
@ eInvLaplacianWithUnityMean
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)