47 namespace LocalRegions
50 StdExpansion(),
Expansion(pGeom), StdExpansion2D()
64 "Routine only set up for two-dimensions");
85 if (edgeExp->GetRightAdjacentElementExp())
87 if (edgeExp->GetRightAdjacentElementExp()->GetGeom2D()
88 ->GetGlobalID() ==
GetGeom2D()->GetGlobalID())
100 int nquad_e = min(EdgeExp->GetNumPoints(0),
101 int(normals[0].num_elements()));
103 int nEdgePts = EdgeExp->GetTotPoints();
105 Vmath::Vmul (nquad_e, normals[0], 1, Fx, 1, edgePhys, 1);
106 Vmath::Vvtvp(nquad_e, normals[1], 1, Fy, 1, edgePhys, 1,
115 else if (locExp->GetRightAdjacentElementEdge() != -1)
117 if (locExp->GetRightAdjacentElementExp()->GetGeom2D()->GetGlobalID()
151 if (edgeExp->GetRightAdjacentElementExp())
153 if (edgeExp->GetRightAdjacentElementExp()->GetGeom2D()
154 ->GetGlobalID() ==
GetGeom2D()->GetGlobalID())
167 StdExpansion::GetIndexMap(ikey);
170 int order_e = map->num_elements();
172 int n_coeffs = EdgeExp->GetNcoeffs();
175 if(n_coeffs!=order_e)
177 EdgeExp->FwdTrans(Fn, edgeCoeffs);
187 LibUtilities::BasisKey bkey_ortho(btype,EdgeExp->GetBasis(0)->GetNumModes(),EdgeExp->GetBasis(0)->GetPointsKey());
188 LibUtilities::BasisKey bkey(EdgeExp->GetBasis(0)->GetBasisType(),EdgeExp->GetBasis(0)->GetNumModes(),EdgeExp->GetBasis(0)->GetPointsKey());
192 for(i = order_e; i < n_coeffs; i++)
200 EdgeExp->MassMatrixOp(edgeCoeffs, edgeCoeffs, masskey);
204 EdgeExp->IProductWRTBase(Fn, edgeCoeffs);
215 if(EdgeExp->GetBasis(0)->GetBasisType() !=
219 for(i = 0; i < order_e; ++i)
221 outarray[(*map)[i].index] +=
222 (*map)[i].sign * edgeCoeffs[i];
227 int nCoeffs0, nCoeffs1;
241 nCoeffs1 =
m_base[1]->GetNumModes();
243 for(i = 0; i < order_e; ++i)
245 for(j = 0; j < nCoeffs1; j++)
247 outarray[(*map)[i].index + j*order_e] +=
248 mat_gauss->GetPtr()[j]*
249 (*map)[i].sign*edgeCoeffs[i];
256 nCoeffs0 =
m_base[0]->GetNumModes();
258 for(i = 0; i < order_e; ++i)
260 for(j = 0; j < nCoeffs0; j++)
262 outarray[(*map)[i].index - j] +=
263 mat_gauss->GetPtr()[order_e - 1 -j]*
264 (*map)[i].sign*edgeCoeffs[i];
271 nCoeffs1 =
m_base[1]->GetNumModes();
273 for(i = 0; i < order_e; ++i)
275 for(j = 0; j < nCoeffs1; j++)
277 outarray[(*map)[i].index - j*order_e] +=
278 mat_gauss->GetPtr()[order_e - 1 - j]*
279 (*map)[i].sign*edgeCoeffs[i];
286 nCoeffs0 =
m_base[0]->GetNumModes();
288 for(i = 0; i < order_e; ++i)
290 for(j = 0; j < nCoeffs0; j++)
292 outarray[(*map)[i].index + j] +=
293 mat_gauss->GetPtr()[j]*
294 (*map)[i].sign*edgeCoeffs[i];
300 ASSERTL0(
false,
"edge value (< 3) is out of range");
314 for(i = 0; i < nedges; ++i)
316 EdgeExp[i]->SetCoeffsToOrientation(
GetEorient(i),
318 e_tmp = inout + cnt);
339 for(e = 0; e < nedges; ++e)
341 order_e = EdgeExp[e]->GetNcoeffs();
342 nquad_e = EdgeExp[e]->GetNumPoints(0);
349 for(i = 0; i < order_e; ++i)
351 edgeCoeffs[i] = inarray[i+cnt];
355 EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
371 Vmath::Vmul(nquad_e, normals[dir], 1, edgePhys, 1, edgePhys, 1);
392 for(e = 0; e < nedges; ++e)
394 nquad_e = EdgeExp[e]->GetNumPoints(0);
400 EdgeExp[e]->BwdTrans(edgeCoeffs[e], edgePhys);
402 Vmath::Vmul(nquad_e, normals[dir], 1, edgePhys, 1, edgePhys, 1);
425 int order_e = EdgeExp->GetNcoeffs();
426 int nquad_e = EdgeExp->GetNumPoints(0);
436 StdRegions::VarCoeffMap::const_iterator x;
439 if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
443 edge, EdgeExp, x->second, work);
444 Vmath::Vmul(nquad_e, work, 1, edgePhys, 1, edgePhys, 1);
447 EdgeExp->IProductWRTBase(edgePhys, coeff);
450 for(i = 0; i < order_e; ++i)
452 outarray[map[i]] += sign[i]*coeff[i];
466 ASSERTL0(&inarray[0] != &outarray[0],
467 "Input and output arrays use the same memory");
469 int e, cnt, order_e, nedges =
GetNedges();
474 for(e = 0; e < nedges; ++e)
476 order_e = EdgeExp[e]->GetNcoeffs();
480 Vmath::Vcopy(order_e, tmp = inarray + cnt, 1, edgeCoeffs, 1);
481 EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
483 tau, e, EdgeExp, edgePhys, dirForcing, outarray);
500 int nquad_e = EdgeExp[edge]->GetNumPoints(0);
501 int order_e = EdgeExp[edge]->GetNcoeffs();
532 StdRegions::VarCoeffMap::const_iterator x;
534 if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
538 edge, EdgeExp[edge], x->second, work);
539 Vmath::Vmul(nquad_e, work, 1, edgePhys, 1, edgePhys, 1);
545 EdgeExp[edge]->IProductWRTBase(edgePhys, outcoeff);
547 for(i = 0; i < order_e; ++i)
549 outarray[emap[i]] += sign[i] * tau * outcoeff[i];
558 for(n = 0; n < coordim; ++n)
560 Vmath::Vmul(nquad_e, normals[n], 1, edgePhys, 1, inval, 1);
576 EdgeExp[edge]->IProductWRTBase(inval, outcoeff);
579 for(i = 0; i < ncoeffs; ++i)
582 for(j = 0; j < order_e; ++j)
584 tmpcoeff[i] += invMass(i,emap[j])*sign[j]*outcoeff[j];
588 if(varcoeffs.find(VarCoeff[n]) != varcoeffs.end())
592 Coeffs = Coeffs + Dmat*Tmpcoeff;
598 Coeffs = Coeffs + Dmat*Tmpcoeff;
624 for (
unsigned int i = 0; i < EdgeExp->GetNcoeffs(); ++i)
626 edgetmp[i] = tmp[emap[i]];
630 EdgeExp->BwdTrans(edgetmp, outarray);
651 "HybridDGHelmholtz matrix not set up "
652 "for non boundary-interior expansions");
671 DNekMat LocMat(ncoeffs,ncoeffs);
681 for(i=0; i < coordim; ++i)
690 Mat = Mat + DmatL*invMass*
Transpose(DmatR);
695 Mat = Mat + Dmat*invMass*
Transpose(Dmat);
702 Mat = Mat + lambdaval*Mass;
705 for(i = 0; i < nedges; ++i)
709 order_e = EdgeExp->GetNcoeffs();
710 int nq = EdgeExp->GetNumPoints(0);
724 for(j = 0; j < order_e; ++j)
726 for(k = 0; k < order_e; ++k)
728 Mat(emap[j],emap[k]) = Mat(emap[j],emap[k]) + tau*sign[j]*sign[k]*
eMass(j,k);
762 for(i = 0; i < nedges; ++i)
770 for(j = 0; j < nbndry; ++j)
787 for(k = 0; k < ncoeffs; ++k)
828 for(i = 0; i < nedges; ++i)
850 ASSERTL0(
false,
"Direction not known");
857 for(j = 0; j < nbndry; ++j)
863 for(k = 0; k < ncoeffs; ++k)
865 Ulam[k] = lamToU(k,j);
882 Vmath::Vcopy(ncoeffs,&ulam[0],1,&(Qmat.GetPtr())[0]+j*ncoeffs,1);
890 int order_e, nquad_e;
931 for(i = 0; i < nedges; ++i)
937 for(i = 0; i < nbndry; ++i)
945 for(e = 0; e < nedges; ++e)
947 order_e = EdgeExp[e]->GetNcoeffs();
948 nquad_e = EdgeExp[e]->GetNumPoints(0);
963 StdRegions::VarCoeffMap::const_iterator x;
968 for(j = 0; j < order_e; ++j)
970 edgeCoeffs[j] = sign[j]*(*LamToQ[0])(emap[j],i);
973 EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
982 Vmath::Vmul(nquad_e, normals[0], 1, edgePhys, 1, work, 1);
985 for(j = 0; j < order_e; ++j)
987 edgeCoeffs[j] = sign[j]*(*LamToQ[1])(emap[j],i);
990 EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1006 for(j = 0; j < order_e; ++j)
1008 edgeCoeffs[j] = sign[j]*(*LamToQ[2])(emap[j],i);
1011 EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1031 for(j = 0; j < order_e; ++j)
1033 edgeCoeffs[j] = sign[j]*LamToU(emap[j],i) - lam[cnt+j];
1036 EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1039 if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
1042 Vmath::Vmul(nquad_e,varcoeff_work,1,edgePhys,1,edgePhys,1);
1048 EdgeExp[e]->IProductWRTBase(work, edgeCoeffs);
1050 EdgeExp[e]->SetCoeffsToOrientation(edgeCoeffs, edgedir);
1052 for(j = 0; j < order_e; ++j)
1054 BndMat(cnt+j,i) = edgeCoeffs[j];
1081 &(lmat->GetPtr())[0],1);
1087 ASSERTL0(
false,
"This matrix type cannot be generated from this class");
1125 Out_d = InvMass*Coeffs;
1138 "Not set up for non boundary-interior expansions");
1139 ASSERTL1(inoutmat->GetRows() == inoutmat->GetColumns(),
1140 "Assuming that input matrix was square");
1144 int order_e = edgeExp->GetNcoeffs();
1153 DNekScalMat &edgemat = *edgeExp->GetLocMatrix(mkey);
1167 int rows = inoutmat->GetRows();
1182 for(i = 0; i < order_e; ++i)
1184 for(j = 0; j < nbndry; ++j)
1186 if(map[i] == bmap[j])
1192 ASSERTL1(j != nbndry,
"Did not find number in map");
1202 for(i = 0; i < edge; ++i)
1207 for(i = 0; i < order_e; ++i)
1214 switch(edgeExp->GetBasis(0)->GetBasisType())
1217 reverse( map.get() , map.get()+order_e);
1220 reverse( map.get() , map.get()+order_e);
1224 swap(map[0],map[1]);
1225 for(i = 3; i < order_e; i+=2)
1232 ASSERTL0(
false,
"Edge boundary type not valid for this method");
1238 ASSERTL0(
false,
"Could not identify matrix type from dimension");
1241 for(i = 0; i < order_e; ++i)
1244 for(j = 0; j < order_e; ++j)
1247 (*inoutmat)(id1,id2) += edgemat(i,j)*sign[i]*sign[j];
1263 "Not set up for non boundary-interior expansions");
1266 int order_e = edgeExp->GetNcoeffs();
1275 DNekScalMat &edgemat = *edgeExp->GetLocMatrix(mkey);
1281 for (i = 0; i < order_e; ++i)
1283 vEdgeCoeffs[i] = coeffs[map[i]]*sign[i];
1287 vEdgeCoeffs = edgemat * vEdgeCoeffs;
1289 for (i = 0; i < order_e; ++i)
1291 coeffs[map[i]] = vEdgeCoeffs[i]*sign[i];
1301 int nVerts, vid1, vid2, vMap1, vMap2;
1308 nVerts, nVerts, 0.0, storage);
1309 DNekMat &VertexMat = (*m_vertexmatrix);
1311 for (vid1 = 0; vid1 < nVerts; ++vid1)
1315 for (vid2 = 0; vid2 < nVerts; ++vid2)
1318 VertexValue = (*r_bnd)(vMap1, vMap2);
1319 VertexMat.SetValue(vid1, vid2, VertexValue);
1323 return m_vertexmatrix;
1338 map<int, int> invmap;
1339 for (j = 0; j < nBndCoeffs; ++j)
1341 invmap[bmap[j]] = j;
1357 for (n = 0; n < nEdgeCoeffs; ++n)
1359 edgemaparray[n] = invmap[maparray[n]];
1362 return edgemaparray;
1372 std::map<int, StdRegions::NormalVector>::const_iterator x;
1375 "Edge normal not computed.");
1406 if (idmap.num_elements() != nq0)
1414 for (
int i = 0; i < nq0; ++i)
1422 for (
int i = 0; i < nq0; ++i)
1429 ASSERTL0(
false,
"Unknown orientation");
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
NekDouble GetConstFactor(const ConstFactorType &factor) const
#define ASSERTL0(condition, msg)
const NormalVector & GetEdgeNormal(const int edge) const
const StdRegions::NormalVector & v_GetEdgeNormal(const int edge) const
const ConstFactorMap & GetConstFactors() const
const VarCoeffMap & GetVarCoeffs() const
boost::shared_ptr< IndexMapValues > IndexMapValuesSharedPtr
MatrixType GetMatrixType() const
const VarCoeffMap GetVarCoeffAsMap(const VarCoeffType &coeff) const
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
bool HasVarCoeff(const StdRegions::VarCoeffType &coeff) const
#define sign(a, b)
return the sign(b)*a
void GetPhysEdgeVarCoeffsFromElement(const int edge, ExpansionSharedPtr &EdgeExp, const Array< OneD, const NekDouble > &varcoeff, Array< OneD, NekDouble > &outarray)
void AddNormTraceInt(const int dir, Array< OneD, ExpansionSharedPtr > &EdgeExp, Array< OneD, Array< OneD, NekDouble > > &edgeCoeffs, Array< OneD, NekDouble > &outarray)
void GetEdgeInteriorMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
void InterpCoeff1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
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...
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
const StdRegions::NormalVector & v_GetSurfaceNormal(const int id) const
const Array< OneD, const NekDouble > & GetVarCoeff(const StdRegions::VarCoeffType &coeff) const
std::map< int, bool > m_negatedNormals
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*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
Principle Modified Functions .
std::map< int, StdRegions::NormalVector > m_edgeNormals
Lagrange Polynomials using the Gauss points .
StdRegions::Orientation GetEorient(int edge)
void SetTraceToGeomOrientation(Array< OneD, ExpansionSharedPtr > &EdgeExp, Array< OneD, NekDouble > &inout)
virtual void v_AddRobinEdgeContribution(const int edgeid, const Array< OneD, const NekDouble > &primCoeffs, Array< OneD, NekDouble > &coeffs)
std::map< ConstFactorType, NekDouble > ConstFactorMap
virtual DNekMatSharedPtr v_BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd)
int GetEdgeNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th edge.
virtual void v_AddRobinMassMatrix(const int edgeid, const Array< OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat)
boost::shared_ptr< DNekMat > DNekMatSharedPtr
void AddEdgeBoundaryInt(const int edge, ExpansionSharedPtr &EdgeExp, Array< OneD, NekDouble > &edgePhys, Array< OneD, NekDouble > &outarray, const StdRegions::VarCoeffMap &varcoeffs=StdRegions::NullVarCoeffMap)
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
int NumBndryCoeffs(void) const
virtual void v_NegateEdgeNormal(const int edge)
bool IsBoundaryInteriorExpansion()
std::vector< bool > m_requireNeg
SpatialDomains::Geometry2DSharedPtr GetGeom2D() const
int GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Expansion1DSharedPtr GetEdgeExp(int edge, bool SetUpNormal=true)
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
virtual void v_AddEdgeNormBoundaryInt(const int edge, const ExpansionSharedPtr &EdgeExp, const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
void Neg(int n, T *x, const int incx)
Negate x = -x.
boost::shared_ptr< Expansion > ExpansionSharedPtr
void AddHDGHelmholtzEdgeTerms(const NekDouble tau, const int edge, Array< OneD, ExpansionSharedPtr > &EdgeExp, Array< OneD, NekDouble > &edgePhys, const StdRegions::VarCoeffMap &dirForcing, Array< OneD, NekDouble > &outarray)
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
virtual bool v_EdgeNormalNegated(const int edge)
int GetNverts() const
This function returns the number of vertices of the expansion domain.
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
void AddHDGHelmholtzTraceTerms(const NekDouble tau, const Array< OneD, const NekDouble > &inarray, Array< OneD, ExpansionSharedPtr > &EdgeExp, const StdRegions::VarCoeffMap &dirForcing, Array< OneD, NekDouble > &outarray)
void GetEdgeToElementMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, int P=-1)
virtual void v_DGDeriv(const int dir, const Array< OneD, const NekDouble > &incoeffs, Array< OneD, ExpansionSharedPtr > &EdgeExp, Array< OneD, Array< OneD, NekDouble > > &edgeCoeffs, Array< OneD, NekDouble > &out_d)
std::vector< Expansion1DWeakPtr > m_edgeExp
virtual StdRegions::Orientation v_GetEorient(int edge)
virtual void v_SetUpPhysNormals(const int edge)
int NumDGBndryCoeffs(void) const
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
LibUtilities::NekManager< StdMatrixKey, DNekMat, StdMatrixKey::opLess > m_stdMatrixManager
void Zero(int n, T *x, const int incx)
Zero vector.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Array< OneD, LibUtilities::BasisSharedPtr > m_base
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
void AddEdgeNormBoundaryInt(const int edge, const boost::shared_ptr< Expansion > &EdgeExp, const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
int GetNedges() const
This function returns the number of edges of the expansion domain.
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Describes the specification for a Basis.
Array< OneD, unsigned int > v_GetEdgeInverseBoundaryMap(int eid)
void ReOrientEdgePhysMap(const int nvert, const StdRegions::Orientation orient, const int nq0, Array< OneD, int > &idmap)
void ComputeEdgeNormal(const int edge)
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.
static ConstFactorMap NullConstFactorMap
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space...