35 #include <boost/core/ignore_unused.hpp>
46 namespace LocalRegions
54 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
57 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
61 std::bind(&
Expansion3D::CreateMatrix, this, std::placeholders::_1),
62 std::string(
"PrismExpMatrix")),
63 m_staticCondMatrixManager(std::bind(&
Expansion::CreateStaticCondMatrix,
64 this, std::placeholders::_1),
65 std::string(
"PrismExpStaticCondMatrix"))
70 : StdExpansion(T), StdExpansion3D(T), StdPrismExp(T),
Expansion(T),
72 m_staticCondMatrixManager(T.m_staticCondMatrixManager)
106 int nquad0 =
m_base[0]->GetNumPoints();
107 int nquad1 =
m_base[1]->GetNumPoints();
108 int nquad2 =
m_base[2]->GetNumPoints();
116 (
NekDouble *)&inarray[0], 1, &tmp[0], 1);
121 (
NekDouble *)&inarray[0], 1, &tmp[0], 1);
125 return StdPrismExp::v_Integral(tmp);
144 StdPrismExp::v_PhysDeriv(inarray, diff0, diff1, diff2);
150 Vmath::Vmul(nqtot, &df[0][0], 1, &diff0[0], 1, &out_d0[0], 1);
151 Vmath::Vvtvp(nqtot, &df[1][0], 1, &diff1[0], 1, &out_d0[0], 1,
153 Vmath::Vvtvp(nqtot, &df[2][0], 1, &diff2[0], 1, &out_d0[0], 1,
159 Vmath::Vmul(nqtot, &df[3][0], 1, &diff0[0], 1, &out_d1[0], 1);
160 Vmath::Vvtvp(nqtot, &df[4][0], 1, &diff1[0], 1, &out_d1[0], 1,
162 Vmath::Vvtvp(nqtot, &df[5][0], 1, &diff2[0], 1, &out_d1[0], 1,
168 Vmath::Vmul(nqtot, &df[6][0], 1, &diff0[0], 1, &out_d2[0], 1);
169 Vmath::Vvtvp(nqtot, &df[7][0], 1, &diff1[0], 1, &out_d2[0], 1,
171 Vmath::Vvtvp(nqtot, &df[8][0], 1, &diff2[0], 1, &out_d2[0], 1,
179 Vmath::Smul(nqtot, df[0][0], &diff0[0], 1, &out_d0[0], 1);
180 Blas::Daxpy(nqtot, df[1][0], &diff1[0], 1, &out_d0[0], 1);
181 Blas::Daxpy(nqtot, df[2][0], &diff2[0], 1, &out_d0[0], 1);
186 Vmath::Smul(nqtot, df[3][0], &diff0[0], 1, &out_d1[0], 1);
187 Blas::Daxpy(nqtot, df[4][0], &diff1[0], 1, &out_d1[0], 1);
188 Blas::Daxpy(nqtot, df[5][0], &diff2[0], 1, &out_d1[0], 1);
193 Vmath::Smul(nqtot, df[6][0], &diff0[0], 1, &out_d2[0], 1);
194 Blas::Daxpy(nqtot, df[7][0], &diff1[0], 1, &out_d2[0], 1);
195 Blas::Daxpy(nqtot, df[8][0], &diff2[0], 1, &out_d2[0], 1);
220 if (
m_base[0]->Collocation() &&
m_base[1]->Collocation() &&
237 out = (*matsys) * in;
279 const int nquad0 =
m_base[0]->GetNumPoints();
280 const int nquad1 =
m_base[1]->GetNumPoints();
281 const int nquad2 =
m_base[2]->GetNumPoints();
282 const int order0 =
m_base[0]->GetNumModes();
283 const int order1 =
m_base[1]->GetNumModes();
287 if (multiplybyweights)
295 tmp, outarray, wsp,
true,
true,
true);
301 inarray, outarray, wsp,
true,
true,
true);
346 const int nquad0 =
m_base[0]->GetNumPoints();
347 const int nquad1 =
m_base[1]->GetNumPoints();
348 const int nquad2 =
m_base[2]->GetNumPoints();
349 const int order0 =
m_base[0]->GetNumModes();
350 const int order1 =
m_base[1]->GetNumModes();
351 const int nqtot = nquad0 * nquad1 * nquad2;
370 m_base[2]->GetBdata(), tmp2, outarray, wsp,
374 m_base[2]->GetBdata(), tmp3, tmp6, wsp,
true,
380 m_base[2]->GetDbdata(), tmp4, tmp6, wsp,
true,
390 const int nquad0 =
m_base[0]->GetNumPoints();
391 const int nquad1 =
m_base[1]->GetNumPoints();
392 const int nquad2 =
m_base[2]->GetNumPoints();
393 const int order0 =
m_base[0]->GetNumModes();
394 const int order1 =
m_base[1]->GetNumModes();
395 const int nqtot = nquad0 * nquad1 * nquad2;
418 Vmath::Vmul(nqtot, &df[3 * dir][0], 1, tmp1.get(), 1, tmp2.get(), 1);
419 Vmath::Vmul(nqtot, &df[3 * dir + 1][0], 1, tmp1.get(), 1, tmp3.get(),
421 Vmath::Vmul(nqtot, &df[3 * dir + 2][0], 1, tmp1.get(), 1, tmp4.get(),
426 Vmath::Smul(nqtot, df[3 * dir][0], tmp1.get(), 1, tmp2.get(), 1);
427 Vmath::Smul(nqtot, df[3 * dir + 1][0], tmp1.get(), 1, tmp3.get(), 1);
428 Vmath::Smul(nqtot, df[3 * dir + 2][0], tmp1.get(), 1, tmp4.get(), 1);
432 for (
int i = 0; i < nquad0; ++i)
434 gfac0[i] = 0.5 * (1 + z0[i]);
438 for (
int i = 0; i < nquad2; ++i)
440 gfac2[i] = 2.0 / (1 - z2[i]);
443 const int nq01 = nquad0 * nquad1;
445 for (
int i = 0; i < nquad2; ++i)
447 Vmath::Smul(nq01, gfac2[i], &tmp2[0] + i * nq01, 1, &tmp2[0] + i * nq01,
449 Vmath::Smul(nq01, gfac2[i], &tmp4[0] + i * nq01, 1, &tmp5[0] + i * nq01,
453 for (
int i = 0; i < nquad1 * nquad2; ++i)
455 Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp5[0] + i * nquad0, 1,
456 &tmp5[0] + i * nquad0, 1);
459 Vmath::Vadd(nqtot, &tmp2[0], 1, &tmp5[0], 1, &tmp2[0], 1);
470 m_base[2]->GetBasisKey());
476 m_base[0]->GetPointsKey());
478 m_base[1]->GetPointsKey());
480 m_base[2]->GetPointsKey());
483 bkey0, bkey1, bkey2);
495 ASSERTL1(Lcoords[0] <= -1.0 && Lcoords[0] >= 1.0 && Lcoords[1] <= -1.0 &&
496 Lcoords[1] >= 1.0 && Lcoords[2] <= -1.0 && Lcoords[2] >= 1.0,
497 "Local coordinates are not in region [-1,1]");
501 for (i = 0; i <
m_geom->GetCoordim(); ++i)
503 coords[i] =
m_geom->GetCoord(i, Lcoords);
524 return StdPrismExp::v_PhysEvaluate(Lcoord, physvals);
534 m_geom->GetLocCoords(coord, Lcoord);
536 return StdPrismExp::v_PhysEvaluate(Lcoord, physvals);
545 return m_geom->GetCoordim();
549 const NekDouble *data,
const std::vector<unsigned int> &nummodes,
550 const int mode_offset,
NekDouble *coeffs,
551 std::vector<LibUtilities::BasisType> &fromType)
553 boost::ignore_unused(fromType);
555 int data_order0 = nummodes[mode_offset];
556 int fillorder0 = min(
m_base[0]->GetNumModes(), data_order0);
557 int data_order1 = nummodes[mode_offset + 1];
558 int order1 =
m_base[1]->GetNumModes();
559 int fillorder1 = min(order1, data_order1);
560 int data_order2 = nummodes[mode_offset + 2];
561 int order2 =
m_base[2]->GetNumModes();
562 int fillorder2 = min(order2, data_order2);
573 "Extraction routine not set up for this basis");
575 "Extraction routine not set up for this basis");
578 for (j = 0; j < fillorder0; ++j)
580 for (i = 0; i < fillorder1; ++i)
582 Vmath::Vcopy(fillorder2 - j, &data[cnt], 1, &coeffs[cnt1],
584 cnt += data_order2 - j;
589 for (i = fillorder1; i < data_order1; ++i)
591 cnt += data_order2 - j;
594 for (i = fillorder1; i < order1; ++i)
602 ASSERTL0(
false,
"basis is either not set up or not "
609 int nquad0 =
m_base[0]->GetNumPoints();
610 int nquad1 =
m_base[1]->GetNumPoints();
611 int nquad2 =
m_base[2]->GetNumPoints();
620 if (outarray.size() != nq0 * nq1)
626 for (
int i = 0; i < nquad0 * nquad1; ++i)
635 if (outarray.size() != nq0 * nq1)
641 for (
int k = 0; k < nquad2; k++)
643 for (
int i = 0; i < nquad0; ++i)
645 outarray[k * nquad0 + i] = (nquad0 * nquad1 * k) + i;
654 if (outarray.size() != nq0 * nq1)
660 for (
int j = 0; j < nquad1 * nquad2; ++j)
662 outarray[j] = nquad0 - 1 + j * nquad0;
668 if (outarray.size() != nq0 * nq1)
674 for (
int k = 0; k < nquad2; k++)
676 for (
int i = 0; i < nquad0; ++i)
678 outarray[k * nquad0 + i] =
679 nquad0 * (nquad1 - 1) + (nquad0 * nquad1 * k) + i;
687 if (outarray.size() != nq0 * nq1)
693 for (
int j = 0; j < nquad1 * nquad2; ++j)
695 outarray[j] = j * nquad0;
699 ASSERTL0(
false,
"face value (> 4) is out of range");
714 for (
int i = 0; i < ptsKeys.size(); ++i)
726 geomFactors->GetDerivFactors(ptsKeys);
729 int nq0 = ptsKeys[0].GetNumPoints();
730 int nq1 = ptsKeys[1].GetNumPoints();
731 int nq2 = ptsKeys[2].GetNumPoints();
732 int nq01 = nq0 * nq1;
746 for (i = 0; i < vCoordDim; ++i)
751 size_t nqb = nq_face;
766 for (i = 0; i < vCoordDim; ++i)
768 normal[i][0] = -df[3 * i + 2][0];
775 for (i = 0; i < vCoordDim; ++i)
777 normal[i][0] = -df[3 * i + 1][0];
783 for (i = 0; i < vCoordDim; ++i)
785 normal[i][0] = df[3 * i][0] + df[3 * i + 2][0];
791 for (i = 0; i < vCoordDim; ++i)
793 normal[i][0] = df[3 * i + 1][0];
799 for (i = 0; i < vCoordDim; ++i)
801 normal[i][0] = -df[3 * i][0];
806 ASSERTL0(
false,
"face is out of range (face < 4)");
811 for (i = 0; i < vCoordDim; ++i)
813 fac += normal[i][0] * normal[i][0];
815 fac = 1.0 /
sqrt(fac);
819 for (i = 0; i < vCoordDim; ++i)
821 Vmath::Fill(nq_face, fac * normal[i][0], normal[i], 1);
834 else if (face == 1 || face == 3)
856 for (j = 0; j < nq01; ++j)
858 normals[j] = -df[2][j] * jac[j];
859 normals[nqtot + j] = -df[5][j] * jac[j];
860 normals[2 * nqtot + j] = -df[8][j] * jac[j];
864 points0 = ptsKeys[0];
865 points1 = ptsKeys[1];
871 for (j = 0; j < nq0; ++j)
873 for (k = 0; k < nq2; ++k)
875 int tmp = j + nq01 * k;
876 normals[j + k * nq0] = -df[1][tmp] * jac[tmp];
877 normals[nqtot + j + k * nq0] = -df[4][tmp] * jac[tmp];
878 normals[2 * nqtot + j + k * nq0] =
879 -df[7][tmp] * jac[tmp];
880 faceJac[j + k * nq0] = jac[tmp];
884 points0 = ptsKeys[0];
885 points1 = ptsKeys[2];
891 for (j = 0; j < nq1; ++j)
893 for (k = 0; k < nq2; ++k)
895 int tmp = nq0 - 1 + nq0 * j + nq01 * k;
896 normals[j + k * nq1] =
897 (df[0][tmp] + df[2][tmp]) * jac[tmp];
898 normals[nqtot + j + k * nq1] =
899 (df[3][tmp] + df[5][tmp]) * jac[tmp];
900 normals[2 * nqtot + j + k * nq1] =
901 (df[6][tmp] + df[8][tmp]) * jac[tmp];
902 faceJac[j + k * nq1] = jac[tmp];
906 points0 = ptsKeys[1];
907 points1 = ptsKeys[2];
913 for (j = 0; j < nq0; ++j)
915 for (k = 0; k < nq2; ++k)
917 int tmp = nq0 * (nq1 - 1) + j + nq01 * k;
918 normals[j + k * nq0] = df[1][tmp] * jac[tmp];
919 normals[nqtot + j + k * nq0] = df[4][tmp] * jac[tmp];
920 normals[2 * nqtot + j + k * nq0] =
921 df[7][tmp] * jac[tmp];
922 faceJac[j + k * nq0] = jac[tmp];
926 points0 = ptsKeys[0];
927 points1 = ptsKeys[2];
933 for (j = 0; j < nq1; ++j)
935 for (k = 0; k < nq2; ++k)
937 int tmp = j * nq0 + nq01 * k;
938 normals[j + k * nq1] = -df[0][tmp] * jac[tmp];
939 normals[nqtot + j + k * nq1] = -df[3][tmp] * jac[tmp];
940 normals[2 * nqtot + j + k * nq1] =
941 -df[6][tmp] * jac[tmp];
942 faceJac[j + k * nq1] = jac[tmp];
946 points0 = ptsKeys[1];
947 points1 = ptsKeys[2];
952 ASSERTL0(
false,
"face is out of range (face < 4)");
960 Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
963 for (i = 0; i < vCoordDim; ++i)
968 Vmath::Vmul(nq_face, work, 1, normal[i], 1, normal[i], 1);
975 Vmath::Vvtvp(nq_face, normal[i], 1, normal[i], 1, work, 1, work, 1);
985 Vmath::Vmul(nq_face, normal[i], 1, work, 1, normal[i], 1);
994 StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
1009 StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
1025 if (inarray.get() == outarray.get())
1031 (mat->GetOwnedMatrix())->GetPtr().get(),
m_ncoeffs,
1032 tmp.get(), 1, 0.0, outarray.get(), 1);
1037 (mat->GetOwnedMatrix())->GetPtr().get(),
m_ncoeffs,
1038 inarray.get(), 1, 0.0, outarray.get(), 1);
1063 StdPrismExp::v_SVVLaplacianFilter(array, mkey);
1089 returnval = StdPrismExp::v_GenMatrix(mkey);
1105 return tmp->GetStdMatrix(mkey);
1146 int nquad0 =
m_base[0]->GetNumPoints();
1147 int nquad1 =
m_base[1]->GetNumPoints();
1148 int nquad2 =
m_base[2]->GetNumPoints();
1149 int nqtot = nquad0 * nquad1 * nquad2;
1183 StdExpansion3D::PhysTensorDeriv(inarray, wsp1, wsp2, wsp3);
1192 for (i = 0; i < nquad2; ++i)
1195 &h0[0] + i * nquad0 * nquad1, 1);
1197 &h1[0] + i * nquad0 * nquad1, 1);
1199 for (i = 0; i < nquad0; i++)
1201 Blas::Dscal(nquad1 * nquad2, 0.5 * (1 + z0[i]), &h1[0] + i, nquad0);
1210 Vmath::Vvtvvtp(nqtot, &df[0][0], 1, &h0[0], 1, &df[2][0], 1, &h1[0], 1,
1213 Vmath::Vvtvvtp(nqtot, &df[3][0], 1, &h0[0], 1, &df[5][0], 1, &h1[0], 1,
1216 Vmath::Vvtvvtp(nqtot, &df[6][0], 1, &h0[0], 1, &df[8][0], 1, &h1[0], 1,
1220 Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0],
1222 Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1225 Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &wsp4[0], 1, &df[4][0], 1, &wsp5[0],
1227 Vmath::Vvtvp(nqtot, &df[7][0], 1, &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1230 Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp4[0], 1, &df[5][0], 1, &wsp5[0],
1232 Vmath::Vvtvp(nqtot, &df[8][0], 1, &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1237 &df[4][0], 1, &g1[0], 1);
1238 Vmath::Vvtvp(nqtot, &df[7][0], 1, &df[7][0], 1, &g1[0], 1, &g1[0], 1);
1242 &df[5][0], 1, &g2[0], 1);
1243 Vmath::Vvtvp(nqtot, &df[8][0], 1, &df[8][0], 1, &g2[0], 1, &g2[0], 1);
1247 &df[5][0], 1, &g5[0], 1);
1248 Vmath::Vvtvp(nqtot, &df[7][0], 1, &df[8][0], 1, &g5[0], 1, &g5[0], 1);
1263 Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0],
1265 Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1268 Vmath::Svtsvtp(nqtot, df[1][0], &wsp4[0], 1, df[4][0], &wsp5[0], 1,
1270 Vmath::Svtvp(nqtot, df[7][0], &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1273 Vmath::Svtsvtp(nqtot, df[2][0], &wsp4[0], 1, df[5][0], &wsp5[0], 1,
1275 Vmath::Svtvp(nqtot, df[8][0], &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1280 df[1][0] * df[1][0] + df[4][0] * df[4][0] +
1281 df[7][0] * df[7][0],
1286 df[2][0] * df[2][0] + df[5][0] * df[5][0] +
1287 df[8][0] * df[8][0],
1292 df[1][0] * df[2][0] + df[4][0] * df[5][0] +
1293 df[7][0] * df[8][0],
1298 Vmath::Vvtvvtp(nqtot, &g0[0], 1, &wsp1[0], 1, &g3[0], 1, &wsp2[0], 1,
1300 Vmath::Vvtvp(nqtot, &g4[0], 1, &wsp3[0], 1, &wsp7[0], 1, &wsp7[0], 1);
1301 Vmath::Vvtvvtp(nqtot, &g1[0], 1, &wsp2[0], 1, &g3[0], 1, &wsp1[0], 1,
1303 Vmath::Vvtvp(nqtot, &g5[0], 1, &wsp3[0], 1, &wsp8[0], 1, &wsp8[0], 1);
1304 Vmath::Vvtvvtp(nqtot, &g2[0], 1, &wsp3[0], 1, &g4[0], 1, &wsp1[0], 1,
1306 Vmath::Vvtvp(nqtot, &g5[0], 1, &wsp2[0], 1, &wsp9[0], 1, &wsp9[0], 1);
1331 boost::ignore_unused(oldstandard);
1333 int np0 =
m_base[0]->GetNumPoints();
1334 int np1 =
m_base[1]->GetNumPoints();
1335 int np2 =
m_base[2]->GetNumPoints();
1336 int np = max(np0, max(np1, np2));
1338 bool standard =
true;
1340 int vid0 =
m_geom->GetVid(0);
1341 int vid1 =
m_geom->GetVid(1);
1342 int vid2 =
m_geom->GetVid(4);
1346 if ((vid2 < vid1) && (vid2 < vid0))
1354 else if ((vid1 < vid2) && (vid1 < vid0))
1362 else if ((vid0 < vid2) && (vid0 < vid1))
1383 rot[0] = (0 + rotate) % 3;
1384 rot[1] = (1 + rotate) % 3;
1385 rot[2] = (2 + rotate) % 3;
1388 for (
int i = 0; i < np - 1; ++i)
1390 planep1 += (np - i) * np;
1395 if (standard ==
false)
1397 for (
int j = 0; j < np - 1; ++j)
1400 row1p1 += np - i - 1;
1401 for (
int k = 0; k < np - i - 2; ++k)
1404 prismpt[rot[0]] = plane + row + k;
1405 prismpt[rot[1]] = plane + row + k + 1;
1406 prismpt[rot[2]] = planep1 + row1 + k;
1408 prismpt[3 + rot[0]] = plane + rowp1 + k;
1409 prismpt[3 + rot[1]] = plane + rowp1 + k + 1;
1410 prismpt[3 + rot[2]] = planep1 + row1p1 + k;
1412 conn[cnt++] = prismpt[0];
1413 conn[cnt++] = prismpt[1];
1414 conn[cnt++] = prismpt[3];
1415 conn[cnt++] = prismpt[2];
1417 conn[cnt++] = prismpt[5];
1418 conn[cnt++] = prismpt[2];
1419 conn[cnt++] = prismpt[3];
1420 conn[cnt++] = prismpt[4];
1422 conn[cnt++] = prismpt[3];
1423 conn[cnt++] = prismpt[1];
1424 conn[cnt++] = prismpt[4];
1425 conn[cnt++] = prismpt[2];
1428 prismpt[rot[0]] = planep1 + row1 + k + 1;
1429 prismpt[rot[1]] = planep1 + row1 + k;
1430 prismpt[rot[2]] = plane + row + k + 1;
1432 prismpt[3 + rot[0]] = planep1 + row1p1 + k + 1;
1433 prismpt[3 + rot[1]] = planep1 + row1p1 + k;
1434 prismpt[3 + rot[2]] = plane + rowp1 + k + 1;
1436 conn[cnt++] = prismpt[0];
1437 conn[cnt++] = prismpt[1];
1438 conn[cnt++] = prismpt[2];
1439 conn[cnt++] = prismpt[5];
1441 conn[cnt++] = prismpt[5];
1442 conn[cnt++] = prismpt[0];
1443 conn[cnt++] = prismpt[4];
1444 conn[cnt++] = prismpt[1];
1446 conn[cnt++] = prismpt[3];
1447 conn[cnt++] = prismpt[4];
1448 conn[cnt++] = prismpt[0];
1449 conn[cnt++] = prismpt[5];
1453 prismpt[rot[0]] = plane + row + np - i - 2;
1454 prismpt[rot[1]] = plane + row + np - i - 1;
1455 prismpt[rot[2]] = planep1 + row1 + np - i - 2;
1457 prismpt[3 + rot[0]] = plane + rowp1 + np - i - 2;
1458 prismpt[3 + rot[1]] = plane + rowp1 + np - i - 1;
1459 prismpt[3 + rot[2]] = planep1 + row1p1 + np - i - 2;
1461 conn[cnt++] = prismpt[0];
1462 conn[cnt++] = prismpt[1];
1463 conn[cnt++] = prismpt[3];
1464 conn[cnt++] = prismpt[2];
1466 conn[cnt++] = prismpt[5];
1467 conn[cnt++] = prismpt[2];
1468 conn[cnt++] = prismpt[3];
1469 conn[cnt++] = prismpt[4];
1471 conn[cnt++] = prismpt[3];
1472 conn[cnt++] = prismpt[1];
1473 conn[cnt++] = prismpt[4];
1474 conn[cnt++] = prismpt[2];
1482 for (
int j = 0; j < np - 1; ++j)
1485 row1p1 += np - i - 1;
1486 for (
int k = 0; k < np - i - 2; ++k)
1489 prismpt[rot[0]] = plane + row + k;
1490 prismpt[rot[1]] = plane + row + k + 1;
1491 prismpt[rot[2]] = planep1 + row1 + k;
1493 prismpt[3 + rot[0]] = plane + rowp1 + k;
1494 prismpt[3 + rot[1]] = plane + rowp1 + k + 1;
1495 prismpt[3 + rot[2]] = planep1 + row1p1 + k;
1497 conn[cnt++] = prismpt[0];
1498 conn[cnt++] = prismpt[1];
1499 conn[cnt++] = prismpt[4];
1500 conn[cnt++] = prismpt[2];
1502 conn[cnt++] = prismpt[4];
1503 conn[cnt++] = prismpt[3];
1504 conn[cnt++] = prismpt[0];
1505 conn[cnt++] = prismpt[2];
1507 conn[cnt++] = prismpt[3];
1508 conn[cnt++] = prismpt[4];
1509 conn[cnt++] = prismpt[5];
1510 conn[cnt++] = prismpt[2];
1513 prismpt[rot[0]] = planep1 + row1 + k + 1;
1514 prismpt[rot[1]] = planep1 + row1 + k;
1515 prismpt[rot[2]] = plane + row + k + 1;
1517 prismpt[3 + rot[0]] = planep1 + row1p1 + k + 1;
1518 prismpt[3 + rot[1]] = planep1 + row1p1 + k;
1519 prismpt[3 + rot[2]] = plane + rowp1 + k + 1;
1521 conn[cnt++] = prismpt[0];
1522 conn[cnt++] = prismpt[2];
1523 conn[cnt++] = prismpt[1];
1524 conn[cnt++] = prismpt[5];
1526 conn[cnt++] = prismpt[3];
1527 conn[cnt++] = prismpt[5];
1528 conn[cnt++] = prismpt[0];
1529 conn[cnt++] = prismpt[1];
1531 conn[cnt++] = prismpt[5];
1532 conn[cnt++] = prismpt[3];
1533 conn[cnt++] = prismpt[4];
1534 conn[cnt++] = prismpt[1];
1538 prismpt[rot[0]] = plane + row + np - i - 2;
1539 prismpt[rot[1]] = plane + row + np - i - 1;
1540 prismpt[rot[2]] = planep1 + row1 + np - i - 2;
1542 prismpt[3 + rot[0]] = plane + rowp1 + np - i - 2;
1543 prismpt[3 + rot[1]] = plane + rowp1 + np - i - 1;
1544 prismpt[3 + rot[2]] = planep1 + row1p1 + np - i - 2;
1546 conn[cnt++] = prismpt[0];
1547 conn[cnt++] = prismpt[1];
1548 conn[cnt++] = prismpt[4];
1549 conn[cnt++] = prismpt[2];
1551 conn[cnt++] = prismpt[4];
1552 conn[cnt++] = prismpt[3];
1553 conn[cnt++] = prismpt[0];
1554 conn[cnt++] = prismpt[2];
1556 conn[cnt++] = prismpt[3];
1557 conn[cnt++] = prismpt[4];
1558 conn[cnt++] = prismpt[5];
1559 conn[cnt++] = prismpt[2];
1565 plane += (np - i) * np;
1586 if (d0factors.size() != 5)
1593 if (d0factors[0].size() != nquad0 * nquad1)
1600 if (d0factors[1].size() != nquad0 * nquad2)
1610 if (d0factors[2].size() != nquad1 * nquad2)
1632 int ncoords = normal_0.size();
1638 for (
int i = 0; i < nquad0 * nquad1; ++i)
1640 d0factors[0][i] = df[0][i] * normal_0[0][i];
1641 d1factors[0][i] = df[1][i] * normal_0[0][i];
1642 d2factors[0][i] = df[2][i] * normal_0[0][i];
1645 for (
int n = 1; n < ncoords; ++n)
1647 for (
int i = 0; i < nquad0 * nquad1; ++i)
1649 d0factors[0][i] += df[3 * n][i] * normal_0[n][i];
1650 d1factors[0][i] += df[3 * n + 1][i] * normal_0[n][i];
1651 d2factors[0][i] += df[3 * n + 2][i] * normal_0[n][i];
1656 for (
int j = 0; j < nquad2; ++j)
1658 for (
int i = 0; i < nquad0; ++i)
1660 d0factors[1][i] = df[0][j * nquad0 * nquad1 + i] *
1661 normal_1[0][j * nquad0 + i];
1663 df[0][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1664 normal_3[0][j * nquad0 + i];
1665 d1factors[1][i] = df[1][j * nquad0 * nquad1 + i] *
1666 normal_1[0][j * nquad0 + i];
1668 df[1][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1669 normal_3[0][j * nquad0 + i];
1670 d2factors[1][i] = df[2][j * nquad0 * nquad1 + i] *
1671 normal_1[0][j * nquad0 + i];
1673 df[2][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1674 normal_3[0][j * nquad0 + i];
1678 for (
int n = 1; n < ncoords; ++n)
1680 for (
int j = 0; j < nquad2; ++j)
1682 for (
int i = 0; i < nquad0; ++i)
1684 d0factors[1][i] = df[3 * n][j * nquad0 * nquad1 + i] *
1685 normal_1[0][j * nquad0 + i];
1687 df[3 * n][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1688 normal_3[0][j * nquad0 + i];
1689 d1factors[1][i] = df[3 * n + 1][j * nquad0 * nquad1 + i] *
1690 normal_1[0][j * nquad0 + i];
1692 df[3 * n + 1][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1693 normal_3[0][j * nquad0 + i];
1694 d2factors[1][i] = df[3 * n + 2][j * nquad0 * nquad1 + i] *
1695 normal_1[0][j * nquad0 + i];
1697 df[3 * n + 2][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1698 normal_3[0][j * nquad0 + i];
1704 for (
int j = 0; j < nquad2; ++j)
1706 for (
int i = 0; i < nquad1; ++i)
1708 d0factors[2][j * nquad1 + i] =
1709 df[0][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1710 normal_2[0][j * nquad1 + i];
1711 d0factors[4][j * nquad1 + i] =
1712 df[0][j * nquad0 * nquad1 + i * nquad0] *
1713 normal_4[0][j * nquad1 + i];
1714 d1factors[2][j * nquad1 + i] =
1715 df[1][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1716 normal_2[0][j * nquad1 + i];
1717 d1factors[4][j * nquad1 + i] =
1718 df[1][j * nquad0 * nquad1 + i * nquad0] *
1719 normal_4[0][j * nquad1 + i];
1720 d2factors[2][j * nquad1 + i] =
1721 df[2][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1722 normal_2[0][j * nquad1 + i];
1723 d2factors[4][j * nquad1 + i] =
1724 df[2][j * nquad0 * nquad1 + i * nquad0] *
1725 normal_4[0][j * nquad1 + i];
1729 for (
int n = 1; n < ncoords; ++n)
1731 for (
int j = 0; j < nquad2; ++j)
1733 for (
int i = 0; i < nquad1; ++i)
1735 d0factors[2][j * nquad1 + i] +=
1736 df[3 * n][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1737 normal_2[n][j * nquad0 + i];
1738 d0factors[4][j * nquad0 + i] +=
1739 df[3 * n][i * nquad0 + j * nquad0 * nquad1] *
1740 normal_4[n][j * nquad0 + i];
1741 d1factors[2][j * nquad1 + i] +=
1743 [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1744 normal_2[n][j * nquad0 + i];
1745 d1factors[4][j * nquad0 + i] +=
1746 df[3 * n + 1][i * nquad0 + j * nquad0 * nquad1] *
1747 normal_4[n][j * nquad0 + i];
1748 d2factors[2][j * nquad1 + i] +=
1750 [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1751 normal_2[n][j * nquad0 + i];
1752 d2factors[4][j * nquad0 + i] +=
1753 df[3 * n + 2][i * nquad0 + j * nquad0 * nquad1] *
1754 normal_4[n][j * nquad0 + i];
1762 for (
int i = 0; i < nquad0 * nquad1; ++i)
1764 d0factors[0][i] = df[0][0] * normal_0[0][i];
1765 d1factors[0][i] = df[1][0] * normal_0[0][i];
1766 d2factors[0][i] = df[2][0] * normal_0[0][i];
1769 for (
int n = 1; n < ncoords; ++n)
1771 for (
int i = 0; i < nquad0 * nquad1; ++i)
1773 d0factors[0][i] += df[3 * n][0] * normal_0[n][i];
1774 d1factors[0][i] += df[3 * n + 1][0] * normal_0[n][i];
1775 d2factors[0][i] += df[3 * n + 2][0] * normal_0[n][i];
1780 for (
int i = 0; i < nquad0 * nquad2; ++i)
1782 d0factors[1][i] = df[0][0] * normal_1[0][i];
1783 d0factors[3][i] = df[0][0] * normal_3[0][i];
1785 d1factors[1][i] = df[1][0] * normal_1[0][i];
1786 d1factors[3][i] = df[1][0] * normal_3[0][i];
1788 d2factors[1][i] = df[2][0] * normal_1[0][i];
1789 d2factors[3][i] = df[2][0] * normal_3[0][i];
1792 for (
int n = 1; n < ncoords; ++n)
1794 for (
int i = 0; i < nquad0 * nquad2; ++i)
1796 d0factors[1][i] += df[3 * n][0] * normal_1[n][i];
1797 d0factors[3][i] += df[3 * n][0] * normal_3[n][i];
1799 d1factors[1][i] += df[3 * n + 1][0] * normal_1[n][i];
1800 d1factors[3][i] += df[3 * n + 1][0] * normal_3[n][i];
1802 d2factors[1][i] += df[3 * n + 2][0] * normal_1[n][i];
1803 d2factors[3][i] += df[3 * n + 2][0] * normal_3[n][i];
1808 for (
int i = 0; i < nquad1 * nquad2; ++i)
1810 d0factors[2][i] = df[0][0] * normal_2[0][i];
1811 d0factors[4][i] = df[0][0] * normal_4[0][i];
1813 d1factors[2][i] = df[1][0] * normal_2[0][i];
1814 d1factors[4][i] = df[1][0] * normal_4[0][i];
1816 d2factors[2][i] = df[2][0] * normal_2[0][i];
1817 d2factors[4][i] = df[2][0] * normal_4[0][i];
1820 for (
int n = 1; n < ncoords; ++n)
1822 for (
int i = 0; i < nquad1 * nquad2; ++i)
1824 d0factors[2][i] += df[3 * n][0] * normal_2[n][i];
1825 d0factors[4][i] += df[3 * n][0] * normal_4[n][i];
1827 d1factors[2][i] += df[3 * n + 1][0] * normal_2[n][i];
1828 d1factors[4][i] += df[3 * n + 1][0] * normal_4[n][i];
1830 d2factors[2][i] += df[3 * n + 2][0] * normal_2[n][i];
1831 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.
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
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
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
const NormalVector & GetTraceNormal(const int id)
void v_ComputeTraceNormal(const int face)
Get the normals along specficied face Get the face normals interplated to a points0 x points 0 type d...
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
virtual void v_NormalTraceDerivFactors(Array< OneD, Array< OneD, NekDouble >> &d0factors, Array< OneD, Array< OneD, NekDouble >> &d1factors, Array< OneD, Array< OneD, NekDouble >> &d2factors)
: This method gets all of the factors which are required as part of the Gradient Jump Penalty stabili...
virtual void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey)
virtual void v_GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
virtual void v_ExtractDataToCoeffs(const NekDouble *data, const std::vector< unsigned int > &nummodes, const int mode_offset, NekDouble *coeffs, std::vector< LibUtilities::BasisType > &fromType)
virtual void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true)
virtual NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coord, const Array< OneD, const NekDouble > &physvals)
This function evaluates the expansion at a single (arbitrary) point of the domain.
virtual NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
virtual void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp)
Calculate the Laplacian multiplication in a matrix-free manner.
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Calculate the inner product of inarray with respect to the basis B=base0*base1*base2 and put into out...
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Calculates the inner product .
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
Calculate the derivative of the physical points.
virtual void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Integrate the physical point list inarray over prismatic region and return the value.
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.
virtual int v_GetCoordim()
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
virtual void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
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)
virtual void v_GetTracePhysMap(const int face, Array< OneD, int > &outarray)
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
virtual void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords)
Get the coordinates #coords at the local coordinates #Lcoords.
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
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)
virtual void v_HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
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
const LibUtilities::BasisKey GetTraceBasisKey(const int i, int k=-1) const
This function returns the basis key belonging to the i-th trace.
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)
Array< OneD, LibUtilities::BasisSharedPtr > m_base
MatrixType GetMatrixType() const
static void Dgemv(const char &trans, const int &m, const int &n, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = A x where A[m x n].
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
static 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
The above copyright notice and this permission notice shall be included.
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)
vvtvvtp (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/y.
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)