42 namespace LocalRegions
50 StdNodalTriExp(Ba,Bb,Ntype),
55 std::string(
"NodalTriExpMatrix")),
56 m_staticCondMatrixManager(
57 boost::bind(&
NodalTriExp::CreateStaticCondMatrix, this, _1),
58 std::string(
"NodalTriExpStaticCondMatrix"))
65 StdRegions::StdNodalTriExp(T),
68 m_matrixManager(T.m_matrixManager),
69 m_staticCondMatrixManager(T.m_staticCondMatrixManager)
100 int nquad0 =
m_base[0]->GetNumPoints();
101 int nquad1 =
m_base[1]->GetNumPoints();
109 Vmath::Vmul(nquad0*nquad1, jac, 1, inarray, 1,tmp, 1);
113 Vmath::Smul(nquad0*nquad1, jac[0], inarray, 1, tmp, 1);
117 ival = StdNodalTriExp::v_Integral(tmp);
124 bool multiplybyweights)
126 int nquad0 =
m_base[0]->GetNumPoints();
127 int nquad1 =
m_base[1]->GetNumPoints();
128 int order1 =
m_base[1]->GetNumModes();
130 if(multiplybyweights)
136 StdTriExp::IProductWRTBase_SumFacKernel(
m_base[0]->GetBdata(),
m_base[1]->GetBdata(),tmp,outarray,wsp);
143 StdTriExp::IProductWRTBase_SumFacKernel(
m_base[0]->GetBdata(),
m_base[1]->GetBdata(),inarray,outarray,wsp);
155 Blas::Dgemv(
'N',
m_ncoeffs,nq,iprodmat->Scale(),(iprodmat->GetOwnedMatrix())->GetPtr().get(),
156 m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
163 ASSERTL1((dir==0)||(dir==1)||(dir==2),
"Invalid direction.");
164 ASSERTL1((dir==2)?(
m_geom->GetCoordim()==3):
true,
"Invalid direction.");
167 int nquad0 =
m_base[0]->GetNumPoints();
168 int nquad1 =
m_base[1]->GetNumPoints();
169 int nqtot = nquad0*nquad1;
186 for(i = 0; i < nquad1; ++i)
188 gfac0[i] = 2.0/(1-z1[i]);
190 for(i = 0; i < nquad0; ++i)
192 gfac1[i] = 0.5*(1+z0[i]);
195 for(i = 0; i < nquad1; ++i)
197 Vmath::Smul(nquad0,gfac0[i],&inarray[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
200 for(i = 0; i < nquad1; ++i)
202 Vmath::Vmul(nquad0,&gfac1[0],1,&tmp0[0]+i*nquad0,1,&tmp1[0]+i*nquad0,1);
207 Vmath::Vmul(nqtot,&df[2*dir][0], 1,&tmp0[0], 1,&tmp0[0],1);
208 Vmath::Vmul(nqtot,&df[2*dir+1][0],1,&tmp1[0], 1,&tmp1[0],1);
209 Vmath::Vmul(nqtot,&df[2*dir+1][0],1,&inarray[0],1,&tmp2[0],1);
213 Vmath::Smul(nqtot, df[2*dir][0], tmp0, 1, tmp0, 1);
214 Vmath::Smul(nqtot, df[2*dir+1][0], tmp1, 1, tmp1, 1);
215 Vmath::Smul(nqtot, df[2*dir+1][0], inarray, 1, tmp2, 1);
255 ASSERTL1(
false,
"input dir is out of range");
263 Blas::Dgemv(
'N',
m_ncoeffs,nq,iprodmat->Scale(),(iprodmat->GetOwnedMatrix())->GetPtr().get(),
264 m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
280 int nquad0 =
m_base[0]->GetNumPoints();
281 int nquad1 =
m_base[1]->GetNumPoints();
282 int nqtot = nquad0*nquad1;
289 StdNodalTriExp::v_PhysDeriv(inarray, diff0, diff1);
293 if(out_d0.num_elements())
296 Vmath::Vvtvp (nqtot,df[1],1,diff1,1, out_d0, 1, out_d0,1);
299 if(out_d1.num_elements())
302 Vmath::Vvtvp (nqtot,df[3],1,diff1,1, out_d1, 1, out_d1,1);
305 if(out_d2.num_elements())
308 Vmath::Vvtvp (nqtot,df[5],1,diff1,1, out_d2, 1, out_d2,1);
313 if(out_d0.num_elements())
315 Vmath::Smul (nqtot, df[0][0], diff0, 1, out_d0, 1);
316 Blas::Daxpy (nqtot, df[1][0], diff1, 1, out_d0, 1);
319 if(out_d1.num_elements())
321 Vmath::Smul (nqtot, df[2][0], diff0, 1, out_d1, 1);
322 Blas::Daxpy (nqtot, df[3][0], diff1, 1, out_d1, 1);
325 if(out_d2.num_elements())
327 Vmath::Smul (nqtot, df[4][0], diff0, 1, out_d2, 1);
328 Blas::Daxpy (nqtot, df[5][0], diff1, 1, out_d2, 1);
368 if(inarray.get() == outarray.get())
373 Blas::Dgemv(
'N',
m_ncoeffs,
m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
374 m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
378 Blas::Dgemv(
'N',
m_ncoeffs,
m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
379 m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
396 ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[1] <= 1.0 &&
397 Lcoords[1] >= -1.0 && Lcoords[1] <=1.0,
398 "Local coordinates are not in region [-1,1]");
402 for(i = 0; i <
m_geom->GetCoordim(); ++i)
404 coords[i] =
m_geom->GetCoord(i,Lcoords);
416 return tmp->GetStdMatrix(mkey);
427 m_geom->GetLocCoords(coord,Lcoord);
429 return StdNodalTriExp::v_PhysEvaluate(Lcoord, physvals);
490 ASSERTL1(
m_geom->GetCoordim() == 2,
"Standard Region Laplacian is only set up for Quads in two-dimensional");
506 int rows = lap00.GetRows();
507 int cols = lap00.GetColumns();
511 (*lap) = gmat[0][0] * lap00 +
512 gmat[1][0] * (lap01 +
Transpose(lap01)) +
529 int rows = LapMat.GetRows();
530 int cols = LapMat.GetColumns();
535 (*helm) = LapMat + factor*MassMat;
556 unsigned int nint = (
unsigned int)(
m_ncoeffs - nbdry);
557 unsigned int exp_size[] = {nbdry,nint};
558 unsigned int nblks = 2;
569 goto UseLocRegionsMatrix;
575 goto UseLocRegionsMatrix;
580 factor = mat->Scale();
581 goto UseStdRegionsMatrix;
614 for(i = 0; i < nbdry; ++i)
616 for(j = 0; j < nbdry; ++j)
618 (*A)(i,j) = mat(bmap[i],bmap[j]);
621 for(j = 0; j < nint; ++j)
623 (*B)(i,j) = mat(bmap[i],imap[j]);
627 for(i = 0; i < nint; ++i)
629 for(j = 0; j < nbdry; ++j)
631 (*C)(i,j) = mat(imap[i],bmap[j]);
634 for(j = 0; j < nint; ++j)
636 (*D)(i,j) = mat(imap[i],imap[j]);
645 (*A) = (*A) - (*B)*(*C);
686 returnval = StdNodalTriExp::v_GenMatrix(mkey);
701 int nqe =
m_base[0]->GetNumPoints();
706 for (i = 0; i < dim; ++i)
727 Vmath::Fill(nqe,df[2*i+1][0] + df[2*i][0],normal[i],1);
737 ASSERTL0(
false,
"Edge is out of range (edge < 3)");
744 fac += normal[i][0]*normal[i][0];
756 int nquad0 = ptsKeys[0].GetNumPoints();
757 int nquad1 = ptsKeys[1].GetNumPoints();
770 for(j = 0; j < nquad0; ++j)
775 normals[i*nquad0+j] = -df[2*i+1][j]*edgejac[j];
778 from_key = ptsKeys[0];
781 for(j = 0; j < nquad1; ++j)
783 edgejac[j] = jac[nquad0*j+nquad0-1];
786 normals[i*nquad1+j] = (df[2*i][nquad0*j + nquad0-1] + df[2*i+1][nquad0*j + nquad0-1])*edgejac[j];
789 from_key = ptsKeys[1];
792 for(j = 0; j < nquad1; ++j)
794 edgejac[j] = jac[nquad0*j];
797 normals[i*nquad1+j] = -df[2*i][nquad0*j]*edgejac[j];
800 from_key = ptsKeys[1];
803 ASSERTL0(
false,
"edge is out of range (edge < 3)");
825 Vmath::Vvtvp(nqe,normal[i],1, normal[i],1,work,1,work,1);
const LibUtilities::PointsKeyVector GetPointsKeys() const
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
NekDouble GetConstFactor(const ConstFactorType &factor) const
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
std::vector< PointsKey > PointsKeyVector
MatrixType GetMatrixType() const
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_ComputeEdgeNormal(const int edge)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords)
void IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
void IProductWRTBase_MatOp(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, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Differentiation Methods.
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
std::map< int, StdRegions::NormalVector > m_edgeNormals
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
LibUtilities::ShapeType GetShapeType() const
SpatialDomains::GeometrySharedPtr m_geom
boost::shared_ptr< DNekMat > DNekMatSharedPtr
void GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
boost::shared_ptr< StdNodalTriExp > StdNodalTriExpSharedPtr
NekDouble PhysEvaluate(const Array< OneD, const NekDouble > &coord, const Array< OneD, const NekDouble > &physvals)
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
void NodalToModalTranspose(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
int NumBndryCoeffs(void) const
DNekBlkMatSharedPtr GetStdStaticCondMatrix(const StdMatrixKey &mkey)
LibUtilities::PointsKey m_nodalPointsKey
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
int getNumberOfCoefficients(int Na)
void IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
void GetInteriorMap(Array< OneD, unsigned int > &outarray)
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Defines a specification for a set of points.
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...
~NodalTriExp()
Destructor.
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
NekDouble Integral(const Array< OneD, const NekDouble > &inarray)
Integrate the physical point list inarray over region.
boost::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
SpatialDomains::GeometrySharedPtr GetGeom() const
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Geometry is straight-sided with constant geometric factors.
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 ...
DNekScalBlkMatSharedPtr CreateStaticCondMatrix(const MatrixKey &mkey)
unsigned int GetNumPoints() const
boost::shared_ptr< TriGeom > TriGeomSharedPtr
void GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3=NullNekDouble1DArray)
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)
GeomType
Indicates the type of element geometry.
void Zero(int n, T *x, const int incx)
Zero vector.
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Array< OneD, LibUtilities::BasisSharedPtr > m_base
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Geometry is curved or has non-constant factors.
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 GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Describes the specification for a Basis.
PointsType GetPointsType() const
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 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.
DNekMatSharedPtr CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
static VarCoeffMap NullVarCoeffMap
static ConstFactorMap NullConstFactorMap