50 : StdExpansion(Ba.GetNumModes() * Bb.GetNumModes(), 2, Ba, Bb),
51 StdExpansion2D(Ba.GetNumModes() * Bb.GetNumModes(), Ba, Bb),
55 std::string(
"QuadExpMatrix")),
56 m_staticCondMatrixManager(
std::bind(&
Expansion::CreateStaticCondMatrix,
57 this,
std::placeholders::_1),
58 std::string(
"QuadExpStaticCondMatrix"))
63 : StdExpansion(T), StdExpansion2D(T), StdQuadExp(T),
Expansion(T),
65 m_staticCondMatrixManager(T.m_staticCondMatrixManager)
71 int nquad0 =
m_base[0]->GetNumPoints();
72 int nquad1 =
m_base[1]->GetNumPoints();
80 Vmath::Vmul(nquad0 * nquad1, jac, 1, inarray, 1, tmp, 1);
84 Vmath::Smul(nquad0 * nquad1, jac[0], inarray, 1, tmp, 1);
88 ival = StdQuadExp::v_Integral(tmp);
97 int nquad0 =
m_base[0]->GetNumPoints();
98 int nquad1 =
m_base[1]->GetNumPoints();
99 int nqtot = nquad0 * nquad1;
105 StdQuadExp::v_PhysDeriv(inarray, diff0, diff1);
112 Vmath::Vvtvp(nqtot, df[1], 1, diff1, 1, out_d0, 1, out_d0, 1);
118 Vmath::Vvtvp(nqtot, df[3], 1, diff1, 1, out_d1, 1, out_d1, 1);
124 Vmath::Vvtvp(nqtot, df[5], 1, diff1, 1, out_d2, 1, out_d2, 1);
175 ASSERTL1(
false,
"input dir is out of range");
185 int nquad0 =
m_base[0]->GetNumPoints();
186 int nquad1 =
m_base[1]->GetNumPoints();
187 int nqtot = nquad0 * nquad1;
195 StdQuadExp::v_PhysDeriv(inarray, diff0, diff1);
202 for (
int i = 0; i < 2; ++i)
205 for (
int k = 0; k < (
m_geom->GetCoordim()); ++k)
207 Vmath::Vvtvp(nqtot, &df[2 * k + i][0], 1, &direction[k * nqtot],
208 1, &tangmat[i][0], 1, &tangmat[i][0], 1);
215 Vmath::Vmul(nqtot, &tangmat[0][0], 1, &diff0[0], 1, &out[0], 1);
216 Vmath::Vvtvp(nqtot, &tangmat[1][0], 1, &diff1[0], 1, &out[0], 1,
230 if ((
m_base[0]->Collocation()) && (
m_base[1]->Collocation()))
246 out = (*matsys) * in;
254 if ((
m_base[0]->Collocation()) && (
m_base[1]->Collocation()))
261 int npoints[2] = {
m_base[0]->GetNumPoints(),
m_base[1]->GetNumPoints()};
262 int nmodes[2] = {
m_base[0]->GetNumModes(),
m_base[1]->GetNumModes()};
264 fill(outarray.data(), outarray.data() +
m_ncoeffs, 0.0);
266 if (nmodes[0] == 1 && nmodes[1] == 1)
268 outarray[0] = inarray[0];
275 for (i = 0; i < 4; i++)
282 for (i = 0; i < npoints[0]; i++)
284 physEdge[0][i] = inarray[i];
285 physEdge[2][i] = inarray[npoints[0] * (npoints[1] - 1) + i];
288 for (i = 0; i < npoints[1]; i++)
290 physEdge[1][i] = inarray[npoints[0] - 1 + i * npoints[0]];
291 physEdge[3][i] = inarray[i * npoints[0]];
294 for (i = 0; i < 4; i++)
298 reverse((physEdge[i]).data(),
299 (physEdge[i]).data() + npoints[i % 2]);
304 for (i = 0; i < 4; i++)
314 for (i = 0; i < 4; i++)
316 segexp[i % 2]->FwdTransBndConstrained(physEdge[i], coeffEdge[i]);
319 for (j = 0; j < nmodes[i % 2]; j++)
322 outarray[mapArray[j]] =
sign * coeffEdge[i][j];
327 int nInteriorDofs =
m_ncoeffs - nBoundaryDofs;
329 if (nInteriorDofs > 0)
353 for (i = 0; i < nInteriorDofs; i++)
355 rhs[i] = tmp1[mapArray[i]];
358 Blas::Dgemv(
'N', nInteriorDofs, nInteriorDofs, matsys->Scale(),
359 &((matsys->GetOwnedMatrix())->GetPtr())[0],
360 nInteriorDofs, rhs.data(), 1, 0.0, result.data(), 1);
362 for (i = 0; i < nInteriorDofs; i++)
364 outarray[mapArray[i]] = result[i];
373 if (
m_base[0]->Collocation() &&
m_base[1]->Collocation())
394 int nquad0 =
m_base[0]->GetNumPoints();
395 int nquad1 =
m_base[1]->GetNumPoints();
396 int order0 =
m_base[0]->GetNumModes();
398 if (multiplybyweights)
404 StdQuadExp::IProductWRTBase_SumFacKernel(
m_base[0]->GetBdata(),
405 m_base[1]->GetBdata(), tmp,
406 outarray, wsp,
true,
true);
412 StdQuadExp::IProductWRTBase_SumFacKernel(
m_base[0]->GetBdata(),
413 m_base[1]->GetBdata(), inarray,
414 outarray, wsp,
true,
true);
422 ASSERTL1((dir == 0) || (dir == 1) || (dir == 2),
"Invalid direction.");
424 "Invalid direction.");
426 int nquad0 =
m_base[0]->GetNumPoints();
427 int nquad1 =
m_base[1]->GetNumPoints();
428 int nqtot = nquad0 * nquad1;
429 int nmodes0 =
m_base[0]->GetNumModes();
446 tmp1, tmp3, tmp4,
false,
true);
448 tmp2, outarray, tmp4,
true,
false);
456 ASSERTL1((dir == 0) || (dir == 1) || (dir == 2),
"Invalid direction.");
458 "Invalid direction.");
460 int nquad0 =
m_base[0]->GetNumPoints();
461 int nquad1 =
m_base[1]->GetNumPoints();
462 int nqtot = nquad0 * nquad1;
463 int nmodes0 =
m_base[0]->GetNumModes();
475 Vmath::Vmul(nqtot, &df[2 * dir][0], 1, inarray.data(), 1, tmp1.data(),
477 Vmath::Vmul(nqtot, &df[2 * dir + 1][0], 1, inarray.data(), 1,
482 Vmath::Smul(nqtot, df[2 * dir][0], inarray.data(), 1, tmp1.data(), 1);
483 Vmath::Smul(nqtot, df[2 * dir + 1][0], inarray.data(), 1, tmp2.data(),
493 int nq =
m_base[0]->GetNumPoints() *
m_base[1]->GetNumPoints();
502 Vmath::Vvtvvtp(nq, &normals[0][0], 1, &Fx[0], 1, &normals[1][0], 1,
503 &Fy[0], 1, &Fn[0], 1);
504 Vmath::Vvtvp(nq, &normals[2][0], 1, &Fz[0], 1, &Fn[0], 1, &Fn[0], 1);
508 Vmath::Svtsvtp(nq, normals[0][0], &Fx[0], 1, normals[1][0], &Fy[0], 1,
510 Vmath::Svtvp(nq, normals[2][0], &Fz[0], 1, &Fn[0], 1, &Fn[0], 1);
532 m_base[0]->GetPointsKey());
534 m_base[1]->GetPointsKey());
552 ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[1] <= 1.0 && Lcoords[1] >= -1.0 &&
554 "Local coordinates are not in region [-1,1]");
557 for (i = 0; i <
m_geom->GetCoordim(); ++i)
559 coords[i] =
m_geom->GetCoord(i, Lcoords);
573 return StdExpansion2D::v_PhysEvaluate(Lcoord, physvals);
582 m_geom->GetLocCoords(coord, Lcoord);
584 return StdExpansion2D::v_PhysEvaluate(Lcoord, physvals);
590 std::array<NekDouble, 3> &firstOrderDerivs)
594 m_geom->GetLocCoords(coord, Lcoord);
595 return StdQuadExp::v_PhysEvalFirstDeriv(Lcoord, inarray, firstOrderDerivs);
603 int nquad0 =
m_base[0]->GetNumPoints();
604 int nquad1 =
m_base[1]->GetNumPoints();
613 Vmath::Vcopy(nquad0, &(inarray[0]), 1, &(outarray[0]), 1);
616 Vmath::Vcopy(nquad1, &(inarray[0]) + (nquad0 - 1), nquad0,
620 Vmath::Vcopy(nquad0, &(inarray[0]) + nquad0 * (nquad1 - 1), 1,
624 Vmath::Vcopy(nquad1, &(inarray[0]), nquad0, &(outarray[0]), 1);
627 ASSERTL0(
false,
"edge value (< 3) is out of range");
637 if (
m_base[edge % 2]->GetPointsKey() !=
638 EdgeExp->GetBasis(0)->GetPointsKey())
645 EdgeExp->GetBasis(0)->GetPointsKey(), outarray);
655 Vmath::Reverse(EdgeExp->GetNumPoints(0), &outarray[0], 1, &outarray[0],
665 int nq0 =
m_base[0]->GetNumPoints();
666 int nq1 =
m_base[1]->GetNumPoints();
680 for (i = 0; i < nq0; i++)
683 nq1, mat_gauss->GetOwnedMatrix()->GetPtr().data(), 1,
690 for (i = 0; i < nq1; i++)
693 nq0, mat_gauss->GetOwnedMatrix()->GetPtr().data(), 1,
694 &inarray[i * nq0], 1);
700 for (i = 0; i < nq0; i++)
703 nq1, mat_gauss->GetOwnedMatrix()->GetPtr().data(), 1,
710 for (i = 0; i < nq1; i++)
713 nq0, mat_gauss->GetOwnedMatrix()->GetPtr().data(), 1,
714 &inarray[i * nq0], 1);
719 ASSERTL0(
false,
"edge value (< 3) is out of range");
726 int nquad0 =
m_base[0]->GetNumPoints();
727 int nquad1 =
m_base[1]->GetNumPoints();
734 for (
int i = 0; i < nquad0; ++i)
741 for (
int i = 0; i < nquad1; ++i)
743 outarray[i] = (nquad0 - 1) + i * nquad0;
748 for (
int i = 0; i < nquad0; ++i)
750 outarray[i] = i + nquad0 * (nquad1 - 1);
755 for (
int i = 0; i < nquad1; ++i)
757 outarray[i] = i * nquad0;
761 ASSERTL0(
false,
"edge value (< 3) is out of range");
770 int nquad0 =
m_base[0]->GetNumPoints();
771 int nquad1 =
m_base[1]->GetNumPoints();
797 for (i = 0; i < nquad0; ++i)
800 j[i] *
sqrt(g1[i] * g1[i] + g3[i] * g3[i]);
804 Vmath::Vcopy(nquad1, &(df[0][0]) + (nquad0 - 1), nquad0,
807 Vmath::Vcopy(nquad1, &(df[2][0]) + (nquad0 - 1), nquad0,
813 for (i = 0; i < nquad1; ++i)
816 j[i] *
sqrt(g0[i] * g0[i] + g2[i] * g2[i]);
821 Vmath::Vcopy(nquad0, &(df[1][0]) + (nquad0 * (nquad1 - 1)),
824 Vmath::Vcopy(nquad0, &(df[3][0]) + (nquad0 * (nquad1 - 1)),
827 Vmath::Vcopy(nquad0, &(jac[0]) + (nquad0 * (nquad1 - 1)), 1,
830 for (i = 0; i < nquad0; ++i)
833 j[i] *
sqrt(g1[i] * g1[i] + g3[i] * g3[i]);
842 for (i = 0; i < nquad1; ++i)
845 j[i] *
sqrt(g0[i] * g0[i] + g2[i] * g2[i]);
849 ASSERTL0(
false,
"edge value (< 3) is out of range");
855 int nqtot = nquad0 * nquad1;
876 for (i = 0; i < nquad0; ++i)
878 outarray[i] =
sqrt(g1_edge[i] * g1_edge[i] +
879 g3_edge[i] * g3_edge[i]);
891 for (i = 0; i < nquad1; ++i)
893 outarray[i] =
sqrt(g0_edge[i] * g0_edge[i] +
894 g2_edge[i] * g2_edge[i]);
907 for (i = 0; i < nquad0; ++i)
909 outarray[i] =
sqrt(g1_edge[i] * g1_edge[i] +
910 g3_edge[i] * g3_edge[i]);
924 for (i = 0; i < nquad1; ++i)
926 outarray[i] =
sqrt(g0_edge[i] * g0_edge[i] +
927 g2_edge[i] * g2_edge[i]);
934 ASSERTL0(
false,
"edge value (< 3) is out of range");
946 for (i = 0; i < nquad0; ++i)
948 outarray[i] = jac[0] *
sqrt(df[1][0] * df[1][0] +
949 df[3][0] * df[3][0]);
953 for (i = 0; i < nquad1; ++i)
955 outarray[i] = jac[0] *
sqrt(df[0][0] * df[0][0] +
956 df[2][0] * df[2][0]);
960 for (i = 0; i < nquad0; ++i)
962 outarray[i] = jac[0] *
sqrt(df[1][0] * df[1][0] +
963 df[3][0] * df[3][0]);
967 for (i = 0; i < nquad1; ++i)
969 outarray[i] = jac[0] *
sqrt(df[0][0] * df[0][0] +
970 df[2][0] * df[2][0]);
974 ASSERTL0(
false,
"edge value (< 3) is out of range");
988 for (i = 0; i < ptsKeys.size(); ++i)
999 geomFactors->GetDerivFactors(ptsKeys);
1010 for (i = 0; i < vCoordDim; ++i)
1029 for (i = 0; i < vCoordDim; ++i)
1031 Vmath::Fill(nqe, -df[2 * i + 1][0], normal[i], 1);
1035 for (i = 0; i < vCoordDim; ++i)
1041 for (i = 0; i < vCoordDim; ++i)
1047 for (i = 0; i < vCoordDim; ++i)
1053 ASSERTL0(
false,
"edge is out of range (edge < 4)");
1058 for (i = 0; i < vCoordDim; ++i)
1060 fac += normal[i][0] * normal[i][0];
1062 fac = 1.0 /
sqrt(fac);
1066 for (i = 0; i < vCoordDim; ++i)
1068 Vmath::Smul(nqe, fac, normal[i], 1, normal[i], 1);
1075 int nquad0 = ptsKeys[0].GetNumPoints();
1076 int nquad1 = ptsKeys[1].GetNumPoints();
1094 for (j = 0; j < nquad0; ++j)
1096 edgejac[j] = jac[j];
1097 for (i = 0; i < vCoordDim; ++i)
1099 normals[i * nquad0 + j] =
1100 -df[2 * i + 1][j] * edgejac[j];
1103 from_key = ptsKeys[0];
1106 for (j = 0; j < nquad1; ++j)
1108 edgejac[j] = jac[nquad0 * j + nquad0 - 1];
1109 for (i = 0; i < vCoordDim; ++i)
1111 normals[i * nquad1 + j] =
1112 df[2 * i][nquad0 * j + nquad0 - 1] * edgejac[j];
1115 from_key = ptsKeys[1];
1118 for (j = 0; j < nquad0; ++j)
1120 edgejac[j] = jac[nquad0 * (nquad1 - 1) + j];
1121 for (i = 0; i < vCoordDim; ++i)
1123 normals[i * nquad0 + j] =
1124 (df[2 * i + 1][nquad0 * (nquad1 - 1) + j]) *
1128 from_key = ptsKeys[0];
1131 for (j = 0; j < nquad1; ++j)
1133 edgejac[j] = jac[nquad0 * j];
1134 for (i = 0; i < vCoordDim; ++i)
1136 normals[i * nquad1 + j] =
1137 -df[2 * i][nquad0 * j] * edgejac[j];
1140 from_key = ptsKeys[1];
1143 ASSERTL0(
false,
"edge is out of range (edge < 3)");
1148 int nqtot = nquad0 * nquad1;
1155 for (j = 0; j < nquad0; ++j)
1157 for (i = 0; i < vCoordDim; ++i)
1159 Vmath::Vmul(nqtot, &(df[2 * i + 1][0]), 1, &jac[0],
1160 1, &(tmp_gmat[0]), 1);
1163 normals[i * nquad0 + j] = -tmp_gmat_edge[j];
1166 from_key = ptsKeys[0];
1169 for (j = 0; j < nquad1; ++j)
1171 for (i = 0; i < vCoordDim; ++i)
1173 Vmath::Vmul(nqtot, &(df[2 * i][0]), 1, &jac[0], 1,
1177 normals[i * nquad1 + j] = tmp_gmat_edge[j];
1180 from_key = ptsKeys[1];
1183 for (j = 0; j < nquad0; ++j)
1185 for (i = 0; i < vCoordDim; ++i)
1187 Vmath::Vmul(nqtot, &(df[2 * i + 1][0]), 1, &jac[0],
1188 1, &(tmp_gmat[0]), 1);
1191 normals[i * nquad0 + j] = tmp_gmat_edge[j];
1194 from_key = ptsKeys[0];
1197 for (j = 0; j < nquad1; ++j)
1199 for (i = 0; i < vCoordDim; ++i)
1201 Vmath::Vmul(nqtot, &(df[2 * i][0]), 1, &jac[0], 1,
1205 normals[i * nquad1 + j] = -tmp_gmat_edge[j];
1208 from_key = ptsKeys[1];
1211 ASSERTL0(
false,
"edge is out of range (edge < 3)");
1227 Vmath::Vmul(nqe, work, 1, normal[i], 1, normal[i], 1);
1234 Vmath::Vvtvp(nqe, normal[i], 1, normal[i], 1, work, 1, work, 1);
1244 Vmath::Vmul(nqe, normal[i], 1, work, 1, normal[i], 1);
1249 for (i = 0; i < vCoordDim; ++i)
1260 const NekDouble *data,
const std::vector<unsigned int> &nummodes,
1262 std::vector<LibUtilities::BasisType> &fromType)
1264 int data_order0 = nummodes[mode_offset];
1265 int fillorder0 = std::min(
m_base[0]->GetNumModes(), data_order0);
1267 int data_order1 = nummodes[mode_offset + 1];
1268 int order1 =
m_base[1]->GetNumModes();
1269 int fillorder1 = min(order1, data_order1);
1280 m_base[0]->GetPointsKey()),
1282 m_base[1]->GetPointsKey()));
1284 m_base[1]->GetBasisKey());
1306 "Extraction routine not set up for this basis");
1309 for (i = 0; i < fillorder0; ++i)
1311 Vmath::Vcopy(fillorder1, data + cnt, 1, coeffs + cnt1, 1);
1345 ASSERTL0(
false,
"basis is either not set up or not hierarchicial");
1351 return m_geom->GetEorient(edge);
1369 returnval = StdQuadExp::v_GenMatrix(mkey);
1381 return tmp->GetStdMatrix(mkey);
1408 StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
1423 StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
1431 StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
1438 StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(inarray, outarray, mkey);
1445 StdExpansion::MassLevelCurvatureMatrixOp_MatFree(inarray, outarray, mkey);
1459 int n_coeffs = inarray.size();
1465 int nmodes0 =
m_base[0]->GetNumModes();
1466 int nmodes1 =
m_base[1]->GetNumModes();
1467 int numMax = nmodes0;
1485 for (
int i = 0; i < numMin + 1; ++i)
1487 Vmath::Vcopy(numMin, tmp = coeff + cnt, 1, tmp2 = coeff_tmp + cnt, 1);
1504 int nquad0 =
m_base[0]->GetNumPoints();
1505 int nquad1 =
m_base[1]->GetNumPoints();
1506 int nqtot = nquad0 * nquad1;
1507 int nmodes0 =
m_base[0]->GetNumModes();
1508 int nmodes1 =
m_base[1]->GetNumModes();
1510 max(max(max(nqtot,
m_ncoeffs), nquad1 * nmodes0), nquad0 * nmodes1);
1512 ASSERTL1(wsp.size() >= 3 * wspsize,
"Workspace is of insufficient size.");
1530 StdExpansion2D::PhysTensorDeriv(inarray, wsp1, wsp2);
1536 Vmath::Vvtvvtp(nqtot, &metric00[0], 1, &wsp1[0], 1, &metric01[0], 1,
1537 &wsp2[0], 1, &wsp0[0], 1);
1538 Vmath::Vvtvvtp(nqtot, &metric01[0], 1, &wsp1[0], 1, &metric11[0], 1,
1539 &wsp2[0], 1, &wsp2[0], 1);
1562 const unsigned int dim = 2;
1570 for (
unsigned int i = 0; i < dim; ++i)
1572 for (
unsigned int j = i; j < dim; ++j)
1611 StdQuadExp::v_SVVLaplacianFilter(array, mkey);
1633 if (d0factors.size() != 4)
1639 if (d0factors[0].size() != nquad0)
1647 if (d0factors[1].size() != nquad1)
1665 int ncoords = normal_0.size();
1672 for (
int i = 0; i < nquad0; ++i)
1674 d0factors[0][i] = df[0][i] * normal_0[0][i];
1675 d0factors[2][i] = df[0][nquad0 * (nquad1 - 1) + i] * normal_2[0][i];
1677 d1factors[0][i] = df[1][i] * normal_0[0][i];
1678 d1factors[2][i] = df[1][nquad0 * (nquad1 - 1) + i] * normal_2[0][i];
1681 for (
int n = 1; n < ncoords; ++n)
1685 for (
int i = 0; i < nquad0; ++i)
1687 d0factors[0][i] += df[2 * n][i] * normal_0[n][i];
1689 df[2 * n][nquad0 * (nquad1 - 1) + i] * normal_2[n][i];
1691 d1factors[0][i] += df[2 * n + 1][i] * normal_0[n][i];
1693 df[2 * n + 1][nquad0 * (nquad1 - 1) + i] * normal_2[n][i];
1698 for (
int i = 0; i < nquad1; ++i)
1700 d0factors[1][i] = df[0][(i + 1) * nquad0 - 1] * normal_1[0][i];
1701 d0factors[3][i] = df[0][i * nquad0] * normal_3[0][i];
1703 d1factors[1][i] = df[1][(i + 1) * nquad0 - 1] * normal_1[0][i];
1704 d1factors[3][i] = df[1][i * nquad0] * normal_3[0][i];
1707 for (
int n = 1; n < ncoords; ++n)
1709 for (
int i = 0; i < nquad1; ++i)
1712 df[2 * n][(i + 1) * nquad0 - 1] * normal_1[n][i];
1713 d0factors[3][i] += df[2 * n][i * nquad0] * normal_3[n][i];
1716 df[2 * n + 1][(i + 1) * nquad0 - 1] * normal_1[n][i];
1717 d1factors[3][i] += df[2 * n + 1][i * nquad0] * normal_3[n][i];
1724 for (
int i = 0; i < nquad0; ++i)
1726 d1factors[0][i] = df[1][0] * normal_0[0][i];
1727 d1factors[2][i] = df[1][0] * normal_2[0][i];
1731 for (
int i = 0; i < nquad1; ++i)
1733 d0factors[1][i] = df[0][0] * normal_1[0][i];
1734 d0factors[3][i] = df[0][0] * normal_3[0][i];
1737 for (
int n = 1; n < ncoords; ++n)
1741 for (
int i = 0; i < nquad0; ++i)
1743 d1factors[0][i] += df[2 * n + 1][0] * normal_0[n][i];
1744 d1factors[2][i] += df[2 * n + 1][0] * normal_2[n][i];
1749 for (
int i = 0; i < nquad1; ++i)
1751 d0factors[1][i] += df[2 * n][0] * normal_1[n][i];
1752 d0factors[3][i] += df[2 * n][0] * normal_3[n][i];
1758 for (
int i = 0; i < nquad0; ++i)
1760 d0factors[0][i] = df[0][0] * normal_0[0][i];
1761 d0factors[2][i] = df[0][0] * normal_2[0][i];
1765 for (
int i = 0; i < nquad1; ++i)
1767 d1factors[1][i] = df[1][0] * normal_1[0][i];
1768 d1factors[3][i] = df[1][0] * normal_3[0][i];
1771 for (
int n = 1; n < ncoords; ++n)
1775 for (
int i = 0; i < nquad0; ++i)
1777 d0factors[0][i] += df[2 * n][0] * normal_0[n][i];
1778 d0factors[2][i] += df[2 * n][0] * normal_2[n][i];
1783 for (
int i = 0; i < nquad1; ++i)
1785 d1factors[1][i] += df[2 * n + 1][0] * normal_1[n][i];
1786 d1factors[3][i] += df[2 * n + 1][0] * normal_3[n][i];
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
#define sign(a, b)
return the sign(b)*a
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.
size_t GetNumPoints() const
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
SpatialDomains::Geometry2DSharedPtr GetGeom2D() 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
ExpansionSharedPtr GetLeftAdjacentElementExp() const
void ComputeLaplacianMetric()
void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
int GetLeftAdjacentElementTrace() const
void ComputeQuadratureMetric()
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
StdRegions::Orientation GetTraceOrient(int trace)
const NormalVector & GetTraceNormal(const int id)
NekDouble v_PhysEvalFirstDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) override
void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
StdRegions::Orientation v_GetTraceOrient(int edge) override
DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey) override
void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Transform a given function from physical quadrature space to coefficient space.
DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey) override
StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const override
void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords) override
void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Calculate the inner product of inarray with respect to the basis B=base0*base1 and put into outarray.
void v_ComputeLaplacianMetric() override
void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp) override
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true) override
void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray) override
void GetEdgeInterpVals(const int edge, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, const Array< OneD, const NekDouble > &Fz, Array< OneD, NekDouble > &outarray) override
void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coord, const Array< OneD, const NekDouble > &physvals) override
This function evaluates the expansion at a single (arbitrary) point of the domain.
void v_GetTracePhysVals(const int edge, const StdRegions::StdExpansionSharedPtr &EdgeExp, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, StdRegions::Orientation orient) override
void v_WeakDirectionalDerivMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray) override
Calculate the derivative of the physical points.
NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray) override
Integrates the specified function over the domain.
void v_DropLocMatrix(const MatrixKey &mkey) override
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
void v_ExtractDataToCoeffs(const NekDouble *data, const std::vector< unsigned int > &nummodes, const int mode_offset, NekDouble *coeffs, std::vector< LibUtilities::BasisType > &fromType) override
void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
void v_ComputeTraceNormal(const int edge) override
void v_PhysDirectionalDeriv(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &direction, Array< OneD, NekDouble > &out) override
Physical derivative along a direction vector.
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const override
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
void v_DropLocStaticCondMatrix(const MatrixKey &mkey) override
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_MassLevelCurvatureMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey) override
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey) override
void v_NormalTraceDerivFactors(Array< OneD, Array< OneD, NekDouble > > &factors, Array< OneD, Array< OneD, NekDouble > > &d0factors, Array< OneD, Array< OneD, NekDouble > > &d1factors) override
: This method gets all of the factors which are required as part of the Gradient Jump Penalty (GJP) s...
NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals) override
void v_FwdTransBndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
void v_GetTracePhysMap(const int edge, Array< OneD, int > &outarray) override
void v_GetTraceQFactors(const int edge, Array< OneD, NekDouble > &outarray) override
void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
void IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0=true, bool doCheckCollDir1=true)
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.
void HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
int NumBndryCoeffs(void) const
const LibUtilities::PointsKeyVector GetPointsKeys() const
void MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, Array< OneD, NekDouble > &outarray)
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...
void LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
void GetTraceToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient=eForwards, int P=-1, int Q=-1)
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
void GetInteriorMap(Array< OneD, unsigned int > &outarray)
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 IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
const LibUtilities::BasisKey GetTraceBasisKey(const int i, int k=-1, bool UseGLL=false) const
This function returns the basis key belonging to the i-th trace.
Array< OneD, LibUtilities::BasisSharedPtr > m_base
MatrixType GetMatrixType() const
static void 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 = alpha A x plus beta y where A[m x n].
static double Ddot(const int &n, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: output = .
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 Interp1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
this function interpolates a 1D function evaluated at the quadrature points of the basis fbasis0 to ...
void InterpCoeff2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, 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
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
@ eGauss_Lagrange
Lagrange Polynomials using the Gauss points.
@ eOrtho_A
Principle Orthogonal Functions .
@ eGLL_Lagrange
Lagrange for SEM basis .
@ eModified_A
Principle Modified Functions .
std::shared_ptr< SegExp > SegExpSharedPtr
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< Geometry2D > Geometry2DSharedPtr
std::shared_ptr< StdQuadExp > StdQuadExpSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
@ eInvLaplacianWithUnityMean
std::map< ConstFactorType, NekDouble > ConstFactorMap
StdRegions::ConstFactorMap factors
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 Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
Svtsvtp (scalar times vector plus scalar times vector):
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Svtvp (scalar times vector plus vector): z = alpha*x + y.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/x.
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
void Zero(int n, T *x, const int incx)
Zero vector.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
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)
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
scalarT< T > sqrt(scalarT< T > in)