Commit 5b3697d4 authored by Dave Moxey's avatar Dave Moxey
Browse files

Removed unneeded metric information from GeomFactors classes.

Fixed an issue with SetUpTangents function.
Removed a few unnecessary typecasts.
parent e4311aa9
......@@ -195,9 +195,8 @@ namespace Nektar
const int mode_offset,
NekDouble *coeffs);
LOCAL_REGIONS_EXPORT virtual void v_SetUpPhysTangents(
const ExpansionSharedPtr &exp2D,
const boost::shared_ptr<Expansion> &exp2D,
const int edge);
LOCAL_REGIONS_EXPORT virtual const
......
......@@ -203,15 +203,7 @@ namespace Nektar
// Loop over all edges of element i
for(j = 0; j < locExpansion->GetNedges(); ++j)
{
if (nDim == 2)
{
meshEdgeId = LocalRegions::Expansion2D::FromStdExp(locExpansion)->GetGeom2D()->GetEid(j);
}
else
{
meshEdgeId = LocalRegions::Expansion3D::FromStdExp(locExpansion)->GetGeom3D()->GetEid(j);
}
meshEdgeId = locExpansion->GetGeom()->GetEid(j);
pIt = perEdges.find(meshEdgeId);
if (pIt != perEdges.end())
{
......@@ -239,11 +231,13 @@ namespace Nektar
// Loop over all faces of element i
for(j = 0; j < locExpansion->GetNfaces(); ++j)
{
faceOrient = LocalRegions::Expansion3D::FromStdExp(locExpansion)->GetGeom3D()->GetFaceOrient(j);
faceOrient = boost::dynamic_pointer_cast<
LocalRegions::Expansion3D>(
locExpansion)->GetGeom3D()->GetFaceOrient(j);
locExpansion->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
dof = locExpansion->GetFaceIntNcoeffs(j);
meshFaceId = LocalRegions::Expansion3D::FromStdExp(locExpansion)->GetGeom3D()->GetFid(j);
meshFaceId = locExpansion->GetGeom()->GetFid(j);
pIt = perFaces.find(meshFaceId);
if (pIt != perFaces.end())
......
......@@ -445,17 +445,15 @@ namespace Nektar
}
// Set up information for periodic boundary conditions.
LocalRegions::Expansion2DSharedPtr exp2d;
boost::unordered_map<int,pair<int,int> > perEdgeToExpMap;
boost::unordered_map<int,pair<int,int> >::iterator it2;
for (cnt = n = 0; n < m_exp->size(); ++n)
{
exp2d = LocalRegions::Expansion2D::FromStdExp((*m_exp)[n]);
for (e = 0; e < exp2d->GetNedges(); ++e, ++cnt)
for (e = 0; e < (*m_exp)[n]->GetNedges(); ++e, ++cnt)
{
PeriodicMap::iterator it = m_periodicEdges.find(
exp2d->GetGeom2D()->GetEid(e));
(*m_exp)[n]->GetGeom()->GetEid(e));
if (it != m_periodicEdges.end())
{
perEdgeToExpMap[it->first] = make_pair(n, e);
......@@ -468,8 +466,7 @@ namespace Nektar
cnt = 0;
for (int i = 0; i < m_exp->size(); ++i)
{
exp2d = LocalRegions::Expansion2D::FromStdExp((*m_exp)[i]);
for (int j = 0; j < exp2d->GetNedges(); ++j, ++cnt)
for (int j = 0; j < (*m_exp)[i]->GetNedges(); ++j, ++cnt)
{
m_leftAdjacentEdges[cnt] = IsLeftAdjacentEdge(i, j);
}
......@@ -479,10 +476,9 @@ namespace Nektar
cnt = 0;
for (int n = 0; n < m_exp->size(); ++n)
{
exp2d = LocalRegions::Expansion2D::FromStdExp((*m_exp)[n]);
for (int e = 0; e < exp2d->GetNedges(); ++e, ++cnt)
for (int e = 0; e < (*m_exp)[n]->GetNedges(); ++e, ++cnt)
{
int edgeGeomId = exp2d->GetGeom2D()->GetEid(e);
int edgeGeomId = (*m_exp)[n]->GetGeom()->GetEid(e);
int offset = m_trace->GetPhys_Offset(
elmtToTrace[n][e]->GetElmtId());
......
......@@ -2794,7 +2794,7 @@ namespace Nektar
}
void ExpList::v_SetUpPhysTangents(
const StdRegions::StdExpansionVector &locexp)
const LocalRegions::ExpansionVector &locexp)
{
ASSERTL0(false,
"This method is not defined or valid for this class type");
......
......@@ -705,7 +705,7 @@ namespace Nektar
inline void SetUpPhysNormals();
inline void SetUpPhysTangents(const StdRegions::StdExpansionVector &locexp);
inline void SetUpPhysTangents(const LocalRegions::ExpansionVector &locexp);
inline void SetUpTangents();
......@@ -1097,7 +1097,7 @@ namespace Nektar
virtual void v_IProductWRTBase_IterPerExp(const Array<OneD,const NekDouble> &inarray, Array<OneD, NekDouble> &outarray);
virtual void v_SetUpPhysTangents(const StdRegions::StdExpansionVector &locexp);
virtual void v_SetUpPhysTangents(const LocalRegions::ExpansionVector &locexp);
virtual void v_GeneralMatrixOp(
const GlobalMatrixKey &gkey,
......@@ -2014,7 +2014,7 @@ namespace Nektar
}
inline void ExpList::SetUpPhysTangents(
const StdRegions::StdExpansionVector &locexp)
const LocalRegions::ExpansionVector &locexp)
{
v_SetUpPhysTangents(locexp);
}
......
......@@ -799,7 +799,7 @@ namespace Nektar
//setup map of all global ids along booundary
for(cnt = i=0; i< (*m_exp).size(); ++i)
{
exp1D = LocalRegions::Expansion1D::FromStdExp((*m_exp)[i]);
exp1D = LocalRegions::Expansion1D::FromStdExp((*m_exp)[i]);
id = exp1D->GetGeom1D()->GetEid();
EdgeGID[id] = cnt++;
}
......@@ -809,13 +809,13 @@ namespace Nektar
{
for(i=0; i < locexp[n]->GetNedges(); ++i)
{
exp2D = LocalRegions::Expansion2D::FromStdExp(locexp[n]);
id = exp2D->GetGeom2D()->GetEid(i);
if(EdgeGID.count(id)> 0)
{
(*m_exp)[EdgeGID.find(id)->second]
->SetUpPhysTangents(locexp[n],i);
}
exp2D = LocalRegions::Expansion2D::FromStdExp(locexp[n]);
id = exp2D->GetGeom2D()->GetEid(i);
if(EdgeGID.count(id)> 0)
{
(*m_exp)[EdgeGID.find(id)->second]
->SetUpPhysTangents(locexp[n],i);
}
}
}
}
......
......@@ -715,7 +715,7 @@ namespace Nektar
map<int, int> ElmtID_to_ExpID;
for(i = 0; i < m_planes[0]->GetExpSize(); ++i)
{
ElmtID_to_ExpID[LocalRegions::Expansion2D::FromStdExp((*m_exp)[i])->GetGeom()->GetGlobalID()] = i;
ElmtID_to_ExpID[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
}
......
......@@ -52,15 +52,11 @@ namespace Nektar
*/
GeomFactors::GeomFactors(const GeomType gtype,
const int expdim,
const int coordim,
const bool UseQuadMet,
const bool UseLaplMet):
const int coordim) :
m_type(gtype),
m_expDim(expdim),
m_coordDim(coordim),
m_coords(Array<OneD, StdRegions::StdExpansionSharedPtr>(m_coordDim)),
m_isUsingQuadMetrics(UseQuadMet),
m_isUsingLaplMetrics(UseLaplMet),
m_pointsKey(expdim)
{
}
......@@ -73,8 +69,6 @@ namespace Nektar
m_expDim(S.m_expDim),
m_coordDim(S.m_coordDim),
m_coords(S.m_coords),
m_isUsingQuadMetrics(S.m_isUsingQuadMetrics),
m_isUsingLaplMetrics(S.m_isUsingLaplMetrics),
m_pointsKey(S.m_pointsKey)
{
}
......@@ -200,17 +194,6 @@ namespace Nektar
ASSERTL0(false, "Cannot compute tangents for this geometry.");
}
/**
* Placeholder function.
*/
void GeomFactors::v_SetUpQuadratureMetrics(
LibUtilities::ShapeType shape,
const Array<OneD, const LibUtilities::BasisSharedPtr>
&tbasis)
{
ASSERTL0(false, "Quadrature Metrics not implemented.");
}
/**
* @brief Establishes if two GeomFactors objects are equal.
*/
......@@ -231,16 +214,6 @@ namespace Nektar
return false;
}
if(!(lhs.m_isUsingQuadMetrics == rhs.m_isUsingQuadMetrics))
{
return false;
}
if(!(lhs.m_isUsingLaplMetrics == rhs.m_isUsingLaplMetrics))
{
return false;
}
for(int i = 0; i < lhs.m_expDim; i++)
{
if(!(lhs.m_pointsKey[i] == rhs.m_pointsKey[i]))
......
......@@ -41,10 +41,9 @@
#include <LibUtilities/Foundations/Basis.h>
#include <SpatialDomains/SpatialDomains.hpp>
#include <SpatialDomains/SpatialDomainsDeclspec.h>
#include <StdRegions/StdExpansion.h> // for StdExpansionSharedPtr
#include <StdRegions/StdExpansion.h>
#include <StdRegions/StdRegions.hpp>
namespace Nektar
{
namespace LibUtilities
......@@ -125,25 +124,12 @@ namespace Nektar
/// Return the number of dimensions of the coordinate system.
inline int GetCoordim() const;
/// Flag indicating if quadrature metrics are in use.
inline bool IsUsingQuadMetrics() const;
/// Flag indicating if Tangents have been computed.
inline bool IsUsingTangents() const;
/// Set up quadrature metrics
inline void SetUpQuadratureMetrics(
LibUtilities::ShapeType shape,
const Array<OneD, const LibUtilities::BasisSharedPtr>
&tbasis);
/// Set up Tangents
inline void SetUpTangents();
/// Retrieve the quadrature metrics.
inline const Array<OneD, const NekDouble> &GetQuadratureMetrics()
const;
/// Computes the edge tangents from 1D element
inline void ComputeEdgeTangents(
const GeometrySharedPtr &geom2D,
......@@ -204,8 +190,6 @@ namespace Nektar
boost::hash_combine(hash, (int)m_type);
boost::hash_combine(hash, m_expDim);
boost::hash_combine(hash, m_coordDim);
boost::hash_combine(hash, m_isUsingQuadMetrics);
boost::hash_combine(hash, m_isUsingLaplMetrics);
boost::hash_range(hash, m_jac.begin(), m_jac.end());
for (int i = 0; i < m_gmat.GetRows(); ++i)
{
......@@ -223,10 +207,6 @@ namespace Nektar
int m_coordDim;
/// Stores information about the expansion.
Array<OneD, StdRegions::StdExpansionSharedPtr> m_coords;
/// Use Quadrature metrics
bool m_isUsingQuadMetrics;
/// Use Laplacian metrics
bool m_isUsingLaplMetrics;
/// Principle tangent direction.
enum GeomTangents m_tangentDir;
/// Principle tangent circular dir coords
......@@ -239,9 +219,6 @@ namespace Nektar
/// at each of those points.
Array<OneD,NekDouble> m_jac;
///
Array<OneD,NekDouble> m_weightedjac;
/// Array of size coordim x nquad which holds the inverse of the
/// derivative of the local map in each direction at each
/// quadrature point.
......@@ -279,9 +256,7 @@ namespace Nektar
/// instantiated externally.
GeomFactors(const GeomType gtype,
const int expdim,
const int coordim,
const bool UseQuadMet,
const bool UseLaplMet);
const int coordim);
/// Copy constructor.
GeomFactors(const GeomFactors &S);
......@@ -298,12 +273,6 @@ namespace Nektar
/// Set up the tangent vectors
virtual void v_ComputeTangents();
/// Set up quadrature metrics
virtual void v_SetUpQuadratureMetrics(
LibUtilities::ShapeType shape,
const Array<OneD, const LibUtilities::BasisSharedPtr>
&tbasis);
};
/// A hash functor for geometric factors. Utilises
......@@ -354,27 +323,12 @@ namespace Nektar
return m_coordDim;
}
/// Flag indicating if quadrature metrics are in use.
inline bool GeomFactors::IsUsingQuadMetrics() const
{
return m_isUsingQuadMetrics;
}
/// Flag indicating if Tangents are in use.
inline bool GeomFactors::IsUsingTangents() const
{
return (m_tangents.num_elements() != 0);
}
/// Set up quadrature metrics
inline void GeomFactors::SetUpQuadratureMetrics(
LibUtilities::ShapeType shape,
const Array<OneD, const LibUtilities::BasisSharedPtr>
&tbasis)
{
v_SetUpQuadratureMetrics(shape, tbasis);
}
/// Set up Tangents
inline void GeomFactors::SetUpTangents()
{
......@@ -382,16 +336,6 @@ namespace Nektar
v_ComputeTangents();
}
/// Retrieve the quadrature metrics.
inline const Array<OneD, const NekDouble>
&GeomFactors::GetQuadratureMetrics() const
{
ASSERTL0(m_isUsingQuadMetrics,
"This metric has not been set up for this type of "
"expansion");
return m_weightedjac;
}
/// Computes the edge tangents from a 1D element
inline void GeomFactors::ComputeEdgeTangents(
const GeometrySharedPtr &geom2D,
......
......@@ -63,37 +63,27 @@ namespace Nektar
* @param coordim Dimension of coordinate system.
* @param Coords ?
* @param tbasis Basis for derivatives
* @param SetUpQuadratureMetrics ?
* @param SetUpLaplacianMetrics ?
*/
GeomFactors1D::GeomFactors1D(const GeomType gtype,
const int coordim,
const Array<OneD, const StdRegions
::StdExpansion1DSharedPtr> &Coords,
const Array<OneD, const LibUtilities::BasisSharedPtr>
&tbasis,
const bool QuadMetrics,
const bool LaplMetrics):
GeomFactors(gtype,1,coordim,QuadMetrics,LaplMetrics)
&tbasis) :
GeomFactors(gtype, 1, coordim)
{
// Perform sanity checks
ASSERTL1((coordim == 1)||(coordim == 2)||(coordim == 3),
ASSERTL1(coordim == 1 || coordim == 2 || coordim == 3,
"The coordinate dimension should be equal to one, two or "
"three for one-dimensional elements");
ASSERTL1(tbasis.num_elements()==1,
ASSERTL1(tbasis.num_elements() == 1,
"tbasis should be an array of size one");
ASSERTL1(LaplMetrics?QuadMetrics:true,
"SetUpQuadratureMetrics should be true if "
"SetUpLaplacianMetrics is true");
for (int i = 0; i < m_coordDim; ++i)
{
m_coords[i] = Coords[i];
}
// Get the shape of the expansion
LibUtilities::ShapeType shape = Coords[0]->DetShapeType();
// The quadrature points of the mapping
// (as specified in Coords)
LibUtilities::PointsKey pkey_map(
......@@ -147,12 +137,6 @@ namespace Nektar
// 1. The (determinant of the) jacobian and the differentation
// metrics
SetUpJacGmat1D();
// 2. the jacobian muliplied with the quadrature weights
if(QuadMetrics)
{
SetUpQuadratureMetrics(shape,tbasis);
}
}
......
......@@ -65,9 +65,7 @@ namespace Nektar
const Array<OneD, const StdRegions
::StdExpansion1DSharedPtr> &Coords,
const Array<OneD, const LibUtilities::BasisSharedPtr>
&tbasis,
const bool QuadMetrics = false,
const bool LaplMetrics = false);
&tbasis);
/// Copy constructor
SPATIAL_DOMAINS_EXPORT GeomFactors1D(const GeomFactors1D& S);
......@@ -84,11 +82,6 @@ namespace Nektar
const GeometrySharedPtr &geom,
const int edge,
const LibUtilities::PointsKey &to_key);
/**
* @todo Implement 1D quadrature metrics.
* @todo Implement 1D laplacian metrics.
*/
};
}
}
......
......@@ -117,8 +117,6 @@ namespace Nektar
* @param coordim Dimension of coordinate system.
* @param Coords Coordinate information.
* @param tbasis Tangent basis.
* @param SetUpQuadratureMetrics Compute quadrature metrics.
* @param SetUpLaplacianMetrics Compute Laplacian metrics.
*/
GeomFactors2D::GeomFactors2D(
const GeomType gtype,
......@@ -127,20 +125,15 @@ namespace Nektar
::StdExpansion2DSharedPtr> &Coords,
const Array<OneD, const LibUtilities::BasisSharedPtr>
&tbasis,
const bool QuadMetrics,
const bool LaplMetrics,
const bool CheckJacPositive ) :
GeomFactors(gtype,2,coordim,QuadMetrics,LaplMetrics)
GeomFactors(gtype, 2, coordim)
{
// Sanity checks.
ASSERTL1((coordim == 2)||(coordim == 3),
ASSERTL1(coordim == 2 || coordim == 3,
"The coordinate dimension should be equal to two or three"
"for two-dimensional elements");
ASSERTL1(tbasis.num_elements()==2,
ASSERTL1(tbasis.num_elements() == 2,
"tbasis should be an array of size two");
ASSERTL1(LaplMetrics?QuadMetrics:true,
"SetUpQuadratureMetrics should be true if "
"SetUpLaplacianMetrics is true");
// Copy shared pointers.
for (int i = 0; i < m_coordDim; ++i)
......@@ -148,8 +141,6 @@ namespace Nektar
m_coords[i] = Coords[i];
}
LibUtilities::ShapeType shape = Coords[0]->DetShapeType();
// The quadrature points of the mapping
// (as specified in Coords)
LibUtilities::PointsKey pkey0_map(
......
......@@ -49,16 +49,15 @@ namespace Nektar
public:
/// One dimensional geometric factors based on one-, two- or three-
/// dimensional coordinate description.
SPATIAL_DOMAINS_EXPORT GeomFactors2D(const GeomType gtype,
const int coordim,
const Array<OneD, const StdRegions
::StdExpansion2DSharedPtr> &Coords,
const Array<OneD, const LibUtilities::BasisSharedPtr>
&tbasis,
const bool QuadMetrics = false,
const bool LaplMetrics = false,
const bool CheckJacPositive = true);
SPATIAL_DOMAINS_EXPORT GeomFactors2D(
const GeomType gtype,
const int coordim,
const Array<OneD, const StdRegions
::StdExpansion2DSharedPtr> &Coords,
const Array<OneD, const LibUtilities::BasisSharedPtr>
&tbasis,
const bool CheckJacPositive = true);
/// Copy constructor.
SPATIAL_DOMAINS_EXPORT GeomFactors2D(const GeomFactors2D& S);
......
......@@ -52,33 +52,26 @@ namespace Nektar
* @param coordim Dimension of coordinate system.
* @param Coords ?
* @param tbasis Basis for tangential vectors.
* @param SetUpQuadratureMetrics ?
* @param SetUpLaplacianMetrics ?
*/
GeomFactors3D::GeomFactors3D(const GeomType gtype,
const int coordim,
const Array<OneD, const StdRegions
::StdExpansion3DSharedPtr> &Coords,
const Array<OneD, const LibUtilities::BasisSharedPtr>
&tbasis,
const bool QuadMetrics,
const bool LaplMetrics) :
GeomFactors(gtype,3,coordim,QuadMetrics,LaplMetrics)
&tbasis) :
GeomFactors(gtype, 3, coordim)
{
ASSERTL1((coordim == 3),
ASSERTL1(coordim == 3,
"The coordinate dimension should be to three"
"for three-dimensional elements");
ASSERTL1(tbasis.num_elements()==3,"tbasis should be an array of size three");
ASSERTL1(LaplMetrics?QuadMetrics:true,
"SetUpQuadratureMetrics should be true if SetUpLaplacianMetrics is true");
ASSERTL1(tbasis.num_elements() == 3,
"tbasis should be an array of size three");
for (int i = 0; i < m_coordDim; ++i)
{
m_coords[i] = Coords[i];
}
LibUtilities::ShapeType shape = Coords[0]->DetShapeType();
// The quadrature points of the mapping
// (as specified in Coords)
LibUtilities::PointsKey pkey0_map(Coords[0]->GetBasis(0)->GetPointsKey());
......@@ -257,208 +250,5 @@ namespace Nektar
}
}
}
/**
* @brief Set up the m_weightedjac array, which holds the Jacobian at
* each quadrature point multipled by the quadrature weight.
*/
/* void GeomFactors3D::v_SetUpQuadratureMetrics(
StdRegions::ExpansionType shape,
const Array<OneD, const LibUtilities::BasisSharedPtr> &tbasis)
{
ASSERTL1(tbasis.num_elements() == m_expDim,
"Inappropriate dimension of tbasis");
int i,j;
int nquad0 = m_pointsKey[0].GetNumPoints();
int nquad1 = m_pointsKey[1].GetNumPoints();
int nquad2 = m_pointsKey[2].GetNumPoints();
int nqtot = nquad0*nquad1*nquad2;
m_weightedjac = Array<OneD, NekDouble>(nqtot);
m_isUsingQuadMetrics = true;
// Fill the array m_weighted jac with the values
// of the (already computed) jacobian (=m_jac)
if (m_type == eRegular || m_type == eMovingRegular)
{
Vmath::Fill(nqtot,m_jac[0],m_weightedjac.get(),1);
}
else
{
Vmath::Vcopy(nqtot,m_jac.get(),1,m_weightedjac.get(),1);
}
// Get hold of the quadrature weights
const Array<OneD, const NekDouble>& w0 = tbasis[0]->GetW();
const Array<OneD, const NekDouble>& w1 = tbasis[1]->GetW();
const Array<OneD, const NekDouble>& w2 = tbasis[2]->GetW();
if (w0.num_elements() == 0 ||
w1.num_elements() == 0 ||
w2.num_elements() == 0)
{
m_isUsingQuadMetrics = false;
m_weightedjac = Array<OneD, NekDouble>();
return;
}
// Multiply the jacobian with the quadrature weights
switch(shape)
{
case LibUtilities::eHexahedron:
{