Skip to content

GitLab

Commit b635336c authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'master' into cardiac-fibre

parents 5f6999a8 966d08c6
Hide whitespace changes
Inline Side-by-side
library/LibUtilities/Interpreter/AnalyticExpressionEvaluator.cpp
View file @ b635336c
......@@ -684,7 +684,9 @@ namespace Nektar
boost::algorithm::trim(valueStr);
const parser_id parserID = location->value.id();
#if defined(NEKTAR_DEBUG) || defined(NEKTAR_FULLDEBUG)
const int num_children = location->children.size();
#endif
if (parserID == AnalyticExpression::constantID)
{
......
library/LocalRegions/Expansion3D.cpp
View file @ b635336c
......@@ -938,7 +938,11 @@ namespace Nektar
* normals must be negated, since the integral being performed
* here requires an outwards facing normal.
*/
if (locExp->GetRightAdjacentElementFace() != -1)
if (m_negatedNormals[face])
{
Vmath::Neg(order_e,FaceExp->UpdateCoeffs(),1);
}
else if (locExp->GetRightAdjacentElementFace() != -1)
{
if (locExp->GetRightAdjacentElementExp()->GetGeom3D()->GetGlobalID()
== GetGeom3D()->GetGlobalID())
......
library/LocalRegions/HexExp.cpp
View file @ b635336c
......@@ -1292,7 +1292,28 @@ namespace Nektar
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 nqe0 = m_base[0]->GetNumPoints();
int nqe1 = m_base[1]->GetNumPoints();
int nqe2 = m_base[2]->GetNumPoints();
int nqe01 = nqe0*nqe1;
int nqe02 = nqe0*nqe2;
int nqe12 = nqe1*nqe2;
int nqe;
if (face == 0 || face == 5)
{
nqe = nqe01;
}
else if (face == 1 || face == 3)
{
nqe = nqe02;
}
else
{
nqe = nqe12;
}
int vCoordDim = GetCoordim();
m_faceNormals[face] = Array<OneD, Array<OneD, NekDouble> >(vCoordDim);
......@@ -1360,176 +1381,90 @@ namespace Nektar
Vmath::Smul(nqe,fac,normal[i],1,normal[i],1);
}
}
}
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;
LibUtilities::PointsKey points0;
LibUtilities::PointsKey points1;
{
int j, k;
int nq;
Array<OneD,NekDouble> work(nqe,0.0);
Array<OneD,NekDouble> normals(vCoordDim*nquad0*nquad1,0.0);
//Array<OneD,NekDouble> normals(vCoordDim*nqe,0.0);
// Extract Jacobian along face and recover local
// derivates (dx/dr) for polynomial interpolation by
// multiplying m_gmat by jacobian
switch(face)
{
case 0:
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);
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)
{
for(k=0; k<nquad2; ++k)
{
case 0:
for(j = 0; j < nqe; ++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);
// 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 2:
for (j=0; j< nquad1; ++j)
{
for(k=0; k<nquad2; ++k)
normal[0][j] = -gmat[2][j]*jac[j];
normal[1][j] = -gmat[5][j]*jac[j];
normal[2][j] = -gmat[8][j]*jac[j];
}
break;
case 1:
for (j = 0; j< nqe0; ++j)
{
normals[j+k*nquad0] = gmat[0][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]*jac[nquad0-1+nquad0*j+nquad0*nquad1*k];
normals[2*nqtot+j+k*nquad0] = gmat[6][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 < nqe2; ++k)
{
normal[0][j+k*nqe0] = -gmat[1][j+nqe01*k]*jac[j+nqe01*k];
normal[1][j+k*nqe0] = -gmat[4][j+nqe01*k]*jac[j+nqe01*k];
normal[2][j+k*nqe0] = -gmat[7][j+nqe01*k]*jac[j+nqe01*k];
}
}
break;
case 2:
for (j=0; j< nqe1; ++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<nqe2; ++k)
{
normal[0][j+k*nqe0] = gmat[0][nqe0-1+nqe0*j+nqe01*k]*jac[nqe0-1+nqe0*j+nqe01*k];
normal[1][j+k*nqe0] = gmat[3][nqe0-1+nqe0*j+nqe01*k]*jac[nqe0-1+nqe0*j+nqe01*k];
normal[2][j+k*nqe0] = gmat[6][nqe0-1+nqe0*j+nqe01*k]*jac[nqe0-1+nqe0*j+nqe01*k];
}
}
break;
case 3:
for (j=0; j< nqe0; ++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);
for(k=0; k<nqe2; ++k)
{
normal[0][j+k*nqe0] = gmat[1][nqe0*(nqe1-1)+j+nqe0*nqe1*k]*jac[nqe0*(nqe1-1)+j+nqe01*k];
normal[1][j+k*nqe0] = gmat[4][nqe0*(nqe1-1)+j+nqe0*nqe1*k]*jac[nqe0*(nqe1-1)+j+nqe01*k];
normal[1][j+k*nqe0] = gmat[7][nqe0*(nqe1-1)+j+nqe0*nqe1*k]*jac[nqe0*(nqe1-1)+j+nqe01*k];
}
}
break;
case 4:
for (j=0; j< nqe0; ++j)
{
for(k=0; k<nqe2; ++k)
{
normal[0][j+k*nqe0] = -gmat[0][j*nqe0+nqe01*k]*jac[j*nqe0+nqe01*k];
normal[1][j+k*nqe0] = -gmat[3][j*nqe0+nqe01*k]*jac[j*nqe0+nqe01*k];
normal[2][j+k*nqe0] = -gmat[6][j*nqe0+nqe01*k]*jac[j*nqe0+nqe01*k];
}
}
break;
case 5:
for (j=0; j< nqe01; ++j)
{
normal[0][j] = gmat[2][j+nqe01*(nqe2-1)]*jac[j+nqe01*(nqe2-1)];
normal[1][j] = gmat[5][j+nqe01*(nqe2-1)]*jac[j+nqe01*(nqe2-1)];
normal[2][j] = gmat[8][j+nqe01*(nqe2-1)]*jac[j+nqe01*(nqe2-1)];
}
break;
default:
ASSERTL0(false,"face is out of range (face < 5)");
}
// 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 5:
for (j=0; j< nquad0*nquad1; ++j)
{
normals[j] = gmat[2][j+nquad0*nquad1*(nquad2-1)]*jac[j+nquad0*nquad1*(nquad2-1)];
normals[nqtot+j] = gmat[5][j+nquad0*nquad1*(nquad2-1)]*jac[j+nquad0*nquad1*(nquad2-1)];
normals[2*nqtot+j] = gmat[8][j+nquad0*nquad1*(nquad2-1)]*jac[j+nquad0*nquad1*(nquad2-1)];
}
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);
Vmath::Sdiv(nqe,1.0,&jac[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;
default:
ASSERTL0(false,"face is out of range (face < 5)");
}
// interpolate
for(i = 0; i < GetCoordim(); ++i)
{
Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
}
//normalise normal vectors
Vmath::Zero(nqe,work,1);
......@@ -1544,7 +1479,7 @@ namespace Nektar
for(i = 0; i < GetCoordim(); ++i)
{
Vmath::Vmul(nqe,normal[i],1,work,1,normal[i],1);
}
}
}
}
......@@ -2362,28 +2297,21 @@ namespace Nektar
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble> &outarray)
{
const int nqtot = m_base[0]->GetNumPoints() *
m_base[1]->GetNumPoints() *
m_base[2]->GetNumPoints();
if(m_metricinfo->IsUsingQuadMetrics())
{
int nqtot = m_base[0]->GetNumPoints()
* m_base[1]->GetNumPoints()
* m_base[2]->GetNumPoints();
const Array<OneD, const NekDouble>& metric
= m_metricinfo->GetQuadratureMetrics();
const Array<OneD, const NekDouble> &metric
= m_metricinfo->GetQuadratureMetrics();
Vmath::Vmul(nqtot, metric, 1, inarray, 1, outarray, 1);
}
else
{
int i;
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nquad2 = m_base[2]->GetNumPoints();
int nqtot = nquad0*nquad1*nquad2;
const Array<OneD, const NekDouble>& jac = m_metricinfo->GetJac();
const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
const Array<OneD, const NekDouble>& w2 = m_base[2]->GetW();
const Array<OneD, const NekDouble> &jac
= m_metricinfo->GetJac();
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
......@@ -2394,101 +2322,9 @@ namespace Nektar
Vmath::Smul(nqtot, jac[0], inarray, 1, outarray, 1);
}
// First coordinate
for(i = 0; i < nquad1*nquad2; ++i)
{
Vmath::Vmul(nquad0, outarray.get()+i*nquad0, 1,
w0.get(), 1, outarray.get()+i*nquad0,1);
}
// Second coordinate
for(i = 0; i < nquad1*nquad2; ++i)
{
Vmath::Smul(nquad0, w1[i%nquad2], outarray.get()+i*nquad0, 1,
outarray.get()+i*nquad0, 1);
}
// Third coordinate
for(i = 0; i < nquad2; ++i)
{
Vmath::Smul(nquad0*nquad1, w2[i], outarray.get()+i*nquad0*nquad1, 1,
outarray.get()+i*nquad0*nquad1, 1);
}
/*
// multiply by integration constants
for(i = 0; i < nquad1; ++i)
{
Vmath::Vmul(nquad0, outarray.get()+i*nquad0,1,
w0.get(),1,outarray.get()+i*nquad0,1);
}
for(i = 0; i < nquad0; ++i)
{
Vmath::Vmul(nquad1,outarray.get()+i,nquad0,w1.get(),1,
outarray.get()+i,nquad0);
}
*/ }
}
// DNekScalMatSharedPtr& HexExp::v_GetLocMatrix(
// const StdRegions::MatrixType mtype,
// NekDouble lambdaval,
// NekDouble tau)
// {
// MatrixKey mkey( mtype,DetExpansionType(),*this,lambdaval,tau );
// return m_matrixManager[mkey];
// }
//
// DNekScalMatSharedPtr& HexExp::v_GetLocMatrix(
// const StdRegions::MatrixType mtype,
// const Array<OneD, NekDouble> &dir1Forcing,
// NekDouble lambdaval,
// NekDouble tau)
// {
// MatrixKey mkey( mtype,DetExpansionType(),*this,lambdaval,tau,
// dir1Forcing);
// return m_matrixManager[mkey];
// }
//
// DNekScalMatSharedPtr& HexExp::v_GetLocMatrix(
// const StdRegions::MatrixType mtype,
// const Array<OneD, Array<OneD, const NekDouble> >&
// dirForcing,
// NekDouble lambdaval,
// NekDouble tau)
// {
// MatrixKey mkey( mtype,DetExpansionType(),*this,lambdaval,tau,
// dirForcing);
// return m_matrixManager[mkey];
// }
/*
void HexExp::v_IProductWRTBase_SumFac(const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble> &outarray)
{
cout << "Hex::IProduct::SumFac" << endl;
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nquad2 = m_base[2]->GetNumPoints();
int order0 = m_base[0]->GetNumModes();
int order1 = m_base[1]->GetNumModes();
Array<OneD,NekDouble> tmp(nquad0*nquad1*nquad2);
Array<OneD,NekDouble> wsp(nquad0*nquad1*(nquad2+order0)
+ order0*order1*nquad2);
MultiplyByQuadratureMetric(inarray, tmp);
StdHexExp::IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
m_base[1]->GetBdata(),
m_base[2]->GetBdata(),
tmp, outarray, wsp,
true, true, true);
StdHexExp::MultiplyByQuadratureMetric(outarray, outarray);
}
}
*/
void HexExp::LaplacianMatrixOp_MatFree_Kernel(
const Array<OneD, const NekDouble> &inarray,
......
library/LocalRegions/HexExp.h
View file @ b635336c
......@@ -162,7 +162,7 @@ namespace Nektar
LOCAL_REGIONS_EXPORT virtual
StdRegions::ExpansionType v_DetExpansionType() const;
LOCAL_REGIONS_EXPORT virtual const
LOCAL_REGIONS_EXPORT virtual const
SpatialDomains::GeomFactorsSharedPtr& v_GetMetricInfo() const;
LOCAL_REGIONS_EXPORT virtual const
......
library/LocalRegions/PrismExp.cpp
View file @ b635336c
......@@ -1758,76 +1758,36 @@ namespace Nektar
return returnval;
}
/// @todo add functionality for IsUsingQuadMetrics
void PrismExp::MultiplyByQuadratureMetric(
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble>& outarray)
{
int i, j;
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nquad2 = m_base[2]->GetNumPoints();
int nqtot = nquad0*nquad1*nquad2;
const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
const Array<OneD, const NekDouble>& w2 = m_base[2]->GetW();
const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
const Array<OneD, const NekDouble>& jac = m_metricinfo->GetJac();
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
const int nqtot = m_base[0]->GetNumPoints() *
m_base[1]->GetNumPoints() *
m_base[2]->GetNumPoints();
if(m_metricinfo->IsUsingQuadMetrics())
{
Vmath::Vmul(nqtot, jac, 1, inarray, 1, outarray, 1);
const Array<OneD, const NekDouble> &metric
= m_metricinfo->GetQuadratureMetrics();
Vmath::Vmul(nqtot, metric, 1, inarray, 1, outarray, 1);
}
else
{
Vmath::Smul(nqtot, jac[0], inarray, 1, outarray, 1);
}
// Multiply by integration constants in x-direction
for(i = 0; i < nquad1*nquad2; ++i)
{
Vmath::Vmul(nquad0, outarray.get()+i*nquad0, 1,
w0.get(), 1, outarray.get()+i*nquad0,1);
}
// Multiply by integration constants in y-direction
for(j = 0; j < nquad2; ++j)
{
for(i = 0; i < nquad1; ++i)
const Array<OneD, const NekDouble> &jac
= m_metricinfo->GetJac();
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
Blas::Dscal(nquad0,w1[i], &outarray[0]+i*nquad0 +
j*nquad0*nquad1,1);
Vmath::Vmul(nqtot, jac, 1, inarray, 1, outarray, 1);
}
}
// Multiply by integration constants in z-direction; need to
// incorporate factor (1-eta_3)/2 into weights, but only if using
// GLL quadrature points.
switch(m_base[2]->GetPointsType())
{
// Legendre inner product.
case LibUtilities::eGaussLobattoLegendre:
for(i = 0; i < nquad2; ++i)
{
Blas::Dscal(nquad0*nquad1,0.25*(1-z2[i])*w2[i],
&outarray[0]+i*nquad0*nquad1,1);
}
break;
// (1,0) Jacobi inner product.
case LibUtilities::eGaussRadauMAlpha1Beta0:
for(i = 0; i < nquad2; ++i)
{
Blas::Dscal(nquad0*nquad1, 0.5*w2[i],
&outarray[0]+i*nquad0*nquad1, 1);
}
break;
default:
ASSERTL0(false, "Quadrature point type not supported for this element.");
break;
else
{
Vmath::Smul(nqtot, jac[0], inarray, 1, outarray, 1);
}
StdPrismExp::MultiplyByQuadratureMetric(outarray, outarray);
}
}
......
library/LocalRegions/QuadExp.cpp
View file @ b635336c
......@@ -2270,24 +2270,21 @@ namespace Nektar
void QuadExp::MultiplyByQuadratureMetric(const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble> &outarray)
{
const int nqtot = m_base[0]->GetNumPoints() *
m_base[1]->GetNumPoints();
if(m_metricinfo->IsUsingQuadMetrics())