diff --git a/library/LibUtilities/BasicUtils/SessionReader.cpp b/library/LibUtilities/BasicUtils/SessionReader.cpp
index 3cc8baf397edc9837f50f7aba6ec8b2a886a6f12..8245c23ed127d43297db0d675ae0c05bfb63d47c 100644
--- a/library/LibUtilities/BasicUtils/SessionReader.cpp
+++ b/library/LibUtilities/BasicUtils/SessionReader.cpp
@@ -842,7 +842,7 @@ namespace Nektar
m_comm = GetCommFactory().CreateInstance(vCommModule,argc,argv);
- // If running in parallel change the default global sys soln type
+ //If running in parallel change the default global sys soln type
if (m_comm->GetSize() > 1)
{
m_solverInfoDefaults["GLOBALSYSSOLN"] = "IterativeStaticCond";
diff --git a/library/LocalRegions/HexExp.cpp b/library/LocalRegions/HexExp.cpp
index cc26474f4af36daecdbcb0ab3cddf055a9f5bfa8..edee5acf897fbe9534a2c8d24e946c8df0e4482a 100644
--- a/library/LocalRegions/HexExp.cpp
+++ b/library/LocalRegions/HexExp.cpp
@@ -521,12 +521,18 @@ namespace Nektar
*/
NekDouble HexExp::v_PhysEvaluate(
const Array<OneD, const NekDouble> &coord)
+ {
+ v_PhysEvaluate(coord,m_phys);
+ }
+
+ NekDouble HexExp::v_PhysEvaluate(
+ const Array<OneD, const NekDouble> &coord, const Array<OneD, const NekDouble> & physvals)
{
Array<OneD,NekDouble> Lcoord = Array<OneD,NekDouble>(3);
ASSERTL0(m_geom,"m_geom not defined");
m_geom->GetLocCoords(coord,Lcoord);
- return StdHexExp::v_PhysEvaluate(Lcoord);
+ return StdHexExp::v_PhysEvaluate(Lcoord, physvals);
}
diff --git a/library/LocalRegions/HexExp.h b/library/LocalRegions/HexExp.h
index e7487e84693f947f56d94869414d097e838cb16a..5a5521de9f975d1bd90dd2c51952830a7167d1fc 100644
--- a/library/LocalRegions/HexExp.h
+++ b/library/LocalRegions/HexExp.h
@@ -135,6 +135,9 @@ namespace Nektar
virtual NekDouble v_PhysEvaluate(
const Array<OneD, const NekDouble> &coords);
+ virtual NekDouble v_PhysEvaluate(
+ const Array<OneD, const NekDouble> &coords, const Array<OneD, const NekDouble> & physvals);
+
/// Retrieve the local coordinates of each quadrature point.
virtual void v_GetCoords( Array<OneD,NekDouble> &coords_1,
Array<OneD,NekDouble> &coords_2,
diff --git a/library/LocalRegions/PrismExp.cpp b/library/LocalRegions/PrismExp.cpp
index b2f196bd4391eb1e8bcbb9d5f5143bdcb7d6d412..36a9fe89e85646e8919a8c2789c165423671225d 100644
--- a/library/LocalRegions/PrismExp.cpp
+++ b/library/LocalRegions/PrismExp.cpp
@@ -443,6 +443,11 @@ namespace Nektar
NekDouble PrismExp::PhysEvaluate(const Array<OneD, const NekDouble> &coord)
+ {
+ PhysEvaluate(coord,m_phys);
+ }
+
+ NekDouble PrismExp::PhysEvaluate(const Array<OneD, const NekDouble> &coord, const Array<OneD, const NekDouble> & physvals)
{
Array<OneD,NekDouble> Lcoord = Array<OneD,NekDouble>(3);
@@ -450,7 +455,7 @@ namespace Nektar
m_geom->GetLocCoords(coord, Lcoord);
- return StdPrismExp::PhysEvaluate(Lcoord);
+ return StdPrismExp::PhysEvaluate(Lcoord,physvals);
}
DNekMatSharedPtr PrismExp::CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
diff --git a/library/LocalRegions/PrismExp.h b/library/LocalRegions/PrismExp.h
index 33aa24d022d0e17a3f52a0c8ca12e39c6147b484..88a487bc746e31e369664710a63ba2f912dc7b40 100644
--- a/library/LocalRegions/PrismExp.h
+++ b/library/LocalRegions/PrismExp.h
@@ -109,6 +109,8 @@ namespace Nektar
LOCAL_REGIONS_EXPORT NekDouble PhysEvaluate(const Array<OneD, const NekDouble> &coord);
+ LOCAL_REGIONS_EXPORT NekDouble PhysEvaluate(const Array<OneD, const NekDouble> &coord, const Array<OneD, const NekDouble> & physvals);
+
//-----------------------------
// Differentiation Methods
//-----------------------------
@@ -271,6 +273,11 @@ namespace Nektar
return PhysEvaluate(coords);
}
+ virtual NekDouble v_PhysEvaluate(const Array<OneD, const NekDouble> &coords, const Array<OneD, const NekDouble> & physvals)
+ {
+ return PhysEvaluate(coords,physvals);
+ }
+
virtual NekDouble v_Linf(const Array<OneD, const NekDouble> &sol)
{
return Linf(sol);
diff --git a/library/LocalRegions/QuadExp.cpp b/library/LocalRegions/QuadExp.cpp
index 0eaf8e1004ec0a02995e3bec4c85f1e06295d620..4a455b5d0a4d957a64e6645c188bc3348802b301 100644
--- a/library/LocalRegions/QuadExp.cpp
+++ b/library/LocalRegions/QuadExp.cpp
@@ -1314,12 +1314,17 @@ namespace Nektar
NekDouble QuadExp::PhysEvaluate(const Array<OneD, const NekDouble> &coord)
{
- Array<OneD,NekDouble> Lcoord = Array<OneD,NekDouble>(2);
+ PhysEvaluate(coord,m_phys);
+ }
+ NekDouble QuadExp::PhysEvaluate(const Array<OneD, const NekDouble> &coord, const Array<OneD, const NekDouble> & physvals)
+ {
+ Array<OneD,NekDouble> Lcoord = Array<OneD,NekDouble>(2);
+
ASSERTL0(m_geom,"m_geom not defined");
m_geom->GetLocCoords(coord,Lcoord);
- return StdQuadExp::PhysEvaluate(Lcoord);
+ return StdQuadExp::PhysEvaluate(Lcoord, physvals);
}
DNekMatSharedPtr QuadExp::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
diff --git a/library/LocalRegions/QuadExp.h b/library/LocalRegions/QuadExp.h
index 2590102276ab7b13ec0c291c0375f15204c4bf6f..9fa0d43dbdad695be6c769ced98e090b750095fc 100644
--- a/library/LocalRegions/QuadExp.h
+++ b/library/LocalRegions/QuadExp.h
@@ -171,14 +171,14 @@ namespace Nektar
/** \brief Forward transform from physical quadrature space
stored in \a inarray and evaluate the expansion coefficients and
store in \a (this)->_coeffs */
- LOCAL_REGIONS_EXPORT void FwdTrans(const Array<OneD, const NekDouble> &inarray,
- Array<OneD, NekDouble> &outarray);
+ LOCAL_REGIONS_EXPORT void FwdTrans(const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &outarray);
- LOCAL_REGIONS_EXPORT void FwdTrans_BndConstrained(const Array<OneD, const NekDouble>& inarray,
- Array<OneD, NekDouble> &outarray);
+ LOCAL_REGIONS_EXPORT void FwdTrans_BndConstrained(const Array<OneD, const NekDouble>& inarray, Array<OneD, NekDouble> &outarray);
LOCAL_REGIONS_EXPORT NekDouble PhysEvaluate(const Array<OneD, const NekDouble> &coord);
+ LOCAL_REGIONS_EXPORT NekDouble PhysEvaluate(const Array<OneD, const NekDouble> &coord, const Array<OneD, const NekDouble> & physvals);
+
LOCAL_REGIONS_EXPORT void GetEdgePhysVals(const int edge, const Array<OneD, const NekDouble> &inarray, Array<OneD,NekDouble> &outarray);
/** \brief Extract the physical values along edge \a edge
@@ -481,6 +481,11 @@ namespace Nektar
{
return PhysEvaluate(coords);
}
+
+ virtual NekDouble v_PhysEvaluate(const Array<OneD, const NekDouble> &coords, const Array<OneD, const NekDouble> & physvals)
+ {
+ return PhysEvaluate(coords, physvals);
+ }
virtual NekDouble v_Linf(const Array<OneD, const NekDouble> &sol)
{
diff --git a/library/LocalRegions/SegExp.cpp b/library/LocalRegions/SegExp.cpp
index f52ad61dbf7d9cc6075ff9ff05bd8bb75424a75d..b18a87de20ec04c8e15e9db2787384c037f05364 100644
--- a/library/LocalRegions/SegExp.cpp
+++ b/library/LocalRegions/SegExp.cpp
@@ -939,13 +939,18 @@ cout<<"deps/dx ="<<inarray_d0[i]<<" deps/dy="<<inarray_d1[i]<<endl;
}
NekDouble SegExp::PhysEvaluate(const Array<OneD, const NekDouble>& coord)
+ {
+ PhysEvaluate(coord,m_phys);
+ }
+
+ NekDouble SegExp::PhysEvaluate(const Array<OneD, const NekDouble>& coord, const Array<OneD, const NekDouble> &physvals)
{
Array<OneD,NekDouble> Lcoord = Array<OneD,NekDouble>(1);
ASSERTL0(m_geom,"_geom not defined");
m_geom->GetLocCoords(coord,Lcoord);
- return StdSegExp::PhysEvaluate(Lcoord);
+ return StdSegExp::PhysEvaluate(Lcoord, physvals);
}
DNekScalMatSharedPtr SegExp::CreateMatrix(const MatrixKey &mkey)
diff --git a/library/LocalRegions/SegExp.h b/library/LocalRegions/SegExp.h
index 23ff07857f52f3120de2d1e4cace07ac16c2a938..51394cce9198033a593b473b38e288c1e467362e 100644
--- a/library/LocalRegions/SegExp.h
+++ b/library/LocalRegions/SegExp.h
@@ -214,6 +214,8 @@ namespace Nektar
LOCAL_REGIONS_EXPORT NekDouble PhysEvaluate(const Array<OneD, const NekDouble>& coord);
+ LOCAL_REGIONS_EXPORT NekDouble PhysEvaluate(const Array<OneD, const NekDouble>& coord, const Array<OneD, const NekDouble> & physvals);
+
LOCAL_REGIONS_EXPORT void LaplacianMatrixOp(const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray);
@@ -419,6 +421,11 @@ namespace Nektar
return PhysEvaluate(coords);
}
+ virtual NekDouble v_PhysEvaluate(const Array<OneD, const NekDouble>& coords, const Array<OneD, const NekDouble> & physvals)
+ {
+ return PhysEvaluate(coords, physvals);
+ }
+
/** \brief Virtual function to evaluate the discrete \f$ L_\infty\f$
error \f$ |\epsilon|_\infty = \max |u - u_{exact}|\f$ where \f$
u_{exact}\f$ is given by the array \a sol.
diff --git a/library/LocalRegions/TetExp.cpp b/library/LocalRegions/TetExp.cpp
index c6c4a8e890fd6b262e880d96b1c190cdfcb987cc..3cb57a09014606f007705c4144d699847bdf5281 100644
--- a/library/LocalRegions/TetExp.cpp
+++ b/library/LocalRegions/TetExp.cpp
@@ -412,6 +412,16 @@ namespace Nektar
*/
NekDouble TetExp::v_PhysEvaluate(
const Array<OneD, const NekDouble> &coord)
+ {
+ v_PhysEvaluate(coord,m_phys);
+ }
+ /**
+ * @param coord Physical space coordinate
+ * @returns Evaluation of expansion at given coordinate.
+ */
+ NekDouble TetExp::v_PhysEvaluate(
+ const Array<OneD, const NekDouble> &coord,
+ const Array<OneD, const NekDouble> & physvals)
{
ASSERTL0(m_geom,"m_geom not defined");
@@ -421,7 +431,7 @@ namespace Nektar
m_geom->GetLocCoords(coord,Lcoord);
// Evaluate point in local (eta) coordinates.
- return StdExpansion3D::v_PhysEvaluate(Lcoord);
+ return StdExpansion3D::v_PhysEvaluate(Lcoord,physvals);
}
diff --git a/library/LocalRegions/TetExp.h b/library/LocalRegions/TetExp.h
index 3b04a7cb194ce08b97d3f48c82fffa7679eea669..c6d50e1c5665775ffc6dc39ba2ceeefa0e015de6 100644
--- a/library/LocalRegions/TetExp.h
+++ b/library/LocalRegions/TetExp.h
@@ -102,6 +102,10 @@ namespace Nektar
virtual NekDouble v_PhysEvaluate(
const Array<OneD, const NekDouble> &coords);
+ virtual NekDouble v_PhysEvaluate(
+ const Array<OneD, const NekDouble> &coords,
+ const Array<OneD, const NekDouble> & physvals);
+
/// Get the x,y,z coordinates of each quadrature point.
virtual void v_GetCoords(Array<OneD,NekDouble> &coords_0,
Array<OneD,NekDouble> &coords_1,
diff --git a/library/LocalRegions/TriExp.cpp b/library/LocalRegions/TriExp.cpp
index 8ecfe317ee2cd97d40f4ec57409a6bb3b896e1a5..2318d6fc9a0aa6e7881626d8286055e421127230 100644
--- a/library/LocalRegions/TriExp.cpp
+++ b/library/LocalRegions/TriExp.cpp
@@ -1400,13 +1400,18 @@ namespace Nektar
}
NekDouble TriExp::PhysEvaluate(const Array<OneD, const NekDouble> &coord)
+ {
+ PhysEvaluate(coord,m_phys);
+ }
+
+ NekDouble TriExp::PhysEvaluate(const Array<OneD, const NekDouble> &coord, const Array<OneD, const NekDouble> & physvals)
{
Array<OneD,NekDouble> Lcoord = Array<OneD,NekDouble>(2);
ASSERTL0(m_geom,"m_geom not defined");
m_geom->GetLocCoords(coord,Lcoord);
-
- return StdTriExp::PhysEvaluate(Lcoord);
+
+ return StdTriExp::PhysEvaluate(Lcoord, physvals);
}
DNekScalMatSharedPtr TriExp::CreateMatrix(const MatrixKey &mkey)
diff --git a/library/LocalRegions/TriExp.h b/library/LocalRegions/TriExp.h
index 6ae88abeccfadc734e5af5590b922c772c75203e..ba7ee3d8a68cc7131d515682fe229ae2e658c1cf 100644
--- a/library/LocalRegions/TriExp.h
+++ b/library/LocalRegions/TriExp.h
@@ -165,6 +165,8 @@ namespace Nektar
Array<OneD, NekDouble> &outarray);
LOCAL_REGIONS_EXPORT NekDouble PhysEvaluate(const Array<OneD, const NekDouble> &coord);
+ LOCAL_REGIONS_EXPORT NekDouble PhysEvaluate(const Array<OneD, const NekDouble> &coord, const Array<OneD, const NekDouble> & physvals);
+
/** \brief Extract the physical values along edge \a edge
from \a inarray into \a outarray following the local
edge orientation and point distribution defined by
@@ -471,12 +473,16 @@ namespace Nektar
return PhysEvaluate(coords);
}
+ virtual NekDouble v_PhysEvaluate(const Array<OneD, const NekDouble> &coords, const Array<OneD, const NekDouble> & physvals)
+ {
+ return PhysEvaluate(coords, physvals);
+ }
+
virtual NekDouble v_Linf(const Array<OneD, const NekDouble> &sol)
{
return Linf(sol);
}
-
virtual NekDouble v_Linf()
{
return Linf();
diff --git a/library/MultiRegions/ExpList.cpp b/library/MultiRegions/ExpList.cpp
index 57da7ef06b478dfb9139fca3269c49f1bdbe39a7..552299e769b2a5c303c0720d7e2563f5c81422ca 100644
--- a/library/MultiRegions/ExpList.cpp
+++ b/library/MultiRegions/ExpList.cpp
@@ -1294,12 +1294,25 @@ namespace Nektar
/* region defined by vertices. The do point search.
**/
int ExpList::GetExpIndex(
- const Array<OneD, const NekDouble> &gloCoord)
+ const Array<OneD, const NekDouble> &gloCoord,
+ NekDouble tol)
{
- for (int i = 0; i < (*m_exp).size(); ++i)
+ static int start = 0;
+ // start search at previous element or 0
+ for (int i = start; i < (*m_exp).size(); ++i)
{
- if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoord))
+ if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoord,tol))
+ {
+ start = i;
+ return i;
+ }
+ }
+
+ for (int i = 0; i < start; ++i)
+ {
+ if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoord,tol))
{
+ start = i;
return i;
}
}
diff --git a/library/MultiRegions/ExpList.h b/library/MultiRegions/ExpList.h
index b6cf0b33be0d0f2f72a1fb18bf6e9eedb772838e..c3c129218e15f55eeff65014d20abf8a0fe3c7a0 100644
--- a/library/MultiRegions/ExpList.h
+++ b/library/MultiRegions/ExpList.h
@@ -524,7 +524,7 @@ namespace Nektar
/// This function returns the index of the local elemental
/// expansion containing the arbitrary point given by \a gloCoord.
- MULTI_REGIONS_EXPORT int GetExpIndex(const Array<OneD, const NekDouble> &gloCoord);
+ MULTI_REGIONS_EXPORT int GetExpIndex(const Array<OneD, const NekDouble> &gloCoord, NekDouble tol = 0.0);
/// Get the start offset position for a global list of #m_coeffs
/// correspoinding to element n.
@@ -1037,54 +1037,54 @@ namespace Nektar
Array<OneD, NekDouble> &coord_1,
Array<OneD, NekDouble> &coord_2 = NullNekDouble1DArray);
- virtual void v_GetCoords(NekDouble &x,NekDouble &y,NekDouble &z);
-
- virtual void v_GetCoord(Array<OneD, NekDouble> &coords);
+ virtual void v_GetCoords(NekDouble &x,NekDouble &y,NekDouble &z);
+
+ virtual void v_GetCoord(Array<OneD, NekDouble> &coords);
- virtual void v_SetCoeff(NekDouble val);
-
- virtual void v_SetPhys(NekDouble val);
-
- virtual const SpatialDomains::VertexComponentSharedPtr &v_GetGeom(void) const;
-
- virtual const SpatialDomains::VertexComponentSharedPtr &v_GetVertex(void) const;
-
- 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, bool UseContCoeffs = false);
-
- virtual void v_PhysDeriv(const int dir,
- const Array<OneD, const NekDouble> &inarray,
- Array<OneD, NekDouble> &out_d, bool UseContCoeffs = false);
-
- virtual void v_PhysDeriv(Direction edir,
- const Array<OneD, const NekDouble> &inarray,
- Array<OneD, NekDouble> &out_d, bool UseContCoeffs = false);
-
- virtual void v_HomogeneousFwdTrans(const Array<OneD, const NekDouble> &inarray,
- Array<OneD, NekDouble> &outarray,
- bool UseContCoeffs = false);
-
- virtual void v_HomogeneousBwdTrans(const Array<OneD, const NekDouble> &inarray,
- Array<OneD, NekDouble> &outarray,
- bool UseContCoeffs = false);
-
- virtual void v_Dealiasing(Array<OneD, NekDouble> &outarray, bool UseContCoeffs = false);
-
- virtual void v_GetBCValues(Array<OneD, NekDouble> &BndVals,
- const Array<OneD, NekDouble> &TotField,
- int BndID);
-
- virtual void v_NormVectorIProductWRTBase(Array<OneD, const NekDouble> &V1,
- Array<OneD, const NekDouble> &V2,
- Array<OneD, NekDouble> &outarray,
- int BndID);
-
+ virtual void v_SetCoeff(NekDouble val);
+
+ virtual void v_SetPhys(NekDouble val);
+
+ virtual const SpatialDomains::VertexComponentSharedPtr &v_GetGeom(void) const;
+
+ virtual const SpatialDomains::VertexComponentSharedPtr &v_GetVertex(void) const;
+
+ 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, bool UseContCoeffs = false);
+
+ virtual void v_PhysDeriv(const int dir,
+ const Array<OneD, const NekDouble> &inarray,
+ Array<OneD, NekDouble> &out_d, bool UseContCoeffs = false);
+
+ virtual void v_PhysDeriv(Direction edir,
+ const Array<OneD, const NekDouble> &inarray,
+ Array<OneD, NekDouble> &out_d, bool UseContCoeffs = false);
+
+ virtual void v_HomogeneousFwdTrans(const Array<OneD, const NekDouble> &inarray,
+ Array<OneD, NekDouble> &outarray,
+ bool UseContCoeffs = false);
+
+ virtual void v_HomogeneousBwdTrans(const Array<OneD, const NekDouble> &inarray,
+ Array<OneD, NekDouble> &outarray,
+ bool UseContCoeffs = false);
+
+ virtual void v_Dealiasing(Array<OneD, NekDouble> &outarray, bool UseContCoeffs = false);
+
+ virtual void v_GetBCValues(Array<OneD, NekDouble> &BndVals,
+ const Array<OneD, NekDouble> &TotField,
+ int BndID);
+
+ virtual void v_NormVectorIProductWRTBase(Array<OneD, const NekDouble> &V1,
+ Array<OneD, const NekDouble> &V2,
+ Array<OneD, NekDouble> &outarray,
+ int BndID);
+
virtual void v_SetUpPhysNormals();
-
+
virtual void v_SetUpTangents();
-
+
virtual void v_GetBoundaryToElmtMap(Array<OneD, int> &ElmtID,
Array<OneD,int> &EdgeID);
diff --git a/library/MultiRegions/LocalToGlobalC0ContMap.cpp b/library/MultiRegions/LocalToGlobalC0ContMap.cpp
index e6f7dfc933380fb582a70f60d3430573495b9ccc..a997064212716ca76b7ff5062b538d19d71fadef 100644
--- a/library/MultiRegions/LocalToGlobalC0ContMap.cpp
+++ b/library/MultiRegions/LocalToGlobalC0ContMap.cpp
@@ -1003,11 +1003,23 @@ namespace Nektar
systemSingular = (s == 1 ? true : false);
if(systemSingular == true && checkIfSystemSingular)
{
- //last region i and j=0 edge
- bndSegExp = boost::dynamic_pointer_cast<LocalRegions::SegExp>(bndCondExp[bndCondExp.num_elements()-1]->GetExp(0));
+ if(m_session->DefinesParameter("SingularElement"))
+ {
+ int s_eid = m_session->GetParameter("SingularElement");
+
+ ASSERTL1(s_eid < locExpVector.size(),"SingularElement Parameter is too large");
+
+ meshVertId = locExpVector[s_eid]->GetGeom2D()->GetVid(0);
+ }
+ else
+ {
+ //last region i and j=0 edge
+ bndSegExp = boost::dynamic_pointer_cast<LocalRegions::SegExp>(bndCondExp[bndCondExp.num_elements()-1]->GetExp(0));
+
+ //first vertex 0 of the edge
+ meshVertId = (bndSegExp->GetGeom1D())->GetVid(0);
+ }
- //first vertex 0 of the edge
- meshVertId = (bndSegExp->GetGeom1D())->GetVid(0);
if(ReorderedGraphVertId[0].count(meshVertId) == 0)
{
ReorderedGraphVertId[0][meshVertId] = graphVertId++;
diff --git a/library/SpatialDomains/Geometry.h b/library/SpatialDomains/Geometry.h
index bdc8b2b865ec44e60d7a2e2f779f3d4cbbcec0ca..b70e22bedb292877eb2375f41b94bd3cb7d0c577 100644
--- a/library/SpatialDomains/Geometry.h
+++ b/library/SpatialDomains/Geometry.h
@@ -166,9 +166,9 @@ namespace Nektar
}
inline bool ContainsPoint(
- const Array<OneD, const NekDouble> &gloCoord)
+ const Array<OneD, const NekDouble> &gloCoord, NekDouble tol = 0.0)
{
- return v_ContainsPoint(gloCoord);
+ return v_ContainsPoint(gloCoord,tol);
}
protected:
@@ -242,7 +242,7 @@ namespace Nektar
}
virtual bool v_ContainsPoint(
- const Array<OneD, const NekDouble> &gloCoord)
+ const Array<OneD, const NekDouble> &gloCoord, NekDouble tol = 0.0)
{
NEKERROR(ErrorUtil::efatal,
"This function has not been defined for this geometry");
diff --git a/library/SpatialDomains/Geometry2D.h b/library/SpatialDomains/Geometry2D.h
index 93fc96d629e7b0052c734a52259498fe361532b4..c4c42079ae38b1af9ecc8724cb1c3a55e2e0487c 100644
--- a/library/SpatialDomains/Geometry2D.h
+++ b/library/SpatialDomains/Geometry2D.h
@@ -353,6 +353,14 @@ namespace Nektar
{
return 2;
}
+
+ virtual bool v_ContainsPoint(
+ const Array<OneD, const NekDouble> &gloCoord, NekDouble tol = 0.0)
+ {
+ NEKERROR(ErrorUtil::efatal,
+ "This function has not been defined for this geometry");
+ return false;
+ }
};
diff --git a/library/SpatialDomains/HexGeom.cpp b/library/SpatialDomains/HexGeom.cpp
index 402349378dbee63484fc684d66fa8ef2bbfb801b..c3580afd197a84e853931ac13c523123b8e787e7 100644
--- a/library/SpatialDomains/HexGeom.cpp
+++ b/library/SpatialDomains/HexGeom.cpp
@@ -792,16 +792,16 @@ namespace Nektar
}
bool HexGeom::v_ContainsPoint(
- const Array<OneD, const NekDouble> &gloCoord)
+ const Array<OneD, const NekDouble> &gloCoord, NekDouble tol)
{
ASSERTL1(gloCoord.num_elements() == 3,
"Three dimensional geometry expects three coordinates.");
Array<OneD,NekDouble> stdCoord(GetCoordim(),0.0);
GetLocCoords(gloCoord, stdCoord);
- if (stdCoord[0] >= -1 && stdCoord[0] <= 1
- && stdCoord[1] >= -1 && stdCoord[1] <= 1
- && stdCoord[2] >= -1 && stdCoord[2] <= 1)
+ if (stdCoord[0] >= -(1+tol) && stdCoord[0] <= 1+tol
+ && stdCoord[1] >= -(1+tol) && stdCoord[1] <= 1+tol
+ && stdCoord[2] >= -(1+tol) && stdCoord[2] <= 1+tol)
{
return true;
}
diff --git a/library/SpatialDomains/HexGeom.h b/library/SpatialDomains/HexGeom.h
index a189f20427dfaca99a2b8f8193818359ad447b72..6185aaa48e0cdf2a20e990fef13e1e5443467416 100644
--- a/library/SpatialDomains/HexGeom.h
+++ b/library/SpatialDomains/HexGeom.h
@@ -318,8 +318,7 @@ namespace Nektar
return kNedges;
}
- virtual bool v_ContainsPoint(
- const Array<OneD, const NekDouble> &gloCoord);
+ virtual bool v_ContainsPoint(const Array<OneD, const NekDouble> &gloCoord, NekDouble tol = 0.0);
};
diff --git a/library/SpatialDomains/QuadGeom.cpp b/library/SpatialDomains/QuadGeom.cpp
index a14de53aaac7955c32fad8c0a82165fadf0c0a5d..61bd74f9519f11b13b4243d80af5b4f4d461b123 100644
--- a/library/SpatialDomains/QuadGeom.cpp
+++ b/library/SpatialDomains/QuadGeom.cpp
@@ -407,15 +407,15 @@ namespace Nektar
}
bool QuadGeom::v_ContainsPoint(
- const Array<OneD, const NekDouble> &gloCoord)
+ const Array<OneD, const NekDouble> &gloCoord, NekDouble tol)
{
ASSERTL1(gloCoord.num_elements() >= 2,
"Two dimensional geometry expects at least two coordinates.");
Array<OneD,NekDouble> stdCoord(GetCoordim(),0.0);
GetLocCoords(gloCoord, stdCoord);
- if (stdCoord[0] >= -1 && stdCoord[1] >= -1
- && stdCoord[0] <= 1 && stdCoord[1] <= 1)
+ if (stdCoord[0] >= -(1+tol) && stdCoord[1] >= -(1+tol)
+ && stdCoord[0] <= (1+tol) && stdCoord[1] <= (1+tol))
{
return true;
}
diff --git a/library/SpatialDomains/QuadGeom.h b/library/SpatialDomains/QuadGeom.h
index e917463e2f0d580c3156d9252d075156c0314866..da4788a3b957f83b1c33dad1f9442930661bc393 100644
--- a/library/SpatialDomains/QuadGeom.h
+++ b/library/SpatialDomains/QuadGeom.h
@@ -319,7 +319,7 @@ namespace Nektar
}
virtual bool v_ContainsPoint(
- const Array<OneD, const NekDouble> &gloCoord);
+ const Array<OneD, const NekDouble> &gloCoord, NekDouble tol = 0.0);
};
}; //end of namespace
diff --git a/library/SpatialDomains/SegGeom.cpp b/library/SpatialDomains/SegGeom.cpp
index 1ffaa57a23fd652e492d1962776a72240639b428..f804f9050154502f09f959b00e5db10b5f3a9c85 100644
--- a/library/SpatialDomains/SegGeom.cpp
+++ b/library/SpatialDomains/SegGeom.cpp
@@ -393,12 +393,11 @@ namespace Nektar
}
}
- bool SegGeom::v_ContainsPoint(
- const Array<OneD, const NekDouble> &gloCoord)
+ bool SegGeom::v_ContainsPoint(const Array<OneD, const NekDouble> &gloCoord, NekDouble tol)
{
Array<OneD,NekDouble> stdCoord(GetCoordim(),0.0);
GetLocCoords(gloCoord, stdCoord);
- if (stdCoord[0] >= -1 && stdCoord[0] <= 1)
+ if (stdCoord[0] >= -(1+tol) && stdCoord[0] <= 1+tol)
{
return true;
}
diff --git a/library/SpatialDomains/SegGeom.h b/library/SpatialDomains/SegGeom.h
index 3167b9a15c5d24f1cb2ac5839013076fc9da4372..67156034cdfc91293e608cc932c214bbc1a02940 100644
--- a/library/SpatialDomains/SegGeom.h
+++ b/library/SpatialDomains/SegGeom.h
@@ -259,7 +259,7 @@ namespace Nektar
}
virtual bool v_ContainsPoint(
- const Array<OneD, const NekDouble> &gloCoord);
+ const Array<OneD, const NekDouble> &gloCoord, NekDouble tol = 0.0);
};
// shorthand for boost pointer
diff --git a/library/SpatialDomains/TetGeom.cpp b/library/SpatialDomains/TetGeom.cpp
index 4e55d0f8e3527e7a8ea3a46a43213fa842eed11a..13c030b8d3d827ceecebeaa618c884651878d259 100644
--- a/library/SpatialDomains/TetGeom.cpp
+++ b/library/SpatialDomains/TetGeom.cpp
@@ -726,7 +726,7 @@ namespace Nektar
* within this tetrahedral geometry.
*/
bool TetGeom::v_ContainsPoint(
- const Array<OneD, const NekDouble> &gloCoord)
+ const Array<OneD, const NekDouble> &gloCoord, NekDouble tol)
{
// Validation checks
ASSERTL1(gloCoord.num_elements() == 3,
@@ -737,8 +737,9 @@ namespace Nektar
GetLocCoords(gloCoord, locCoord);
// Check local coordinate is within [-1,1]^3 bounds.
- if (locCoord[0] >= -1 && locCoord[1] >= -1 && locCoord[2] >= -1
- && locCoord[0] <= 1 && locCoord[1] <= 1 && locCoord[2] <= 1)
+ if (locCoord[0] >= -(1+tol) && locCoord[1] >= -(1+tol)
+ && locCoord[2] >= -(1+tol)
+ && locCoord[0] + locCoord[1] + locCoord[2] <= tol)
{
return true;
}
diff --git a/library/SpatialDomains/TetGeom.h b/library/SpatialDomains/TetGeom.h
index baf659564205aa98cb061e1a984b438e0a5ed220..5dd30dacb91f44833d320582cae496f2cc1d7c95 100644
--- a/library/SpatialDomains/TetGeom.h
+++ b/library/SpatialDomains/TetGeom.h
@@ -326,7 +326,7 @@ namespace Nektar
}
virtual bool v_ContainsPoint(
- const Array<OneD, const NekDouble> &gloCoord);
+ const Array<OneD, const NekDouble> &gloCoord, NekDouble tol = 0.0);
};
diff --git a/library/SpatialDomains/TriGeom.cpp b/library/SpatialDomains/TriGeom.cpp
index 557daeadee6f4f1ff51135b4c19a1ea38da2c9c8..08f8e43b921fe651b5cd1b2be3e101bdbcacea62 100644
--- a/library/SpatialDomains/TriGeom.cpp
+++ b/library/SpatialDomains/TriGeom.cpp
@@ -361,16 +361,15 @@ namespace Nektar
return returnval;
}
- bool TriGeom::v_ContainsPoint(
- const Array<OneD, const NekDouble> &gloCoord)
+ bool TriGeom::v_ContainsPoint(const Array<OneD, const NekDouble> &gloCoord, NekDouble tol)
{
ASSERTL1(gloCoord.num_elements() >= 2,
"Two dimensional geometry expects at least two coordinates.");
Array<OneD,NekDouble> stdCoord(GetCoordim(),0.0);
GetLocCoords(gloCoord, stdCoord);
- if (stdCoord[0] >= -1 && stdCoord[1] >= -1
- && stdCoord[0] + stdCoord[1] <= 0)
+ if (stdCoord[0] >= -(1+tol) && stdCoord[1] >= -(1+tol)
+ && stdCoord[0] + stdCoord[1] <= tol)
{
return true;
}
diff --git a/library/SpatialDomains/TriGeom.h b/library/SpatialDomains/TriGeom.h
index cdcf9714031c80a6bb045a3b97e836fe2221ba37..7e6ec3b54bb1fd4ec944996c6e529ba6bbfc14e2 100644
--- a/library/SpatialDomains/TriGeom.h
+++ b/library/SpatialDomains/TriGeom.h
@@ -315,8 +315,7 @@ namespace Nektar
return kNedges;
}
- virtual bool v_ContainsPoint(
- const Array<OneD, const NekDouble> &gloCoord);
+ virtual bool v_ContainsPoint(const Array<OneD, const NekDouble> &gloCoord, NekDouble tol = 0.0);
};
}; //end of namespace SpatialDomains
}; //end of namespace Nektar
diff --git a/library/StdRegions/StdExpansion.cpp b/library/StdRegions/StdExpansion.cpp
index c4e5dc537134d41541b5428a7b753d89cbfe5123..2cdb7102fc674484d66102f2fbf43d3411a56132 100644
--- a/library/StdRegions/StdExpansion.cpp
+++ b/library/StdRegions/StdExpansion.cpp
@@ -1373,6 +1373,13 @@ namespace Nektar
}
+ NekDouble StdExpansion::v_PhysEvaluate(const Array<OneD, const NekDouble>& coords, const Array<OneD, const NekDouble>& physvals)
+ {
+ NEKERROR(ErrorUtil::efatal, "Method does not exist for this shape");
+ return 0;
+ }
+
+
void StdExpansion::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
{
NEKERROR(ErrorUtil::efatal, "This function has not "
diff --git a/library/StdRegions/StdExpansion.h b/library/StdRegions/StdExpansion.h
index e50678c72e3bf7c7b1f5e8ecd6e4494949efb248..6703699d79f9759327ae5caee1ed8a553c7e6fa2 100644
--- a/library/StdRegions/StdExpansion.h
+++ b/library/StdRegions/StdExpansion.h
@@ -1137,6 +1137,36 @@ namespace Nektar
return v_PhysEvaluate(coords);
}
+
+
+ /** \brief This function evaluates the expansion at a single
+ * (arbitrary) point of the domain
+ *
+ * This function is a wrapper around the virtual function
+ * \a v_PhysEvaluate()
+ *
+ * Based on the value of the expansion at the quadrature
+ * points provided in \a physvals, this function
+ * calculates the value of the expansion at an arbitrary
+ * single points (with coordinates \f$ \mathbf{x_c}\f$
+ * given by the pointer \a coords). This operation,
+ * equivalent to \f[ u(\mathbf{x_c}) = \sum_p
+ * \phi_p(\mathbf{x_c}) \hat{u}_p \f] is evaluated using
+ * Lagrangian interpolants through the quadrature points:
+ * \f[ u(\mathbf{x_c}) = \sum_p h_p(\mathbf{x_c}) u_p\f]
+ *
+ * \param coords the coordinates of the single point
+ * \param physvals the interpolated field at the quadrature points
+ *
+ * \return returns the value of the expansion at the
+ * single point
+ */
+ NekDouble PhysEvaluate(const Array<OneD, const NekDouble>& coords,
+ const Array<OneD, const NekDouble>& physvals)
+ {
+ return v_PhysEvaluate(coords,physvals);
+ }
+
const boost::shared_ptr<SpatialDomains::GeomFactors>& GetMetricInfo(void) const
{
return v_GetMetricInfo();
@@ -1531,6 +1561,8 @@ namespace Nektar
STD_REGIONS_EXPORT virtual NekDouble v_PhysEvaluate(const Array<OneD, const NekDouble>& coords);
+ STD_REGIONS_EXPORT virtual NekDouble v_PhysEvaluate(const Array<OneD, const NekDouble>& coords, const Array<OneD, const NekDouble> & physvals);
+
STD_REGIONS_EXPORT virtual void v_FillMode(const int mode, Array<OneD, NekDouble> &outarray);
STD_REGIONS_EXPORT virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey);
diff --git a/library/StdRegions/StdExpansion1D.cpp b/library/StdRegions/StdExpansion1D.cpp
index 94132020f2589c60777d648f51778b346bb02139..4aad708c9a6d605f7ede28a7fe8394152177b554 100644
--- a/library/StdRegions/StdExpansion1D.cpp
+++ b/library/StdRegions/StdExpansion1D.cpp
@@ -103,8 +103,13 @@ namespace Nektar
#endif //NEKTAR_USING_DIRECT_BLAS_CALLS
}
- NekDouble StdExpansion1D::PhysEvaluate(const Array<OneD, const NekDouble>& Lcoord)
- {
+ NekDouble StdExpansion1D::PhysEvaluate(const Array<OneD, const NekDouble>& Lcoord)
+ {
+ PhysEvaluate(Lcoord,m_phys);
+ }
+
+ NekDouble StdExpansion1D::PhysEvaluate(const Array<OneD, const NekDouble>& Lcoord, const Array<OneD, const NekDouble> & physvals)
+ {
int nquad = GetTotPoints();
NekDouble val;
DNekMatSharedPtr I = m_base[0]->GetI(Lcoord);
@@ -112,7 +117,7 @@ namespace Nektar
ASSERTL2(Lcoord[0] >= -1,"Lcoord[0] < -1");
ASSERTL2(Lcoord[0] <= 1,"Lcoord[0] > 1");
- val = Blas::Ddot(nquad, I->GetPtr(), 1, m_phys, 1);
+ val = Blas::Ddot(nquad, I->GetPtr(), 1, physvals, 1);
return val;
}
diff --git a/library/StdRegions/StdExpansion1D.h b/library/StdRegions/StdExpansion1D.h
index 8e86510d6616e428dd9819e7d8bfe48f1d749166..1c8c5fee3b47c118ac8034fd320eaf5904e93be6 100644
--- a/library/StdRegions/StdExpansion1D.h
+++ b/library/StdRegions/StdExpansion1D.h
@@ -77,6 +77,8 @@ namespace Nektar
*/
STD_REGIONS_EXPORT NekDouble PhysEvaluate(const Array<OneD, const NekDouble>& coords);
+ STD_REGIONS_EXPORT NekDouble PhysEvaluate(const Array<OneD, const NekDouble>& coords, const Array<OneD, const NekDouble> & physvals);
+
void SetCoeffsToOrientation(StdRegions::EdgeOrientation dir)
{
v_SetCoeffsToOrientation(dir);
@@ -120,6 +122,11 @@ namespace Nektar
return PhysEvaluate(coords);
}
+ virtual NekDouble v_PhysEvaluate(const Array<OneD, const NekDouble>& coords, const Array<OneD, const NekDouble> & physvals)
+ {
+ return PhysEvaluate(coords,physvals);
+ }
+
virtual void v_SetCoeffsToOrientation(StdRegions::EdgeOrientation dir,
Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
diff --git a/library/StdRegions/StdExpansion2D.cpp b/library/StdRegions/StdExpansion2D.cpp
index 9439861169483f5f036cfb6c7910cf7ff5ae9c31..b969d94351d35c37999cfb81d32b7bfbf40530ce 100644
--- a/library/StdRegions/StdExpansion2D.cpp
+++ b/library/StdRegions/StdExpansion2D.cpp
@@ -115,6 +115,12 @@ namespace Nektar
}
NekDouble StdExpansion2D::PhysEvaluate(const Array<OneD, const NekDouble>& coords)
+ {
+ PhysEvaluate(coords,m_phys);
+ }
+
+
+ NekDouble StdExpansion2D::PhysEvaluate(const Array<OneD, const NekDouble>& coords, const Array<OneD, const NekDouble> & physvals)
{
NekDouble val;
int i;
@@ -133,7 +139,7 @@ namespace Nektar
for (i = 0; i < nq1;++i)
{
wsp1[i] = Blas::Ddot(nq0, &(I->GetPtr())[0], 1,
- &m_phys[i * nq0], 1);
+ &physvals[i * nq0], 1);
}
// interpolate in second coordinate direction
diff --git a/library/StdRegions/StdExpansion2D.h b/library/StdRegions/StdExpansion2D.h
index 656b0ece0e82051edb951f9a36bcf662f59e104c..f3bd1f7d5be1242742cc7f13ea9c10a60c8cf6cb 100644
--- a/library/StdRegions/StdExpansion2D.h
+++ b/library/StdRegions/StdExpansion2D.h
@@ -118,6 +118,8 @@ namespace Nektar
*/
STD_REGIONS_EXPORT NekDouble PhysEvaluate(const Array<OneD, const NekDouble>& coords);
+ STD_REGIONS_EXPORT NekDouble PhysEvaluate(const Array<OneD, const NekDouble>& coords, const Array<OneD, const NekDouble> & physvals);
+
STD_REGIONS_EXPORT NekDouble Integral(const Array<OneD, const NekDouble>& inarray,
const Array<OneD, const NekDouble>& w0,
const Array<OneD, const NekDouble>& w1);
@@ -150,6 +152,11 @@ namespace Nektar
return PhysEvaluate(coords);
}
+ virtual NekDouble v_PhysEvaluate(const Array<OneD, const NekDouble>& coords, const Array<OneD, const NekDouble> & physvals)
+ {
+ return PhysEvaluate(coords,physvals);
+ }
+
virtual void v_SetUpPhysNormals(const int edge);
const NormalVector & v_GetEdgeNormal(const int edge) const;
diff --git a/library/StdRegions/StdExpansion3D.cpp b/library/StdRegions/StdExpansion3D.cpp
index d2280fe82034fdba4d632dac34d5bbba8bc7eebb..b0369052299505bd878ab8a9e0de2fe5d74a3a87 100644
--- a/library/StdRegions/StdExpansion3D.cpp
+++ b/library/StdRegions/StdExpansion3D.cpp
@@ -201,7 +201,12 @@ namespace Nektar
}
- NekDouble StdExpansion3D::v_PhysEvaluate(const Array<OneD, const NekDouble> &coords)
+ NekDouble StdExpansion3D::v_PhysEvaluate(const Array<OneD, const NekDouble> &coords )
+ {
+ v_PhysEvaluate(coords,m_phys);
+ }
+
+ NekDouble StdExpansion3D::v_PhysEvaluate(const Array<OneD, const NekDouble> &coords, const Array<OneD, const NekDouble> & physvals)
{
NekDouble value;
ASSERTL2(coords[0] >= -1,"coord[0] < -1");
@@ -227,7 +232,7 @@ namespace Nektar
interpolatingNodes = &I->GetPtr()[0];
for(int i = 0; i < Qy*Qz;++i)
{
- sumFactorization_qr[i] = Blas::Ddot(Qx, interpolatingNodes, 1, &m_phys[ i*Qx ], 1);
+ sumFactorization_qr[i] = Blas::Ddot(Qx, interpolatingNodes, 1, &physvals[ i*Qx ], 1);
}
// Interpolate in second coordinate direction
diff --git a/library/StdRegions/StdExpansion3D.h b/library/StdRegions/StdExpansion3D.h
index 1e588a5793563746004102360a1884fa99570268..752c0f116770c93b783c49e013643995cf3abd8d 100644
--- a/library/StdRegions/StdExpansion3D.h
+++ b/library/StdRegions/StdExpansion3D.h
@@ -130,10 +130,18 @@ namespace Nektar
return v_PhysEvaluate(coords);
}
+ NekDouble PhysEvaluate(const Array<OneD, const NekDouble>& coords,
+ const Array<OneD, const NekDouble>& physvals)
+ {
+ return v_PhysEvaluate(coords, physvals);
+ }
+
protected:
STD_REGIONS_EXPORT virtual NekDouble v_PhysEvaluate(const Array<OneD, const NekDouble>& coords);
+ STD_REGIONS_EXPORT virtual NekDouble v_PhysEvaluate(const Array<OneD, const NekDouble>& coords, const Array<OneD, const NekDouble> & physvals);
+
private:
virtual int v_GetShapeDimension() const
diff --git a/library/StdRegions/StdHexExp.cpp b/library/StdRegions/StdHexExp.cpp
index 451171e201cf5d5b18e356cbcad9a8fe654a55a5..df895002982e63382b9201fda488437c3b4dd04b 100644
--- a/library/StdRegions/StdHexExp.cpp
+++ b/library/StdRegions/StdHexExp.cpp
@@ -2540,6 +2540,12 @@ namespace Nektar
return StdExpansion3D::v_PhysEvaluate(Lcoords);
}
+ NekDouble StdHexExp::v_PhysEvaluate(
+ Array<OneD, const NekDouble>& Lcoords, const Array<OneD, const NekDouble> & physvals)
+ {
+ return StdExpansion3D::v_PhysEvaluate(Lcoords, physvals);
+ }
+
int StdHexExp::v_GetEdgeNcoeffs(const int i) const
{
ASSERTL2((i >= 0)&&(i <= 11),"edge id is out of range");
diff --git a/library/StdRegions/StdHexExp.h b/library/StdRegions/StdHexExp.h
index d973c64a680c92e90f8135ae896067cf96a7d9c4..31e6d38b9bbe8d2c027c4b78a3f2cc22dd9f8a0c 100644
--- a/library/StdRegions/StdHexExp.h
+++ b/library/StdRegions/StdHexExp.h
@@ -240,6 +240,9 @@ namespace Nektar
STD_REGIONS_EXPORT virtual NekDouble v_PhysEvaluate(Array<OneD,
const NekDouble>& Lcoords);
+
+ STD_REGIONS_EXPORT virtual NekDouble v_PhysEvaluate(Array<OneD, const NekDouble>& Lcoords, const Array<OneD, const NekDouble> & physvals);
+
STD_REGIONS_EXPORT virtual void v_MassMatrixOp(
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
diff --git a/library/StdRegions/StdPrismExp.cpp b/library/StdRegions/StdPrismExp.cpp
index 38136a6e3c8371dda45e87fd90324c4cd03f3232..b8a5c275c9e7eb4115fe2d0cf7a83ec17028c578 100644
--- a/library/StdRegions/StdPrismExp.cpp
+++ b/library/StdRegions/StdPrismExp.cpp
@@ -658,6 +658,11 @@ namespace Nektar
NekDouble StdPrismExp::PhysEvaluate(const Array<OneD, const NekDouble>& xi)
+ {
+ PhysEvaluate(xi,m_phys);
+ }
+
+ NekDouble StdPrismExp::PhysEvaluate(const Array<OneD, const NekDouble>& xi, const Array<OneD, const NekDouble> & physvals)
{
Array<OneD, NekDouble> eta = Array<OneD, NekDouble>(3);
@@ -682,7 +687,7 @@ namespace Nektar
}
- return StdExpansion3D::PhysEvaluate(eta);
+ return StdExpansion3D::PhysEvaluate(eta, physvals);
}
diff --git a/library/StdRegions/StdPrismExp.h b/library/StdRegions/StdPrismExp.h
index d8a25e352f1bd51e545b21137700b43cb7afda87..0f9a348f26d022aa5ab4832b83b91e69a5db5da9 100644
--- a/library/StdRegions/StdPrismExp.h
+++ b/library/StdRegions/StdPrismExp.h
@@ -205,7 +205,9 @@ namespace Nektar
/** \brief Single Point Evaluation */
STD_REGIONS_EXPORT NekDouble PhysEvaluate(const Array<OneD, const NekDouble>& xi);
-
+
+ STD_REGIONS_EXPORT NekDouble PhysEvaluate(const Array<OneD, const NekDouble>& xi, const Array<OneD, const NekDouble> & physvals);
+
STD_REGIONS_EXPORT void GetCoords( Array<OneD, NekDouble> & xi_x, Array<OneD, NekDouble> & xi_y, Array<OneD, NekDouble> & xi_z);
STD_REGIONS_EXPORT void FillMode(const int mode, Array<OneD, NekDouble> &outarray);
STD_REGIONS_EXPORT void WriteToFile(std::ofstream &outfile, OutputFormat format, const bool dumpVar = true, std::string var = "v");
@@ -404,6 +406,11 @@ namespace Nektar
return PhysEvaluate(Lcoords);
}
+ virtual NekDouble v_PhysEvaluate(const Array<OneD, const NekDouble>& Lcoords, const Array<OneD, const NekDouble> & physvals)
+ {
+ return PhysEvaluate(Lcoords, physvals);
+ }
+
virtual void v_WriteToFile(std::ofstream &outfile, OutputFormat format, const bool dumpVar = true, std::string var = "v")
{
WriteToFile(outfile,format,dumpVar,var);
diff --git a/library/StdRegions/StdPyrExp.cpp b/library/StdRegions/StdPyrExp.cpp
index 4f584754141e50fb6d40fbe4aee4b1a2913a783e..970634a837de0f5e499241538db6dad4b308d844 100644
--- a/library/StdRegions/StdPyrExp.cpp
+++ b/library/StdRegions/StdPyrExp.cpp
@@ -570,6 +570,11 @@ namespace Nektar
}
NekDouble StdPyrExp::PhysEvaluate(const Array<OneD, const NekDouble>& xi)
+ {
+ PhysEvaluate(xi,m_phys);
+ }
+
+ NekDouble StdPyrExp::PhysEvaluate(const Array<OneD, const NekDouble>& xi, const Array<OneD, const NekDouble> & physvals)
{
Array<OneD, NekDouble> eta = Array<OneD, NekDouble>(3);
@@ -592,7 +597,7 @@ namespace Nektar
// cout << "xi_x = " << xi[0] << ", xi_y = " << xi[1] << ", xi_z = " << xi[2] << endl;
// cout << "eta_x = " << eta[0] << ", eta_y = " << eta[1] << ", eta_z = " << eta[2] << endl;
- return StdExpansion3D::PhysEvaluate(eta);
+ return StdExpansion3D::PhysEvaluate(eta, physvals);
}
void StdPyrExp::GetCoords( Array<OneD, NekDouble> & xi_x, Array<OneD, NekDouble> & xi_y, Array<OneD, NekDouble> & xi_z)
diff --git a/library/StdRegions/StdPyrExp.h b/library/StdRegions/StdPyrExp.h
index 89a73fb1aeececf06ccd405917fd87e7d78da0c8..a5fcc0171019404bf96bf6910e479b5ee703c103 100644
--- a/library/StdRegions/StdPyrExp.h
+++ b/library/StdRegions/StdPyrExp.h
@@ -215,6 +215,7 @@ namespace Nektar
/** \brief Single Point Evaluation */
STD_REGIONS_EXPORT NekDouble PhysEvaluate(const Array<OneD, const NekDouble>& xi);
+ STD_REGIONS_EXPORT NekDouble PhysEvaluate(const Array<OneD, const NekDouble>& xi, const Array<OneD, const NekDouble> & physvals);
STD_REGIONS_EXPORT void GetCoords( Array<OneD, NekDouble> & xi_x, Array<OneD, NekDouble> & xi_y, Array<OneD, NekDouble> & xi_z);
STD_REGIONS_EXPORT void WriteToFile(std::ofstream &outfile, OutputFormat format, const bool dumpVar = true, std::string var = "v");
@@ -392,6 +393,11 @@ namespace Nektar
{
return PhysEvaluate(Lcoords);
}
+
+ virtual NekDouble v_PhysEvaluate(const Array<OneD, const NekDouble>& Lcoords, const Array<OneD, const NekDouble> & physvals)
+ {
+ return PhysEvaluate(Lcoords, physvals);
+ }
virtual int v_GetEdgeNcoeffs(const int i) const
{
diff --git a/library/StdRegions/StdQuadExp.h b/library/StdRegions/StdQuadExp.h
index e58dd6b1588d721b6961eba1eb4501726692c998..f84aa7aa3e36db28ca0f94d7e73682c4041f7d89 100644
--- a/library/StdRegions/StdQuadExp.h
+++ b/library/StdRegions/StdQuadExp.h
@@ -213,7 +213,13 @@ namespace Nektar
NekDouble PhysEvaluate(Array<OneD, const NekDouble>& coords)
{
- return StdExpansion2D::PhysEvaluate(coords);
+ return StdExpansion2D::PhysEvaluate(coords, m_phys);
+ }
+
+ NekDouble PhysEvaluate(Array<OneD, const NekDouble>& coords,
+ const Array<OneD, const NekDouble> & physvals)
+ {
+ return StdExpansion2D::PhysEvaluate(coords, physvals);
}
STD_REGIONS_EXPORT void GetBoundaryMap(Array<OneD, unsigned int>& outarray);
@@ -629,6 +635,12 @@ namespace Nektar
return PhysEvaluate(Lcoords);
}
+
+ virtual NekDouble v_PhysEvaluate(Array<OneD, const NekDouble>& Lcoords, const Array<OneD, const NekDouble> &physvals)
+ {
+ return PhysEvaluate(Lcoords, physvals);
+ }
+
virtual void v_GetBoundaryMap(Array<OneD, unsigned int>& outarray)
{
GetBoundaryMap(outarray);
diff --git a/library/StdRegions/StdTetExp.cpp b/library/StdRegions/StdTetExp.cpp
index b7c4ccd8b133f1c6631e76b2e25fc1e26a413c7f..bdb114dcc7e5a108d3192ad890feb7c8c3112a47 100644
--- a/library/StdRegions/StdTetExp.cpp
+++ b/library/StdRegions/StdTetExp.cpp
@@ -1079,6 +1079,11 @@ namespace Nektar
///////////////////////////////
NekDouble StdTetExp::v_PhysEvaluate(const Array<OneD, const NekDouble>& xi)
+ {
+ v_PhysEvaluate(xi,m_phys);
+ }
+
+ NekDouble StdTetExp::v_PhysEvaluate(const Array<OneD, const NekDouble>& xi, const Array<OneD, const NekDouble> & physvals)
{
// Validation checks
ASSERTL0(xi[0] + xi[1] + xi[2] <= -1,
@@ -1124,7 +1129,7 @@ namespace Nektar
|| (eta[2] - NekConstants::kNekZeroTol <= 1),
"Eta Coordinate outside bounds of tetrahedron.");
- return StdExpansion3D::v_PhysEvaluate(eta);
+ return StdExpansion3D::v_PhysEvaluate(eta, physvals);
}
diff --git a/library/StdRegions/StdTetExp.h b/library/StdRegions/StdTetExp.h
index 9d3f55e069f50dbc96283f457aac1752b3e110e0..7622ca7a963a2e036bcc2dffe5a0e39b0c461309 100644
--- a/library/StdRegions/StdTetExp.h
+++ b/library/StdRegions/StdTetExp.h
@@ -143,6 +143,8 @@ namespace Nektar
/** \brief Single Point Evaluation */
STD_REGIONS_EXPORT NekDouble PhysEvaluate3D(const Array<OneD, const NekDouble>& coords);
+ STD_REGIONS_EXPORT NekDouble PhysEvaluate3D(const Array<OneD, const NekDouble>& coords, const Array<OneD, const NekDouble> & physvals);
+
//----------------------------------
// Generate Matrix Routine
@@ -236,6 +238,8 @@ namespace Nektar
STD_REGIONS_EXPORT virtual NekDouble v_PhysEvaluate(const Array<OneD, const NekDouble>& coords);
+ STD_REGIONS_EXPORT virtual NekDouble v_PhysEvaluate(const Array<OneD, const NekDouble>& coords, const Array<OneD, const NekDouble> & physvals);
+
private:
STD_REGIONS_EXPORT virtual void v_BwdTrans_SumFac(const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble> &outarray);
diff --git a/library/StdRegions/StdTriExp.h b/library/StdRegions/StdTriExp.h
index 03efbc7e533c42c50ea000b1f2ae004c05d5f51f..285a70732b17191ec10960a1b95dc644298ab724 100644
--- a/library/StdRegions/StdTriExp.h
+++ b/library/StdRegions/StdTriExp.h
@@ -232,6 +232,8 @@ namespace Nektar
/** \brief Single Point Evaluation */
STD_REGIONS_EXPORT NekDouble PhysEvaluate(const Array<OneD, const NekDouble>& coords);
+
+ STD_REGIONS_EXPORT NekDouble PhysEvaluate(const Array<OneD, const NekDouble>& coords,const Array<OneD, const NekDouble> & physvals);
STD_REGIONS_EXPORT void GetBoundaryMap(Array<OneD, unsigned int>& outarray);
@@ -641,6 +643,11 @@ namespace Nektar
return PhysEvaluate(coords);
}
+ virtual NekDouble v_PhysEvaluate(const Array<OneD, const NekDouble>& coords, const Array<OneD, const NekDouble> &physvals)
+ {
+ return PhysEvaluate(coords, physvals);
+ }
+
virtual void v_GetBoundaryMap(Array<OneD, unsigned int>& outarray)
{
GetBoundaryMap(outarray);