47 Ba.GetNumModes(), (Bb.GetNumModes())),
50 Ba.GetNumModes(), (Bb.GetNumModes())),
55 std::string(
"NodalTriExpMatrix")),
56 m_staticCondMatrixManager(
std::bind(&
Expansion::CreateStaticCondMatrix,
57 this,
std::placeholders::_1),
58 std::string(
"NodalTriExpStaticCondMatrix"))
63 : StdExpansion(T), StdExpansion2D(T), StdRegions::StdTriExp(T),
65 m_matrixManager(T.m_matrixManager),
66 m_staticCondMatrixManager(T.m_staticCondMatrixManager)
72 int nquad0 =
m_base[0]->GetNumPoints();
73 int nquad1 =
m_base[1]->GetNumPoints();
81 Vmath::Vmul(nquad0 * nquad1, jac, 1, inarray, 1, tmp, 1);
85 Vmath::Smul(nquad0 * nquad1, jac[0], inarray, 1, tmp, 1);
89 ival = StdNodalTriExp::v_Integral(tmp);
98 int nquad0 =
m_base[0]->GetNumPoints();
99 int nquad1 =
m_base[1]->GetNumPoints();
100 int nqtot = nquad0 * nquad1;
107 StdNodalTriExp::v_PhysDeriv(inarray, diff0, diff1);
114 Vmath::Vvtvp(nqtot, df[1], 1, diff1, 1, out_d0, 1, out_d0, 1);
120 Vmath::Vvtvp(nqtot, df[3], 1, diff1, 1, out_d1, 1, out_d1, 1);
126 Vmath::Vvtvp(nqtot, df[5], 1, diff1, 1, out_d2, 1, out_d2, 1);
170 ASSERTL1(dir >= 0 && dir < 2,
"input dir is out of range");
192 out = (*matsys) * in;
212 int nquad0 =
m_base[0]->GetNumPoints();
213 int nquad1 =
m_base[1]->GetNumPoints();
214 int order1 =
m_base[1]->GetNumModes();
216 if (multiplybyweights)
222 StdTriExp::IProductWRTBase_SumFacKernel(
223 m_base[0]->GetBdata(),
m_base[1]->GetBdata(), tmp, outarray, wsp);
230 StdTriExp::IProductWRTBase_SumFacKernel(
m_base[0]->GetBdata(),
231 m_base[1]->GetBdata(), inarray,
241 int nquad0 =
m_base[0]->GetNumPoints();
242 int nquad1 =
m_base[1]->GetNumPoints();
243 int nqtot = nquad0 * nquad1;
263 tmp2, outarray, tmp0);
273 ASSERTL1((dir == 0) || (dir == 1) || (dir == 2),
"Invalid direction.");
275 "Invalid direction.");
277 int nquad0 =
m_base[0]->GetNumPoints();
278 int nquad1 =
m_base[1]->GetNumPoints();
279 int nqtot = nquad0 * nquad1;
297 for (
int i = 0; i < nquad1; ++i)
299 gfac0[i] = 2.0 / (1 - z1[i]);
301 for (
int i = 0; i < nquad0; ++i)
303 gfac1[i] = 0.5 * (1 + z0[i]);
306 for (
int i = 0; i < nquad1; ++i)
308 Vmath::Smul(nquad0, gfac0[i], &inarray[0] + i * nquad0, 1,
309 &tmp0[0] + i * nquad0, 1);
312 for (
int i = 0; i < nquad1; ++i)
314 Vmath::Vmul(nquad0, &gfac1[0], 1, &tmp0[0] + i * nquad0, 1,
315 &tmp1[0] + i * nquad0, 1);
320 Vmath::Vmul(nqtot, &df[2 * dir][0], 1, &tmp0[0], 1, &tmp0[0], 1);
321 Vmath::Vmul(nqtot, &df[2 * dir + 1][0], 1, &tmp1[0], 1, &tmp1[0], 1);
322 Vmath::Vmul(nqtot, &df[2 * dir + 1][0], 1, &inarray[0], 1, &tmp2[0], 1);
326 Vmath::Smul(nqtot, df[2 * dir][0], tmp0, 1, tmp0, 1);
327 Vmath::Smul(nqtot, df[2 * dir + 1][0], tmp1, 1, tmp1, 1);
328 Vmath::Smul(nqtot, df[2 * dir + 1][0], inarray, 1, tmp2, 1);
344 m_base[0]->GetPointsKey());
346 m_base[1]->GetPointsKey());
364 ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[1] <= 1.0 && Lcoords[1] >= -1.0 &&
366 "Local coordinates are not in region [-1,1]");
370 for (i = 0; i <
m_geom->GetCoordim(); ++i)
372 coords[i] =
m_geom->GetCoord(i, Lcoords);
384 m_geom->GetLocCoords(coord, Lcoord);
386 return StdExpansion2D::v_PhysEvaluate(Lcoord, physvals);
398 geomFactors->GetDerivFactors(ptsKeys);
400 int nqe =
m_base[0]->GetNumPoints();
405 for (i = 0; i < dim; ++i)
432 Vmath::Fill(nqe, df[2 * i + 1][0] + df[2 * i][0], normal[i],
443 ASSERTL0(
false,
"Edge is out of range (edge < 3)");
450 fac += normal[i][0] * normal[i][0];
452 fac = 1.0 /
sqrt(fac);
465 int nquad0 = ptsKeys[0].GetNumPoints();
466 int nquad1 = ptsKeys[1].GetNumPoints();
479 for (j = 0; j < nquad0; ++j)
484 normals[i * nquad0 + j] =
485 -df[2 * i + 1][j] * edgejac[j];
488 from_key = ptsKeys[0];
491 for (j = 0; j < nquad1; ++j)
493 edgejac[j] = jac[nquad0 * j + nquad0 - 1];
496 normals[i * nquad1 + j] =
497 (df[2 * i][nquad0 * j + nquad0 - 1] +
498 df[2 * i + 1][nquad0 * j + nquad0 - 1]) *
502 from_key = ptsKeys[1];
505 for (j = 0; j < nquad1; ++j)
507 edgejac[j] = jac[nquad0 * j];
510 normals[i * nquad1 + j] =
511 -df[2 * i][nquad0 * j] * edgejac[j];
514 from_key = ptsKeys[1];
517 ASSERTL0(
false,
"edge is out of range (edge < 3)");
531 m_base[0]->GetPointsKey(), &normal[i][0]);
532 Vmath::Vmul(nqe, work, 1, normal[i], 1, normal[i], 1);
539 Vmath::Vvtvp(nqe, normal[i], 1, normal[i], 1, work, 1, work, 1);
549 Vmath::Vmul(nqe, normal[i], 1, work, 1, normal[i], 1);
565 const NekDouble *data,
const std::vector<unsigned int> &nummodes,
566 const int mode_offset,
NekDouble *coeffs,
567 [[maybe_unused]] std::vector<LibUtilities::BasisType> &fromType)
582 int nquad0 =
m_base[0]->GetNumPoints();
583 int nquad1 =
m_base[1]->GetNumPoints();
590 Vmath::Vcopy(nquad0, &(inarray[0]), 1, &(outarray[0]), 1);
594 Vmath::Vcopy(nquad1, &(inarray[0]) + (nquad0 - 1), nquad0,
599 Vmath::Vcopy(nquad1, &(inarray[0]), nquad0, &(outarray[0]), 1);
603 ASSERTL0(
false,
"edge value (< 3) is out of range");
607 ASSERTL1(EdgeExp->GetBasis(0)->GetPointsType() ==
609 "Edge expansion should be GLL");
612 if (
m_base[edge ? 1 : 0]->GetPointsKey() !=
613 EdgeExp->GetBasis(0)->GetPointsKey())
620 EdgeExp->GetBasis(0)->GetPointsKey(), outarray);
631 Vmath::Reverse(EdgeExp->GetNumPoints(0), &outarray[0], 1, &outarray[0],
651 returnval = StdNodalTriExp::v_GenMatrix(mkey);
666 return tmp->GetStdMatrix(mkey);
687 StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
694 StdExpansion::LaplacianMatrixOp_MatFree_GenericImpl(inarray, outarray,
702 StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
709 StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
716 StdExpansion::HelmholtzMatrixOp_MatFree_GenericImpl(inarray, outarray,
#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.
Defines a specification for a set of points.
PointsType GetPointsType() const
size_t GetNumPoints() const
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 AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
void ExtractDataToCoeffs(const NekDouble *data, const std::vector< unsigned int > &nummodes, const int nmodes_offset, NekDouble *coeffs, std::vector< LibUtilities::BasisType > &fromType)
StdRegions::Orientation GetTraceOrient(int trace)
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) 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_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_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Calculates the inner product of a given function f with the different modes of the expansion.
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray) override
NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray) override
Integrates the specified function over the domain.
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
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_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) 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
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.
StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const override
void v_GetCoords(Array< OneD, NekDouble > &coords_0, Array< OneD, NekDouble > &coords_1=NullNekDouble1DArray, Array< OneD, NekDouble > &coords_2=NullNekDouble1DArray) override
void v_DropLocMatrix(const MatrixKey &mkey) override
void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey) override
StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const override
NodalTriExp(const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const LibUtilities::PointsType Ntype, const SpatialDomains::TriGeomSharedPtr &geom)
Constructor using BasisKey class for quadrature points and order definition.
DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey) override
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true) override
DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey) override
void v_ComputeTraceNormal(const int edge) override
void v_GetCoord(const Array< OneD, const NekDouble > &lcoord, Array< OneD, NekDouble > &coord) 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)
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::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Array< OneD, LibUtilities::BasisSharedPtr > m_base
MatrixType GetMatrixType() const
void NodalToModalTranspose(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void ModalToNodal(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
LibUtilities::PointsKey m_nodalPointsKey
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 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 ...
std::vector< PointsKey > PointsKeyVector
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
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< TriGeom > TriGeomSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
std::shared_ptr< StdNodalTriExp > StdNodalTriExpSharedPtr
static ConstFactorMap NullConstFactorMap
static VarCoeffMap NullVarCoeffMap
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 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 Vcopy(int n, const T *x, const int incx, T *y, const int incy)
scalarT< T > sqrt(scalarT< T > in)