Commit 14e3e2f1 authored by Spencer Sherwin's avatar Spencer Sherwin
Browse files

Added collection routien for IProductWRTDerivBase and associated regression tests

parent 53ae6c4d
......@@ -1453,5 +1453,236 @@ OperatorKey IProductWRTDerivBase_SumFac_Prism::m_type = GetOperatorFactory().
IProductWRTDerivBase_SumFac_Prism::create,
"IProductWRTDerivBase_SumFac_Prism");
/**
* @brief Inner product WRT deriv base operator using sum-factorisation (Pyr)
*/
class IProductWRTDerivBase_SumFac_Pyr : public Operator
{
public:
OPERATOR_CREATE(IProductWRTDerivBase_SumFac_Pyr)
virtual ~IProductWRTDerivBase_SumFac_Pyr()
{
}
virtual void operator()(
const Array<OneD, const NekDouble> &entry0,
Array<OneD, NekDouble> &entry1,
Array<OneD, NekDouble> &entry2,
Array<OneD, NekDouble> &entry3,
Array<OneD, NekDouble> &wsp)
{
unsigned int nPhys = m_stdExp->GetTotPoints();
unsigned int ntot = m_numElmt*nPhys;
unsigned int nmodes = m_stdExp->GetNcoeffs();
unsigned int nmax = max(ntot,m_numElmt*nmodes);
Array<OneD, Array<OneD, const NekDouble> > in(3);
Array<OneD, NekDouble> output, wsp1;
Array<OneD, Array<OneD, NekDouble> > tmp(3);
in[0] = entry0; in[1] = entry1;
in[2] = entry2;
output = entry3;
for(int i = 0; i < 3; ++i)
{
tmp[i] = wsp + i*nmax;
}
// calculate (dphi/dx,in[0]) = ((dphi/dxi_0 dxi_0/dx +
// dphi/dxi_1 dxi_1/dx),in[0])
// + (dphi/dy,in[1]) = ((dphi/dxi_0 dxi_0/dy +
// dphi/dxi_1 dxi_1/dy),in[1])
// + (dphi/dz,in[2]) = ((dphi/dxi_0 dxi_0/dz +
// dphi/dxi_1 dxi_1/dz),in[2])
//
// Note dphi/dxi_0 =
// dphi/deta_0 deta_0/dxi_0 = dphi/deta_0 2/(1-eta_2)
//
// Note dphi/dxi_1 =
// dphi/deta_1 deta_1/dxi_1 = dphi/deta_1 2/(1-eta_2)
//
// dphi/dxi_2 =
// dphi/deta_0 deta_0/dxi_2 + dphi/deta_1 deta_1/dxi_2 +
// + dphi/deta_2 deta_2/dxi_2 =
// dphi/deta_0 (1+eta_0)/(1-eta_2) +
// dphi/deta_1 (1+eta_1)/(1-eta_2) + dphi/deta_2
//
// and so the full inner products are
//
// (dphi/dx,in[0]) + (dphi/dy,in[1]) + (dphi/dz,in[2])
// = (dphi/deta_0, ((2/(1-eta_2) (dxi_0/dx in[0] + dxi_0/dy in[1]
// + dxi_0/dz in[2])
// + (1_eta_0)/(1-eta_2) (dxi_2/dx in[0] + dxi_2/dy in[1]
// + dxi_2/dz in[2] ))
// + (dphi/deta_1, ((2/(1-eta_2) (dxi_1/dx in[0] + dxi_0/dy in[1]
// + dxi_0/dz in[2])
// + (1_eta_1)/(1-eta_2) (dxi_2/dx in[0] + dxi_2/dy in[1]
// + dxi_2/dz in[2] ))
// + (dphi/deta_2, (dxi_2/dx in[0] + dxi_2/dy in[1]
// + dxi_2/dz in[2]))
for(int i = 0; i < 3; ++i)
{
Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1, tmp[i],1);
for(int j = 1; j < 3; ++j)
{
Vmath::Vvtvp (ntot,m_derivFac[i+3*j],1,
in[j],1, tmp[i], 1, tmp[i],1);
}
}
wsp1 = wsp + 3*nmax;
// Sort into eta factors
for (int i = 0; i < m_numElmt; ++i)
{
// scale tmp[0] by fac0
Vmath::Vmul(nPhys,&m_fac0[0],1,tmp[0].get()+i*nPhys,1,
tmp[0].get()+i*nPhys,1);
// scale tmp[2] by fac1 and add to tmp0
Vmath::Vvtvp(nPhys,&m_fac1[0],1,tmp[2].get()+i*nPhys,1,
tmp[0].get()+i*nPhys,1,tmp[0].get()+i*nPhys,1);
// scale tmp[1] by fac0
Vmath::Vmul(nPhys,&m_fac0[0],1,tmp[1].get()+i*nPhys,1,
tmp[1].get()+i*nPhys,1);
// scale tmp[2] by fac2 and add to tmp1
Vmath::Vvtvp(nPhys,&m_fac2[0],1,tmp[2].get()+i*nPhys,1,
tmp[1].get()+i*nPhys,1,tmp[1].get()+i*nPhys,1);
}
// calculate Iproduct WRT Std Deriv
PyrIProduct(m_sortTopVertex, m_numElmt,
m_nquad0, m_nquad1, m_nquad2,
m_nmodes0, m_nmodes1, m_nmodes2,
m_derbase0, m_base1, m_base2,
m_jac,tmp[0],output,wsp1);
PyrIProduct(m_sortTopVertex, m_numElmt,
m_nquad0, m_nquad1, m_nquad2,
m_nmodes0, m_nmodes1, m_nmodes2,
m_base0, m_derbase1, m_base2,
m_jac,tmp[1],tmp[0],wsp1);
Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
PyrIProduct(m_sortTopVertex, m_numElmt,
m_nquad0, m_nquad1, m_nquad2,
m_nmodes0, m_nmodes1, m_nmodes2,
m_base0, m_base1, m_derbase2,
m_jac,tmp[2],tmp[0],wsp1);
Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
}
virtual void operator()(
int dir,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
ASSERTL0(false, "Not valid for this operator.");
}
protected:
const int m_nquad0;
const int m_nquad1;
const int m_nquad2;
const int m_nmodes0;
const int m_nmodes1;
const int m_nmodes2;
Array<OneD, const NekDouble> m_jac;
Array<OneD, const NekDouble> m_base0;
Array<OneD, const NekDouble> m_base1;
Array<OneD, const NekDouble> m_base2;
Array<OneD, const NekDouble> m_derbase0;
Array<OneD, const NekDouble> m_derbase1;
Array<OneD, const NekDouble> m_derbase2;
Array<TwoD, const NekDouble> m_derivFac;
Array<OneD, NekDouble> m_fac0;
Array<OneD, NekDouble> m_fac1;
Array<OneD, NekDouble> m_fac2;
bool m_sortTopVertex;
private:
IProductWRTDerivBase_SumFac_Pyr(
vector<StdRegions::StdExpansionSharedPtr> pCollExp,
CoalescedGeomDataSharedPtr pGeomData)
: Operator (pCollExp, pGeomData),
m_nquad0 (m_stdExp->GetNumPoints(0)),
m_nquad1 (m_stdExp->GetNumPoints(1)),
m_nquad2 (m_stdExp->GetNumPoints(2)),
m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
m_nmodes2 (m_stdExp->GetBasisNumModes(2)),
m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
m_base2 (m_stdExp->GetBasis(2)->GetBdata()),
m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
m_derbase1(m_stdExp->GetBasis(1)->GetDbdata()),
m_derbase2(m_stdExp->GetBasis(2)->GetDbdata())
{
m_jac = pGeomData->GetJacWithStdWeights(pCollExp);
m_wspSize = 6 * m_numElmt * (max(m_nquad0*m_nquad1*m_nquad2,
m_nmodes0*m_nmodes1*m_nmodes2));
m_derivFac = pGeomData->GetDerivFactors(pCollExp);
if(m_stdExp->GetBasis(0)->GetBasisType()
== LibUtilities::eModified_A)
{
m_sortTopVertex = true;
}
else
{
m_sortTopVertex = false;
}
const Array<OneD, const NekDouble>& z0
= m_stdExp->GetBasis(0)->GetZ();
const Array<OneD, const NekDouble>& z1
= m_stdExp->GetBasis(1)->GetZ();
const Array<OneD, const NekDouble>& z2
= m_stdExp->GetBasis(2)->GetZ();
m_fac0 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
m_fac1 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
m_fac2 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
for (int i = 0; i < m_nquad0; ++i)
{
for(int j = 0; j < m_nquad1; ++j)
{
for(int k = 0; k < m_nquad2; ++k)
{
// set up geometric factor: 2/(1-z2)
m_fac0[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
= 2.0/(1-z2[k]);
// set up geometric factor: (1+z0)/(1-z2)
m_fac1[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
= (1+z0[i])/(1-z2[k]);
// set up geometric factor: (1+z1)/(1-z2)
m_fac2[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
= (1+z1[j])/(1-z2[k]);
}
}
}
}
};
/// Factory initialisation for the IProductWRTDerivBase_SumFac_Pyr operator
OperatorKey IProductWRTDerivBase_SumFac_Pyr::m_type = GetOperatorFactory().
RegisterCreatorFunction(
OperatorKey(ePyramid, eIProductWRTDerivBase, eSumFac, false),
IProductWRTDerivBase_SumFac_Pyr::create,
"IProductWRTDerivBase_SumFac_Pyr");
}
}
......@@ -308,7 +308,6 @@ namespace Nektar
LibUtilities::SessionReaderSharedPtr dummySession;
Collections::CollectionOptimisation colOpt(dummySession, Collections::eSumFac);
//Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
Collections::OperatorImpMap impTypes;
impTypes[Collections::eBwdTrans] = Collections::eSumFac;
Collections::Collection c(CollExp, impTypes);
......@@ -497,7 +496,6 @@ namespace Nektar
LibUtilities::SessionReaderSharedPtr dummySession;
Collections::CollectionOptimisation colOpt(dummySession, Collections::eSumFac);
//Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
Collections::OperatorImpMap impTypes;
impTypes[Collections::eBwdTrans] = Collections::eSumFac;
Collections::Collection c(CollExp, impTypes);
......@@ -727,7 +725,6 @@ namespace Nektar
LibUtilities::SessionReaderSharedPtr dummySession;
Collections::CollectionOptimisation colOpt(dummySession, Collections::eSumFac);
//Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
Collections::OperatorImpMap impTypes;
impTypes[Collections::ePhysDeriv] = Collections::eSumFac;
Collections::Collection c(CollExp, impTypes);
......@@ -896,7 +893,6 @@ namespace Nektar
LibUtilities::SessionReaderSharedPtr dummySession;
Collections::CollectionOptimisation colOpt(dummySession, Collections::eSumFac);
//Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
Collections::OperatorImpMap impTypes;
impTypes[Collections::ePhysDeriv] = Collections::eSumFac;
......@@ -1138,7 +1134,6 @@ namespace Nektar
LibUtilities::SessionReaderSharedPtr dummySession;
Collections::CollectionOptimisation colOpt(dummySession, Collections::eSumFac);
//Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
Collections::OperatorImpMap impTypes;
impTypes[Collections::eIProductWRTBase] = Collections::eSumFac;
Collections::Collection c(CollExp, impTypes);
......@@ -1380,7 +1375,6 @@ namespace Nektar
LibUtilities::SessionReaderSharedPtr dummySession;
Collections::CollectionOptimisation colOpt(dummySession, Collections::eSumFac);
//Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
Collections::OperatorImpMap impTypes;
impTypes[Collections::eIProductWRTBase] = Collections::eSumFac;
Collections::Collection c(CollExp, impTypes);
......@@ -1418,7 +1412,6 @@ namespace Nektar
}
}
#if 0
BOOST_AUTO_TEST_CASE(TestPyrIProductWRTDerivBase_IterPerExp_UniformP_MultiElmt)
{
......@@ -1427,9 +1420,8 @@ namespace Nektar
SpatialDomains::PointGeomSharedPtr v2(new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
SpatialDomains::PointGeomSharedPtr v3(new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
SpatialDomains::PointGeomSharedPtr v4(new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
SpatialDomains::PointGeomSharedPtr v5(new SpatialDomains::PointGeom(3u, 5u, -1.0, 1.0, 1.0));
SpatialDomains::PyrGeomSharedPtr pyrGeom = CreatePyr(v0, v1, v2, v3, v4, v5);
SpatialDomains::PyrGeomSharedPtr pyrGeom = CreatePyr(v0, v1, v2, v3, v4);
Nektar::LibUtilities::PointsType PointsTypeDir1 = Nektar::LibUtilities::eGaussLobattoLegendre;
const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
......@@ -1520,9 +1512,8 @@ namespace Nektar
SpatialDomains::PointGeomSharedPtr v2(new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
SpatialDomains::PointGeomSharedPtr v3(new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
SpatialDomains::PointGeomSharedPtr v4(new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
SpatialDomains::PointGeomSharedPtr v5(new SpatialDomains::PointGeom(3u, 5u, -1.0, 1.0, 1.0));
SpatialDomains::PyrGeomSharedPtr pyrGeom = CreatePyr(v0, v1, v2, v3, v4, v5);
SpatialDomains::PyrGeomSharedPtr pyrGeom = CreatePyr(v0, v1, v2, v3, v4);
Nektar::LibUtilities::PointsType PointsTypeDir1 = Nektar::LibUtilities::eGaussLobattoLegendre;
const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
......@@ -1612,9 +1603,8 @@ namespace Nektar
SpatialDomains::PointGeomSharedPtr v2(new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
SpatialDomains::PointGeomSharedPtr v3(new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
SpatialDomains::PointGeomSharedPtr v4(new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
SpatialDomains::PointGeomSharedPtr v5(new SpatialDomains::PointGeom(3u, 5u, -1.0, 1.0, 1.0));
SpatialDomains::PyrGeomSharedPtr pyrGeom = CreatePyr(v0, v1, v2, v3, v4, v5);
SpatialDomains::PyrGeomSharedPtr pyrGeom = CreatePyr(v0, v1, v2, v3, v4);
Nektar::LibUtilities::PointsType PointsTypeDir1 = Nektar::LibUtilities::eGaussLobattoLegendre;
const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
......@@ -1707,9 +1697,8 @@ namespace Nektar
SpatialDomains::PointGeomSharedPtr v2(new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
SpatialDomains::PointGeomSharedPtr v3(new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
SpatialDomains::PointGeomSharedPtr v4(new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
SpatialDomains::PointGeomSharedPtr v5(new SpatialDomains::PointGeom(3u, 5u, -1.0, 1.0, 1.0));
SpatialDomains::PyrGeomSharedPtr pyrGeom = CreatePyr(v0, v1, v2, v3, v4, v5);
SpatialDomains::PyrGeomSharedPtr pyrGeom = CreatePyr(v0, v1, v2, v3, v4);
Nektar::LibUtilities::PointsType PointsTypeDir1 = Nektar::LibUtilities::eGaussLobattoLegendre;
const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
......@@ -1802,9 +1791,8 @@ namespace Nektar
SpatialDomains::PointGeomSharedPtr v2(new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
SpatialDomains::PointGeomSharedPtr v3(new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
SpatialDomains::PointGeomSharedPtr v4(new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
SpatialDomains::PointGeomSharedPtr v5(new SpatialDomains::PointGeom(3u, 5u, -1.0, 1.0, 1.0));
SpatialDomains::PyrGeomSharedPtr pyrGeom = CreatePyr(v0, v1, v2, v3, v4, v5);
SpatialDomains::PyrGeomSharedPtr pyrGeom = CreatePyr(v0, v1, v2, v3, v4);
Nektar::LibUtilities::PointsType PointsTypeDir1 = Nektar::LibUtilities::eGaussLobattoLegendre;
const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
......@@ -1896,9 +1884,8 @@ namespace Nektar
SpatialDomains::PointGeomSharedPtr v2(new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
SpatialDomains::PointGeomSharedPtr v3(new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
SpatialDomains::PointGeomSharedPtr v4(new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
SpatialDomains::PointGeomSharedPtr v5(new SpatialDomains::PointGeom(3u, 5u, -1.0, 1.0, 1.0));
SpatialDomains::PyrGeomSharedPtr pyrGeom = CreatePyr(v0, v1, v2, v3, v4, v5);
SpatialDomains::PyrGeomSharedPtr pyrGeom = CreatePyr(v0, v1, v2, v3, v4);
Nektar::LibUtilities::PointsType PointsTypeDir1 = Nektar::LibUtilities::eGaussLobattoLegendre;
const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
......@@ -1982,7 +1969,5 @@ namespace Nektar
BOOST_CHECK_CLOSE(coeffs1[i],coeffs2[i], epsilon);
}
}
#endif
}
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment