47 Ba.GetNumModes(), (Bb.GetNumModes())),
50 Ba.GetNumModes(), (Bb.GetNumModes())),
54 std::bind(&
Expansion2D::CreateMatrix, this, std::placeholders::_1),
55 std::string(
"NodalTriExpMatrix")),
56 m_staticCondMatrixManager(std::bind(&
Expansion::CreateStaticCondMatrix,
57 this, std::placeholders::_1),
58 std::string(
"NodalTriExpStaticCondMatrix"))
64 StdExpansion2D(T), StdRegions::StdTriExp(T), StdRegions::StdNodalTriExp(
67 m_staticCondMatrixManager(T.m_staticCondMatrixManager)
93 int nquad0 =
m_base[0]->GetNumPoints();
94 int nquad1 =
m_base[1]->GetNumPoints();
102 Vmath::Vmul(nquad0 * nquad1, jac, 1, inarray, 1, tmp, 1);
106 Vmath::Smul(nquad0 * nquad1, jac[0], inarray, 1, tmp, 1);
110 ival = StdNodalTriExp::v_Integral(tmp);
118 int nquad0 =
m_base[0]->GetNumPoints();
119 int nquad1 =
m_base[1]->GetNumPoints();
120 int order1 =
m_base[1]->GetNumModes();
122 if (multiplybyweights)
128 StdTriExp::IProductWRTBase_SumFacKernel(
129 m_base[0]->GetBdata(),
m_base[1]->GetBdata(), tmp, outarray, wsp);
136 StdTriExp::IProductWRTBase_SumFacKernel(
m_base[0]->GetBdata(),
137 m_base[1]->GetBdata(), inarray,
152 (iprodmat->GetOwnedMatrix())->GetPtr().get(),
m_ncoeffs,
153 inarray.get(), 1, 0.0, outarray.get(), 1);
160 int nquad0 =
m_base[0]->GetNumPoints();
161 int nquad1 =
m_base[1]->GetNumPoints();
162 int nqtot = nquad0 * nquad1;
182 tmp2, outarray, tmp0);
192 ASSERTL1((dir == 0) || (dir == 1) || (dir == 2),
"Invalid direction.");
194 "Invalid direction.");
196 int nquad0 =
m_base[0]->GetNumPoints();
197 int nquad1 =
m_base[1]->GetNumPoints();
198 int nqtot = nquad0 * nquad1;
216 for (
int i = 0; i < nquad1; ++i)
218 gfac0[i] = 2.0 / (1 - z1[i]);
220 for (
int i = 0; i < nquad0; ++i)
222 gfac1[i] = 0.5 * (1 + z0[i]);
225 for (
int i = 0; i < nquad1; ++i)
227 Vmath::Smul(nquad0, gfac0[i], &inarray[0] + i * nquad0, 1,
228 &tmp0[0] + i * nquad0, 1);
231 for (
int i = 0; i < nquad1; ++i)
233 Vmath::Vmul(nquad0, &gfac1[0], 1, &tmp0[0] + i * nquad0, 1,
234 &tmp1[0] + i * nquad0, 1);
239 Vmath::Vmul(nqtot, &df[2 * dir][0], 1, &tmp0[0], 1, &tmp0[0], 1);
240 Vmath::Vmul(nqtot, &df[2 * dir + 1][0], 1, &tmp1[0], 1, &tmp1[0], 1);
241 Vmath::Vmul(nqtot, &df[2 * dir + 1][0], 1, &inarray[0], 1, &tmp2[0], 1);
245 Vmath::Smul(nqtot, df[2 * dir][0], tmp0, 1, tmp0, 1);
246 Vmath::Smul(nqtot, df[2 * dir + 1][0], tmp1, 1, tmp1, 1);
247 Vmath::Smul(nqtot, df[2 * dir + 1][0], inarray, 1, tmp2, 1);
278 ASSERTL1(
false,
"input dir is out of range");
287 (iprodmat->GetOwnedMatrix())->GetPtr().get(),
m_ncoeffs,
288 inarray.get(), 1, 0.0, outarray.get(), 1);
303 int nquad0 =
m_base[0]->GetNumPoints();
304 int nquad1 =
m_base[1]->GetNumPoints();
305 int nqtot = nquad0 * nquad1;
312 StdNodalTriExp::v_PhysDeriv(inarray, diff0, diff1);
319 Vmath::Vvtvp(nqtot, df[1], 1, diff1, 1, out_d0, 1, out_d0, 1);
325 Vmath::Vvtvp(nqtot, df[3], 1, diff1, 1, out_d1, 1, out_d1, 1);
331 Vmath::Vvtvp(nqtot, df[5], 1, diff1, 1, out_d2, 1, out_d2, 1);
385 out = (*matsys) * in;
394 if (inarray.get() == outarray.get())
400 (mat->GetOwnedMatrix())->GetPtr().get(),
m_ncoeffs,
401 tmp.get(), 1, 0.0, outarray.get(), 1);
406 (mat->GetOwnedMatrix())->GetPtr().get(),
m_ncoeffs,
407 inarray.get(), 1, 0.0, outarray.get(), 1);
424 ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[1] <= 1.0 && Lcoords[1] >= -1.0 &&
426 "Local coordinates are not in region [-1,1]");
430 for (i = 0; i <
m_geom->GetCoordim(); ++i)
432 coords[i] =
m_geom->GetCoord(i, Lcoords);
445 return tmp->GetStdMatrix(mkey);
456 m_geom->GetLocCoords(coord, Lcoord);
458 return StdExpansion2D::v_PhysEvaluate(Lcoord, physvals);
464 std::array<NekDouble, 3> &firstOrderDerivs)
468 m_geom->GetLocCoords(coord, Lcoord);
469 return StdExpansion2D::v_PhysEvaluate(Lcoord, inarray, firstOrderDerivs);
483 m_base[0]->GetPointsKey());
485 m_base[1]->GetPointsKey());
506 returnval = StdNodalTriExp::v_GenMatrix(mkey);
521 geomFactors->GetDerivFactors(ptsKeys);
523 int nqe =
m_base[0]->GetNumPoints();
528 for (i = 0; i < dim; ++i)
555 Vmath::Fill(nqe, df[2 * i + 1][0] + df[2 * i][0], normal[i],
566 ASSERTL0(
false,
"Edge is out of range (edge < 3)");
573 fac += normal[i][0] * normal[i][0];
575 fac = 1.0 /
sqrt(fac);
588 int nquad0 = ptsKeys[0].GetNumPoints();
589 int nquad1 = ptsKeys[1].GetNumPoints();
602 for (j = 0; j < nquad0; ++j)
607 normals[i * nquad0 + j] =
608 -df[2 * i + 1][j] * edgejac[j];
611 from_key = ptsKeys[0];
614 for (j = 0; j < nquad1; ++j)
616 edgejac[j] = jac[nquad0 * j + nquad0 - 1];
619 normals[i * nquad1 + j] =
620 (df[2 * i][nquad0 * j + nquad0 - 1] +
621 df[2 * i + 1][nquad0 * j + nquad0 - 1]) *
625 from_key = ptsKeys[1];
628 for (j = 0; j < nquad1; ++j)
630 edgejac[j] = jac[nquad0 * j];
633 normals[i * nquad1 + j] =
634 -df[2 * i][nquad0 * j] * edgejac[j];
637 from_key = ptsKeys[1];
640 ASSERTL0(
false,
"edge is out of range (edge < 3)");
654 m_base[0]->GetPointsKey(), &normal[i][0]);
655 Vmath::Vmul(nqe, work, 1, normal[i], 1, normal[i], 1);
662 Vmath::Vvtvp(nqe, normal[i], 1, normal[i], 1, work, 1, work, 1);
672 Vmath::Vmul(nqe, normal[i], 1, work, 1, normal[i], 1);
#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
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
void GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords)
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Differentiation Methods.
DNekMatSharedPtr CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
NekDouble v_PhysEvaluate(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) final
void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray) override
void GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3=NullNekDouble1DArray)
void GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Inner product of inarray over region with respect to the expansion basis (this)->_Base[0] and return ...
void IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void 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...
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
NekDouble PhysEvaluate(const Array< OneD, const NekDouble > &coord, const Array< OneD, const NekDouble > &physvals)
StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const override
NekDouble Integral(const Array< OneD, const NekDouble > &inarray)
Integrate the physical point list inarray over region.
void IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
void IProductWRTBase_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
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.
void v_ComputeTraceNormal(const int edge) 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 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
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)
Array< OneD, LibUtilities::BasisSharedPtr > m_base
MatrixType GetMatrixType() const
void NodalToModalTranspose(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
LibUtilities::PointsKey m_nodalPointsKey
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 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
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< 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)