42 namespace LocalRegions
50 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
53 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
57 std::bind(&
Expansion3D::CreateMatrix, this, std::placeholders::_1),
58 std::string(
"PyrExpMatrix")),
59 m_staticCondMatrixManager(std::bind(&
Expansion::CreateStaticCondMatrix,
60 this, std::placeholders::_1),
61 std::string(
"PyrExpStaticCondMatrix"))
66 : StdExpansion(T), StdExpansion3D(T), StdPyrExp(T),
Expansion(T),
68 m_staticCondMatrixManager(T.m_staticCondMatrixManager)
102 int nquad0 =
m_base[0]->GetNumPoints();
103 int nquad1 =
m_base[1]->GetNumPoints();
104 int nquad2 =
m_base[2]->GetNumPoints();
112 (
NekDouble *)&inarray[0], 1, &tmp[0], 1);
117 (
NekDouble *)&inarray[0], 1, &tmp[0], 1);
121 return StdPyrExp::v_Integral(tmp);
133 int nquad0 =
m_base[0]->GetNumPoints();
134 int nquad1 =
m_base[1]->GetNumPoints();
135 int nquad2 =
m_base[2]->GetNumPoints();
142 StdPyrExp::v_PhysDeriv(inarray, diff0, diff1, diff2);
148 Vmath::Vmul(nquad0 * nquad1 * nquad2, &gmat[0][0], 1, &diff0[0], 1,
150 Vmath::Vvtvp(nquad0 * nquad1 * nquad2, &gmat[1][0], 1, &diff1[0], 1,
151 &out_d0[0], 1, &out_d0[0], 1);
152 Vmath::Vvtvp(nquad0 * nquad1 * nquad2, &gmat[2][0], 1, &diff2[0], 1,
153 &out_d0[0], 1, &out_d0[0], 1);
158 Vmath::Vmul(nquad0 * nquad1 * nquad2, &gmat[3][0], 1, &diff0[0], 1,
160 Vmath::Vvtvp(nquad0 * nquad1 * nquad2, &gmat[4][0], 1, &diff1[0], 1,
161 &out_d1[0], 1, &out_d1[0], 1);
162 Vmath::Vvtvp(nquad0 * nquad1 * nquad2, &gmat[5][0], 1, &diff2[0], 1,
163 &out_d1[0], 1, &out_d1[0], 1);
168 Vmath::Vmul(nquad0 * nquad1 * nquad2, &gmat[6][0], 1, &diff0[0], 1,
170 Vmath::Vvtvp(nquad0 * nquad1 * nquad2, &gmat[7][0], 1, &diff1[0], 1,
171 &out_d2[0], 1, &out_d2[0], 1);
172 Vmath::Vvtvp(nquad0 * nquad1 * nquad2, &gmat[8][0], 1, &diff2[0], 1,
173 &out_d2[0], 1, &out_d2[0], 1);
180 Vmath::Smul(nquad0 * nquad1 * nquad2, gmat[0][0], &diff0[0], 1,
182 Blas::Daxpy(nquad0 * nquad1 * nquad2, gmat[1][0], &diff1[0], 1,
184 Blas::Daxpy(nquad0 * nquad1 * nquad2, gmat[2][0], &diff2[0], 1,
190 Vmath::Smul(nquad0 * nquad1 * nquad2, gmat[3][0], &diff0[0], 1,
192 Blas::Daxpy(nquad0 * nquad1 * nquad2, gmat[4][0], &diff1[0], 1,
194 Blas::Daxpy(nquad0 * nquad1 * nquad2, gmat[5][0], &diff2[0], 1,
200 Vmath::Smul(nquad0 * nquad1 * nquad2, gmat[6][0], &diff0[0], 1,
202 Blas::Daxpy(nquad0 * nquad1 * nquad2, gmat[7][0], &diff1[0], 1,
204 Blas::Daxpy(nquad0 * nquad1 * nquad2, gmat[8][0], &diff2[0], 1,
230 if (
m_base[0]->Collocation() &&
m_base[1]->Collocation() &&
247 out = (*matsys) * in;
292 const int nquad0 =
m_base[0]->GetNumPoints();
293 const int nquad1 =
m_base[1]->GetNumPoints();
294 const int nquad2 =
m_base[2]->GetNumPoints();
295 const int order0 =
m_base[0]->GetNumModes();
296 const int order1 =
m_base[1]->GetNumModes();
300 if (multiplybyweights)
308 tmp, outarray, wsp,
true,
true,
true);
314 inarray, outarray, wsp,
true,
true,
true);
359 const int nquad0 =
m_base[0]->GetNumPoints();
360 const int nquad1 =
m_base[1]->GetNumPoints();
361 const int nquad2 =
m_base[2]->GetNumPoints();
362 const int order0 =
m_base[0]->GetNumModes();
363 const int order1 =
m_base[1]->GetNumModes();
364 const int nqtot = nquad0 * nquad1 * nquad2;
372 std::max(nqtot, order0 * nquad2 * (nquad1 + order1)));
384 m_base[2]->GetBdata(), tmp2, outarray, wsp,
388 m_base[2]->GetBdata(), tmp3, tmp6, wsp,
true,
394 m_base[2]->GetDbdata(), tmp4, tmp6, wsp,
true,
404 const int nquad0 =
m_base[0]->GetNumPoints();
405 const int nquad1 =
m_base[1]->GetNumPoints();
406 const int nquad2 =
m_base[2]->GetNumPoints();
407 const int order0 =
m_base[0]->GetNumModes();
408 const int order1 =
m_base[1]->GetNumModes();
409 const int nqtot = nquad0 * nquad1 * nquad2;
420 std::max(nqtot, order0 * nquad2 * (nquad1 + order1)));
434 Vmath::Vmul(nqtot, &df[3 * dir][0], 1, tmp1.get(), 1, tmp2.get(), 1);
435 Vmath::Vmul(nqtot, &df[3 * dir + 1][0], 1, tmp1.get(), 1, tmp3.get(),
437 Vmath::Vmul(nqtot, &df[3 * dir + 2][0], 1, tmp1.get(), 1, tmp4.get(),
442 Vmath::Smul(nqtot, df[3 * dir][0], tmp1.get(), 1, tmp2.get(), 1);
443 Vmath::Smul(nqtot, df[3 * dir + 1][0], tmp1.get(), 1, tmp3.get(), 1);
444 Vmath::Smul(nqtot, df[3 * dir + 2][0], tmp1.get(), 1, tmp4.get(), 1);
448 for (
int i = 0; i < nquad0; ++i)
450 gfac0[i] = 0.5 * (1 + z0[i]);
454 for (
int i = 0; i < nquad1; ++i)
456 gfac1[i] = 0.5 * (1 + z1[i]);
460 for (
int i = 0; i < nquad2; ++i)
462 gfac2[i] = 2.0 / (1 - z2[i]);
465 const int nq01 = nquad0 * nquad1;
467 for (
int i = 0; i < nquad2; ++i)
469 Vmath::Smul(nq01, gfac2[i], &tmp2[0] + i * nq01, 1, &tmp2[0] + i * nq01,
471 Vmath::Smul(nq01, gfac2[i], &tmp3[0] + i * nq01, 1, &tmp3[0] + i * nq01,
473 Vmath::Smul(nq01, gfac2[i], &tmp4[0] + i * nq01, 1, &tmp5[0] + i * nq01,
478 for (
int i = 0; i < nquad1 * nquad2; ++i)
480 Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp5[0] + i * nquad0, 1,
481 &wsp[0] + i * nquad0, 1);
484 Vmath::Vadd(nqtot, &tmp2[0], 1, &wsp[0], 1, &tmp2[0], 1);
487 for (
int i = 0; i < nquad1 * nquad2; ++i)
489 Vmath::Smul(nquad0, gfac1[i % nquad1], &tmp5[0] + i * nquad0, 1,
490 &tmp5[0] + i * nquad0, 1);
492 Vmath::Vadd(nqtot, &tmp3[0], 1, &tmp5[0], 1, &tmp3[0], 1);
503 m_base[2]->GetBasisKey());
509 m_base[0]->GetPointsKey());
511 m_base[1]->GetPointsKey());
513 m_base[2]->GetPointsKey());
528 ASSERTL1(Lcoords[0] <= -1.0 && Lcoords[0] >= 1.0 && Lcoords[1] <= -1.0 &&
529 Lcoords[1] >= 1.0 && Lcoords[2] <= -1.0 && Lcoords[2] >= 1.0,
530 "Local coordinates are not in region [-1,1]");
534 for (i = 0; i <
m_geom->GetCoordim(); ++i)
536 coords[i] =
m_geom->GetCoord(i, Lcoords);
548 const NekDouble *data,
const std::vector<unsigned int> &nummodes,
549 const int mode_offset,
NekDouble *coeffs,
550 std::vector<LibUtilities::BasisType> &fromType)
552 int data_order0 = nummodes[mode_offset];
553 int fillorder0 = min(
m_base[0]->GetNumModes(), data_order0);
554 int data_order1 = nummodes[mode_offset + 1];
555 int order1 =
m_base[1]->GetNumModes();
556 int fillorder1 = min(order1, data_order1);
557 int data_order2 = nummodes[mode_offset + 2];
558 int order2 =
m_base[2]->GetNumModes();
559 int fillorder2 = min(order2, data_order2);
566 data_order1 != fillorder1 || data_order2 != fillorder2)
573 m_base[0]->GetPointsKey()),
575 m_base[1]->GetPointsKey()),
577 m_base[2]->GetPointsKey()));
581 m_base[2]->GetBasisKey());
607 return StdPyrExp::v_PhysEvaluate(Lcoord, physvals);
618 m_geom->GetLocCoords(coord, Lcoord);
620 return StdPyrExp::v_PhysEvaluate(Lcoord, physvals);
629 return m_geom->GetCoordim();
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);
1100 const unsigned int dim = 3;
1106 for (
unsigned int i = 0; i < dim; ++i)
1108 for (
unsigned int j = i; j < dim; ++j)
1139 const unsigned int nquad0 =
m_base[0]->GetNumPoints();
1140 const unsigned int nquad1 =
m_base[1]->GetNumPoints();
1141 const unsigned int nquad2 =
m_base[2]->GetNumPoints();
1144 for (j = 0; j < nquad2; ++j)
1146 for (i = 0; i < nquad1; ++i)
1149 &h0[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1151 &h1[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1152 Vmath::Fill(nquad0, (1.0 + z1[i]) / (1.0 - z2[j]),
1153 &h2[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1156 for (i = 0; i < nquad0; i++)
1158 Blas::Dscal(nquad1 * nquad2, 1 + z0[i], &h1[0] + i, nquad0);
1167 Vmath::Vvtvvtp(nqtot, &df[0][0], 1, &h0[0], 1, &df[2][0], 1, &h1[0], 1,
1169 Vmath::Vvtvvtp(nqtot, &df[3][0], 1, &h0[0], 1, &df[5][0], 1, &h1[0], 1,
1171 Vmath::Vvtvvtp(nqtot, &df[6][0], 1, &h0[0], 1, &df[8][0], 1, &h1[0], 1,
1175 Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp1[0], 1, &wsp2[0], 1, &wsp2[0],
1177 Vmath::Vvtvp(nqtot, &wsp3[0], 1, &wsp3[0], 1, &g0[0], 1, &g0[0], 1);
1180 Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp1[0], 1, &df[5][0], 1, &wsp2[0],
1182 Vmath::Vvtvp(nqtot, &df[8][0], 1, &wsp3[0], 1, &g4[0], 1, &g4[0], 1);
1185 Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &h0[0], 1, &df[2][0], 1, &h2[0], 1,
1187 Vmath::Vvtvvtp(nqtot, &df[4][0], 1, &h0[0], 1, &df[5][0], 1, &h2[0], 1,
1189 Vmath::Vvtvvtp(nqtot, &df[7][0], 1, &h0[0], 1, &df[8][0], 1, &h2[0], 1,
1193 Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0],
1195 Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp6[0], 1, &g1[0], 1, &g1[0], 1);
1198 Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp4[0], 1, &wsp2[0], 1, &wsp5[0],
1200 Vmath::Vvtvp(nqtot, &wsp3[0], 1, &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1203 Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp4[0], 1, &df[5][0], 1, &wsp5[0],
1205 Vmath::Vvtvp(nqtot, &df[8][0], 1, &wsp6[0], 1, &g5[0], 1, &g5[0], 1);
1209 &df[5][0], 1, &g2[0], 1);
1210 Vmath::Vvtvp(nqtot, &df[8][0], 1, &df[8][0], 1, &g2[0], 1, &g2[0], 1);
1223 Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp1[0], 1, &wsp2[0], 1, &wsp2[0],
1225 Vmath::Vvtvp(nqtot, &wsp3[0], 1, &wsp3[0], 1, &g0[0], 1, &g0[0], 1);
1228 Vmath::Svtsvtp(nqtot, df[2][0], &wsp1[0], 1, df[5][0], &wsp2[0], 1,
1230 Vmath::Svtvp(nqtot, df[8][0], &wsp3[0], 1, &g4[0], 1, &g4[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, &g1[0], 1, &g1[0], 1);
1246 Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp4[0], 1, &wsp2[0], 1, &wsp5[0],
1248 Vmath::Vvtvp(nqtot, &wsp3[0], 1, &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, &g5[0], 1, &g5[0], 1);
1257 df[2][0] * df[2][0] + df[5][0] * df[5][0] +
1258 df[8][0] * df[8][0],
1262 for (
unsigned int i = 0; i < dim; ++i)
1264 for (
unsigned int j = i; j < dim; ++j)
1282 int nquad0 =
m_base[0]->GetNumPoints();
1283 int nquad1 =
m_base[1]->GetNumPoints();
1284 int nq2 =
m_base[2]->GetNumPoints();
1285 int nqtot = nquad0 * nquad1 * nq2;
1287 ASSERTL1(wsp.size() >= 6 * nqtot,
"Insufficient workspace size.");
1288 ASSERTL1(m_ncoeffs <= nqtot, "Workspace not set up for ncoeffs > nqtot
");
1290 const Array<OneD, const NekDouble> &base0 = m_base[0]->GetBdata();
1291 const Array<OneD, const NekDouble> &base1 = m_base[1]->GetBdata();
1292 const Array<OneD, const NekDouble> &base2 = m_base[2]->GetBdata();
1293 const Array<OneD, const NekDouble> &dbase0 = m_base[0]->GetDbdata();
1294 const Array<OneD, const NekDouble> &dbase1 = m_base[1]->GetDbdata();
1295 const Array<OneD, const NekDouble> &dbase2 = m_base[2]->GetDbdata();
1296 const Array<OneD, const NekDouble> &metric00 =
1297 m_metrics[eMetricLaplacian00];
1298 const Array<OneD, const NekDouble> &metric01 =
1299 m_metrics[eMetricLaplacian01];
1300 const Array<OneD, const NekDouble> &metric02 =
1301 m_metrics[eMetricLaplacian02];
1302 const Array<OneD, const NekDouble> &metric11 =
1303 m_metrics[eMetricLaplacian11];
1304 const Array<OneD, const NekDouble> &metric12 =
1305 m_metrics[eMetricLaplacian12];
1306 const Array<OneD, const NekDouble> &metric22 =
1307 m_metrics[eMetricLaplacian22];
1309 // Allocate temporary storage
1310 Array<OneD, NekDouble> wsp0(2 * nqtot, wsp);
1311 Array<OneD, NekDouble> wsp1(nqtot, wsp + 1 * nqtot);
1312 Array<OneD, NekDouble> wsp2(nqtot, wsp + 2 * nqtot);
1313 Array<OneD, NekDouble> wsp3(nqtot, wsp + 3 * nqtot);
1314 Array<OneD, NekDouble> wsp4(nqtot, wsp + 4 * nqtot);
1315 Array<OneD, NekDouble> wsp5(nqtot, wsp + 5 * nqtot);
1317 // LAPLACIAN MATRIX OPERATION
1318 // wsp1 = du_dxi1 = D_xi1 * inarray = D_xi1 * u
1319 // wsp2 = du_dxi2 = D_xi2 * inarray = D_xi2 * u
1320 // wsp2 = du_dxi3 = D_xi3 * inarray = D_xi3 * u
1321 StdExpansion3D::PhysTensorDeriv(inarray, wsp0, wsp1, wsp2);
1323 // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1324 // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1325 // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
1326 // especially for this purpose
1327 Vmath::Vvtvvtp(nqtot, &metric00[0], 1, &wsp0[0], 1, &metric01[0], 1,
1328 &wsp1[0], 1, &wsp3[0], 1);
1329 Vmath::Vvtvp(nqtot, &metric02[0], 1, &wsp2[0], 1, &wsp3[0], 1, &wsp3[0], 1);
1330 Vmath::Vvtvvtp(nqtot, &metric01[0], 1, &wsp0[0], 1, &metric11[0], 1,
1331 &wsp1[0], 1, &wsp4[0], 1);
1332 Vmath::Vvtvp(nqtot, &metric12[0], 1, &wsp2[0], 1, &wsp4[0], 1, &wsp4[0], 1);
1333 Vmath::Vvtvvtp(nqtot, &metric02[0], 1, &wsp0[0], 1, &metric12[0], 1,
1334 &wsp1[0], 1, &wsp5[0], 1);
1335 Vmath::Vvtvp(nqtot, &metric22[0], 1, &wsp2[0], 1, &wsp5[0], 1, &wsp5[0], 1);
1337 // outarray = m = (D_xi1 * B)^T * k
1338 // wsp1 = n = (D_xi2 * B)^T * l
1339 IProductWRTBase_SumFacKernel(dbase0, base1, base2, wsp3, outarray, wsp0,
1341 IProductWRTBase_SumFacKernel(base0, dbase1, base2, wsp4, wsp2, wsp0, true,
1343 Vmath::Vadd(m_ncoeffs, wsp2.get(), 1, outarray.get(), 1, outarray.get(), 1);
1344 IProductWRTBase_SumFacKernel(base0, base1, dbase2, wsp5, wsp2, wsp0, true,
1346 Vmath::Vadd(m_ncoeffs, wsp2.get(), 1, outarray.get(), 1, outarray.get(), 1);
1354 void PyrExp::v_NormalTraceDerivFactors(
1355 Array<OneD, Array<OneD, NekDouble>> &d0factors,
1356 Array<OneD, Array<OneD, NekDouble>> &d1factors,
1357 Array<OneD, Array<OneD, NekDouble>> &d2factors)
1359 int nquad0 = GetNumPoints(0);
1360 int nquad1 = GetNumPoints(1);
1361 int nquad2 = GetNumPoints(2);
1363 const Array<TwoD, const NekDouble> &df =
1364 m_metricinfo->GetDerivFactors(GetPointsKeys());
1366 if (d0factors.size() != 5)
1368 d0factors = Array<OneD, Array<OneD, NekDouble>>(5);
1369 d1factors = Array<OneD, Array<OneD, NekDouble>>(5);
1370 d2factors = Array<OneD, Array<OneD, NekDouble>>(5);
1373 if (d0factors[0].size() != nquad0 * nquad1)
1375 d0factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1376 d1factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1377 d2factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1380 if (d0factors[1].size() != nquad0 * nquad2)
1382 d0factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1383 d0factors[3] = Array<OneD, NekDouble>(nquad0 * nquad2);
1384 d1factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1385 d1factors[3] = Array<OneD, NekDouble>(nquad0 * nquad2);
1386 d2factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1387 d2factors[3] = Array<OneD, NekDouble>(nquad0 * nquad2);
1390 if (d0factors[2].size() != nquad1 * nquad2)
1392 d0factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1393 d0factors[4] = Array<OneD, NekDouble>(nquad1 * nquad2);
1394 d1factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1395 d1factors[4] = Array<OneD, NekDouble>(nquad1 * nquad2);
1396 d2factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1397 d2factors[4] = Array<OneD, NekDouble>(nquad1 * nquad2);
1401 const Array<OneD, const Array<OneD, NekDouble>> &normal_0 =
1403 const Array<OneD, const Array<OneD, NekDouble>> &normal_1 =
1405 const Array<OneD, const Array<OneD, NekDouble>> &normal_2 =
1407 const Array<OneD, const Array<OneD, NekDouble>> &normal_3 =
1409 const Array<OneD, const Array<OneD, NekDouble>> &normal_4 =
1412 int ncoords = normal_0.size();
1414 // first gather together standard cartesian inner products
1415 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1418 for (int i = 0; i < nquad0 * nquad1; ++i)
1420 d0factors[0][i] = df[0][i] * normal_0[0][i];
1421 d1factors[0][i] = df[1][i] * normal_0[0][i];
1422 d2factors[0][i] = df[2][i] * normal_0[0][i];
1425 for (int n = 1; n < ncoords; ++n)
1427 for (int i = 0; i < nquad0 * nquad1; ++i)
1429 d0factors[0][i] += df[3 * n][i] * normal_0[n][i];
1430 d1factors[0][i] += df[3 * n + 1][i] * normal_0[n][i];
1431 d2factors[0][i] += df[3 * n + 2][i] * normal_0[n][i];
1436 for (int j = 0; j < nquad2; ++j)
1438 for (int i = 0; i < nquad0; ++i)
1440 d0factors[1][i] = df[0][j * nquad0 * nquad1 + i] *
1441 normal_1[0][j * nquad0 + i];
1443 df[0][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1444 normal_3[0][j * nquad0 + i];
1445 d1factors[1][i] = df[1][j * nquad0 * nquad1 + i] *
1446 normal_1[0][j * nquad0 + i];
1448 df[1][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1449 normal_3[0][j * nquad0 + i];
1450 d2factors[1][i] = df[2][j * nquad0 * nquad1 + i] *
1451 normal_1[0][j * nquad0 + i];
1453 df[2][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1454 normal_3[0][j * nquad0 + i];
1458 for (int n = 1; n < ncoords; ++n)
1460 for (int j = 0; j < nquad2; ++j)
1462 for (int i = 0; i < nquad0; ++i)
1464 d0factors[1][i] = df[3 * n][j * nquad0 * nquad1 + i] *
1465 normal_1[0][j * nquad0 + i];
1467 df[3 * n][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1468 normal_3[0][j * nquad0 + i];
1469 d1factors[1][i] = df[3 * n + 1][j * nquad0 * nquad1 + i] *
1470 normal_1[0][j * nquad0 + i];
1472 df[3 * n + 1][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1473 normal_3[0][j * nquad0 + i];
1474 d2factors[1][i] = df[3 * n + 2][j * nquad0 * nquad1 + i] *
1475 normal_1[0][j * nquad0 + i];
1477 df[3 * n + 2][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1478 normal_3[0][j * nquad0 + i];
1484 for (int j = 0; j < nquad2; ++j)
1486 for (int i = 0; i < nquad1; ++i)
1488 d0factors[2][j * nquad1 + i] =
1489 df[0][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1490 normal_2[0][j * nquad1 + i];
1491 d0factors[4][j * nquad1 + i] =
1492 df[0][j * nquad0 * nquad1 + i * nquad0] *
1493 normal_4[0][j * nquad1 + i];
1494 d1factors[2][j * nquad1 + i] =
1495 df[1][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1496 normal_2[0][j * nquad1 + i];
1497 d1factors[4][j * nquad1 + i] =
1498 df[1][j * nquad0 * nquad1 + i * nquad0] *
1499 normal_4[0][j * nquad1 + i];
1500 d2factors[2][j * nquad1 + i] =
1501 df[2][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1502 normal_2[0][j * nquad1 + i];
1503 d2factors[4][j * nquad1 + i] =
1504 df[2][j * nquad0 * nquad1 + i * nquad0] *
1505 normal_4[0][j * nquad1 + i];
1509 for (int n = 1; n < ncoords; ++n)
1511 for (int j = 0; j < nquad2; ++j)
1513 for (int i = 0; i < nquad1; ++i)
1515 d0factors[2][j * nquad1 + i] +=
1516 df[3 * n][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1517 normal_2[n][j * nquad0 + i];
1518 d0factors[4][j * nquad0 + i] +=
1519 df[3 * n][i * nquad0 + j * nquad0 * nquad1] *
1520 normal_4[n][j * nquad0 + i];
1521 d1factors[2][j * nquad1 + i] +=
1523 [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1524 normal_2[n][j * nquad0 + i];
1525 d1factors[4][j * nquad0 + i] +=
1526 df[3 * n + 1][i * nquad0 + j * nquad0 * nquad1] *
1527 normal_4[n][j * nquad0 + i];
1528 d2factors[2][j * nquad1 + i] +=
1530 [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1531 normal_2[n][j * nquad0 + i];
1532 d2factors[4][j * nquad0 + i] +=
1533 df[3 * n + 2][i * nquad0 + j * nquad0 * nquad1] *
1534 normal_4[n][j * nquad0 + i];
1542 for (int i = 0; i < nquad0 * nquad1; ++i)
1544 d0factors[0][i] = df[0][0] * normal_0[0][i];
1545 d1factors[0][i] = df[1][0] * normal_0[0][i];
1546 d2factors[0][i] = df[2][0] * normal_0[0][i];
1549 for (int n = 1; n < ncoords; ++n)
1551 for (int i = 0; i < nquad0 * nquad1; ++i)
1553 d0factors[0][i] += df[3 * n][0] * normal_0[n][i];
1554 d1factors[0][i] += df[3 * n + 1][0] * normal_0[n][i];
1555 d2factors[0][i] += df[3 * n + 2][0] * normal_0[n][i];
1560 for (int i = 0; i < nquad0 * nquad2; ++i)
1562 d0factors[1][i] = df[0][0] * normal_1[0][i];
1563 d0factors[3][i] = df[0][0] * normal_3[0][i];
1565 d1factors[1][i] = df[1][0] * normal_1[0][i];
1566 d1factors[3][i] = df[1][0] * normal_3[0][i];
1568 d2factors[1][i] = df[2][0] * normal_1[0][i];
1569 d2factors[3][i] = df[2][0] * normal_3[0][i];
1572 for (int n = 1; n < ncoords; ++n)
1574 for (int i = 0; i < nquad0 * nquad2; ++i)
1576 d0factors[1][i] += df[3 * n][0] * normal_1[n][i];
1577 d0factors[3][i] += df[3 * n][0] * normal_3[n][i];
1579 d1factors[1][i] += df[3 * n + 1][0] * normal_1[n][i];
1580 d1factors[3][i] += df[3 * n + 1][0] * normal_3[n][i];
1582 d2factors[1][i] += df[3 * n + 2][0] * normal_1[n][i];
1583 d2factors[3][i] += df[3 * n + 2][0] * normal_3[n][i];
1588 for (int i = 0; i < nquad1 * nquad2; ++i)
1590 d0factors[2][i] = df[0][0] * normal_2[0][i];
1591 d0factors[4][i] = df[0][0] * normal_4[0][i];
1593 d1factors[2][i] = df[1][0] * normal_2[0][i];
1594 d1factors[4][i] = df[1][0] * normal_4[0][i];
1596 d2factors[2][i] = df[2][0] * normal_2[0][i];
1597 d2factors[4][i] = df[2][0] * normal_4[0][i];
1600 for (int n = 1; n < ncoords; ++n)
1602 for (int i = 0; i < nquad1 * nquad2; ++i)
1604 d0factors[2][i] += df[3 * n][0] * normal_2[n][i];
1605 d0factors[4][i] += df[3 * n][0] * normal_4[n][i];
1607 d1factors[2][i] += df[3 * n + 1][0] * normal_2[n][i];
1608 d1factors[4][i] += df[3 * n + 1][0] * normal_4[n][i];
1610 d2factors[2][i] += df[3 * n + 2][0] * normal_2[n][i];
1611 d2factors[4][i] += df[3 * n + 2][0] * normal_4[n][i];
1617 } // namespace LocalRegions
1618 } // namespace Nektar
#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
void ComputeLaplacianMetric()
void ComputeQuadratureMetric()
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
virtual int v_GetCoordim()
virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey)
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
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_GetTracePhysMap(const int face, Array< OneD, int > &outarray)
virtual DNekMatSharedPtr v_CreateStdMatrix(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...
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
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.
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
void v_ComputeTraceNormal(const int face)
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
virtual void v_ComputeLaplacianMetric()
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
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.
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
virtual void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords)
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Integrate the physical point list inarray over pyramidic region and return the value.
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...
virtual void v_ExtractDataToCoeffs(const NekDouble *data, const std::vector< unsigned int > &nummodes, const int mode_offset, NekDouble *coeffs, std::vector< LibUtilities::BasisType > &fromType)
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Calculates the inner product .
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
virtual void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp)
virtual void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
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
const LibUtilities::BasisKey GetTraceBasisKey(const int i, int k=-1) const
This function returns the basis key belonging to the i-th trace.
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)
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
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)