62 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
65 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
70 std::string(
"TetExpMatrix")),
71 m_staticCondMatrixManager(
std::bind(&
Expansion::CreateStaticCondMatrix,
72 this,
std::placeholders::_1),
73 std::string(
"TetExpStaticCondMatrix"))
81 : StdRegions::StdExpansion(T), StdRegions::StdExpansion3D(T),
83 m_matrixManager(T.m_matrixManager),
84 m_staticCondMatrixManager(T.m_staticCondMatrixManager)
104 int nquad0 =
m_base[0]->GetNumPoints();
105 int nquad1 =
m_base[1]->GetNumPoints();
106 int nquad2 =
m_base[2]->GetNumPoints();
115 (
NekDouble *)&inarray[0], 1, &tmp[0], 1);
120 (
NekDouble *)&inarray[0], 1, &tmp[0], 1);
124 retrunVal = StdTetExp::v_Integral(tmp);
146 int TotPts =
m_base[0]->GetNumPoints() *
m_base[1]->GetNumPoints() *
147 m_base[2]->GetNumPoints();
155 StdTetExp::v_PhysDeriv(inarray, Diff0, Diff1, Diff2);
161 Vmath::Vmul(TotPts, &df[0][0], 1, &Diff0[0], 1, &out_d0[0], 1);
162 Vmath::Vvtvp(TotPts, &df[1][0], 1, &Diff1[0], 1, &out_d0[0], 1,
164 Vmath::Vvtvp(TotPts, &df[2][0], 1, &Diff2[0], 1, &out_d0[0], 1,
170 Vmath::Vmul(TotPts, &df[3][0], 1, &Diff0[0], 1, &out_d1[0], 1);
171 Vmath::Vvtvp(TotPts, &df[4][0], 1, &Diff1[0], 1, &out_d1[0], 1,
173 Vmath::Vvtvp(TotPts, &df[5][0], 1, &Diff2[0], 1, &out_d1[0], 1,
179 Vmath::Vmul(TotPts, &df[6][0], 1, &Diff0[0], 1, &out_d2[0], 1);
180 Vmath::Vvtvp(TotPts, &df[7][0], 1, &Diff1[0], 1, &out_d2[0], 1,
182 Vmath::Vvtvp(TotPts, &df[8][0], 1, &Diff2[0], 1, &out_d2[0], 1,
190 Vmath::Smul(TotPts, df[0][0], &Diff0[0], 1, &out_d0[0], 1);
191 Blas::Daxpy(TotPts, df[1][0], &Diff1[0], 1, &out_d0[0], 1);
192 Blas::Daxpy(TotPts, df[2][0], &Diff2[0], 1, &out_d0[0], 1);
197 Vmath::Smul(TotPts, df[3][0], &Diff0[0], 1, &out_d1[0], 1);
198 Blas::Daxpy(TotPts, df[4][0], &Diff1[0], 1, &out_d1[0], 1);
199 Blas::Daxpy(TotPts, df[5][0], &Diff2[0], 1, &out_d1[0], 1);
204 Vmath::Smul(TotPts, df[6][0], &Diff0[0], 1, &out_d2[0], 1);
205 Blas::Daxpy(TotPts, df[7][0], &Diff1[0], 1, &out_d2[0], 1);
206 Blas::Daxpy(TotPts, df[8][0], &Diff2[0], 1, &out_d2[0], 1);
226 if ((
m_base[0]->Collocation()) && (
m_base[1]->Collocation()) &&
227 (
m_base[2]->Collocation()))
243 out = (*matsys) * in;
286 const int nquad0 =
m_base[0]->GetNumPoints();
287 const int nquad1 =
m_base[1]->GetNumPoints();
288 const int nquad2 =
m_base[2]->GetNumPoints();
289 const int order0 =
m_base[0]->GetNumModes();
290 const int order1 =
m_base[1]->GetNumModes();
292 nquad2 * order0 * (order1 + 1) / 2);
294 if (multiplybyweights)
301 tmp, outarray, wsp,
true,
true,
true);
307 inarray, outarray, wsp,
true,
true,
true);
345 const int nquad0 =
m_base[0]->GetNumPoints();
346 const int nquad1 =
m_base[1]->GetNumPoints();
347 const int nquad2 =
m_base[2]->GetNumPoints();
348 const int order0 =
m_base[0]->GetNumModes();
349 const int order1 =
m_base[1]->GetNumModes();
350 const int nqtot = nquad0 * nquad1 * nquad2;
358 nquad2 * order0 * (order1 + 1) / 2);
370 m_base[2]->GetBdata(), tmp2, outarray, wsp,
374 m_base[2]->GetBdata(), tmp3, tmp6, wsp,
true,
380 m_base[2]->GetDbdata(), tmp4, tmp6, wsp,
true,
392 const int nquad0 =
m_base[0]->GetNumPoints();
393 const int nquad1 =
m_base[1]->GetNumPoints();
394 const int nquad2 =
m_base[2]->GetNumPoints();
395 const int nqtot = nquad0 * nquad1 * nquad2;
409 Vmath::Vmul(nqtot, &df[3 * dir][0], 1, inarray.data(), 1, tmp2.data(),
411 Vmath::Vmul(nqtot, &df[3 * dir + 1][0], 1, inarray.data(), 1,
413 Vmath::Vmul(nqtot, &df[3 * dir + 2][0], 1, inarray.data(), 1,
414 outarray[2].data(), 1);
418 Vmath::Smul(nqtot, df[3 * dir][0], inarray.data(), 1, tmp2.data(), 1);
419 Vmath::Smul(nqtot, df[3 * dir + 1][0], inarray.data(), 1, tmp3.data(),
421 Vmath::Smul(nqtot, df[3 * dir + 2][0], inarray.data(), 1,
422 outarray[2].data(), 1);
428 for (cnt = 0, k = 0; k < nquad2; ++k)
430 for (j = 0; j < nquad1; ++j)
432 g2 = 2.0 / (1.0 - z2[k]);
433 g1 = g2 / (1.0 - z1[j]);
435 g3 = (1.0 + z1[j]) * g2 * 0.5;
437 for (i = 0; i < nquad0; ++i, ++cnt)
439 g1a = g1 * (1 + z0[i]);
442 g0 * tmp2[cnt] + g1a * (tmp3[cnt] + outarray[2][cnt]);
444 outarray[1][cnt] = g2 * tmp3[cnt] + g3 * outarray[2][cnt];
464 return StdExpansion3D::v_PhysEvaluate(Lcoord, physvals);
479 m_geom->GetLocCoords(coord, Lcoord);
482 return StdExpansion3D::v_PhysEvaluate(Lcoord, physvals);
488 std::array<NekDouble, 3> &firstOrderDerivs)
492 m_geom->GetLocCoords(coord, Lcoord);
493 return StdTetExp::v_PhysEvalFirstDeriv(Lcoord, inarray, firstOrderDerivs);
504 ASSERTL1(Lcoords[0] <= -1.0 && Lcoords[0] >= 1.0 && Lcoords[1] <= -1.0 &&
505 Lcoords[1] >= 1.0 && Lcoords[2] <= -1.0 && Lcoords[2] >= 1.0,
506 "Local coordinates are not in region [-1,1]");
510 for (i = 0; i <
m_geom->GetCoordim(); ++i)
512 coords[i] =
m_geom->GetCoord(i, Lcoords);
539 m_base[2]->GetBasisKey());
545 m_base[0]->GetPointsKey());
547 m_base[1]->GetPointsKey());
549 m_base[2]->GetPointsKey());
556 const NekDouble *data,
const std::vector<unsigned int> &nummodes,
557 const int mode_offset,
NekDouble *coeffs,
558 [[maybe_unused]] std::vector<LibUtilities::BasisType> &fromType)
560 int data_order0 = nummodes[mode_offset];
561 int fillorder0 = min(
m_base[0]->GetNumModes(), data_order0);
562 int data_order1 = nummodes[mode_offset + 1];
563 int order1 =
m_base[1]->GetNumModes();
564 int fillorder1 = min(order1, data_order1);
565 int data_order2 = nummodes[mode_offset + 2];
566 int order2 =
m_base[2]->GetNumModes();
567 int fillorder2 = min(order2, data_order2);
578 "Extraction routine not set up for this basis");
580 "Extraction routine not set up for this basis");
583 for (j = 0; j < fillorder0; ++j)
585 for (i = 0; i < fillorder1 - j; ++i)
589 cnt += data_order2 - j - i;
590 cnt1 += order2 - j - i;
594 for (i = fillorder1 - j; i < data_order1 - j; ++i)
596 cnt += data_order2 - j - i;
599 for (i = fillorder1 - j; i < order1 - j; ++i)
601 cnt1 += order2 - j - i;
607 ASSERTL0(
false,
"basis is either not set up or not "
617 int nquad0 =
m_base[0]->GetNumPoints();
618 int nquad1 =
m_base[1]->GetNumPoints();
619 int nquad2 =
m_base[2]->GetNumPoints();
631 if (outarray.size() != nq0 * nq1)
636 for (
int i = 0; i < nquad0 * nquad1; ++i)
647 if (outarray.size() != nq0 * nq1)
653 for (
int k = 0; k < nquad2; k++)
655 for (
int i = 0; i < nquad0; ++i)
657 outarray[k * nquad0 + i] = (nquad0 * nquad1 * k) + i;
666 if (outarray.size() != nq0 * nq1)
672 for (
int j = 0; j < nquad1 * nquad2; ++j)
674 outarray[j] = nquad0 - 1 + j * nquad0;
682 if (outarray.size() != nq0 * nq1)
688 for (
int j = 0; j < nquad1 * nquad2; ++j)
690 outarray[j] = j * nquad0;
695 ASSERTL0(
false,
"face value (> 3) is out of range");
710 for (
int i = 0; i < ptsKeys.size(); ++i)
722 geomFactors->GetDerivFactors(ptsKeys);
735 for (i = 0; i < vCoordDim; ++i)
740 size_t nqb = nq_face;
756 for (i = 0; i < vCoordDim; ++i)
758 normal[i][0] = -df[3 * i + 2][0];
765 for (i = 0; i < vCoordDim; ++i)
767 normal[i][0] = -df[3 * i + 1][0];
774 for (i = 0; i < vCoordDim; ++i)
777 df[3 * i][0] + df[3 * i + 1][0] + df[3 * i + 2][0];
784 for (i = 0; i < vCoordDim; ++i)
786 normal[i][0] = -df[3 * i][0];
791 ASSERTL0(
false,
"face is out of range (edge < 3)");
796 for (i = 0; i < vCoordDim; ++i)
798 fac += normal[i][0] * normal[i][0];
800 fac = 1.0 /
sqrt(fac);
803 for (i = 0; i < vCoordDim; ++i)
805 Vmath::Fill(nq_face, fac * normal[i][0], normal[i], 1);
813 int nq0 = ptsKeys[0].GetNumPoints();
814 int nq1 = ptsKeys[1].GetNumPoints();
815 int nq2 = ptsKeys[2].GetNumPoints();
817 int nq01 = nq0 * nq1;
846 for (j = 0; j < nq01; ++j)
848 normals[j] = -df[2][j] * jac[j];
849 normals[nqtot + j] = -df[5][j] * jac[j];
850 normals[2 * nqtot + j] = -df[8][j] * jac[j];
854 points0 = ptsKeys[0];
855 points1 = ptsKeys[1];
861 for (j = 0; j < nq0; ++j)
863 for (k = 0; k < nq2; ++k)
865 int tmp = j + nq01 * k;
866 normals[j + k * nq0] = -df[1][tmp] * jac[tmp];
867 normals[nqtot + j + k * nq0] = -df[4][tmp] * jac[tmp];
868 normals[2 * nqtot + j + k * nq0] =
869 -df[7][tmp] * jac[tmp];
870 faceJac[j + k * nq0] = jac[tmp];
874 points0 = ptsKeys[0];
875 points1 = ptsKeys[2];
881 for (j = 0; j < nq1; ++j)
883 for (k = 0; k < nq2; ++k)
885 int tmp = nq0 - 1 + nq0 * j + nq01 * k;
886 normals[j + k * nq1] =
887 (df[0][tmp] + df[1][tmp] + df[2][tmp]) * jac[tmp];
888 normals[nqtot + j + k * nq1] =
889 (df[3][tmp] + df[4][tmp] + df[5][tmp]) * jac[tmp];
890 normals[2 * nqtot + j + k * nq1] =
891 (df[6][tmp] + df[7][tmp] + df[8][tmp]) * jac[tmp];
892 faceJac[j + k * nq1] = jac[tmp];
896 points0 = ptsKeys[1];
897 points1 = ptsKeys[2];
903 for (j = 0; j < nq1; ++j)
905 for (k = 0; k < nq2; ++k)
907 int tmp = j * nq0 + nq01 * k;
908 normals[j + k * nq1] = -df[0][tmp] * jac[tmp];
909 normals[nqtot + j + k * nq1] = -df[3][tmp] * jac[tmp];
910 normals[2 * nqtot + j + k * nq1] =
911 -df[6][tmp] * jac[tmp];
912 faceJac[j + k * nq1] = jac[tmp];
916 points0 = ptsKeys[1];
917 points1 = ptsKeys[2];
922 ASSERTL0(
false,
"face is out of range (face < 3)");
930 Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
933 for (i = 0; i < vCoordDim; ++i)
938 Vmath::Vmul(nq_face, work, 1, normal[i], 1, normal[i], 1);
945 Vmath::Vvtvp(nq_face, normal[i], 1, normal[i], 1, work, 1, work, 1);
955 Vmath::Vmul(nq_face, normal[i], 1, work, 1, normal[i], 1);
982 StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
1006 StdTetExp::v_SVVLaplacianFilter(array, mkey);
1031 returnval = StdTetExp::v_GenMatrix(mkey);
1045 return tmp->GetStdMatrix(mkey);
1074 if (inarray.data() == outarray.data())
1080 (mat->GetOwnedMatrix())->GetPtr().data(),
m_ncoeffs,
1081 tmp.data(), 1, 0.0, outarray.data(), 1);
1086 (mat->GetOwnedMatrix())->GetPtr().data(),
m_ncoeffs,
1087 inarray.data(), 1, 0.0, outarray.data(), 1);
1102 int nquad0 =
m_base[0]->GetNumPoints();
1103 int nquad1 =
m_base[1]->GetNumPoints();
1104 int nquad2 =
m_base[2]->GetNumPoints();
1105 int nqtot = nquad0 * nquad1 * nquad2;
1107 ASSERTL1(wsp.size() >= 6 * nqtot,
"Insufficient workspace size.");
1108 ASSERTL1(m_ncoeffs <= nqtot, "Workspace not set up for ncoeffs > nqtot
");
1110 const Array<OneD, const NekDouble> &base0 = m_base[0]->GetBdata();
1111 const Array<OneD, const NekDouble> &base1 = m_base[1]->GetBdata();
1112 const Array<OneD, const NekDouble> &base2 = m_base[2]->GetBdata();
1113 const Array<OneD, const NekDouble> &dbase0 = m_base[0]->GetDbdata();
1114 const Array<OneD, const NekDouble> &dbase1 = m_base[1]->GetDbdata();
1115 const Array<OneD, const NekDouble> &dbase2 = m_base[2]->GetDbdata();
1116 const Array<OneD, const NekDouble> &metric00 =
1117 m_metrics[eMetricLaplacian00];
1118 const Array<OneD, const NekDouble> &metric01 =
1119 m_metrics[eMetricLaplacian01];
1120 const Array<OneD, const NekDouble> &metric02 =
1121 m_metrics[eMetricLaplacian02];
1122 const Array<OneD, const NekDouble> &metric11 =
1123 m_metrics[eMetricLaplacian11];
1124 const Array<OneD, const NekDouble> &metric12 =
1125 m_metrics[eMetricLaplacian12];
1126 const Array<OneD, const NekDouble> &metric22 =
1127 m_metrics[eMetricLaplacian22];
1129 // Allocate temporary storage
1130 Array<OneD, NekDouble> wsp0(2 * nqtot, wsp);
1131 Array<OneD, NekDouble> wsp1(nqtot, wsp + 1 * nqtot);
1132 Array<OneD, NekDouble> wsp2(nqtot, wsp + 2 * nqtot);
1133 Array<OneD, NekDouble> wsp3(nqtot, wsp + 3 * nqtot);
1134 Array<OneD, NekDouble> wsp4(nqtot, wsp + 4 * nqtot);
1135 Array<OneD, NekDouble> wsp5(nqtot, wsp + 5 * nqtot);
1137 // LAPLACIAN MATRIX OPERATION
1138 // wsp1 = du_dxi1 = D_xi1 * inarray = D_xi1 * u
1139 // wsp2 = du_dxi2 = D_xi2 * inarray = D_xi2 * u
1140 // wsp2 = du_dxi3 = D_xi3 * inarray = D_xi3 * u
1141 StdExpansion3D::PhysTensorDeriv(inarray, wsp0, wsp1, wsp2);
1143 // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1144 // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1145 // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
1146 // especially for this purpose
1147 Vmath::Vvtvvtp(nqtot, &metric00[0], 1, &wsp0[0], 1, &metric01[0], 1,
1148 &wsp1[0], 1, &wsp3[0], 1);
1149 Vmath::Vvtvp(nqtot, &metric02[0], 1, &wsp2[0], 1, &wsp3[0], 1, &wsp3[0], 1);
1150 Vmath::Vvtvvtp(nqtot, &metric01[0], 1, &wsp0[0], 1, &metric11[0], 1,
1151 &wsp1[0], 1, &wsp4[0], 1);
1152 Vmath::Vvtvp(nqtot, &metric12[0], 1, &wsp2[0], 1, &wsp4[0], 1, &wsp4[0], 1);
1153 Vmath::Vvtvvtp(nqtot, &metric02[0], 1, &wsp0[0], 1, &metric12[0], 1,
1154 &wsp1[0], 1, &wsp5[0], 1);
1155 Vmath::Vvtvp(nqtot, &metric22[0], 1, &wsp2[0], 1, &wsp5[0], 1, &wsp5[0], 1);
1157 // outarray = m = (D_xi1 * B)^T * k
1158 // wsp1 = n = (D_xi2 * B)^T * l
1159 IProductWRTBase_SumFacKernel(dbase0, base1, base2, wsp3, outarray, wsp0,
1161 IProductWRTBase_SumFacKernel(base0, dbase1, base2, wsp4, wsp2, wsp0, true,
1163 Vmath::Vadd(m_ncoeffs, wsp2.data(), 1, outarray.data(), 1, outarray.data(),
1165 IProductWRTBase_SumFacKernel(base0, base1, dbase2, wsp5, wsp2, wsp0, true,
1167 Vmath::Vadd(m_ncoeffs, wsp2.data(), 1, outarray.data(), 1, outarray.data(),
1171void TetExp::v_ComputeLaplacianMetric()
1173 if (m_metrics.count(eMetricQuadrature) == 0)
1175 ComputeQuadratureMetric();
1179 const unsigned int nqtot = GetTotPoints();
1180 const unsigned int dim = 3;
1181 const MetricType m[3][3] = {
1182 {eMetricLaplacian00, eMetricLaplacian01, eMetricLaplacian02},
1183 {eMetricLaplacian01, eMetricLaplacian11, eMetricLaplacian12},
1184 {eMetricLaplacian02, eMetricLaplacian12, eMetricLaplacian22}};
1186 for (unsigned int i = 0; i < dim; ++i)
1188 for (unsigned int j = i; j < dim; ++j)
1190 m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
1194 // Define shorthand synonyms for m_metrics storage
1195 Array<OneD, NekDouble> g0(m_metrics[m[0][0]]);
1196 Array<OneD, NekDouble> g1(m_metrics[m[1][1]]);
1197 Array<OneD, NekDouble> g2(m_metrics[m[2][2]]);
1198 Array<OneD, NekDouble> g3(m_metrics[m[0][1]]);
1199 Array<OneD, NekDouble> g4(m_metrics[m[0][2]]);
1200 Array<OneD, NekDouble> g5(m_metrics[m[1][2]]);
1202 // Allocate temporary storage
1203 Array<OneD, NekDouble> alloc(7 * nqtot, 0.0);
1204 Array<OneD, NekDouble> h0(alloc); // h0
1205 Array<OneD, NekDouble> h1(alloc + 1 * nqtot); // h1
1206 Array<OneD, NekDouble> h2(alloc + 2 * nqtot); // h2
1207 Array<OneD, NekDouble> h3(alloc + 3 * nqtot); // h3
1208 Array<OneD, NekDouble> wsp4(alloc + 4 * nqtot); // wsp4
1209 Array<OneD, NekDouble> wsp5(alloc + 5 * nqtot); // wsp5
1210 Array<OneD, NekDouble> wsp6(alloc + 6 * nqtot); // wsp6
1211 // Reuse some of the storage as workspace
1212 Array<OneD, NekDouble> wsp7(alloc); // wsp7
1213 Array<OneD, NekDouble> wsp8(alloc + 1 * nqtot); // wsp8
1214 Array<OneD, NekDouble> wsp9(alloc + 2 * nqtot); // wsp9
1216 const Array<TwoD, const NekDouble> &df =
1217 m_metricinfo->GetDerivFactors(GetPointsKeys());
1218 const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
1219 const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
1220 const Array<OneD, const NekDouble> &z2 = m_base[2]->GetZ();
1221 const unsigned int nquad0 = m_base[0]->GetNumPoints();
1222 const unsigned int nquad1 = m_base[1]->GetNumPoints();
1223 const unsigned int nquad2 = m_base[2]->GetNumPoints();
1225 for (j = 0; j < nquad2; ++j)
1227 for (i = 0; i < nquad1; ++i)
1229 Vmath::Fill(nquad0, 4.0 / (1.0 - z1[i]) / (1.0 - z2[j]),
1230 &h0[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1231 Vmath::Fill(nquad0, 2.0 / (1.0 - z1[i]) / (1.0 - z2[j]),
1232 &h1[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1233 Vmath::Fill(nquad0, 2.0 / (1.0 - z2[j]),
1234 &h2[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1235 Vmath::Fill(nquad0, (1.0 + z1[i]) / (1.0 - z2[j]),
1236 &h3[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1239 for (i = 0; i < nquad0; i++)
1241 Blas::Dscal(nquad1 * nquad2, 1 + z0[i], &h1[0] + i, nquad0);
1244 // Step 3. Construct combined metric terms for physical space to
1245 // collapsed coordinate system.
1246 // Order of construction optimised to minimise temporary storage
1247 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1250 Vmath::Vadd(nqtot, &df[1][0], 1, &df[2][0], 1, &wsp4[0], 1);
1251 Vmath::Vvtvvtp(nqtot, &df[0][0], 1, &h0[0], 1, &wsp4[0], 1, &h1[0], 1,
1254 Vmath::Vadd(nqtot, &df[4][0], 1, &df[5][0], 1, &wsp5[0], 1);
1255 Vmath::Vvtvvtp(nqtot, &df[3][0], 1, &h0[0], 1, &wsp5[0], 1, &h1[0], 1,
1258 Vmath::Vadd(nqtot, &df[7][0], 1, &df[8][0], 1, &wsp6[0], 1);
1259 Vmath::Vvtvvtp(nqtot, &df[6][0], 1, &h0[0], 1, &wsp6[0], 1, &h1[0], 1,
1263 Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0],
1265 Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1268 Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp4[0], 1, &df[5][0], 1, &wsp5[0],
1270 Vmath::Vvtvp(nqtot, &df[8][0], 1, &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1272 // overwrite h0, h1, h2
1273 // wsp7 (h2f1 + h3f2)
1274 Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &h2[0], 1, &df[2][0], 1, &h3[0], 1,
1276 // wsp8 (h2f4 + h3f5)
1277 Vmath::Vvtvvtp(nqtot, &df[4][0], 1, &h2[0], 1, &df[5][0], 1, &h3[0], 1,
1279 // wsp9 (h2f7 + h3f8)
1280 Vmath::Vvtvvtp(nqtot, &df[7][0], 1, &h2[0], 1, &df[8][0], 1, &h3[0], 1,
1284 Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp7[0], 1, &wsp5[0], 1, &wsp8[0],
1286 Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp9[0], 1, &g3[0], 1, &g3[0], 1);
1288 // overwrite wsp4, wsp5, wsp6
1290 Vmath::Vvtvvtp(nqtot, &wsp7[0], 1, &wsp7[0], 1, &wsp8[0], 1, &wsp8[0],
1292 Vmath::Vvtvp(nqtot, &wsp9[0], 1, &wsp9[0], 1, &g1[0], 1, &g1[0], 1);
1295 Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp7[0], 1, &df[5][0], 1, &wsp8[0],
1297 Vmath::Vvtvp(nqtot, &df[8][0], 1, &wsp9[0], 1, &g5[0], 1, &g5[0], 1);
1300 Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &df[2][0], 1, &df[5][0], 1,
1301 &df[5][0], 1, &g2[0], 1);
1302 Vmath::Vvtvp(nqtot, &df[8][0], 1, &df[8][0], 1, &g2[0], 1, &g2[0], 1);
1307 Vmath::Svtsvtp(nqtot, df[0][0], &h0[0], 1, df[1][0] + df[2][0], &h1[0],
1310 Vmath::Svtsvtp(nqtot, df[3][0], &h0[0], 1, df[4][0] + df[5][0], &h1[0],
1313 Vmath::Svtsvtp(nqtot, df[6][0], &h0[0], 1, df[7][0] + df[8][0], &h1[0],
1317 Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0],
1319 Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1322 Vmath::Svtsvtp(nqtot, df[2][0], &wsp4[0], 1, df[5][0], &wsp5[0], 1,
1324 Vmath::Svtvp(nqtot, df[8][0], &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1326 // overwrite h0, h1, h2
1327 // wsp7 (h2f1 + h3f2)
1328 Vmath::Svtsvtp(nqtot, df[1][0], &h2[0], 1, df[2][0], &h3[0], 1,
1330 // wsp8 (h2f4 + h3f5)
1331 Vmath::Svtsvtp(nqtot, df[4][0], &h2[0], 1, df[5][0], &h3[0], 1,
1333 // wsp9 (h2f7 + h3f8)
1334 Vmath::Svtsvtp(nqtot, df[7][0], &h2[0], 1, df[8][0], &h3[0], 1,
1338 Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp7[0], 1, &wsp5[0], 1, &wsp8[0],
1340 Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp9[0], 1, &g3[0], 1, &g3[0], 1);
1342 // overwrite wsp4, wsp5, wsp6
1344 Vmath::Vvtvvtp(nqtot, &wsp7[0], 1, &wsp7[0], 1, &wsp8[0], 1, &wsp8[0],
1346 Vmath::Vvtvp(nqtot, &wsp9[0], 1, &wsp9[0], 1, &g1[0], 1, &g1[0], 1);
1349 Vmath::Svtsvtp(nqtot, df[2][0], &wsp7[0], 1, df[5][0], &wsp8[0], 1,
1351 Vmath::Svtvp(nqtot, df[8][0], &wsp9[0], 1, &g5[0], 1, &g5[0], 1);
1355 df[2][0] * df[2][0] + df[5][0] * df[5][0] +
1356 df[8][0] * df[8][0],
1360 for (unsigned int i = 0; i < dim; ++i)
1362 for (unsigned int j = i; j < dim; ++j)
1364 MultiplyByQuadratureMetric(m_metrics[m[i][j]], m_metrics[m[i][j]]);
1374void TetExp::v_NormalTraceDerivFactors(
1375 Array<OneD, Array<OneD, NekDouble>> &d0factors,
1376 Array<OneD, Array<OneD, NekDouble>> &d1factors,
1377 Array<OneD, Array<OneD, NekDouble>> &d2factors)
1379 int nquad0 = GetNumPoints(0);
1380 int nquad1 = GetNumPoints(1);
1381 int nquad2 = GetNumPoints(2);
1383 const Array<TwoD, const NekDouble> &df =
1384 m_metricinfo->GetDerivFactors(GetPointsKeys());
1386 if (d0factors.size() != 4)
1388 d0factors = Array<OneD, Array<OneD, NekDouble>>(4);
1389 d1factors = Array<OneD, Array<OneD, NekDouble>>(4);
1390 d2factors = Array<OneD, Array<OneD, NekDouble>>(4);
1393 if (d0factors[0].size() != nquad0 * nquad1)
1395 d0factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1396 d1factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1397 d2factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1400 if (d0factors[1].size() != nquad0 * nquad2)
1402 d0factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1403 d1factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1404 d2factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1407 if (d0factors[2].size() != nquad1 * nquad2)
1409 d0factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1410 d0factors[3] = Array<OneD, NekDouble>(nquad1 * nquad2);
1411 d1factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1412 d1factors[3] = Array<OneD, NekDouble>(nquad1 * nquad2);
1413 d2factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1414 d2factors[3] = Array<OneD, NekDouble>(nquad1 * nquad2);
1418 const Array<OneD, const Array<OneD, NekDouble>> &normal_0 =
1420 const Array<OneD, const Array<OneD, NekDouble>> &normal_1 =
1422 const Array<OneD, const Array<OneD, NekDouble>> &normal_2 =
1424 const Array<OneD, const Array<OneD, NekDouble>> &normal_3 =
1427 int ncoords = normal_0.size();
1429 // first gather together standard cartesian inner products
1430 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1433 for (int i = 0; i < nquad0 * nquad1; ++i)
1435 d0factors[0][i] = df[0][i] * normal_0[0][i];
1436 d1factors[0][i] = df[1][i] * normal_0[0][i];
1437 d2factors[0][i] = df[2][i] * normal_0[0][i];
1440 for (int n = 1; n < ncoords; ++n)
1442 for (int i = 0; i < nquad0 * nquad1; ++i)
1444 d0factors[0][i] += df[3 * n][i] * normal_0[n][i];
1445 d1factors[0][i] += df[3 * n + 1][i] * normal_0[n][i];
1446 d2factors[0][i] += df[3 * n + 2][i] * normal_0[n][i];
1451 for (int j = 0; j < nquad2; ++j)
1453 for (int i = 0; i < nquad0; ++i)
1455 d0factors[1][j * nquad0 + i] = df[0][j * nquad0 * nquad1 + i] *
1456 normal_1[0][j * nquad0 + i];
1457 d1factors[1][j * nquad0 + i] = df[1][j * nquad0 * nquad1 + i] *
1458 normal_1[0][j * nquad0 + i];
1459 d2factors[1][j * nquad0 + i] = df[2][j * nquad0 * nquad1 + i] *
1460 normal_1[0][j * nquad0 + i];
1464 for (int n = 1; n < ncoords; ++n)
1466 for (int j = 0; j < nquad2; ++j)
1468 for (int i = 0; i < nquad0; ++i)
1470 d0factors[1][j * nquad0 + i] +=
1471 df[3 * n][j * nquad0 * nquad1 + i] *
1472 normal_1[0][j * nquad0 + i];
1473 d1factors[1][j * nquad0 + i] +=
1474 df[3 * n + 1][j * nquad0 * nquad1 + i] *
1475 normal_1[0][j * nquad0 + i];
1476 d2factors[1][j * nquad0 + i] +=
1477 df[3 * n + 2][j * nquad0 * nquad1 + i] *
1478 normal_1[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 * nquad1 + i + 1) * nquad0 - 1] *
1490 normal_2[0][j * nquad1 + i];
1491 d1factors[2][j * nquad1 + i] =
1492 df[1][(j * nquad1 + i + 1) * nquad0 - 1] *
1493 normal_2[0][j * nquad1 + i];
1494 d2factors[2][j * nquad1 + i] =
1495 df[2][(j * nquad1 + i + 1) * nquad0 - 1] *
1496 normal_2[0][j * nquad1 + i];
1498 d0factors[3][j * nquad1 + i] =
1499 df[0][(j * nquad1 + i) * nquad0] *
1500 normal_3[0][j * nquad1 + i];
1501 d1factors[3][j * nquad1 + i] =
1502 df[1][(j * nquad1 + i) * nquad0] *
1503 normal_3[0][j * nquad1 + i];
1504 d2factors[3][j * nquad1 + i] =
1505 df[2][(j * nquad1 + i) * nquad0] *
1506 normal_3[0][j * nquad1 + i];
1510 for (int n = 1; n < ncoords; ++n)
1512 for (int j = 0; j < nquad2; ++j)
1514 for (int i = 0; i < nquad1; ++i)
1516 d0factors[2][j * nquad1 + i] +=
1517 df[3 * n][(j * nquad1 + i + 1) * nquad0 - 1] *
1518 normal_2[n][j * nquad1 + i];
1519 d1factors[2][j * nquad1 + i] +=
1520 df[3 * n + 1][(j * nquad1 + i + 1) * nquad0 - 1] *
1521 normal_2[n][j * nquad1 + i];
1522 d2factors[2][j * nquad1 + i] +=
1523 df[3 * n + 2][(j * nquad1 + i + 1) * nquad0 - 1] *
1524 normal_2[n][j * nquad1 + i];
1526 d0factors[3][j * nquad1 + i] +=
1527 df[3 * n][(j * nquad1 + i) * nquad0] *
1528 normal_3[n][j * nquad1 + i];
1529 d1factors[3][j * nquad1 + i] +=
1530 df[3 * n + 1][(j * nquad1 + i) * nquad0] *
1531 normal_3[n][j * nquad1 + i];
1532 d2factors[3][j * nquad1 + i] +=
1533 df[3 * n + 2][(j * nquad1 + i) * nquad0] *
1534 normal_3[n][j * nquad1 + 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 d1factors[1][i] = df[1][0] * normal_1[0][i];
1564 d2factors[1][i] = df[2][0] * normal_1[0][i];
1567 for (int n = 1; n < ncoords; ++n)
1569 for (int i = 0; i < nquad0 * nquad2; ++i)
1571 d0factors[1][i] += df[3 * n][0] * normal_1[n][i];
1572 d1factors[1][i] += df[3 * n + 1][0] * normal_1[n][i];
1573 d2factors[1][i] += df[3 * n + 2][0] * normal_1[n][i];
1578 for (int i = 0; i < nquad1 * nquad2; ++i)
1580 d0factors[2][i] = df[0][0] * normal_2[0][i];
1581 d0factors[3][i] = df[0][0] * normal_3[0][i];
1583 d1factors[2][i] = df[1][0] * normal_2[0][i];
1584 d1factors[3][i] = df[1][0] * normal_3[0][i];
1586 d2factors[2][i] = df[2][0] * normal_2[0][i];
1587 d2factors[3][i] = df[2][0] * normal_3[0][i];
1590 for (int n = 1; n < ncoords; ++n)
1592 for (int i = 0; i < nquad1 * nquad2; ++i)
1594 d0factors[2][i] += df[3 * n][0] * normal_2[n][i];
1595 d0factors[3][i] += df[3 * n][0] * normal_3[n][i];
1597 d1factors[2][i] += df[3 * n + 1][0] * normal_2[n][i];
1598 d1factors[3][i] += df[3 * n + 1][0] * normal_3[n][i];
1600 d2factors[2][i] += df[3 * n + 2][0] * normal_2[n][i];
1601 d2factors[3][i] += df[3 * n + 2][0] * normal_3[n][i];
1606} // namespace Nektar::LocalRegions
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Describes the specification for a Basis.
int GetNumPoints() const
Return points order at which basis is defined.
PointsKey GetPointsKey() const
Return distribution of points.
Defines a specification for a set of points.
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
std::map< int, NormalVector > m_traceNormals
std::map< int, Array< OneD, NekDouble > > m_elmtBndNormDirElmtLen
the element length in each element boundary(Vertex, edge or face) normal direction calculated based o...
SpatialDomains::GeometrySharedPtr GetGeom() const
SpatialDomains::GeometrySharedPtr m_geom
void ComputeLaplacianMetric()
void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
void v_ComputeTraceNormal(const int face) override
Compute the normal of a triangular face.
NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &physvals) override
NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray) override
Integrate the physical point list inarray over region.
DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey) override
NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals) 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) override
Differentiate inarray in the three coordinate directions.
void v_DropLocMatrix(const MatrixKey &mkey) override
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true) override
NekDouble v_PhysEvalFirstDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) override
DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey) override
void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords) override
Get the coordinates "coords" at the local coordinates "Lcoords".
void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Calculate the inner product of inarray with respect to the basis B=m_base0*m_base1*m_base2 and put in...
void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Forward transform from physical quadrature space stored in inarray and evaluate the expansion coeffic...
StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const override
LibUtilities::ShapeType v_DetShapeType() const override
Return Shape of region, using ShapeType enum list.
DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey) override
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Calculates the inner product .
StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const override
TetExp(const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const LibUtilities::BasisKey &Bc, const SpatialDomains::TetGeomSharedPtr &geom)
Constructor using BasisKey class for quadrature points and order definition.
void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
void GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
void v_DropLocStaticCondMatrix(const MatrixKey &mkey) override
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey) override
void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray) override
void v_GetTracePhysMap(const int face, Array< OneD, int > &outarray) override
Returns the physical values at the quadrature points of a face.
void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp) override
void v_ExtractDataToCoeffs(const NekDouble *data, const std::vector< unsigned int > &nummodes, const int mode_offset, NekDouble *coeffs, std::vector< LibUtilities::BasisType > &fromType) override
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
void IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
void v_HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
void v_LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
const LibUtilities::PointsKeyVector GetPointsKeys() const
void 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::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
const LibUtilities::BasisKey GetTraceBasisKey(const int i, int k=-1, bool UseGLL=false) const
This function returns the basis key belonging to the i-th trace.
Array< OneD, LibUtilities::BasisSharedPtr > m_base
MatrixType GetMatrixType() const
LibUtilities::ShapeType DetShapeType() 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 void Daxpy(const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: y = alpha x plus y.
int getNumberOfCoefficients(int Na)
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis,...
std::vector< PointsKey > PointsKeyVector
@ eModified_B
Principle Modified Functions .
@ eModified_C
Principle Modified Functions .
@ 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< TetGeom > TetGeomSharedPtr
std::shared_ptr< StdTetExp > StdTetExpSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
@ eInvLaplacianWithUnityMean
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
void 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/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 Vcopy(int n, const T *x, const int incx, T *y, const int incy)
scalarT< T > sqrt(scalarT< T > in)