35 #include <boost/core/ignore_unused.hpp>
46 namespace LocalRegions
66 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
69 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
73 std::bind(&
Expansion3D::CreateMatrix, this, std::placeholders::_1),
74 std::string(
"TetExpMatrix")),
75 m_staticCondMatrixManager(std::bind(&
Expansion::CreateStaticCondMatrix,
76 this, std::placeholders::_1),
77 std::string(
"TetExpStaticCondMatrix"))
85 : StdRegions::StdExpansion(T), StdRegions::StdExpansion3D(T),
87 m_matrixManager(T.m_matrixManager),
88 m_staticCondMatrixManager(T.m_staticCondMatrixManager)
115 int nquad0 =
m_base[0]->GetNumPoints();
116 int nquad1 =
m_base[1]->GetNumPoints();
117 int nquad2 =
m_base[2]->GetNumPoints();
126 (
NekDouble *)&inarray[0], 1, &tmp[0], 1);
131 (
NekDouble *)&inarray[0], 1, &tmp[0], 1);
135 retrunVal = StdTetExp::v_Integral(tmp);
157 int TotPts =
m_base[0]->GetNumPoints() *
m_base[1]->GetNumPoints() *
158 m_base[2]->GetNumPoints();
166 StdTetExp::v_PhysDeriv(inarray, Diff0, Diff1, Diff2);
172 Vmath::Vmul(TotPts, &df[0][0], 1, &Diff0[0], 1, &out_d0[0], 1);
173 Vmath::Vvtvp(TotPts, &df[1][0], 1, &Diff1[0], 1, &out_d0[0], 1,
175 Vmath::Vvtvp(TotPts, &df[2][0], 1, &Diff2[0], 1, &out_d0[0], 1,
181 Vmath::Vmul(TotPts, &df[3][0], 1, &Diff0[0], 1, &out_d1[0], 1);
182 Vmath::Vvtvp(TotPts, &df[4][0], 1, &Diff1[0], 1, &out_d1[0], 1,
184 Vmath::Vvtvp(TotPts, &df[5][0], 1, &Diff2[0], 1, &out_d1[0], 1,
190 Vmath::Vmul(TotPts, &df[6][0], 1, &Diff0[0], 1, &out_d2[0], 1);
191 Vmath::Vvtvp(TotPts, &df[7][0], 1, &Diff1[0], 1, &out_d2[0], 1,
193 Vmath::Vvtvp(TotPts, &df[8][0], 1, &Diff2[0], 1, &out_d2[0], 1,
201 Vmath::Smul(TotPts, df[0][0], &Diff0[0], 1, &out_d0[0], 1);
202 Blas::Daxpy(TotPts, df[1][0], &Diff1[0], 1, &out_d0[0], 1);
203 Blas::Daxpy(TotPts, df[2][0], &Diff2[0], 1, &out_d0[0], 1);
208 Vmath::Smul(TotPts, df[3][0], &Diff0[0], 1, &out_d1[0], 1);
209 Blas::Daxpy(TotPts, df[4][0], &Diff1[0], 1, &out_d1[0], 1);
210 Blas::Daxpy(TotPts, df[5][0], &Diff2[0], 1, &out_d1[0], 1);
215 Vmath::Smul(TotPts, df[6][0], &Diff0[0], 1, &out_d2[0], 1);
216 Blas::Daxpy(TotPts, df[7][0], &Diff1[0], 1, &out_d2[0], 1);
217 Blas::Daxpy(TotPts, df[8][0], &Diff2[0], 1, &out_d2[0], 1);
237 if ((
m_base[0]->Collocation()) && (
m_base[1]->Collocation()) &&
238 (
m_base[2]->Collocation()))
254 out = (*matsys) * in;
297 const int nquad0 =
m_base[0]->GetNumPoints();
298 const int nquad1 =
m_base[1]->GetNumPoints();
299 const int nquad2 =
m_base[2]->GetNumPoints();
300 const int order0 =
m_base[0]->GetNumModes();
301 const int order1 =
m_base[1]->GetNumModes();
303 nquad2 * order0 * (order1 + 1) / 2);
305 if (multiplybyweights)
312 tmp, outarray, wsp,
true,
true,
true);
318 inarray, outarray, wsp,
true,
true,
true);
356 const int nquad0 =
m_base[0]->GetNumPoints();
357 const int nquad1 =
m_base[1]->GetNumPoints();
358 const int nquad2 =
m_base[2]->GetNumPoints();
359 const int order0 =
m_base[0]->GetNumModes();
360 const int order1 =
m_base[1]->GetNumModes();
361 const int nqtot = nquad0 * nquad1 * nquad2;
369 nquad2 * order0 * (order1 + 1) / 2);
381 m_base[2]->GetBdata(), tmp2, outarray, wsp,
385 m_base[2]->GetBdata(), tmp3, tmp6, wsp,
true,
391 m_base[2]->GetDbdata(), tmp4, tmp6, wsp,
true,
403 const int nquad0 =
m_base[0]->GetNumPoints();
404 const int nquad1 =
m_base[1]->GetNumPoints();
405 const int nquad2 =
m_base[2]->GetNumPoints();
406 const int nqtot = nquad0 * nquad1 * nquad2;
420 Vmath::Vmul(nqtot, &df[3 * dir][0], 1, inarray.get(), 1, tmp2.get(), 1);
421 Vmath::Vmul(nqtot, &df[3 * dir + 1][0], 1, inarray.get(), 1, tmp3.get(),
423 Vmath::Vmul(nqtot, &df[3 * dir + 2][0], 1, inarray.get(), 1,
424 outarray[2].get(), 1);
428 Vmath::Smul(nqtot, df[3 * dir][0], inarray.get(), 1, tmp2.get(), 1);
429 Vmath::Smul(nqtot, df[3 * dir + 1][0], inarray.get(), 1, tmp3.get(), 1);
430 Vmath::Smul(nqtot, df[3 * dir + 2][0], inarray.get(), 1,
431 outarray[2].get(), 1);
437 for (cnt = 0, k = 0; k < nquad2; ++k)
439 for (j = 0; j < nquad1; ++j)
441 g2 = 2.0 / (1.0 - z2[k]);
442 g1 = g2 / (1.0 - z1[j]);
444 g3 = (1.0 + z1[j]) * g2 * 0.5;
446 for (i = 0; i < nquad0; ++i, ++cnt)
448 g1a = g1 * (1 + z0[i]);
451 g0 * tmp2[cnt] + g1a * (tmp3[cnt] + outarray[2][cnt]);
453 outarray[1][cnt] = g2 * tmp3[cnt] + g3 * outarray[2][cnt];
473 return StdTetExp::v_PhysEvaluate(Lcoord, physvals);
488 m_geom->GetLocCoords(coord, Lcoord);
491 return StdTetExp::v_PhysEvaluate(Lcoord, physvals);
502 ASSERTL1(Lcoords[0] <= -1.0 && Lcoords[0] >= 1.0 && Lcoords[1] <= -1.0 &&
503 Lcoords[1] >= 1.0 && Lcoords[2] <= -1.0 && Lcoords[2] >= 1.0,
504 "Local coordinates are not in region [-1,1]");
508 for (i = 0; i <
m_geom->GetCoordim(); ++i)
510 coords[i] =
m_geom->GetCoord(i, Lcoords);
537 m_base[2]->GetBasisKey());
543 m_base[0]->GetPointsKey());
545 m_base[1]->GetPointsKey());
547 m_base[2]->GetPointsKey());
555 return m_geom->GetCoordim();
559 const NekDouble *data,
const std::vector<unsigned int> &nummodes,
560 const int mode_offset,
NekDouble *coeffs,
561 std::vector<LibUtilities::BasisType> &fromType)
563 boost::ignore_unused(fromType);
565 int data_order0 = nummodes[mode_offset];
566 int fillorder0 = min(
m_base[0]->GetNumModes(), data_order0);
567 int data_order1 = nummodes[mode_offset + 1];
568 int order1 =
m_base[1]->GetNumModes();
569 int fillorder1 = min(order1, data_order1);
570 int data_order2 = nummodes[mode_offset + 2];
571 int order2 =
m_base[2]->GetNumModes();
572 int fillorder2 = min(order2, data_order2);
583 "Extraction routine not set up for this basis");
585 "Extraction routine not set up for this basis");
588 for (j = 0; j < fillorder0; ++j)
590 for (i = 0; i < fillorder1 - j; ++i)
594 cnt += data_order2 - j - i;
595 cnt1 += order2 - j - i;
599 for (i = fillorder1 - j; i < data_order1 - j; ++i)
601 cnt += data_order2 - j - i;
604 for (i = fillorder1 - j; i < order1 - j; ++i)
606 cnt1 += order2 - j - i;
612 ASSERTL0(
false,
"basis is either not set up or not "
622 int nquad0 =
m_base[0]->GetNumPoints();
623 int nquad1 =
m_base[1]->GetNumPoints();
624 int nquad2 =
m_base[2]->GetNumPoints();
636 if (outarray.size() != nq0 * nq1)
641 for (
int i = 0; i < nquad0 * nquad1; ++i)
652 if (outarray.size() != nq0 * nq1)
658 for (
int k = 0; k < nquad2; k++)
660 for (
int i = 0; i < nquad0; ++i)
662 outarray[k * nquad0 + i] = (nquad0 * nquad1 * k) + i;
671 if (outarray.size() != nq0 * nq1)
677 for (
int j = 0; j < nquad1 * nquad2; ++j)
679 outarray[j] = nquad0 - 1 + j * nquad0;
687 if (outarray.size() != nq0 * nq1)
693 for (
int j = 0; j < nquad1 * nquad2; ++j)
695 outarray[j] = j * nquad0;
700 ASSERTL0(
false,
"face value (> 3) is out of range");
715 for (
int i = 0; i < ptsKeys.size(); ++i)
727 geomFactors->GetDerivFactors(ptsKeys);
740 for (i = 0; i < vCoordDim; ++i)
745 size_t nqb = nq_face;
761 for (i = 0; i < vCoordDim; ++i)
763 normal[i][0] = -df[3 * i + 2][0];
770 for (i = 0; i < vCoordDim; ++i)
772 normal[i][0] = -df[3 * i + 1][0];
779 for (i = 0; i < vCoordDim; ++i)
782 df[3 * i][0] + df[3 * i + 1][0] + df[3 * i + 2][0];
789 for (i = 0; i < vCoordDim; ++i)
791 normal[i][0] = -df[3 * i][0];
796 ASSERTL0(
false,
"face is out of range (edge < 3)");
801 for (i = 0; i < vCoordDim; ++i)
803 fac += normal[i][0] * normal[i][0];
805 fac = 1.0 /
sqrt(fac);
808 for (i = 0; i < vCoordDim; ++i)
810 Vmath::Fill(nq_face, fac * normal[i][0], normal[i], 1);
818 int nq0 = ptsKeys[0].GetNumPoints();
819 int nq1 = ptsKeys[1].GetNumPoints();
820 int nq2 = ptsKeys[2].GetNumPoints();
822 int nq01 = nq0 * nq1;
851 for (j = 0; j < nq01; ++j)
853 normals[j] = -df[2][j] * jac[j];
854 normals[nqtot + j] = -df[5][j] * jac[j];
855 normals[2 * nqtot + j] = -df[8][j] * jac[j];
859 points0 = ptsKeys[0];
860 points1 = ptsKeys[1];
866 for (j = 0; j < nq0; ++j)
868 for (k = 0; k < nq2; ++k)
870 int tmp = j + nq01 * k;
871 normals[j + k * nq0] = -df[1][tmp] * jac[tmp];
872 normals[nqtot + j + k * nq0] = -df[4][tmp] * jac[tmp];
873 normals[2 * nqtot + j + k * nq0] =
874 -df[7][tmp] * jac[tmp];
875 faceJac[j + k * nq0] = jac[tmp];
879 points0 = ptsKeys[0];
880 points1 = ptsKeys[2];
886 for (j = 0; j < nq1; ++j)
888 for (k = 0; k < nq2; ++k)
890 int tmp = nq0 - 1 + nq0 * j + nq01 * k;
891 normals[j + k * nq1] =
892 (df[0][tmp] + df[1][tmp] + df[2][tmp]) * jac[tmp];
893 normals[nqtot + j + k * nq1] =
894 (df[3][tmp] + df[4][tmp] + df[5][tmp]) * jac[tmp];
895 normals[2 * nqtot + j + k * nq1] =
896 (df[6][tmp] + df[7][tmp] + df[8][tmp]) * jac[tmp];
897 faceJac[j + k * nq1] = jac[tmp];
901 points0 = ptsKeys[1];
902 points1 = ptsKeys[2];
908 for (j = 0; j < nq1; ++j)
910 for (k = 0; k < nq2; ++k)
912 int tmp = j * nq0 + nq01 * k;
913 normals[j + k * nq1] = -df[0][tmp] * jac[tmp];
914 normals[nqtot + j + k * nq1] = -df[3][tmp] * jac[tmp];
915 normals[2 * nqtot + j + k * nq1] =
916 -df[6][tmp] * jac[tmp];
917 faceJac[j + k * nq1] = jac[tmp];
921 points0 = ptsKeys[1];
922 points1 = ptsKeys[2];
927 ASSERTL0(
false,
"face is out of range (face < 3)");
935 Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
938 for (i = 0; i < vCoordDim; ++i)
943 Vmath::Vmul(nq_face, work, 1, normal[i], 1, normal[i], 1);
950 Vmath::Vvtvp(nq_face, normal[i], 1, normal[i], 1, work, 1, work, 1);
960 Vmath::Vmul(nq_face, normal[i], 1, work, 1, normal[i], 1);
987 StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
1011 StdTetExp::v_SVVLaplacianFilter(array, mkey);
1036 returnval = StdTetExp::v_GenMatrix(mkey);
1050 return tmp->GetStdMatrix(mkey);
1074 if (inarray.get() == outarray.get())
1080 (mat->GetOwnedMatrix())->GetPtr().get(),
m_ncoeffs,
1081 tmp.get(), 1, 0.0, outarray.get(), 1);
1086 (mat->GetOwnedMatrix())->GetPtr().get(),
m_ncoeffs,
1087 inarray.get(), 1, 0.0, outarray.get(), 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.get(), 1, outarray.get(), 1, outarray.get(), 1);
1164 IProductWRTBase_SumFacKernel(base0, base1, dbase2, wsp5, wsp2, wsp0, true,
1166 Vmath::Vadd(m_ncoeffs, wsp2.get(), 1, outarray.get(), 1, outarray.get(), 1);
1169 void TetExp::v_ComputeLaplacianMetric()
1171 if (m_metrics.count(eMetricQuadrature) == 0)
1173 ComputeQuadratureMetric();
1177 const unsigned int nqtot = GetTotPoints();
1178 const unsigned int dim = 3;
1179 const MetricType m[3][3] = {
1180 {eMetricLaplacian00, eMetricLaplacian01, eMetricLaplacian02},
1181 {eMetricLaplacian01, eMetricLaplacian11, eMetricLaplacian12},
1182 {eMetricLaplacian02, eMetricLaplacian12, eMetricLaplacian22}};
1184 for (unsigned int i = 0; i < dim; ++i)
1186 for (unsigned int j = i; j < dim; ++j)
1188 m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
1192 // Define shorthand synonyms for m_metrics storage
1193 Array<OneD, NekDouble> g0(m_metrics[m[0][0]]);
1194 Array<OneD, NekDouble> g1(m_metrics[m[1][1]]);
1195 Array<OneD, NekDouble> g2(m_metrics[m[2][2]]);
1196 Array<OneD, NekDouble> g3(m_metrics[m[0][1]]);
1197 Array<OneD, NekDouble> g4(m_metrics[m[0][2]]);
1198 Array<OneD, NekDouble> g5(m_metrics[m[1][2]]);
1200 // Allocate temporary storage
1201 Array<OneD, NekDouble> alloc(7 * nqtot, 0.0);
1202 Array<OneD, NekDouble> h0(alloc); // h0
1203 Array<OneD, NekDouble> h1(alloc + 1 * nqtot); // h1
1204 Array<OneD, NekDouble> h2(alloc + 2 * nqtot); // h2
1205 Array<OneD, NekDouble> h3(alloc + 3 * nqtot); // h3
1206 Array<OneD, NekDouble> wsp4(alloc + 4 * nqtot); // wsp4
1207 Array<OneD, NekDouble> wsp5(alloc + 5 * nqtot); // wsp5
1208 Array<OneD, NekDouble> wsp6(alloc + 6 * nqtot); // wsp6
1209 // Reuse some of the storage as workspace
1210 Array<OneD, NekDouble> wsp7(alloc); // wsp7
1211 Array<OneD, NekDouble> wsp8(alloc + 1 * nqtot); // wsp8
1212 Array<OneD, NekDouble> wsp9(alloc + 2 * nqtot); // wsp9
1214 const Array<TwoD, const NekDouble> &df =
1215 m_metricinfo->GetDerivFactors(GetPointsKeys());
1216 const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
1217 const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
1218 const Array<OneD, const NekDouble> &z2 = m_base[2]->GetZ();
1219 const unsigned int nquad0 = m_base[0]->GetNumPoints();
1220 const unsigned int nquad1 = m_base[1]->GetNumPoints();
1221 const unsigned int nquad2 = m_base[2]->GetNumPoints();
1223 for (j = 0; j < nquad2; ++j)
1225 for (i = 0; i < nquad1; ++i)
1227 Vmath::Fill(nquad0, 4.0 / (1.0 - z1[i]) / (1.0 - z2[j]),
1228 &h0[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1229 Vmath::Fill(nquad0, 2.0 / (1.0 - z1[i]) / (1.0 - z2[j]),
1230 &h1[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1231 Vmath::Fill(nquad0, 2.0 / (1.0 - z2[j]),
1232 &h2[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1233 Vmath::Fill(nquad0, (1.0 + z1[i]) / (1.0 - z2[j]),
1234 &h3[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1237 for (i = 0; i < nquad0; i++)
1239 Blas::Dscal(nquad1 * nquad2, 1 + z0[i], &h1[0] + i, nquad0);
1242 // Step 3. Construct combined metric terms for physical space to
1243 // collapsed coordinate system.
1244 // Order of construction optimised to minimise temporary storage
1245 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1248 Vmath::Vadd(nqtot, &df[1][0], 1, &df[2][0], 1, &wsp4[0], 1);
1249 Vmath::Vvtvvtp(nqtot, &df[0][0], 1, &h0[0], 1, &wsp4[0], 1, &h1[0], 1,
1252 Vmath::Vadd(nqtot, &df[4][0], 1, &df[5][0], 1, &wsp5[0], 1);
1253 Vmath::Vvtvvtp(nqtot, &df[3][0], 1, &h0[0], 1, &wsp5[0], 1, &h1[0], 1,
1256 Vmath::Vadd(nqtot, &df[7][0], 1, &df[8][0], 1, &wsp6[0], 1);
1257 Vmath::Vvtvvtp(nqtot, &df[6][0], 1, &h0[0], 1, &wsp6[0], 1, &h1[0], 1,
1261 Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0],
1263 Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1266 Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp4[0], 1, &df[5][0], 1, &wsp5[0],
1268 Vmath::Vvtvp(nqtot, &df[8][0], 1, &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1270 // overwrite h0, h1, h2
1271 // wsp7 (h2f1 + h3f2)
1272 Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &h2[0], 1, &df[2][0], 1, &h3[0], 1,
1274 // wsp8 (h2f4 + h3f5)
1275 Vmath::Vvtvvtp(nqtot, &df[4][0], 1, &h2[0], 1, &df[5][0], 1, &h3[0], 1,
1277 // wsp9 (h2f7 + h3f8)
1278 Vmath::Vvtvvtp(nqtot, &df[7][0], 1, &h2[0], 1, &df[8][0], 1, &h3[0], 1,
1282 Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp7[0], 1, &wsp5[0], 1, &wsp8[0],
1284 Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp9[0], 1, &g3[0], 1, &g3[0], 1);
1286 // overwrite wsp4, wsp5, wsp6
1288 Vmath::Vvtvvtp(nqtot, &wsp7[0], 1, &wsp7[0], 1, &wsp8[0], 1, &wsp8[0],
1290 Vmath::Vvtvp(nqtot, &wsp9[0], 1, &wsp9[0], 1, &g1[0], 1, &g1[0], 1);
1293 Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp7[0], 1, &df[5][0], 1, &wsp8[0],
1295 Vmath::Vvtvp(nqtot, &df[8][0], 1, &wsp9[0], 1, &g5[0], 1, &g5[0], 1);
1298 Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &df[2][0], 1, &df[5][0], 1,
1299 &df[5][0], 1, &g2[0], 1);
1300 Vmath::Vvtvp(nqtot, &df[8][0], 1, &df[8][0], 1, &g2[0], 1, &g2[0], 1);
1305 Vmath::Svtsvtp(nqtot, df[0][0], &h0[0], 1, df[1][0] + df[2][0], &h1[0],
1308 Vmath::Svtsvtp(nqtot, df[3][0], &h0[0], 1, df[4][0] + df[5][0], &h1[0],
1311 Vmath::Svtsvtp(nqtot, df[6][0], &h0[0], 1, df[7][0] + df[8][0], &h1[0],
1315 Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0],
1317 Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1320 Vmath::Svtsvtp(nqtot, df[2][0], &wsp4[0], 1, df[5][0], &wsp5[0], 1,
1322 Vmath::Svtvp(nqtot, df[8][0], &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1324 // overwrite h0, h1, h2
1325 // wsp7 (h2f1 + h3f2)
1326 Vmath::Svtsvtp(nqtot, df[1][0], &h2[0], 1, df[2][0], &h3[0], 1,
1328 // wsp8 (h2f4 + h3f5)
1329 Vmath::Svtsvtp(nqtot, df[4][0], &h2[0], 1, df[5][0], &h3[0], 1,
1331 // wsp9 (h2f7 + h3f8)
1332 Vmath::Svtsvtp(nqtot, df[7][0], &h2[0], 1, df[8][0], &h3[0], 1,
1336 Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp7[0], 1, &wsp5[0], 1, &wsp8[0],
1338 Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp9[0], 1, &g3[0], 1, &g3[0], 1);
1340 // overwrite wsp4, wsp5, wsp6
1342 Vmath::Vvtvvtp(nqtot, &wsp7[0], 1, &wsp7[0], 1, &wsp8[0], 1, &wsp8[0],
1344 Vmath::Vvtvp(nqtot, &wsp9[0], 1, &wsp9[0], 1, &g1[0], 1, &g1[0], 1);
1347 Vmath::Svtsvtp(nqtot, df[2][0], &wsp7[0], 1, df[5][0], &wsp8[0], 1,
1349 Vmath::Svtvp(nqtot, df[8][0], &wsp9[0], 1, &g5[0], 1, &g5[0], 1);
1353 df[2][0] * df[2][0] + df[5][0] * df[5][0] +
1354 df[8][0] * df[8][0],
1358 for (unsigned int i = 0; i < dim; ++i)
1360 for (unsigned int j = i; j < dim; ++j)
1362 MultiplyByQuadratureMetric(m_metrics[m[i][j]], m_metrics[m[i][j]]);
1372 void TetExp::v_NormalTraceDerivFactors(
1373 Array<OneD, Array<OneD, NekDouble>> &d0factors,
1374 Array<OneD, Array<OneD, NekDouble>> &d1factors,
1375 Array<OneD, Array<OneD, NekDouble>> &d2factors)
1377 int nquad0 = GetNumPoints(0);
1378 int nquad1 = GetNumPoints(1);
1379 int nquad2 = GetNumPoints(2);
1381 const Array<TwoD, const NekDouble> &df =
1382 m_metricinfo->GetDerivFactors(GetPointsKeys());
1384 if (d0factors.size() != 4)
1386 d0factors = Array<OneD, Array<OneD, NekDouble>>(4);
1387 d1factors = Array<OneD, Array<OneD, NekDouble>>(4);
1388 d2factors = Array<OneD, Array<OneD, NekDouble>>(4);
1391 if (d0factors[0].size() != nquad0 * nquad1)
1393 d0factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1394 d1factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1395 d2factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1398 if (d0factors[1].size() != nquad0 * nquad2)
1400 d0factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1401 d1factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1402 d2factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1405 if (d0factors[2].size() != nquad1 * nquad2)
1407 d0factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1408 d0factors[3] = Array<OneD, NekDouble>(nquad1 * nquad2);
1409 d1factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1410 d1factors[3] = Array<OneD, NekDouble>(nquad1 * nquad2);
1411 d2factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1412 d2factors[3] = Array<OneD, NekDouble>(nquad1 * nquad2);
1416 const Array<OneD, const Array<OneD, NekDouble>> &normal_0 =
1418 const Array<OneD, const Array<OneD, NekDouble>> &normal_1 =
1420 const Array<OneD, const Array<OneD, NekDouble>> &normal_2 =
1422 const Array<OneD, const Array<OneD, NekDouble>> &normal_3 =
1425 int ncoords = normal_0.size();
1427 // first gather together standard cartesian inner products
1428 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1431 for (int i = 0; i < nquad0 * nquad1; ++i)
1433 d0factors[0][i] = df[0][i] * normal_0[0][i];
1434 d1factors[0][i] = df[1][i] * normal_0[0][i];
1435 d2factors[0][i] = df[2][i] * normal_0[0][i];
1438 for (int n = 1; n < ncoords; ++n)
1440 for (int i = 0; i < nquad0 * nquad1; ++i)
1442 d0factors[0][i] += df[3 * n][i] * normal_0[n][i];
1443 d1factors[0][i] += df[3 * n + 1][i] * normal_0[n][i];
1444 d2factors[0][i] += df[3 * n + 2][i] * normal_0[n][i];
1449 for (int j = 0; j < nquad2; ++j)
1451 for (int i = 0; i < nquad0; ++i)
1453 d0factors[1][i] = df[0][j * nquad0 * nquad1 + i] *
1454 normal_1[0][j * nquad0 + i];
1455 d1factors[1][i] = df[1][j * nquad0 * nquad1 + i] *
1456 normal_1[0][j * nquad0 + i];
1457 d2factors[1][i] = df[2][j * nquad0 * nquad1 + i] *
1458 normal_1[0][j * nquad0 + i];
1462 for (int n = 1; n < ncoords; ++n)
1464 for (int j = 0; j < nquad2; ++j)
1466 for (int i = 0; i < nquad0; ++i)
1468 d0factors[1][i] = df[3 * n][j * nquad0 * nquad1 + i] *
1469 normal_1[0][j * nquad0 + i];
1470 d1factors[1][i] = df[3 * n + 1][j * nquad0 * nquad1 + i] *
1471 normal_1[0][j * nquad0 + i];
1472 d2factors[1][i] = df[3 * n + 2][j * nquad0 * nquad1 + i] *
1473 normal_1[0][j * nquad0 + i];
1479 for (int j = 0; j < nquad2; ++j)
1481 for (int i = 0; i < nquad1; ++i)
1483 d0factors[2][j * nquad1 + i] =
1484 df[0][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1485 normal_2[0][j * nquad1 + i];
1486 d0factors[3][j * nquad1 + i] =
1487 df[0][j * nquad0 * nquad1 + i * nquad0] *
1488 normal_3[0][j * nquad1 + i];
1489 d1factors[2][j * nquad1 + i] =
1490 df[1][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1491 normal_2[0][j * nquad1 + i];
1492 d1factors[3][j * nquad1 + i] =
1493 df[1][j * nquad0 * nquad1 + i * nquad0] *
1494 normal_3[0][j * nquad1 + i];
1495 d2factors[2][j * nquad1 + i] =
1496 df[2][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1497 normal_2[0][j * nquad1 + i];
1498 d2factors[3][j * nquad1 + i] =
1499 df[2][j * nquad0 * nquad1 + i * nquad0] *
1500 normal_3[0][j * nquad1 + i];
1504 for (int n = 1; n < ncoords; ++n)
1506 for (int j = 0; j < nquad2; ++j)
1508 for (int i = 0; i < nquad1; ++i)
1510 d0factors[2][j * nquad1 + i] +=
1511 df[3 * n][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1512 normal_2[n][j * nquad0 + i];
1513 d0factors[3][j * nquad0 + i] +=
1514 df[3 * n][i * nquad0 + j * nquad0 * nquad1] *
1515 normal_3[n][j * nquad0 + i];
1516 d1factors[2][j * nquad1 + i] +=
1518 [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1519 normal_2[n][j * nquad0 + i];
1520 d1factors[3][j * nquad0 + i] +=
1521 df[3 * n + 1][i * nquad0 + j * nquad0 * nquad1] *
1522 normal_3[n][j * nquad0 + i];
1523 d2factors[2][j * nquad1 + i] +=
1525 [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1526 normal_2[n][j * nquad0 + i];
1527 d2factors[3][j * nquad0 + i] +=
1528 df[3 * n + 2][i * nquad0 + j * nquad0 * nquad1] *
1529 normal_3[n][j * nquad0 + i];
1537 for (int i = 0; i < nquad0 * nquad1; ++i)
1539 d0factors[0][i] = df[0][0] * normal_0[0][i];
1540 d1factors[0][i] = df[1][0] * normal_0[0][i];
1541 d2factors[0][i] = df[2][0] * normal_0[0][i];
1544 for (int n = 1; n < ncoords; ++n)
1546 for (int i = 0; i < nquad0 * nquad1; ++i)
1548 d0factors[0][i] += df[3 * n][0] * normal_0[n][i];
1549 d1factors[0][i] += df[3 * n + 1][0] * normal_0[n][i];
1550 d2factors[0][i] += df[3 * n + 2][0] * normal_0[n][i];
1555 for (int i = 0; i < nquad0 * nquad2; ++i)
1557 d0factors[1][i] = df[0][0] * normal_1[0][i];
1558 d1factors[1][i] = df[1][0] * normal_1[0][i];
1559 d2factors[1][i] = df[2][0] * normal_1[0][i];
1562 for (int n = 1; n < ncoords; ++n)
1564 for (int i = 0; i < nquad0 * nquad2; ++i)
1566 d0factors[1][i] += df[3 * n][0] * normal_1[n][i];
1567 d1factors[1][i] += df[3 * n + 1][0] * normal_1[n][i];
1568 d2factors[1][i] += df[3 * n + 2][0] * normal_1[n][i];
1573 for (int i = 0; i < nquad1 * nquad2; ++i)
1575 d0factors[2][i] = df[0][0] * normal_2[0][i];
1576 d0factors[3][i] = df[0][0] * normal_3[0][i];
1578 d1factors[2][i] = df[1][0] * normal_2[0][i];
1579 d1factors[3][i] = df[1][0] * normal_3[0][i];
1581 d2factors[2][i] = df[2][0] * normal_2[0][i];
1582 d2factors[3][i] = df[2][0] * normal_3[0][i];
1585 for (int n = 1; n < ncoords; ++n)
1587 for (int i = 0; i < nquad1 * nquad2; ++i)
1589 d0factors[2][i] += df[3 * n][0] * normal_2[n][i];
1590 d0factors[3][i] += df[3 * n][0] * normal_3[n][i];
1592 d1factors[2][i] += df[3 * n + 1][0] * normal_2[n][i];
1593 d1factors[3][i] += df[3 * n + 1][0] * normal_3[n][i];
1595 d2factors[2][i] += df[3 * n + 2][0] * normal_2[n][i];
1596 d2factors[3][i] += df[3 * n + 2][0] * normal_3[n][i];
1601 } // namespace LocalRegions
1602 } // namespace Nektar
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Describes the specification for a Basis.
int GetNumPoints() const
Return points order at which basis is defined.
PointsKey GetPointsKey() const
Return distribution of points.
Defines a specification for a set of points.
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
std::map< int, NormalVector > m_traceNormals
std::map< int, Array< OneD, NekDouble > > m_elmtBndNormDirElmtLen
the element length in each element boundary(Vertex, edge or face) normal direction calculated based o...
SpatialDomains::GeometrySharedPtr GetGeom() const
SpatialDomains::GeometrySharedPtr m_geom
void ComputeLaplacianMetric()
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)
virtual void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords)
Get the coordinates "coords" at the local coordinates "Lcoords".
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_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
Differentiate inarray in the three coordinate directions.
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
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 int v_GetCoordim()
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Integrate the physical point list inarray over region.
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Calculate the inner product of inarray with respect to the basis B=m_base0*m_base1*m_base2 and put in...
virtual NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
void v_ComputeTraceNormal(const int face)
Compute the normal of a triangular face.
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey)
virtual LibUtilities::ShapeType v_DetShapeType() const
Return Shape of region, using ShapeType enum list.
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
void GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Calculates the inner product .
virtual void v_GetTracePhysMap(const int face, Array< OneD, int > &outarray)
Returns the physical values at the quadrature points of a face.
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
virtual NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &physvals)
virtual void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp)
virtual void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
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::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Array< OneD, LibUtilities::BasisSharedPtr > m_base
MatrixType GetMatrixType() const
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 = 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.
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
The above copyright notice and this permission notice shall be included.
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
void 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 Vcopy(int n, const T *x, const int incx, T *y, const int incy)
scalarT< T > sqrt(scalarT< T > in)