Skip to content

GitLab

Commit f412f430 authored by Dave Moxey's avatar Dave Moxey
Browse files

Fixed prism/tet deformed normals. 


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@4037 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent dabd9f4a
Hide whitespace changes
Inline Side-by-side
library/LocalRegions/PrismExp.cpp
View file @ f412f430
......@@ -894,255 +894,239 @@ namespace Nektar
void PrismExp::v_ComputeFaceNormal(const int face)
{
int i;
const SpatialDomains::GeomFactorsSharedPtr & geomFactors = GetGeom()->GetMetricInfo();
SpatialDomains::GeomType type = geomFactors->GetGtype();
const Array<TwoD, const NekDouble> & gmat = geomFactors->GetGmat();
const Array<OneD, const NekDouble> & jac = geomFactors->GetJac();
int nqe = m_base[0]->GetNumPoints()*m_base[1]->GetNumPoints();
//int nqe;
const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
GetGeom()->GetMetricInfo();
SpatialDomains::GeomType type = geomFactors->GetGtype();
const Array<TwoD, const NekDouble> &gmat = geomFactors->GetGmat();
const Array<OneD, const NekDouble> &jac = geomFactors->GetJac();
// Number of quadrature points in face expansion.
int nq = m_base[0]->GetNumPoints()*m_base[0]->GetNumPoints();
int vCoordDim = GetCoordim();
/*
switch(face)
{
case 0:
nqe = m_base[0]->GetNumPoints()*m_base[1]->GetNumPoints();
break;
case 1:
case 3:
nqe = m_base[0]->GetNumPoints()*m_base[2]->GetNumPoints();
break;
case 2:
case 4:
nqe = m_base[1]->GetNumPoints()*m_base[2]->GetNumPoints();
break;
}
*/
int i;
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>(nqe);
normal[i] = Array<OneD, NekDouble>(nq);
}
// Regular geometry case
if((type == SpatialDomains::eRegular)||(type == SpatialDomains::eMovingRegular))
if (type == SpatialDomains::eRegular ||
type == SpatialDomains::eMovingRegular)
{
NekDouble fac;
// Set up normals
switch(face)
{
case 0:
for(i = 0; i < vCoordDim; ++i)
case 0:
{
Vmath::Fill(nqe,-gmat[3*i+2][0],normal[i],1);
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-gmat[3*i+2][0],normal[i],1);
}
break;
}
break;
case 1:
for(i = 0; i < vCoordDim; ++i)
case 1:
{
Vmath::Fill(nqe,-gmat[3*i+1][0],normal[i],1);
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-gmat[3*i+1][0],normal[i],1);
}
break;
}
break;
case 2:
for(i = 0; i < vCoordDim; ++i)
case 2:
{
Vmath::Fill(nqe,gmat[3*i][0]+gmat[3*i+2][0],normal[i],1);
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,gmat[3*i][0]+gmat[3*i+2][0],normal[i],1);
}
break;
}
break;
case 3:
for(i = 0; i < vCoordDim; ++i)
case 3:
{
Vmath::Fill(nqe,gmat[3*i+1][0],normal[i],1);
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,gmat[3*i+1][0],normal[i],1);
}
break;
}
break;
case 4:
for(i = 0; i < vCoordDim; ++i)
case 4:
{
Vmath::Fill(nqe,-gmat[3*i][0],normal[i],1);
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-gmat[3*i][0],normal[i],1);
}
break;
}
break;
default:
ASSERTL0(false,"face is out of range (face < 4)");
default:
ASSERTL0(false,"face is out of range (face < 4)");
}
// normalise
// Normalise resulting vector.
fac = 0.0;
for(i =0 ; i < vCoordDim; ++i)
for(i = 0; i < vCoordDim; ++i)
{
fac += normal[i][0]*normal[i][0];
}
fac = 1.0/sqrt(fac);
for (i = 0; i < vCoordDim; ++i)
{
Vmath::Smul(nqe,fac,normal[i],1,normal[i],1);
Vmath::Smul(nq,fac,normal[i],1,normal[i],1);
}
}
else // Set up deformed normals
{
else
{
// Set up deformed normals.
int j, k;
int nquad0 = geomFactors->GetPointsKey(0).GetNumPoints();
int nquad1 = geomFactors->GetPointsKey(1).GetNumPoints();
int nquad2 = geomFactors->GetPointsKey(2).GetNumPoints();
int nqtot = nquad0*nquad1;
int nqtot;
// Determine number of quadrature points on the face.
if (face == 0)
{
nqtot = nquad0*nquad1;
}
else if (face == 1 || face == 3)
{
nqtot = nquad0*nquad2;
}
else
{
nqtot = nquad1*nquad2;
}
LibUtilities::PointsKey points0;
LibUtilities::PointsKey points1;
int nq;
Array<OneD,NekDouble> work(nqe,0.0);
Array<OneD,NekDouble> normals(vCoordDim*nquad0*nquad1,0.0);
// Extract Jacobian along face and recover local
// derivates (dx/dr) for polynomial interpolation by
// multiplying m_gmat by jacobian
Array<OneD, NekDouble> work (nq, 0.0);
Array<OneD, NekDouble> normals(vCoordDim*nqtot,0.0);
// Extract Jacobian along face and recover local derivatives
// (dx/dr) for polynomial interpolation by multiplying m_gmat by
// jacobian
switch(face)
{
case 0:
for(j = 0; j < nquad0*nquad1; ++j)
case 0:
{
normals[j] = -gmat[2][j]*jac[j];
normals[nqtot+j] = -gmat[5][j]*jac[j];
normals[2*nqtot+j] = -gmat[8][j]*jac[j];
for(j = 0; j < nquad0*nquad1; ++j)
{
normals[j] = -gmat[2][j]*jac[j];
normals[nqtot+j] = -gmat[5][j]*jac[j];
normals[2*nqtot+j] = -gmat[8][j]*jac[j];
}
points0 = geomFactors->GetPointsKey(0);
points1 = geomFactors->GetPointsKey(1);
break;
}
points0 = geomFactors->GetPointsKey(0);
points1 = geomFactors->GetPointsKey(1);
nq=points0.GetNumPoints()*points1.GetNumPoints();
// interpolate Jacobian and invert
LibUtilities::Interp2D(points0,points1,jac,m_base[0]->GetPointsKey(),m_base[1]->GetPointsKey(),work);
Vmath::Sdiv(nq,1.0,&work[0],1,&work[0],1);
// interpolate
for(i = 0; i < GetCoordim(); ++i)
{
LibUtilities::Interp2D(points0,points1,&normals[i*nq],m_base[0]->GetPointsKey(),m_base[1]->GetPointsKey(),&normal[i][0]);
Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
}
break;
case 1:
for (j=0; j< nquad0; ++j)
case 1:
{
for(k=0; k<nquad2; ++k)
int tmp = nquad0*nquad1;
for (j = 0; j < nquad0; ++j)
{
normals[j+k*nquad0] = -gmat[1][j+nquad0*nquad1*k]*jac[j+nquad0*nquad1*k];
normals[nqtot+j+k*nquad0] = -gmat[4][j+nquad0*nquad1*k]*jac[j+nquad0*nquad1*k];
normals[2*nqtot+j+k*nquad0] = -gmat[7][j+nquad0*nquad1*k]*jac[j+nquad0*nquad1*k];
}
}
points0 = geomFactors->GetPointsKey(0);
points1 = geomFactors->GetPointsKey(2);
nq=points0.GetNumPoints()*points1.GetNumPoints();
// interpolate Jacobian and invert
LibUtilities::Interp2D(points0,points1,jac,m_base[0]->GetPointsKey(),m_base[2]->GetPointsKey(),work);
Vmath::Sdiv(nq,1.0,&work[0],1,&work[0],1);
for(k = 0; k < nquad2; ++k)
{
normals[j+k*nquad0] =
-gmat[1][j+tmp*k]*jac[j+tmp*k];
normals[nqtot+j+k*nquad0] =
-gmat[4][j+tmp*k]*jac[j+tmp*k];
normals[2*nqtot+j+k*nquad0] =
-gmat[7][j+tmp*k]*jac[j+tmp*k];
}
}
// interpolate
for(i = 0; i < GetCoordim(); ++i)
{
LibUtilities::Interp2D(points0,points1,&normals[i*nq],m_base[0]->GetPointsKey(),m_base[2]->GetPointsKey(),&normal[i][0]);
Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
points0 = geomFactors->GetPointsKey(0);
points1 = geomFactors->GetPointsKey(2);
break;
}
break;
case 2:
for (j=0; j< nquad1; ++j)
{
for(k=0; k<nquad2; ++k)
case 2:
for (j=0; j< nquad1; ++j)
{
normals[j+k*nquad0] = (gmat[0][nquad0-1+nquad0*j+nquad0*nquad1*k]
+gmat[2][nquad0-1+nquad0*j+nquad0*nquad1*k])*jac[nquad0-1+nquad0*j+nquad0*nquad1*k];
normals[nqtot+j+k*nquad0] = (gmat[3][nquad0-1+nquad0*j+nquad0*nquad1*k]
+gmat[5][nquad0-1+nquad0*j+nquad0*nquad1*k])*jac[nquad0-1+nquad0*j+nquad0*nquad1*k];
normals[2*nqtot+j+k*nquad0] = (gmat[6][nquad0-1+nquad0*j+nquad0*nquad1*k]
+gmat[8][nquad0-1+nquad0*j+nquad0*nquad1*k])*jac[nquad0-1+nquad0*j+nquad0*nquad1*k];
}
}
points0 = geomFactors->GetPointsKey(1);
points1 = geomFactors->GetPointsKey(2);
nq=points0.GetNumPoints()*points1.GetNumPoints();
// interpolate Jacobian and invert
LibUtilities::Interp2D(points0,points1,jac,m_base[1]->GetPointsKey(),m_base[2]->GetPointsKey(),work);
Vmath::Sdiv(nq,1.0,&work[0],1,&work[0],1);
// interpolate
for(i = 0; i < GetCoordim(); ++i)
{
LibUtilities::Interp2D(points0,points1,&normals[i*nq],m_base[1]->GetPointsKey(),m_base[2]->GetPointsKey(),&normal[i][0]);
Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
}
break;
case 3:
for (j=0; j< nquad0; ++j)
{
for(k=0; k<nquad2; ++k)
for(k=0; k<nquad2; ++k)
{
normals[j+k*nquad0] = (gmat[0][nquad0-1+nquad0*j+nquad0*nquad1*k]
+gmat[2][nquad0-1+nquad0*j+nquad0*nquad1*k])*jac[nquad0-1+nquad0*j+nquad0*nquad1*k];
normals[nqtot+j+k*nquad0] = (gmat[3][nquad0-1+nquad0*j+nquad0*nquad1*k]
+gmat[5][nquad0-1+nquad0*j+nquad0*nquad1*k])*jac[nquad0-1+nquad0*j+nquad0*nquad1*k];
normals[2*nqtot+j+k*nquad0] = (gmat[6][nquad0-1+nquad0*j+nquad0*nquad1*k]
+gmat[8][nquad0-1+nquad0*j+nquad0*nquad1*k])*jac[nquad0-1+nquad0*j+nquad0*nquad1*k];
}
}
points0 = geomFactors->GetPointsKey(1);
points1 = geomFactors->GetPointsKey(2);
break;
case 3:
for (j=0; j< nquad0; ++j)
{
normals[j+k*nquad0] = gmat[1][nquad0*(nquad1-1)+j+nquad0*nquad1*k]*jac[nquad0*(nquad1-1)+j+nquad0*nquad1*k];
normals[nqtot+j+k*nquad0] = gmat[4][nquad0*(nquad1-1)+j+nquad0*nquad1*k]*jac[nquad0*(nquad1-1)+j+nquad0*nquad1*k];
normals[2*nqtot+j+k*nquad0] = gmat[7][nquad0*(nquad1-1)+j+nquad0*nquad1*k]*jac[nquad0*(nquad1-1)+j+nquad0*nquad1*k];
}
}
points0 = geomFactors->GetPointsKey(0);
points1 = geomFactors->GetPointsKey(2);
nq=points0.GetNumPoints()*points1.GetNumPoints();
// interpolate Jacobian and invert
LibUtilities::Interp2D(points0,points1,jac,m_base[0]->GetPointsKey(),m_base[2]->GetPointsKey(),work);
Vmath::Sdiv(nq,1.0,&work[0],1,&work[0],1);
// interpolate
for(i = 0; i < GetCoordim(); ++i)
{
LibUtilities::Interp2D(points0,points1,&normals[i*nq],m_base[0]->GetPointsKey(),m_base[2]->GetPointsKey(),&normal[i][0]);
Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
}
break;
case 4:
for (j=0; j< nquad0; ++j)
{
for(k=0; k<nquad2; ++k)
for(k=0; k<nquad2; ++k)
{
normals[j+k*nquad0] = gmat[1][nquad0*(nquad1-1)+j+nquad0*nquad1*k]*jac[nquad0*(nquad1-1)+j+nquad0*nquad1*k];
normals[nqtot+j+k*nquad0] = gmat[4][nquad0*(nquad1-1)+j+nquad0*nquad1*k]*jac[nquad0*(nquad1-1)+j+nquad0*nquad1*k];
normals[2*nqtot+j+k*nquad0] = gmat[7][nquad0*(nquad1-1)+j+nquad0*nquad1*k]*jac[nquad0*(nquad1-1)+j+nquad0*nquad1*k];
}
}
points0 = geomFactors->GetPointsKey(0);
points1 = geomFactors->GetPointsKey(2);
break;
case 4:
for (j=0; j< nquad0; ++j)
{
normals[j+k*nquad0] = -gmat[0][j*nquad0+nquad0*nquad1*k]*jac[j*nquad0+nquad0*nquad1*k];
normals[nqtot+j+k*nquad0] = -gmat[3][j*nquad0+nquad0*nquad1*k]*jac[j*nquad0+nquad0*nquad1*k];
normals[2*nqtot+j+k*nquad0] = -gmat[6][j*nquad0+nquad0*nquad1*k]*jac[j*nquad0+nquad0*nquad1*k];
}
}
points0 = geomFactors->GetPointsKey(1);
points1 = geomFactors->GetPointsKey(2);
nq=points0.GetNumPoints()*points1.GetNumPoints();
// interpolate Jacobian and invert
LibUtilities::Interp2D(points0,points1,jac,m_base[1]->GetPointsKey(),m_base[2]->GetPointsKey(),work);
Vmath::Sdiv(nq,1.0,&work[0],1,&work[0],1);
// interpolate
for(i = 0; i < GetCoordim(); ++i)
{
LibUtilities::Interp2D(points0,points1,&normals[i*nq],m_base[1]->GetPointsKey(),m_base[2]->GetPointsKey(),&normal[i][0]);
Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
}
break;
default:
ASSERTL0(false,"face is out of range (face < 4)");
for(k=0; k<nquad2; ++k)
{
normals[j+k*nquad0] = -gmat[0][j*nquad0+nquad0*nquad1*k]*jac[j*nquad0+nquad0*nquad1*k];
normals[nqtot+j+k*nquad0] = -gmat[3][j*nquad0+nquad0*nquad1*k]*jac[j*nquad0+nquad0*nquad1*k];
normals[2*nqtot+j+k*nquad0] = -gmat[6][j*nquad0+nquad0*nquad1*k]*jac[j*nquad0+nquad0*nquad1*k];
}
}
points0 = geomFactors->GetPointsKey(1);
points1 = geomFactors->GetPointsKey(2);
break;
default:
ASSERTL0(false,"face is out of range (face < 4)");
}
//normalise normal vectors
Vmath::Zero(nqe,work,1);
for(i = 0; i < GetCoordim(); ++i)
// Interpolate Jacobian and invert
LibUtilities::Interp2D(points0, points1, jac,
m_base[0]->GetPointsKey(),
m_base[0]->GetPointsKey(),
work);
Vmath::Sdiv(nq, 1.0, &work[0], 1, &work[0], 1);
// Interpolate normal and multiply by inverse Jacobian.
for(i = 0; i < vCoordDim; ++i)
{
Vmath::Vvtvp(nqe,normal[i],1, normal[i],1,work,1,work,1);
LibUtilities::Interp2D(points0, points1,
&normals[i*nqtot],
m_base[0]->GetPointsKey(),
m_base[0]->GetPointsKey(),
&normal[i][0]);
Vmath::Vmul(nq,work,1,normal[i],1,normal[i],1);
}
Vmath::Vsqrt(nqe,work,1,work,1);
Vmath::Sdiv(nqe,1.0,work,1,work,1);
// Normalise to obtain unit normals.
Vmath::Zero(nq,work,1);
for(i = 0; i < GetCoordim(); ++i)
{
Vmath::Vvtvp(nq,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);
for(i = 0; i < GetCoordim(); ++i)
{
Vmath::Vmul(nqe,normal[i],1,work,1,normal[i],1);
Vmath::Vmul(nq,normal[i],1,work,1,normal[i],1);
}
}
}
......
library/LocalRegions/QuadExp.cpp
View file @ f412f430
......@@ -509,19 +509,30 @@ namespace Nektar
void QuadExp::v_NormVectorIProductWRTBase(
const Array<OneD, const NekDouble> &Fx,
const Array<OneD, const NekDouble> &Fy,
const Array<OneD, const NekDouble> &Fz,
Array< OneD, NekDouble> &outarray)
const Array<OneD, const NekDouble> &Fx,
const Array<OneD, const NekDouble> &Fy,
const Array<OneD, const NekDouble> &Fz,
Array<OneD, NekDouble> &outarray)
{
int nq = m_base[0]->GetNumPoints()*m_base[1]->GetNumPoints();
Array<OneD, NekDouble > Fn(nq);
const Array<OneD, const Array<OneD, NekDouble> > &normals = GetLeftAdjacentElementExp()->GetFaceNormal(GetLeftAdjacentElementFace());
Vmath::Vmul (nq,&Fx[0],1,&normals[0][0], 1,&Fn[0],1);
Vmath::Vvtvp(nq,&Fy[0],1,&normals[1][0],1,&Fn[0],1,&Fn[0],1);
Vmath::Vvtvp(nq,&Fz[0],1,&normals[2][0],1,&Fn[0],1,&Fn[0],1);
Array<OneD, NekDouble> Fn(nq);
const Array<OneD, const Array<OneD, NekDouble> > &normals =
GetLeftAdjacentElementExp()->GetFaceNormal(
GetLeftAdjacentElementFace());
if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
Vmath::Vmul (nq,&normals[0][0],1,&Fx[0],1,&Fn[0],1);
Vmath::Vvtvvtp(nq,&normals[1][0],1,&Fy[0],1,
&normals[2][0],1,&Fz[0],1,&Fn[0],1);
}
else
{
Vmath::Smul (nq,normals[0][0],&Fx[0],1,&Fn[0],1);
Vmath::Svtsvtp(nq,normals[1][0],&Fy[0],1,
normals[2][0],&Fz[0],1,&Fn[0],1);
}
IProductWRTBase(Fn,outarray);
}
......
library/LocalRegions/TetExp.cpp
View file @ f412f430
......@@ -790,223 +790,221 @@ namespace Nektar
void TetExp::v_ComputeFaceNormal(const int face)
{
int i;
const SpatialDomains::GeomFactorsSharedPtr & geomFactors = GetGeom()->GetMetricInfo();
SpatialDomains::GeomType type = geomFactors->GetGtype();
const Array<TwoD, const NekDouble> & gmat = geomFactors->GetGmat();
const Array<OneD, const NekDouble> & jac = geomFactors->GetJac();
int nqe;
//int nqe = m_base[0]->GetNumPoints()*m_base[1]->GetNumPoints();
const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
GetGeom()->GetMetricInfo();
SpatialDomains::GeomType type = geomFactors->GetGtype();
const Array<TwoD, const NekDouble> &gmat = geomFactors->GetGmat();
const Array<OneD, const NekDouble> &jac = geomFactors->GetJac();
int nq = m_base[0]->GetNumPoints()*m_base[0]->GetNumPoints();
int vCoordDim = GetCoordim();
switch(face)
{
case 0:
nqe = m_base[0]->GetNumPoints()*m_base[1]->GetNumPoints();
break;
case 1:
nqe = m_base[0]->GetNumPoints()*m_base[2]->GetNumPoints();
break;
case 2:
case 3:
nqe = m_base[1]->GetNumPoints()*m_base[2]->GetNumPoints();
break;
}
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>(nqe);
normal[i] = Array<OneD, NekDouble>(nq);
}
// Regular geometry case
if((type == SpatialDomains::eRegular)||(type == SpatialDomains::eMovingRegular))
if (type == SpatialDomains::eRegular ||
type == SpatialDomains::eMovingRegular)
{
NekDouble fac;
// Set up normals
switch(face)
switch (face)
{
case 0:
for(i = 0; i < vCoordDim; ++i)
case 0:
{
Vmath::Fill(nqe,-gmat[3*i+2][0],normal[i],1);
for (i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-gmat[3*i+2][0],normal[i],1);
}
break;
}
break;
case 1:
for(i = 0; i < vCoordDim; ++i)
case 1:
{
Vmath::Fill(nqe,-gmat[3*i+1][0],normal[i],1);
for (i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,-gmat[3*i+1][0],normal[i],1);
}
break;
}
break;
case 2:
for(i = 0; i < vCoordDim; ++i)
case 2:
{
Vmath::Fill(nqe,gmat[3*i][0]+gmat[3*i+1][0]+gmat[3*i+2][0],normal[i],1);
for (i = 0; i < vCoordDim; ++i)
{
Vmath::Fill(nq,gmat[3*i][0]+gmat[3*i+1][0]+
gmat[3*i+2][0],normal[i],1);
}
break;