44 namespace LocalRegions
63 : StdExpansion(Ba.GetNumModes() * Bb.GetNumModes() * Bc.GetNumModes(), 3,
65 StdExpansion3D(Ba.GetNumModes() * Bb.GetNumModes() * Bc.GetNumModes(), Ba,
69 std::bind(&
Expansion3D::CreateMatrix, this, std::placeholders::_1),
70 std::string(
"HexExpMatrix")),
71 m_staticCondMatrixManager(std::bind(&
Expansion::CreateStaticCondMatrix,
72 this, std::placeholders::_1),
73 std::string(
"HexExpStaticCondMatrix"))
83 : StdExpansion(T), StdExpansion3D(T), StdHexExp(T),
Expansion(T),
85 m_staticCondMatrixManager(T.m_staticCondMatrixManager)
112 int nquad0 =
m_base[0]->GetNumPoints();
113 int nquad1 =
m_base[1]->GetNumPoints();
114 int nquad2 =
m_base[2]->GetNumPoints();
124 (
NekDouble *)&inarray[0], 1, &tmp[0], 1);
129 (
NekDouble *)&inarray[0], 1, &tmp[0], 1);
133 returnVal = StdHexExp::v_Integral(tmp);
156 int nquad0 =
m_base[0]->GetNumPoints();
157 int nquad1 =
m_base[1]->GetNumPoints();
158 int nquad2 =
m_base[2]->GetNumPoints();
159 int ntot = nquad0 * nquad1 * nquad2;
167 StdHexExp::v_PhysDeriv(inarray, Diff0, Diff1, Diff2);
173 Vmath::Vmul(ntot, &df[0][0], 1, &Diff0[0], 1, &out_d0[0], 1);
174 Vmath::Vvtvp(ntot, &df[1][0], 1, &Diff1[0], 1, &out_d0[0], 1,
176 Vmath::Vvtvp(ntot, &df[2][0], 1, &Diff2[0], 1, &out_d0[0], 1,
182 Vmath::Vmul(ntot, &df[3][0], 1, &Diff0[0], 1, &out_d1[0], 1);
183 Vmath::Vvtvp(ntot, &df[4][0], 1, &Diff1[0], 1, &out_d1[0], 1,
185 Vmath::Vvtvp(ntot, &df[5][0], 1, &Diff2[0], 1, &out_d1[0], 1,
191 Vmath::Vmul(ntot, &df[6][0], 1, &Diff0[0], 1, &out_d2[0], 1);
192 Vmath::Vvtvp(ntot, &df[7][0], 1, &Diff1[0], 1, &out_d2[0], 1,
194 Vmath::Vvtvp(ntot, &df[8][0], 1, &Diff2[0], 1, &out_d2[0], 1,
202 Vmath::Smul(ntot, df[0][0], &Diff0[0], 1, &out_d0[0], 1);
203 Blas::Daxpy(ntot, df[1][0], &Diff1[0], 1, &out_d0[0], 1);
204 Blas::Daxpy(ntot, df[2][0], &Diff2[0], 1, &out_d0[0], 1);
209 Vmath::Smul(ntot, df[3][0], &Diff0[0], 1, &out_d1[0], 1);
210 Blas::Daxpy(ntot, df[4][0], &Diff1[0], 1, &out_d1[0], 1);
211 Blas::Daxpy(ntot, df[5][0], &Diff2[0], 1, &out_d1[0], 1);
216 Vmath::Smul(ntot, df[6][0], &Diff0[0], 1, &out_d2[0], 1);
217 Blas::Daxpy(ntot, df[7][0], &Diff1[0], 1, &out_d2[0], 1);
218 Blas::Daxpy(ntot, df[8][0], &Diff2[0], 1, &out_d2[0], 1);
258 ASSERTL1(
false,
"input dir is out of range");
271 int nquad0 =
m_base[0]->GetNumPoints();
272 int nquad1 =
m_base[1]->GetNumPoints();
273 int nquad2 =
m_base[2]->GetNumPoints();
274 int ntot = nquad0 * nquad1 * nquad2;
282 StdHexExp::v_PhysDeriv(inarray, Diff0, Diff1, Diff2);
287 Vmath::Vmul(ntot, &dfdir[0][0], 1, &Diff0[0], 1, &outarray[0], 1);
288 Vmath::Vvtvp(ntot, &dfdir[1][0], 1, &Diff1[0], 1, &outarray[0], 1,
290 Vmath::Vvtvp(ntot, &dfdir[2][0], 1, &Diff2[0], 1, &outarray[0], 1,
309 if (
m_base[0]->Collocation() &&
m_base[1]->Collocation() &&
326 out = (*matsys) * in;
383 int nquad0 =
m_base[0]->GetNumPoints();
384 int nquad1 =
m_base[1]->GetNumPoints();
385 int nquad2 =
m_base[2]->GetNumPoints();
386 int order0 =
m_base[0]->GetNumModes();
387 int order1 =
m_base[1]->GetNumModes();
390 order0 * order1 * nquad2);
392 if (multiplybyweights)
399 tmp, outarray, wsp,
true,
true,
true);
405 inarray, outarray, wsp,
true,
true,
true);
440 ASSERTL1((dir == 0) || (dir == 1) || (dir == 2),
"Invalid direction.");
442 const int nq0 =
m_base[0]->GetNumPoints();
443 const int nq1 =
m_base[1]->GetNumPoints();
444 const int nq2 =
m_base[2]->GetNumPoints();
445 const int nq = nq0 * nq1 * nq2;
446 const int nm0 =
m_base[0]->GetNumModes();
447 const int nm1 =
m_base[1]->GetNumModes();
467 m_base[2]->GetBdata(), tmp2, outarray, wsp,
471 m_base[2]->GetBdata(), tmp3, tmp5, wsp,
true,
476 m_base[2]->GetDbdata(), tmp4, tmp5, wsp,
true,
485 ASSERTL1((dir == 0) || (dir == 1) || (dir == 2),
"Invalid direction.");
487 const int nq0 =
m_base[0]->GetNumPoints();
488 const int nq1 =
m_base[1]->GetNumPoints();
489 const int nq2 =
m_base[2]->GetNumPoints();
490 const int nq = nq0 * nq1 * nq2;
505 Vmath::Vmul(nq, &df[3 * dir][0], 1, tmp1.get(), 1, tmp2.get(), 1);
506 Vmath::Vmul(nq, &df[3 * dir + 1][0], 1, tmp1.get(), 1, tmp3.get(), 1);
507 Vmath::Vmul(nq, &df[3 * dir + 2][0], 1, tmp1.get(), 1, tmp4.get(), 1);
511 Vmath::Smul(nq, df[3 * dir][0], tmp1.get(), 1, tmp2.get(), 1);
512 Vmath::Smul(nq, df[3 * dir + 1][0], tmp1.get(), 1, tmp3.get(), 1);
513 Vmath::Smul(nq, df[3 * dir + 2][0], tmp1.get(), 1, tmp4.get(), 1);
543 ASSERTL1(
false,
"input dir is out of range");
552 (iprodmat->GetOwnedMatrix())->GetPtr().get(),
m_ncoeffs,
553 inarray.get(), 1, 0.0, outarray.get(), 1);
568 const int nq0 =
m_base[0]->GetNumPoints();
569 const int nq1 =
m_base[1]->GetNumPoints();
570 const int nq2 =
m_base[2]->GetNumPoints();
571 const int nq = nq0 * nq1 * nq2;
572 const int nm0 =
m_base[0]->GetNumModes();
573 const int nm1 =
m_base[1]->GetNumModes();
591 Vmath::Vmul(nq, &dfdir[0][0], 1, tmp1.get(), 1, tmp2.get(), 1);
592 Vmath::Vmul(nq, &dfdir[1][0], 1, tmp1.get(), 1, tmp3.get(), 1);
593 Vmath::Vmul(nq, &dfdir[2][0], 1, tmp1.get(), 1, tmp4.get(), 1);
596 m_base[2]->GetBdata(), tmp2, outarray, wsp,
600 m_base[2]->GetBdata(), tmp3, tmp5, wsp,
true,
606 m_base[2]->GetDbdata(), tmp4, tmp5, wsp,
true,
626 return StdHexExp::v_PhysEvaluate(Lcoord, physvals);
635 m_geom->GetLocCoords(coord, Lcoord);
636 return StdHexExp::v_PhysEvaluate(Lcoord, physvals);
643 m_base[2]->GetBasisKey());
649 m_base[0]->GetPointsKey());
651 m_base[1]->GetPointsKey());
653 m_base[2]->GetPointsKey());
671 ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[0] <= 1.0 && Lcoords[1] >= -1.0 &&
672 Lcoords[1] <= 1.0 && Lcoords[2] >= -1.0 && Lcoords[2] <= 1.0,
673 "Local coordinates are not in region [-1,1]");
677 for (i = 0; i <
m_geom->GetCoordim(); ++i)
679 coords[i] =
m_geom->GetCoord(i, Lcoords);
701 const NekDouble *data,
const std::vector<unsigned int> &nummodes,
702 const int mode_offset,
NekDouble *coeffs,
703 std::vector<LibUtilities::BasisType> &fromType)
705 int data_order0 = nummodes[mode_offset];
706 int fillorder0 = min(
m_base[0]->GetNumModes(), data_order0);
707 int data_order1 = nummodes[mode_offset + 1];
708 int order1 =
m_base[1]->GetNumModes();
709 int fillorder1 = min(order1, data_order1);
710 int data_order2 = nummodes[mode_offset + 2];
711 int order2 =
m_base[2]->GetNumModes();
712 int fillorder2 = min(order2, data_order2);
724 m_base[0]->GetPointsKey()),
726 m_base[1]->GetPointsKey()),
728 m_base[2]->GetPointsKey()));
731 m_base[2]->GetBasisKey());
753 "Extraction routine not set up for this basis");
755 "Extraction routine not set up for this basis");
758 for (j = 0; j < fillorder0; ++j)
760 for (i = 0; i < fillorder1; ++i)
762 Vmath::Vcopy(fillorder2, &data[cnt], 1, &coeffs[cnt1], 1);
768 for (i = fillorder1; i < data_order1; ++i)
773 for (i = fillorder1; i < order1; ++i)
798 ASSERTL0(
false,
"basis is either not set up or not "
815 int nquad0 =
m_base[0]->GetNumPoints();
816 int nquad1 =
m_base[1]->GetNumPoints();
817 int nquad2 =
m_base[2]->GetNumPoints();
829 if (outarray.size() != nq0 * nq1)
834 for (
int i = 0; i < nquad0 * nquad1; ++i)
844 if (outarray.size() != nq0 * nq1)
850 for (
int k = 0; k < nquad2; k++)
852 for (
int i = 0; i < nquad0; ++i)
854 outarray[k * nquad0 + i] = nquad0 * nquad1 * k + i;
863 if (outarray.size() != nq0 * nq1)
868 for (
int i = 0; i < nquad1 * nquad2; i++)
870 outarray[i] = nquad0 - 1 + i * nquad0;
878 if (outarray.size() != nq0 * nq1)
883 for (
int k = 0; k < nquad2; k++)
885 for (
int i = 0; i < nquad0; i++)
887 outarray[k * nquad0 + i] =
888 (nquad0 * (nquad1 - 1)) + (k * nquad0 * nquad1) + i;
897 if (outarray.size() != nq0 * nq1)
902 for (
int i = 0; i < nquad1 * nquad2; i++)
904 outarray[i] = i * nquad0;
911 if (outarray.size() != nq0 * nq1)
916 for (
int i = 0; i < nquad0 * nquad1; i++)
918 outarray[i] = nquad0 * nquad1 * (nquad2 - 1) + i;
923 ASSERTL0(
false,
"face value (> 5) is out of range");
936 for (i = 0; i < ptsKeys.size(); ++i)
947 geomFactors->GetDerivFactors(ptsKeys);
960 for (i = 0; i < vCoordDim; ++i)
965 size_t nqb = nq_face;
979 for (i = 0; i < vCoordDim; ++i)
981 normal[i][0] = -df[3 * i + 2][0];
985 for (i = 0; i < vCoordDim; ++i)
987 normal[i][0] = -df[3 * i + 1][0];
991 for (i = 0; i < vCoordDim; ++i)
993 normal[i][0] = df[3 * i][0];
997 for (i = 0; i < vCoordDim; ++i)
999 normal[i][0] = df[3 * i + 1][0];
1003 for (i = 0; i < vCoordDim; ++i)
1005 normal[i][0] = -df[3 * i][0];
1009 for (i = 0; i < vCoordDim; ++i)
1011 normal[i][0] = df[3 * i + 2][0];
1015 ASSERTL0(
false,
"face is out of range (edge < 5)");
1020 for (i = 0; i < vCoordDim; ++i)
1022 fac += normal[i][0] * normal[i][0];
1024 fac = 1.0 /
sqrt(fac);
1027 for (i = 0; i < vCoordDim; ++i)
1029 Vmath::Fill(nq_face, fac * normal[i][0], normal[i], 1);
1036 int nqe0 = ptsKeys[0].GetNumPoints();
1037 int nqe1 = ptsKeys[0].GetNumPoints();
1038 int nqe2 = ptsKeys[0].GetNumPoints();
1039 int nqe01 = nqe0 * nqe1;
1040 int nqe02 = nqe0 * nqe2;
1041 int nqe12 = nqe1 * nqe2;
1044 if (face == 0 || face == 5)
1048 else if (face == 1 || face == 3)
1069 for (j = 0; j < nqe; ++j)
1071 normals[j] = -df[2][j] * jac[j];
1072 normals[nqe + j] = -df[5][j] * jac[j];
1073 normals[2 * nqe + j] = -df[8][j] * jac[j];
1074 faceJac[j] = jac[j];
1077 points0 = ptsKeys[0];
1078 points1 = ptsKeys[1];
1081 for (j = 0; j < nqe0; ++j)
1083 for (k = 0; k < nqe2; ++k)
1085 int idx = j + nqe01 * k;
1086 normals[j + k * nqe0] = -df[1][idx] * jac[idx];
1087 normals[nqe + j + k * nqe0] = -df[4][idx] * jac[idx];
1088 normals[2 * nqe + j + k * nqe0] =
1089 -df[7][idx] * jac[idx];
1090 faceJac[j + k * nqe0] = jac[idx];
1093 points0 = ptsKeys[0];
1094 points1 = ptsKeys[2];
1097 for (j = 0; j < nqe1; ++j)
1099 for (k = 0; k < nqe2; ++k)
1101 int idx = nqe0 - 1 + nqe0 * j + nqe01 * k;
1102 normals[j + k * nqe1] = df[0][idx] * jac[idx];
1103 normals[nqe + j + k * nqe1] = df[3][idx] * jac[idx];
1104 normals[2 * nqe + j + k * nqe1] = df[6][idx] * jac[idx];
1105 faceJac[j + k * nqe1] = jac[idx];
1108 points0 = ptsKeys[1];
1109 points1 = ptsKeys[2];
1112 for (j = 0; j < nqe0; ++j)
1114 for (k = 0; k < nqe2; ++k)
1116 int idx = nqe0 * (nqe1 - 1) + j + nqe01 * k;
1117 normals[j + k * nqe0] = df[1][idx] * jac[idx];
1118 normals[nqe + j + k * nqe0] = df[4][idx] * jac[idx];
1119 normals[2 * nqe + j + k * nqe0] = df[7][idx] * jac[idx];
1120 faceJac[j + k * nqe0] = jac[idx];
1123 points0 = ptsKeys[0];
1124 points1 = ptsKeys[2];
1127 for (j = 0; j < nqe1; ++j)
1129 for (k = 0; k < nqe2; ++k)
1131 int idx = j * nqe0 + nqe01 * k;
1132 normals[j + k * nqe1] = -df[0][idx] * jac[idx];
1133 normals[nqe + j + k * nqe1] = -df[3][idx] * jac[idx];
1134 normals[2 * nqe + j + k * nqe1] =
1135 -df[6][idx] * jac[idx];
1136 faceJac[j + k * nqe1] = jac[idx];
1139 points0 = ptsKeys[1];
1140 points1 = ptsKeys[2];
1143 for (j = 0; j < nqe01; ++j)
1145 int idx = j + nqe01 * (nqe2 - 1);
1146 normals[j] = df[2][idx] * jac[idx];
1147 normals[nqe + j] = df[5][idx] * jac[idx];
1148 normals[2 * nqe + j] = df[8][idx] * jac[idx];
1149 faceJac[j] = jac[idx];
1151 points0 = ptsKeys[0];
1152 points1 = ptsKeys[1];
1155 ASSERTL0(
false,
"face is out of range (face < 5)");
1164 Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
1172 Vmath::Vmul(nq_face, work, 1, normal[i], 1, normal[i], 1);
1179 Vmath::Vvtvp(nq_face, normal[i], 1, normal[i], 1, work, 1, work, 1);
1189 Vmath::Vmul(nq_face, normal[i], 1, work, 1, normal[i], 1);
1201 StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
1216 StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
1224 StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
1231 StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(inarray, outarray, mkey);
1238 StdExpansion::MassLevelCurvatureMatrixOp_MatFree(inarray, outarray, mkey);
1283 if (inarray.get() == outarray.get())
1289 (mat->GetOwnedMatrix())->GetPtr().get(),
m_ncoeffs,
1290 tmp.get(), 1, 0.0, outarray.get(), 1);
1295 (mat->GetOwnedMatrix())->GetPtr().get(),
m_ncoeffs,
1296 inarray.get(), 1, 0.0, outarray.get(), 1);
1315 int n_coeffs = inarray.size();
1316 int nmodes0 =
m_base[0]->GetNumModes();
1317 int nmodes1 =
m_base[1]->GetNumModes();
1318 int nmodes2 =
m_base[2]->GetNumModes();
1319 int numMax = nmodes0;
1347 int cnt = 0, cnt2 = 0;
1349 for (
int u = 0; u < numMin + 1; ++u)
1351 for (
int i = 0; i < numMin; ++i)
1354 tmp2 = coeff_tmp1 + cnt, 1);
1360 tmp4 = coeff_tmp2 + cnt2, 1);
1362 cnt2 = u * nmodes0 * nmodes1;
1390 StdHexExp::v_SVVLaplacianFilter(array, mkey);
1415 returnval = StdHexExp::v_GenMatrix(mkey);
1430 return tmp->GetStdMatrix(mkey);
1459 int nquad0 =
m_base[0]->GetNumPoints();
1460 int nquad1 =
m_base[1]->GetNumPoints();
1461 int nquad2 =
m_base[2]->GetNumPoints();
1462 int nqtot = nquad0 * nquad1 * nquad2;
1464 ASSERTL1(wsp.size() >= 6 * nqtot,
"Insufficient workspace size.");
1493 StdExpansion3D::PhysTensorDeriv(inarray, wsp0, wsp1, wsp2);
1499 Vmath::Vvtvvtp(nqtot, &metric00[0], 1, &wsp0[0], 1, &metric01[0], 1,
1500 &wsp1[0], 1, &wsp3[0], 1);
1501 Vmath::Vvtvp(nqtot, &metric02[0], 1, &wsp2[0], 1, &wsp3[0], 1, &wsp3[0], 1);
1502 Vmath::Vvtvvtp(nqtot, &metric01[0], 1, &wsp0[0], 1, &metric11[0], 1,
1503 &wsp1[0], 1, &wsp4[0], 1);
1504 Vmath::Vvtvp(nqtot, &metric12[0], 1, &wsp2[0], 1, &wsp4[0], 1, &wsp4[0], 1);
1505 Vmath::Vvtvvtp(nqtot, &metric02[0], 1, &wsp0[0], 1, &metric12[0], 1,
1506 &wsp1[0], 1, &wsp5[0], 1);
1507 Vmath::Vvtvp(nqtot, &metric22[0], 1, &wsp2[0], 1, &wsp5[0], 1, &wsp5[0], 1);
1530 const unsigned int dim = 3;
1536 for (
unsigned int i = 0; i < dim; ++i)
1538 for (
unsigned int j = i; j < dim; ++j)
1575 if (d0factors.size() != 6)
1582 if (d0factors[0].size() != nquad0 * nquad1)
1592 if (d0factors[1].size() != nquad0 * nquad2)
1602 if (d0factors[2].size() != nquad1 * nquad2)
1626 int ncoords = normal_0.size();
1631 for (
int i = 0; i < nquad0 * nquad1; ++i)
1633 d0factors[0][i] = df[0][i] * normal_0[0][i];
1635 df[0][nquad0 * nquad1 * (nquad2 - 1) + i] * normal_5[0][i];
1636 d1factors[0][i] = df[1][i] * normal_0[0][i];
1638 df[1][nquad0 * nquad1 * (nquad2 - 1) + i] * normal_5[0][i];
1639 d2factors[0][i] = df[2][i] * normal_0[0][i];
1641 df[2][nquad0 * nquad1 * (nquad2 - 1) + i] * normal_5[0][i];
1644 for (
int n = 1; n < ncoords; ++n)
1646 for (
int i = 0; i < nquad0 * nquad1; ++i)
1648 d0factors[0][i] += df[3 * n][i] * normal_0[n][i];
1650 df[3 * n][nquad0 * nquad1 * (nquad2 - 1) + i] *
1652 d1factors[0][i] += df[3 * n + 1][i] * normal_0[n][i];
1654 df[3 * n + 1][nquad0 * nquad1 * (nquad2 - 1) + i] *
1656 d2factors[0][i] += df[3 * n + 2][i] * normal_0[n][i];
1658 df[3 * n + 2][nquad0 * nquad1 * (nquad2 - 1) + i] *
1664 for (
int j = 0; j < nquad2; ++j)
1666 for (
int i = 0; i < nquad0; ++i)
1668 d0factors[1][i] = df[0][j * nquad0 * nquad1 + i] *
1669 normal_1[0][j * nquad0 + i];
1671 df[0][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1672 normal_3[0][j * nquad0 + i];
1673 d1factors[1][i] = df[1][j * nquad0 * nquad1 + i] *
1674 normal_1[0][j * nquad0 + i];
1676 df[1][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1677 normal_3[0][j * nquad0 + i];
1678 d2factors[1][i] = df[2][j * nquad0 * nquad1 + i] *
1679 normal_1[0][j * nquad0 + i];
1681 df[2][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1682 normal_3[0][j * nquad0 + i];
1686 for (
int n = 1; n < ncoords; ++n)
1688 for (
int j = 0; j < nquad2; ++j)
1690 for (
int i = 0; i < nquad0; ++i)
1692 d0factors[1][i] = df[3 * n][j * nquad0 * nquad1 + i] *
1693 normal_1[0][j * nquad0 + i];
1695 df[3 * n][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1696 normal_3[0][j * nquad0 + i];
1697 d1factors[1][i] = df[3 * n + 1][j * nquad0 * nquad1 + i] *
1698 normal_1[0][j * nquad0 + i];
1700 df[3 * n + 1][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1701 normal_3[0][j * nquad0 + i];
1702 d2factors[1][i] = df[3 * n + 2][j * nquad0 * nquad1 + i] *
1703 normal_1[0][j * nquad0 + i];
1705 df[3 * n + 2][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1706 normal_3[0][j * nquad0 + i];
1712 for (
int j = 0; j < nquad2; ++j)
1714 for (
int i = 0; i < nquad1; ++i)
1716 d0factors[2][j * nquad1 + i] =
1717 df[0][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1718 normal_2[0][j * nquad1 + i];
1719 d0factors[4][j * nquad1 + i] =
1720 df[0][j * nquad0 * nquad1 + i * nquad0] *
1721 normal_4[0][j * nquad1 + i];
1722 d1factors[2][j * nquad1 + i] =
1723 df[1][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1724 normal_2[0][j * nquad1 + i];
1725 d1factors[4][j * nquad1 + i] =
1726 df[1][j * nquad0 * nquad1 + i * nquad0] *
1727 normal_4[0][j * nquad1 + i];
1728 d2factors[2][j * nquad1 + i] =
1729 df[2][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1730 normal_2[0][j * nquad1 + i];
1731 d2factors[4][j * nquad1 + i] =
1732 df[2][j * nquad0 * nquad1 + i * nquad0] *
1733 normal_4[0][j * nquad1 + i];
1737 for (
int n = 1; n < ncoords; ++n)
1739 for (
int j = 0; j < nquad2; ++j)
1741 for (
int i = 0; i < nquad1; ++i)
1743 d0factors[2][j * nquad1 + i] +=
1744 df[3 * n][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1745 normal_2[n][j * nquad0 + i];
1746 d0factors[4][j * nquad0 + i] +=
1747 df[3 * n][i * nquad0 + j * nquad0 * nquad1] *
1748 normal_4[n][j * nquad0 + i];
1749 d1factors[2][j * nquad1 + i] +=
1751 [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1752 normal_2[n][j * nquad0 + i];
1753 d1factors[4][j * nquad0 + i] +=
1754 df[3 * n + 1][i * nquad0 + j * nquad0 * nquad1] *
1755 normal_4[n][j * nquad0 + i];
1756 d2factors[2][j * nquad1 + i] +=
1758 [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1759 normal_2[n][j * nquad0 + i];
1760 d2factors[4][j * nquad0 + i] +=
1761 df[3 * n + 2][i * nquad0 + j * nquad0 * nquad1] *
1762 normal_4[n][j * nquad0 + i];
1770 for (
int i = 0; i < nquad0 * nquad1; ++i)
1772 d0factors[0][i] = df[0][0] * normal_0[0][i];
1773 d0factors[5][i] = df[0][0] * normal_5[0][i];
1775 d1factors[0][i] = df[1][0] * normal_0[0][i];
1776 d1factors[5][i] = df[1][0] * normal_5[0][i];
1778 d2factors[0][i] = df[2][0] * normal_0[0][i];
1779 d2factors[5][i] = df[2][0] * normal_5[0][i];
1782 for (
int n = 1; n < ncoords; ++n)
1784 for (
int i = 0; i < nquad0 * nquad1; ++i)
1786 d0factors[0][i] += df[3 * n][0] * normal_0[n][i];
1787 d0factors[5][i] += df[3 * n][0] * normal_5[n][i];
1789 d1factors[0][i] += df[3 * n + 1][0] * normal_0[n][i];
1790 d1factors[5][i] += df[3 * n + 1][0] * normal_5[n][i];
1792 d2factors[0][i] += df[3 * n + 2][0] * normal_0[n][i];
1793 d2factors[5][i] += df[3 * n + 2][0] * normal_5[n][i];
1798 for (
int i = 0; i < nquad0 * nquad2; ++i)
1800 d0factors[1][i] = df[0][0] * normal_1[0][i];
1801 d0factors[3][i] = df[0][0] * normal_3[0][i];
1803 d1factors[1][i] = df[1][0] * normal_1[0][i];
1804 d1factors[3][i] = df[1][0] * normal_3[0][i];
1806 d2factors[1][i] = df[2][0] * normal_1[0][i];
1807 d2factors[3][i] = df[2][0] * normal_3[0][i];
1810 for (
int n = 1; n < ncoords; ++n)
1812 for (
int i = 0; i < nquad0 * nquad2; ++i)
1814 d0factors[1][i] += df[3 * n][0] * normal_1[n][i];
1815 d0factors[3][i] += df[3 * n][0] * normal_3[n][i];
1817 d1factors[1][i] += df[3 * n + 1][0] * normal_1[n][i];
1818 d1factors[3][i] += df[3 * n + 1][0] * normal_3[n][i];
1820 d2factors[1][i] += df[3 * n + 2][0] * normal_1[n][i];
1821 d2factors[3][i] += df[3 * n + 2][0] * normal_3[n][i];
1826 for (
int i = 0; i < nquad1 * nquad2; ++i)
1828 d0factors[2][i] = df[0][0] * normal_2[0][i];
1829 d0factors[4][i] = df[0][0] * normal_4[0][i];
1831 d1factors[2][i] = df[1][0] * normal_2[0][i];
1832 d1factors[4][i] = df[1][0] * normal_4[0][i];
1834 d2factors[2][i] = df[2][0] * normal_2[0][i];
1835 d2factors[4][i] = df[2][0] * normal_4[0][i];
1838 for (
int n = 1; n < ncoords; ++n)
1840 for (
int i = 0; i < nquad1 * nquad2; ++i)
1842 d0factors[2][i] += df[3 * n][0] * normal_2[n][i];
1843 d0factors[4][i] += df[3 * n][0] * normal_4[n][i];
1845 d1factors[2][i] += df[3 * n + 1][0] * normal_2[n][i];
1846 d1factors[4][i] += df[3 * n + 1][0] * normal_4[n][i];
1848 d2factors[2][i] += df[3 * n + 2][0] * normal_2[n][i];
1849 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)
SpatialDomains::Geometry3DSharedPtr GetGeom3D() const
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 ComputeGmatcdotMF(const Array< TwoD, const NekDouble > &df, const Array< OneD, const NekDouble > &direction, Array< OneD, Array< OneD, NekDouble >> &dfdir)
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)
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
const NormalVector & GetTraceNormal(const int id)
void IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Calculates the inner product .
virtual NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &physvals)
This function evaluates the expansion at a single (arbitrary) point of the domain.
virtual void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp)
virtual NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
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_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
virtual void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
virtual void v_WeakDirectionalDerivMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Calculate the inner product of inarray with respect to the given basis B = base0 * base1 * base2.
virtual void v_GetTracePhysMap(const int face, Array< OneD, int > &outarray)
virtual LibUtilities::ShapeType v_DetShapeType() const
Return the region shape using the enum-list of ShapeType.
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void IProductWRTDirectionalDerivBase_SumFac(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
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_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey)
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_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
virtual void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
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 StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const
virtual void v_ComputeLaplacianMetric()
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
void v_GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
virtual void v_MassLevelCurvatureMatrixOp(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 elements basis.
virtual void v_NormalTraceDerivFactors(Array< OneD, Array< OneD, NekDouble >> &factors, Array< OneD, Array< OneD, NekDouble >> &d0factors, Array< OneD, Array< OneD, NekDouble >> &d1factors)
: This method gets all of the factors which are required as part of the Gradient Jump Penalty stabili...
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
void v_PhysDirectionalDeriv(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &direction, Array< OneD, NekDouble > &out)
Physical derivative along a direction vector.
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Integrate the physical point list inarray over region.
virtual bool v_GetFaceDGForwards(const int i) const
virtual void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords)
Retrieves the physical coordinates of a given set of reference coordinates.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
virtual void v_LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
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 IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
this function calculates the inner product of a given function f with the different modes of the expa...
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space.
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space.
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Class representing a hexehedral element in reference space.
MatrixType GetMatrixType() const
static void 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 Daxpy(const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: y = alpha x plus y.
void Interp3D(const BasisKey &fbasis0, const BasisKey &fbasis1, const BasisKey &fbasis2, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, const BasisKey &tbasis2, Array< OneD, NekDouble > &to)
this function interpolates a 3D function evaluated at the quadrature points of the 3D basis,...
void InterpCoeff3D(const BasisKey &fbasis0, const BasisKey &fbasis1, const BasisKey &fbasis2, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, const BasisKey &tbasis2, Array< OneD, NekDouble > &to)
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis,...
std::vector< PointsKey > PointsKeyVector
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
@ eOrtho_A
Principle Orthogonal Functions .
@ eGLL_Lagrange
Lagrange for SEM basis .
@ eModified_A
Principle Modified Functions .
std::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
GeomType
Indicates the type of element geometry.
@ eRegular
Geometry is straight-sided with constant geometric factors.
@ eMovingRegular
Currently unused.
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< HexGeom > HexGeomSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
@ eInvLaplacianWithUnityMean
std::shared_ptr< StdHexExp > StdHexExpSharedPtr
@ eDir1FwdDir1_Dir2FwdDir2
@ eDir1BwdDir1_Dir2BwdDir2
@ eDir1BwdDir2_Dir2FwdDir1
@ eDir1FwdDir2_Dir2BwdDir1
The above copyright notice and this permission notice shall be included.
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
static Array< OneD, NekDouble > NullNekDouble1DArray
std::shared_ptr< DNekMat > DNekMatSharedPtr
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/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)