Skip to content

GitLab

Commit 8fdd0e80 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'feature/cfs-optimise' of localhost:nektar

parents bebcf155 89dc3fef
Hide whitespace changes
Inline Side-by-side
library/LocalRegions/Expansion2D.cpp
View file @ 8fdd0e80
......@@ -62,6 +62,32 @@ namespace Nektar
const Array<OneD, const Array<OneD, NekDouble> > normals
= GetEdgeNormal(edge);
if (m_requireNeg.size() == 0)
{
m_requireNeg.resize(GetNedges());
for (int i = 0; i < GetNedges(); ++i)
{
m_requireNeg[i] = false;
if (m_negatedNormals[i])
{
m_requireNeg[i] = true;
}
Expansion1DSharedPtr edgeExp = boost::dynamic_pointer_cast<
Expansion1D>(m_edgeExp[i].lock());
if (edgeExp->GetRightAdjacentElementExp())
{
if (edgeExp->GetRightAdjacentElementExp()->GetGeom2D()
->GetGlobalID() == GetGeom2D()->GetGlobalID())
{
m_requireNeg[i] = true;
}
}
}
}
// We allow the case of mixed polynomial order by supporting only
// those modes on the edge common to both adjoining elements. This
// is enforced here by taking the minimum size and padding with
......@@ -79,6 +105,7 @@ namespace Nektar
boost::dynamic_pointer_cast<
LocalRegions::Expansion1D>(EdgeExp);
if (m_negatedNormals[edge])
{
Vmath::Neg(nquad_e,EdgeExp->UpdatePhys(),1);
......@@ -104,6 +131,32 @@ namespace Nektar
{
int i;
if (m_requireNeg.size() == 0)
{
m_requireNeg.resize(GetNedges());
for (i = 0; i < GetNedges(); ++i)
{
m_requireNeg[i] = false;
if (m_negatedNormals[i])
{
m_requireNeg[i] = true;
}
Expansion1DSharedPtr edgeExp = boost::dynamic_pointer_cast<
Expansion1D>(m_edgeExp[i].lock());
if (edgeExp->GetRightAdjacentElementExp())
{
if (edgeExp->GetRightAdjacentElementExp()->GetGeom2D()
->GetGlobalID() == GetGeom2D()->GetGlobalID())
{
m_requireNeg[i] = true;
}
}
}
}
StdRegions::Orientation edgedir = GetEorient(edge);
unsigned short num_mod0 = EdgeExp->GetBasis(0)->GetNumModes();
......@@ -126,18 +179,10 @@ namespace Nektar
boost::dynamic_pointer_cast<
LocalRegions::Expansion1D>(EdgeExp);
if (m_negatedNormals[edge])
if (m_requireNeg[edge])
{
Vmath::Neg(n_coeffs,EdgeExp->UpdateCoeffs(),1);
}
else if (locExp->GetRightAdjacentElementEdge() != -1)
{
if (locExp->GetRightAdjacentElementExp()->GetGeom2D()->GetGlobalID()
== GetGeom2D()->GetGlobalID())
{
Vmath::Neg(n_coeffs,EdgeExp->UpdateCoeffs(),1);
}
}
Array<OneD, NekDouble> coeff(n_coeffs,0.0);
LibUtilities::BasisType btype = ((LibUtilities::BasisType) 1); //1-->Ortho_A
......@@ -162,18 +207,10 @@ namespace Nektar
boost::dynamic_pointer_cast<
LocalRegions::Expansion1D>(EdgeExp);
if (m_negatedNormals[edge])
if (m_requireNeg[edge])
{
Vmath::Neg(n_coeffs,EdgeExp->UpdateCoeffs(),1);
}
else if (locExp->GetRightAdjacentElementEdge() != -1)
{
if (locExp->GetRightAdjacentElementExp()->GetGeom2D()->GetGlobalID()
== GetGeom2D()->GetGlobalID())
{
Vmath::Neg(order_e,EdgeExp->UpdateCoeffs(),1);
}
}
}
// add data to outarray if forward edge normal is outwards
......
library/LocalRegions/Expansion2D.h
View file @ 8fdd0e80
......@@ -134,6 +134,7 @@ namespace Nektar
private:
std::vector<ExpansionWeakPtr> m_edgeExp;
std::vector<bool> m_requireNeg;
Expansion3DWeakPtr m_elementLeft;
Expansion3DWeakPtr m_elementRight;
......
library/LocalRegions/Expansion3D.cpp
View file @ 8fdd0e80
......@@ -42,7 +42,7 @@ namespace Nektar
{
namespace LocalRegions
{
Expansion3D::Expansion3D(){}
Expansion3D::Expansion3D() : m_requireNeg() {}
void Expansion3D::AddHDGHelmholtzTraceTerms(
const NekDouble tau,
......@@ -914,6 +914,40 @@ namespace Nektar
StdRegions::IndexMapValuesSharedPtr map =
StdExpansion::GetIndexMap(ikey);
/*
* Coming into this routine, the velocity V will have been
* multiplied by the trace normals to give the input vector Vn. By
* convention, these normals are inwards facing for elements which
* have FaceExp as their right-adjacent face. This conditional
* statement therefore determines whether the normals must be
* negated, since the integral being performed here requires an
* outwards facing normal.
*/
if (m_requireNeg.size() == 0)
{
m_requireNeg.resize(GetNfaces());
for (i = 0; i < GetNfaces(); ++i)
{
m_requireNeg[i] = false;
if (m_negatedNormals[i])
{
m_requireNeg[i] = true;
}
Expansion2DSharedPtr faceExp = m_faceExp[i].lock();
if (faceExp->GetRightAdjacentElementExp())
{
if (faceExp->GetRightAdjacentElementExp()->GetGeom3D()->GetGlobalID()
== GetGeom3D()->GetGlobalID())
{
m_requireNeg[i] = true;
}
}
}
}
int order_e = (*map).num_elements(); // Order of the element
int n_coeffs = FaceExp->GetCoeffs().num_elements(); // Order of the trace
......@@ -924,37 +958,21 @@ namespace Nektar
else
{
FaceExp->IProductWRTBase(Fn,FaceExp->UpdateCoeffs());
LocalRegions::Expansion2DSharedPtr locExp =
boost::dynamic_pointer_cast<
LocalRegions::Expansion2D>(FaceExp);
/*
* Coming into this routine, the velocity V will have been
* multiplied by the trace normals to give the input vector
* Vn. By convention, these normals are inwards facing for
* elements which have FaceExp as their right-adjacent face.
* This conditional statement therefore determines whether the
* normals must be negated, since the integral being performed
* here requires an outwards facing normal.
*/
if (m_negatedNormals[face])
{
Vmath::Neg(order_e,FaceExp->UpdateCoeffs(),1);
}
else if (locExp->GetRightAdjacentElementFace() != -1)
}
if (m_requireNeg[face])
{
for(i = 0; i < order_e; ++i)
{
if (locExp->GetRightAdjacentElementExp()->GetGeom3D()->GetGlobalID()
== GetGeom3D()->GetGlobalID())
{
Vmath::Neg(order_e,FaceExp->UpdateCoeffs(),1);
}
outarray[(*map)[i].index] -= (*map)[i].sign*FaceExp->GetCoeff(i);
}
}
for(i = 0; i < order_e; ++i)
else
{
outarray[(*map)[i].index] += (*map)[i].sign*FaceExp->GetCoeff(i);
for(i = 0; i < order_e; ++i)
{
outarray[(*map)[i].index] += (*map)[i].sign*FaceExp->GetCoeff(i);
}
}
}
......
library/LocalRegions/Expansion3D.h
View file @ 8fdd0e80
......@@ -106,6 +106,7 @@ namespace Nektar
// Only use this class for member functions
std::vector<Expansion2DWeakPtr> m_faceExp;
std::vector<bool> m_requireNeg;
};
// type defines for use of PrismExp in a boost vector
......
library/LocalRegions/HexExp.cpp
View file @ 8fdd0e80
......@@ -368,24 +368,22 @@ namespace Nektar
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
{
const int nqtot = GetTotPoints();
Array<OneD, const NekDouble> jac = m_metricinfo->GetJac();
Array<OneD, NekDouble> tmp(nqtot);
// multiply inarray with Jacobian
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
Vmath::Vmul(nqtot, &jac[0], 1, (NekDouble*)&inarray[0], 1,
&tmp[0], 1);
}
else
{
Vmath::Smul(nqtot, jac[0], (NekDouble*)&inarray[0], 1,
&tmp[0], 1);
}
StdHexExp::v_IProductWRTBase(tmp, outarray);
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(inarray.num_elements());
Array<OneD, NekDouble> wsp(nquad0*nquad1*(nquad2+order0) +
order0*order1*nquad2);
MultiplyByQuadratureMetric(inarray, tmp);
IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
m_base[1]->GetBdata(),
m_base[2]->GetBdata(),
tmp,outarray,wsp,
true,true,true);
}
void HexExp::v_IProductWRTDerivBase(
......@@ -433,50 +431,49 @@ namespace Nektar
const Array<TwoD, const NekDouble>& gmat = m_metricinfo->GetGmat();
Array<OneD, NekDouble> alloc(3*nqtot + 2*m_ncoeffs +
Array<OneD, NekDouble> alloc(4*nqtot + 2*m_ncoeffs +
nmodes0*nquad2*(nquad1+nmodes1));
Array<OneD, NekDouble> tmp1(alloc); // Dir1 metric
Array<OneD, NekDouble> tmp2(alloc + nqtot); // Dir2 metric
Array<OneD, NekDouble> tmp3(alloc + 2*nqtot); // Dir3 metric
Array<OneD, NekDouble> tmp4(alloc + 3*nqtot); // Dir1 iprod
Array<OneD, NekDouble> tmp5(tmp4 + m_ncoeffs); // Dir2 iprod
Array<OneD, NekDouble> wsp (tmp4 + 2*m_ncoeffs); // Wsp
Array<OneD, NekDouble> tmp1(alloc); // Quad metric
Array<OneD, NekDouble> tmp2(alloc + nqtot); // Dir1 metric
Array<OneD, NekDouble> tmp3(alloc + 2*nqtot); // Dir2 metric
Array<OneD, NekDouble> tmp4(alloc + 3*nqtot); // Dir3 metric
Array<OneD, NekDouble> tmp5(alloc + 4*nqtot); // Dir1 iprod
Array<OneD, NekDouble> tmp6(tmp5 + m_ncoeffs); // Dir2 iprod
Array<OneD, NekDouble> wsp (tmp5 + 2*m_ncoeffs); // Wsp
MultiplyByQuadratureMetric(inarray, tmp1);
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
Vmath::Vmul(nqtot,&gmat[3*dir][0], 1,inarray.get(),1,tmp1.get(),1);
Vmath::Vmul(nqtot,&gmat[3*dir+1][0],1,inarray.get(),1,tmp2.get(),1);
Vmath::Vmul(nqtot,&gmat[3*dir+2][0],1,inarray.get(),1,tmp3.get(),1);
Vmath::Vmul(nqtot,&gmat[3*dir][0], 1,tmp1.get(),1,tmp2.get(),1);
Vmath::Vmul(nqtot,&gmat[3*dir+1][0],1,tmp1.get(),1,tmp3.get(),1);
Vmath::Vmul(nqtot,&gmat[3*dir+2][0],1,tmp1.get(),1,tmp4.get(),1);
}
else
{
Vmath::Smul(nqtot, gmat[3*dir][0], inarray.get(),1,tmp1.get(), 1);
Vmath::Smul(nqtot, gmat[3*dir+1][0],inarray.get(),1,tmp2.get(), 1);
Vmath::Smul(nqtot, gmat[3*dir+2][0],inarray.get(),1,tmp3.get(), 1);
Vmath::Smul(nqtot, gmat[3*dir][0], tmp1.get(),1,tmp2.get(), 1);
Vmath::Smul(nqtot, gmat[3*dir+1][0],tmp1.get(),1,tmp3.get(), 1);
Vmath::Smul(nqtot, gmat[3*dir+2][0],tmp1.get(),1,tmp4.get(), 1);
}
MultiplyByQuadratureMetric(tmp1,tmp1);
MultiplyByQuadratureMetric(tmp2,tmp2);
MultiplyByQuadratureMetric(tmp3,tmp3);
IProductWRTBase_SumFacKernel( m_base[0]->GetDbdata(),
m_base[1]->GetBdata(),
m_base[2]->GetBdata(),
tmp1,tmp4,wsp,
tmp2,tmp5,wsp,
false,true,true);
IProductWRTBase_SumFacKernel( m_base[0]->GetBdata(),
m_base[1]->GetDbdata(),
m_base[2]->GetBdata(),
tmp2,tmp5,wsp,
tmp3,tmp6,wsp,
true,false,true);
IProductWRTBase_SumFacKernel( m_base[0]->GetBdata(),
m_base[1]->GetBdata(),
m_base[2]->GetDbdata(),
tmp3,outarray,wsp,
tmp4,outarray,wsp,
true,true,false);
Vmath::Vadd(GetNcoeffs(), tmp4, 1, outarray, 1, outarray, 1);
Vmath::Vadd(GetNcoeffs(), tmp5, 1, outarray, 1, outarray, 1);
Vmath::Vadd(GetNcoeffs(), tmp6, 1, outarray, 1, outarray, 1);
}
......@@ -2297,19 +2294,18 @@ 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())
{
const Array<OneD, const NekDouble> &metric
= m_metricinfo->GetQuadratureMetrics();
Vmath::Vmul(nqtot, metric, 1, inarray, 1, outarray, 1);
Vmath::Vmul(metric.num_elements(), metric, 1, inarray, 1, outarray, 1);
}
else
{
const int nqtot = m_base[0]->GetNumPoints() *
m_base[1]->GetNumPoints() *
m_base[2]->GetNumPoints();
const Array<OneD, const NekDouble> &jac
= m_metricinfo->GetJac();
......
library/LocalRegions/PrismExp.cpp
View file @ 8fdd0e80
......@@ -279,25 +279,22 @@ namespace Nektar
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble>& outarray)
{
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nquad2 = m_base[2]->GetNumPoints();
Array<OneD, const NekDouble> jac = m_metricinfo->GetJac();
Array<OneD,NekDouble> tmp(nquad0*nquad1*nquad2);
const int nquad0 = m_base[0]->GetNumPoints();
const int nquad1 = m_base[1]->GetNumPoints();
const int nquad2 = m_base[2]->GetNumPoints();
const int order0 = m_base[0]->GetNumModes();
const int order1 = m_base[1]->GetNumModes();
// multiply inarray with Jacobian
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
Vmath::Vmul(nquad0*nquad1*nquad2,&jac[0],1,
(NekDouble*)&inarray[0],1,&tmp[0],1);
}
else
{
Vmath::Smul(nquad0*nquad1*nquad2,jac[0],
(NekDouble*)&inarray[0],1,&tmp[0],1);
}
StdPrismExp::v_IProductWRTBase(tmp,outarray);
Array<OneD, NekDouble> tmp(nquad0*nquad1*nquad2);
Array<OneD, NekDouble> wsp(order0*nquad2*(nquad1+order1));
MultiplyByQuadratureMetric(inarray, tmp);
IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
m_base[1]->GetBdata(),
m_base[2]->GetBdata(),
tmp,outarray,wsp,
true,true,true);
}
/**
......@@ -343,12 +340,12 @@ namespace Nektar
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
{
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 ();
int nqtot = nquad0*nquad1*nquad2;
const int nquad0 = m_base[0]->GetNumPoints();
const int nquad1 = m_base[1]->GetNumPoints();
const int nquad2 = m_base[2]->GetNumPoints();
const int order0 = m_base[0]->GetNumModes ();
const int order1 = m_base[1]->GetNumModes ();
const int nqtot = nquad0*nquad1*nquad2;
int i;
const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
......@@ -360,29 +357,27 @@ namespace Nektar
Array<OneD, NekDouble> tmp2 (nqtot );
Array<OneD, NekDouble> tmp3 (nqtot );
Array<OneD, NekDouble> tmp4 (nqtot );
Array<OneD, NekDouble> tmp5 (m_ncoeffs);
Array<OneD, NekDouble> tmp5 (nqtot );
Array<OneD, NekDouble> tmp6 (m_ncoeffs);
Array<OneD, NekDouble> wsp (order0*nquad2*(nquad1+order1));
const Array<TwoD, const NekDouble>& gmat = m_metricinfo->GetGmat();
MultiplyByQuadratureMetric(inarray, tmp1);
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
Vmath::Vmul(nqtot,&gmat[3*dir][0], 1,inarray.get(),1,tmp1.get(),1);
Vmath::Vmul(nqtot,&gmat[3*dir+1][0],1,inarray.get(),1,tmp2.get(),1);
Vmath::Vmul(nqtot,&gmat[3*dir+2][0],1,inarray.get(),1,tmp3.get(),1);
Vmath::Vmul(nqtot,&gmat[3*dir][0], 1,tmp1.get(),1,tmp2.get(),1);
Vmath::Vmul(nqtot,&gmat[3*dir+1][0],1,tmp1.get(),1,tmp3.get(),1);
Vmath::Vmul(nqtot,&gmat[3*dir+2][0],1,tmp1.get(),1,tmp4.get(),1);
}
else
{
Vmath::Smul(nqtot, gmat[3*dir][0], inarray.get(),1,tmp1.get(), 1);
Vmath::Smul(nqtot, gmat[3*dir+1][0],inarray.get(),1,tmp2.get(), 1);
Vmath::Smul(nqtot, gmat[3*dir+2][0],inarray.get(),1,tmp3.get(), 1);
Vmath::Smul(nqtot, gmat[3*dir][0], tmp1.get(),1,tmp2.get(), 1);
Vmath::Smul(nqtot, gmat[3*dir+1][0],tmp1.get(),1,tmp3.get(), 1);
Vmath::Smul(nqtot, gmat[3*dir+2][0],tmp1.get(),1,tmp4.get(), 1);
}
MultiplyByQuadratureMetric(tmp1,tmp1);
MultiplyByQuadratureMetric(tmp2,tmp2);
MultiplyByQuadratureMetric(tmp3,tmp3);
// set up geometric factor: (1+z0)/2
for (i = 0; i < nquad0; ++i)
{
......@@ -395,43 +390,42 @@ namespace Nektar
gfac2[i] = 2.0/(1-z2[i]);
}
const int nq01 = nquad0*nquad1;
for (i = 0; i < nquad2; ++i)
{
Vmath::Smul(nquad0*nquad1,gfac2[i],
&tmp1[0]+i*nquad0*nquad1,1,
&tmp1[0]+i*nquad0*nquad1,1);
Vmath::Smul(nquad0*nquad1,gfac2[i],
&tmp3[0]+i*nquad0*nquad1,1,
&tmp4[0]+i*nquad0*nquad1,1);
Vmath::Smul(nq01,gfac2[i],&tmp2[0]+i*nq01,1,&tmp2[0]+i*nq01,1);
Vmath::Smul(nq01,gfac2[i],&tmp4[0]+i*nq01,1,&tmp5[0]+i*nq01,1);
}
for(i = 0; i < nquad1*nquad2; ++i)
{
Vmath::Vmul(nquad0,&gfac0[0],1,&tmp4[0]+i*nquad0,1,
&tmp4[0]+i*nquad0,1);
Vmath::Vmul(nquad0,&gfac0[0],1,&tmp5[0]+i*nquad0,1,
&tmp5[0]+i*nquad0,1);
}
Vmath::Vadd(nqtot, &tmp1[0], 1, &tmp4[0], 1, &tmp1[0], 1);
Vmath::Vadd(nqtot, &tmp2[0], 1, &tmp5[0], 1, &tmp2[0], 1);
IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
m_base[1]->GetBdata (),
m_base[2]->GetBdata (),
tmp1,tmp5,wsp,
tmp2,outarray,wsp,
true,true,true);
IProductWRTBase_SumFacKernel(m_base[0]->GetBdata (),
m_base[1]->GetDbdata(),
m_base[2]->GetBdata (),
tmp2,tmp6,wsp,
tmp3,tmp6,wsp,
true,true,true);
Vmath::Vadd(m_ncoeffs, tmp6, 1, outarray, 1, outarray, 1);
IProductWRTBase_SumFacKernel(m_base[0]->GetBdata (),
m_base[1]->GetBdata (),
m_base[2]->GetDbdata(),
tmp3,outarray,wsp,
tmp4,tmp6,wsp,
true,true,true);
Vmath::Vadd(m_ncoeffs, tmp5, 1, outarray, 1, outarray, 1);
Vmath::Vadd(m_ncoeffs, tmp6, 1, outarray, 1, outarray, 1);
}
......@@ -1762,19 +1756,20 @@ 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())
{
const Array<OneD, const NekDouble> &metric
= m_metricinfo->GetQuadratureMetrics();
Vmath::Vmul(nqtot, metric, 1, inarray, 1, outarray, 1);
Vmath::Vmul(metric.num_elements(), metric, 1, inarray, 1,
outarray, 1);
}
else
{
const int nqtot = m_base[0]->GetNumPoints() *
m_base[1]->GetNumPoints() *
m_base[2]->GetNumPoints();
const Array<OneD, const NekDouble> &jac
= m_metricinfo->GetJac();
......
library/LocalRegions/TetExp.cpp
View file @ 8fdd0e80
......@@ -299,25 +299,21 @@ namespace Nektar
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
{
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nquad2 = m_base[2]->GetNumPoints();
Array<OneD, const NekDouble> jac = m_metricinfo->GetJac();
Array<OneD,NekDouble> tmp(nquad0*nquad1*nquad2);
// multiply inarray with Jacobian
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
Vmath::Vmul(nquad0*nquad1*nquad2,&jac[0],1,
(NekDouble*)&inarray[0],1,&tmp[0],1);
}
else
{
Vmath::Smul(nquad0*nquad1*nquad2,jac[0],
(NekDouble*)&inarray[0],1,&tmp[0],1);
}
StdTetExp::v_IProductWRTBase(tmp,outarray);
const int nquad0 = m_base[0]->GetNumPoints();
const int nquad1 = m_base[1]->GetNumPoints();
const int nquad2 = m_base[2]->GetNumPoints();
const int order0 = m_base[0]->GetNumModes();
const int order1 = m_base[1]->GetNumModes();