Commit 52599e7d authored by Spencer Sherwin's avatar Spencer Sherwin

Updated DetFaceBasisKey and SpatialDomains::GetFaceBasisKey to use a unique...

Updated DetFaceBasisKey and SpatialDomains::GetFaceBasisKey to use a unique definition in StdExpansion3D.cpp calling EvaluateTriFaceBasisKey and EvaluateQuadFaceBasisKey. Also updated ComputeFaceNormals to use same definition for interpolation and hopefully be more robust to variable quadrature orders
parent 499785df
This diff is collapsed.
......@@ -869,8 +869,19 @@ namespace Nektar
const Array<TwoD, const NekDouble> &df = geomFactors->GetDerivFactors(ptsKeys);
const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
int nq0 = ptsKeys[0].GetNumPoints();
int nq1 = ptsKeys[1].GetNumPoints();
int nq2 = ptsKeys[2].GetNumPoints();
int nq01 = nq0*nq1;
int nqtot;
LibUtilities::BasisKey tobasis0 = DetFaceBasisKey(face,0);
LibUtilities::BasisKey tobasis1 = DetFaceBasisKey(face,1);
// Number of quadrature points in face expansion.
int nq = m_base[0]->GetNumPoints()*m_base[0]->GetNumPoints();
int nq_face = tobasis0.GetNumPoints()*tobasis1.GetNumPoints();
int vCoordDim = GetCoordim();
int i;
......@@ -878,7 +889,7 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> > &normal = m_faceNormals[face];
for (i = 0; i < vCoordDim; ++i)
{
normal[i] = Array<OneD, NekDouble>(nq);
normal[i] = Array<OneD, NekDouble>(nq_face);
}
// Regular geometry case
......@@ -893,7 +904,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-df[3*i+2][0],normal[i],1);
normal[i][0] = -df[3*i+2][0];;
}
break;
}
......@@ -901,7 +912,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-df[3*i+1][0],normal[i],1);
normal[i][0] = -df[3*i+1][0];
}
break;
}
......@@ -909,7 +920,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,df[3*i][0]+df[3*i+2][0],normal[i],1);
normal[i][0] = df[3*i][0]+df[3*i+2][0];
}
break;
}
......@@ -917,7 +928,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,df[3*i+1][0],normal[i],1);
normal[i][0] = df[3*i+1][0];
}
break;
}
......@@ -925,7 +936,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-df[3*i][0],normal[i],1);
normal[i][0] = -df[3*i][0];
}
break;
}
......@@ -942,21 +953,15 @@ namespace Nektar
fac = 1.0/sqrt(fac);
for (i = 0; i < vCoordDim; ++i)
{
Vmath::Smul(nq,fac,normal[i],1,normal[i],1);
Vmath::Fill(nq_face,fac*normal[i][0],normal[i],1);
}
}
else
{
// Set up deformed normals.
int j, k;
int nq0 = ptsKeys[0].GetNumPoints();
int nq1 = ptsKeys[1].GetNumPoints();
int nq2 = ptsKeys[2].GetNumPoints();
int nq01 = nq0*nq1;
int nqtot;
// Determine number of quadrature points on the face.
// Determine number of quadrature points on the face of 3D elmt
if (face == 0)
{
nqtot = nq0*nq1;
......@@ -973,7 +978,6 @@ namespace Nektar
LibUtilities::PointsKey points0;
LibUtilities::PointsKey points1;
Array<OneD, NekDouble> work (nq, 0.0);
Array<OneD, NekDouble> normals(vCoordDim*nqtot,0.0);
// Extract Jacobian along face and recover local derivatives
......@@ -1083,37 +1087,39 @@ namespace Nektar
ASSERTL0(false,"face is out of range (face < 4)");
}
Array<OneD, NekDouble> work (nq_face, 0.0);
// Interpolate Jacobian and invert
LibUtilities::Interp2D(points0, points1, jac,
m_base[0]->GetPointsKey(),
m_base[0]->GetPointsKey(),
tobasis0.GetPointsKey(),
tobasis1.GetPointsKey(),
work);
Vmath::Sdiv(nq, 1.0, &work[0], 1, &work[0], 1);
Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
// Interpolate normal and multiply by inverse Jacobian.
for(i = 0; i < vCoordDim; ++i)
{
LibUtilities::Interp2D(points0, points1,
&normals[i*nqtot],
m_base[0]->GetPointsKey(),
m_base[0]->GetPointsKey(),
tobasis0.GetPointsKey(),
tobasis1.GetPointsKey(),
&normal[i][0]);
Vmath::Vmul(nq,work,1,normal[i],1,normal[i],1);
Vmath::Vmul(nq_face,work,1,normal[i],1,normal[i],1);
}
// Normalise to obtain unit normals.
Vmath::Zero(nq,work,1);
Vmath::Zero(nq_face,work,1);
for(i = 0; i < GetCoordim(); ++i)
{
Vmath::Vvtvp(nq,normal[i],1,normal[i],1,work,1,work,1);
Vmath::Vvtvp(nq_face,normal[i],1,normal[i],1,work,1,work,1);
}
Vmath::Vsqrt(nq,work,1,work,1);
Vmath::Sdiv (nq,1.0,work,1,work,1);
Vmath::Vsqrt(nq_face,work,1,work,1);
Vmath::Sdiv (nq_face,1.0,work,1,work,1);
for(i = 0; i < GetCoordim(); ++i)
{
Vmath::Vmul(nq,normal[i],1,work,1,normal[i],1);
Vmath::Vmul(nq_face,normal[i],1,work,1,normal[i],1);
}
}
}
......
......@@ -515,8 +515,12 @@ namespace Nektar
const Array<TwoD, const NekDouble> &df = geomFactors->GetDerivFactors(ptsKeys);
const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
LibUtilities::BasisKey tobasis0 = DetFaceBasisKey(face,0);
LibUtilities::BasisKey tobasis1 = DetFaceBasisKey(face,1);
// Number of quadrature points in face expansion.
int nq = m_base[0]->GetNumPoints()*m_base[0]->GetNumPoints();
int nq_face = tobasis0.GetNumPoints()*tobasis1.GetNumPoints();
int vCoordDim = GetCoordim();
int i;
......@@ -524,7 +528,7 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> > &normal = m_faceNormals[face];
for (i = 0; i < vCoordDim; ++i)
{
normal[i] = Array<OneD, NekDouble>(nq);
normal[i] = Array<OneD, NekDouble>(nq_face);
}
// Regular geometry case
......@@ -539,7 +543,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-df[3*i+2][0],normal[i],1);
normal[i][0] = -df[3*i+2][0];
}
break;
}
......@@ -547,7 +551,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-df[3*i+1][0],normal[i],1);
normal[i][0] = -df[3*i+1][0];
}
break;
}
......@@ -555,7 +559,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,df[3*i][0]+df[3*i+2][0],normal[i],1);
normal[i][0] = df[3*i][0]+df[3*i+2][0];
}
break;
}
......@@ -563,7 +567,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,df[3*i+1][0]+df[3*i+2][0],normal[i],1);
normal[i][0] = df[3*i+1][0]+df[3*i+2][0];
}
break;
}
......@@ -571,7 +575,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-df[3*i][0],normal[i],1);
normal[i][0] = -df[3*i][0];
}
break;
}
......@@ -588,7 +592,7 @@ namespace Nektar
fac = 1.0/sqrt(fac);
for (i = 0; i < vCoordDim; ++i)
{
Vmath::Smul(nq,fac,normal[i],1,normal[i],1);
Vmath::Fill(nq_face,fac*normal[i][0],normal[i],1);
}
}
else
......@@ -619,7 +623,6 @@ namespace Nektar
LibUtilities::PointsKey points0;
LibUtilities::PointsKey points1;
Array<OneD, NekDouble> work (nq, 0.0);
Array<OneD, NekDouble> normals(vCoordDim*nqtot,0.0);
// Extract Jacobian along face and recover local derivatives
......@@ -729,37 +732,38 @@ namespace Nektar
ASSERTL0(false,"face is out of range (face < 4)");
}
Array<OneD, NekDouble> work (nq_face, 0.0);
// Interpolate Jacobian and invert
LibUtilities::Interp2D(points0, points1, jac,
m_base[0]->GetPointsKey(),
m_base[0]->GetPointsKey(),
tobasis0.GetPointsKey(),
tobasis1.GetPointsKey(),
work);
Vmath::Sdiv(nq, 1.0, &work[0], 1, &work[0], 1);
Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
// Interpolate normal and multiply by inverse Jacobian.
for(i = 0; i < vCoordDim; ++i)
{
LibUtilities::Interp2D(points0, points1,
&normals[i*nqtot],
m_base[0]->GetPointsKey(),
m_base[0]->GetPointsKey(),
tobasis0.GetPointsKey(),
tobasis1.GetPointsKey(),
&normal[i][0]);
Vmath::Vmul(nq,work,1,normal[i],1,normal[i],1);
Vmath::Vmul(nq_face,work,1,normal[i],1,normal[i],1);
}
// Normalise to obtain unit normals.
Vmath::Zero(nq,work,1);
Vmath::Zero(nq_face,work,1);
for(i = 0; i < GetCoordim(); ++i)
{
Vmath::Vvtvp(nq,normal[i],1,normal[i],1,work,1,work,1);
Vmath::Vvtvp(nq_face,normal[i],1,normal[i],1,work,1,work,1);
}
Vmath::Vsqrt(nq,work,1,work,1);
Vmath::Sdiv (nq,1.0,work,1,work,1);
Vmath::Vsqrt(nq_face,work,1,work,1);
Vmath::Sdiv (nq_face,1.0,work,1,work,1);
for(i = 0; i < GetCoordim(); ++i)
{
Vmath::Vmul(nq,normal[i],1,work,1,normal[i],1);
Vmath::Vmul(nq_face,normal[i],1,work,1,normal[i],1);
}
}
}
......
......@@ -648,8 +648,12 @@ namespace Nektar
//Directions A and B positive
Vmath::Vcopy(nquad0*nquad1,inarray.get(),1,o_tmp.get(),1);
//interpolate
LibUtilities::Interp2D(m_base[0]->GetPointsKey(), m_base[1]->GetPointsKey(), o_tmp.get(),
FaceExp->GetBasis(0)->GetPointsKey(),FaceExp->GetBasis(1)->GetPointsKey(),o_tmp2.get());
LibUtilities::Interp2D(m_base[0]->GetPointsKey(),
m_base[1]->GetPointsKey(),
o_tmp.get(),
FaceExp->GetBasis(0)->GetPointsKey(),
FaceExp->GetBasis(1)->GetPointsKey(),
o_tmp2.get());
break;
}
case 1:
......@@ -660,8 +664,12 @@ namespace Nektar
Vmath::Vcopy(nquad0,inarray.get()+(nquad0*nquad1*k),1,o_tmp.get()+(k*nquad0),1);
}
//interpolate
LibUtilities::Interp2D(m_base[0]->GetPointsKey(), m_base[2]->GetPointsKey(), o_tmp.get(),
FaceExp->GetBasis(0)->GetPointsKey(),FaceExp->GetBasis(1)->GetPointsKey(),o_tmp2.get());
LibUtilities::Interp2D(m_base[0]->GetPointsKey(),
m_base[2]->GetPointsKey(),
o_tmp.get(),
FaceExp->GetBasis(0)->GetPointsKey(),
FaceExp->GetBasis(1)->GetPointsKey(),
o_tmp2.get());
break;
}
case 2:
......@@ -717,14 +725,20 @@ namespace Nektar
const Array<TwoD, const NekDouble> &df = geomFactors->GetDerivFactors(ptsKeys);
const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
int nq = m_base[0]->GetNumPoints()*m_base[0]->GetNumPoints();
LibUtilities::BasisKey tobasis0 = DetFaceBasisKey(face,0);
LibUtilities::BasisKey tobasis1 = DetFaceBasisKey(face,1);
// number of face quadrature points
int nq_face = tobasis0.GetNumPoints()*tobasis1.GetNumPoints();
int vCoordDim = GetCoordim();
m_faceNormals[face] = Array<OneD, Array<OneD, NekDouble> >(vCoordDim);
Array<OneD, Array<OneD, NekDouble> > &normal = m_faceNormals[face];
for (i = 0; i < vCoordDim; ++i)
{
normal[i] = Array<OneD, NekDouble>(nq);
normal[i] = Array<OneD, NekDouble>(nq_face);
}
// Regular geometry case
......@@ -740,7 +754,7 @@ namespace Nektar
{
for (i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-df[3*i+2][0],normal[i],1);
normal[i][0] = -df[3*i+2][0];
}
break;
......@@ -749,7 +763,7 @@ namespace Nektar
{
for (i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-df[3*i+1][0],normal[i],1);
normal[i][0] = -df[3*i+1][0];
}
break;
......@@ -758,8 +772,8 @@ namespace Nektar
{
for (i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,df[3*i][0]+df[3*i+1][0]+
df[3*i+2][0],normal[i],1);
normal[i][0] = df[3*i][0]+df[3*i+1][0]+
df[3*i+2][0];
}
break;
......@@ -768,7 +782,7 @@ namespace Nektar
{
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-df[3*i][0],normal[i],1);
normal[i][0] = -df[3*i][0];
}
break;
}
......@@ -783,9 +797,10 @@ namespace Nektar
fac += normal[i][0]*normal[i][0];
}
fac = 1.0/sqrt(fac);
for (i = 0; i < vCoordDim; ++i)
{
Vmath::Smul(nq,fac,normal[i],1,normal[i],1);
Vmath::Fill(nq_face,fac*normal[i][0],normal[i],1);
}
}
else
......@@ -798,7 +813,8 @@ namespace Nektar
int nq2 = ptsKeys[2].GetNumPoints();
int nqtot;
int nq01 =nq0*nq1;
// number of elemental quad points
if (face == 0)
{
nqtot = nq01;
......@@ -811,11 +827,10 @@ namespace Nektar
{
nqtot = nq1*nq2;
}
LibUtilities::PointsKey points0;
LibUtilities::PointsKey points1;
Array<OneD,NekDouble> work (nq, 0.0);
Array<OneD,NekDouble> normals(vCoordDim*nqtot, 0.0);
// Extract Jacobian along face and recover local derivates
......@@ -907,37 +922,38 @@ namespace Nektar
ASSERTL0(false,"face is out of range (face < 3)");
}
Array<OneD,NekDouble> work (nq_face, 0.0);
// Interpolate Jacobian and invert
LibUtilities::Interp2D(points0, points1, jac,
m_base[0]->GetPointsKey(),
m_base[0]->GetPointsKey(),
tobasis0.GetPointsKey(),
tobasis1.GetPointsKey(),
work);
Vmath::Sdiv(nq, 1.0, &work[0], 1, &work[0], 1);
Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
// Interpolate normal and multiply by inverse Jacobian.
for(i = 0; i < vCoordDim; ++i)
{
LibUtilities::Interp2D(points0, points1,
&normals[i*nqtot],
m_base[0]->GetPointsKey(),
m_base[0]->GetPointsKey(),
tobasis0.GetPointsKey(),
tobasis1.GetPointsKey(),
&normal[i][0]);
Vmath::Vmul(nq,work,1,normal[i],1,normal[i],1);
Vmath::Vmul(nq_face,work,1,normal[i],1,normal[i],1);
}
// Normalise to obtain unit normals.
Vmath::Zero(nq,work,1);
Vmath::Zero(nq_face,work,1);
for(i = 0; i < GetCoordim(); ++i)
{
Vmath::Vvtvp(nq,normal[i],1,normal[i],1,work,1,work,1);
Vmath::Vvtvp(nq_face,normal[i],1,normal[i],1,work,1,work,1);
}
Vmath::Vsqrt(nq,work,1,work,1);
Vmath::Sdiv (nq,1.0,work,1,work,1);
Vmath::Vsqrt(nq_face,work,1,work,1);
Vmath::Sdiv (nq_face,1.0,work,1,work,1);
for(i = 0; i < GetCoordim(); ++i)
{
Vmath::Vmul(nq,normal[i],1,work,1,normal[i],1);
Vmath::Vmul(nq_face,normal[i],1,work,1,normal[i],1);
}
}
}
......
......@@ -993,6 +993,7 @@ namespace Nektar
return it->second;
}
/**
* Retrieve the basis key for a given face direction.
*/
......@@ -1027,127 +1028,24 @@ namespace Nektar
// direction of the element which corresponds to the requested
// coordinate direction of the given face.
int dir = geom3d->GetDir((*elements)[0]->m_FaceIndx, facedir);
;
// Obtain the number of modes for the element basis key in this
// direction.
int nummodes = (int) expansion->m_basisKeyVector[dir].GetNumModes();
int numpoints = (int) expansion->m_basisKeyVector[dir].GetNumPoints();
switch(expansion->m_basisKeyVector[dir].GetBasisType())
if(face->GetNumVerts() == 3)
{
case LibUtilities::eModified_A:
{
const LibUtilities::PointsKey pkey(numpoints, LibUtilities::eGaussLobattoLegendre);
return LibUtilities::BasisKey(LibUtilities::eModified_A,nummodes,pkey);
}
case LibUtilities::eModified_B:
case LibUtilities::eModified_C:
{
switch (facedir)
{
case 0:
{
const LibUtilities::PointsKey pkey(numpoints+1, LibUtilities::eGaussLobattoLegendre);
return LibUtilities::BasisKey(LibUtilities::eModified_A,nummodes,pkey);
}
break;
case 1:
{
if (face->GetNumVerts() == 3)
{
const LibUtilities::PointsKey pkey(numpoints+1,
LibUtilities::eGaussLobattoLegendre);
// Triangle
return LibUtilities::BasisKey(LibUtilities::eModified_B,nummodes,pkey);
}
else
{
const LibUtilities::PointsKey pkey(numpoints+1,
LibUtilities::eGaussLobattoLegendre);
// Quadrilateral
return LibUtilities::BasisKey(LibUtilities::eModified_A,nummodes,pkey);
}
}
break;
default:
ASSERTL0(false,"invalid value to flag");
break;
}
}
break;
case LibUtilities::eGLL_Lagrange:
{
TriGeomSharedPtr triangle = boost::dynamic_pointer_cast<TriGeom>(face);
QuadGeomSharedPtr quadrilateral = boost::dynamic_pointer_cast<QuadGeom>(face);
if(quadrilateral)
{
const LibUtilities::PointsKey pkey(numpoints,LibUtilities::eGaussLobattoLegendre);
return LibUtilities::BasisKey(LibUtilities::eGLL_Lagrange,nummodes,pkey);
}
else if(triangle)
{
switch (facedir)
{
case 0:
{
const LibUtilities::PointsKey pkey(numpoints,LibUtilities::eGaussLobattoLegendre);
return LibUtilities::BasisKey(LibUtilities::eOrtho_A,nummodes,pkey);
}
break;
case 1:
{
const LibUtilities::PointsKey pkey(numpoints,LibUtilities::eGaussRadauMAlpha1Beta0);
return LibUtilities::BasisKey(LibUtilities::eOrtho_B,nummodes,pkey);
}
break;
default:
ASSERTL0(false,"invalid value to flag");
break;
}
}
else
{
ASSERTL0(false,"dynamic cast to a proper Geometry2D failed");
}
}
break;
case LibUtilities::eOrtho_A:
{
switch (facedir)
{
case 0:
{
const LibUtilities::PointsKey pkey(numpoints,LibUtilities::eGaussLobattoLegendre);
return LibUtilities::BasisKey(LibUtilities::eOrtho_A,nummodes,pkey);
}
break;
case 1:
{
const LibUtilities::PointsKey pkey(numpoints,LibUtilities::eGaussRadauMAlpha1Beta0);
return LibUtilities::BasisKey(LibUtilities::eOrtho_B,nummodes,pkey);
}
break;
default:
ASSERTL0(false,"invalid value to flag");
break;
}
}
break;
// case eGLL_Lagrange_SEM:
// {
// const LibUtilities::PointsKey pkey(nummodes,LibUtilities::eGaussLobattoLegendre);
// return LibUtilities::BasisKey(LibUtilities::eGLL_Lagrange,nummodes,pkey);
// }
// break;
default:
ASSERTL0(false,"expansion type unknown");
break;
return StdRegions::EvaluateTriFaceBasisKey(facedir,
expansion->m_basisKeyVector[dir].GetBasisType(),
expansion->m_basisKeyVector[dir].GetNumPoints(),
expansion->m_basisKeyVector[dir].GetNumModes());
}
return LibUtilities::NullBasisKey; // Keep things happy by returning a value.
else
{
return StdRegions::EvaluateQuadFaceBasisKey(facedir,
expansion->m_basisKeyVector[dir].GetBasisType(),
expansion->m_basisKeyVector[dir].GetNumPoints(),
expansion->m_basisKeyVector[dir].GetNumModes());
}
// Keep things happy by returning a value.
return LibUtilities::NullBasisKey;
}
......
......@@ -52,6 +52,28 @@ namespace Nektar
class PrismGeom;
class HexGeom;
/// Rule defining bassis expansion given information about 3D expansion
LibUtilities::BasisKey EvaluateFaceBasisKey(const int facedir,
const LibUtilities::BasisType
faceDirBasisType,
const int numpoints,
const int nummodes,
const int nverts);
LibUtilities::BasisKey EvaluateQuadFaceBasisKey(const int facedir,
const LibUtilities::BasisType
faceDirBasisType,
const int numpoints,
const int nummodes);
LibUtilities::BasisKey EvaluateTriFaceBasisKey(const int facedir,
const LibUtilities::BasisType
faceDirBasisType,
const int numpoints,
const int nummodes);
class MeshGraph3D: public MeshGraph
{
public:
......
......@@ -67,6 +67,8 @@ namespace Nektar
eNotFilled, ///< Geometric information has not been generated.
ePtsFilled ///< Geometric information has been generated.
};
}; // end of namespace
}; // end of namespace
......
......@@ -46,6 +46,7 @@ namespace Nektar
{
namespace StdRegions
{
StdExpansion3D::StdExpansion3D()
{
}
......@@ -341,5 +342,155 @@ namespace Nektar
m_faceNormals[face][i], 1);
}
}
LibUtilities::BasisKey EvaluateQuadFaceBasisKey(const int facedir,
const LibUtilities::BasisType
faceDirBasisType,
const int numpoints,
const int nummodes)
{
switch(faceDirBasisType)
{
case LibUtilities::eModified_A:
{
const LibUtilities::PointsKey pkey(numpoints, LibUtilities::eGaussLobattoLegendre);
return LibUtilities::BasisKey(LibUtilities::eModified_A,nummodes,pkey);
}
break;
case LibUtilities::eModified_B:
case LibUtilities::eModified_C:
{
const LibUtilities::PointsKey pkey(numpoints+1, LibUtilities::eGaussLobattoLegendre);