Commit 2bde118c authored by Spencer Sherwin's avatar Spencer Sherwin
Browse files

Added IproductWRTBase into Collections and associated UnitTests

parent 74d54580
......@@ -422,6 +422,101 @@ void PrismIProduct(bool sortTopVertex, int numElmt,
}
/**
*
*/
void PyrIProduct(bool sortTopVertex, int numElmt,
int nquad0, int nquad1, int nquad2,
int nmodes0, int nmodes1, int nmodes2,
const Array<OneD, const NekDouble> &base0,
const Array<OneD, const NekDouble> &base1,
const Array<OneD, const NekDouble> &base2,
const Array<OneD, const NekDouble> &jac,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp)
{
int totmodes = LibUtilities::StdPyrData::getNumberOfCoefficients(
nmodes0,nmodes1,nmodes2);
int totpoints = nquad0*nquad1*nquad2;
int cnt;
int mode, mode1;
ASSERTL1(wsp.num_elements() >= numElmt*(nquad1*nquad2*nmodes0 +
nquad2*max(nquad0*nquad1,nmodes0*nmodes1)),
"Insufficient workspace size");
Vmath::Vmul(numElmt*totpoints,jac,1,input,1,wsp,1);
Array<OneD, NekDouble> wsp1 = wsp + numElmt * nquad2
* (max(nquad0*nquad1,
nmodes0*nmodes1));
// Perform iproduct with respect to the '0' direction
Blas::Dgemm('T', 'N', nquad1*nquad2*numElmt, nmodes0, nquad0,
1.0, wsp.get(), nquad0, base0.get(),
nquad0, 0.0, wsp1.get(), nquad1*nquad2*numElmt);
// Inner product with respect to the '1' direction
mode = 0;
for(int i=0; i < nmodes0; ++i)
{
Blas::Dgemm('T', 'N', nquad2*numElmt, nmodes1, nquad1,
1.0, wsp1.get()+ i*nquad1*nquad2*numElmt, nquad1,
base1.get(), nquad1,
0.0, wsp.get() + mode*nquad2*numElmt,nquad2*numElmt);
mode += nmodes1;
}
// Inner product with respect to the '2' direction
mode = mode1 = cnt = 0;
for(int i = 0; i < nmodes0; ++i)
{
for(int j = 0; j < nmodes1; ++j, ++cnt)
{
int ijmax = max(i,j);
Blas::Dgemm('T', 'N', nmodes2-ijmax, numElmt, nquad2,
1.0, base2.get()+mode*nquad2, nquad2,
wsp.get()+cnt*nquad2*numElmt, nquad2,
0.0, output.get()+mode1, totmodes);
mode += nmodes2-ijmax;
mode1 += nmodes2-ijmax;
}
//increment mode in case order1!=order2
for(int j = nmodes1; j < nmodes2; ++j)
{
int ijmax = max(i,j);
mode += nmodes2-ijmax;
}
}
// fix for modified basis for top singular vertex component
// Already have evaluated (1+c)/2 (1-b)/2 (1-a)/2
if(sortTopVertex)
{
for(int n = 0; n < numElmt; ++n)
{
// add in (1+c)/2 (1+b)/2 component
output[1+n*totmodes] += Blas::Ddot(nquad2,
base2.get()+nquad2,1,
&wsp[nquad2*numElmt + n*nquad2],1);
// add in (1+c)/2 (1-b)/2 (1+a)/2 component
output[1+n*totmodes] += Blas::Ddot(nquad2,
base2.get()+nquad2,1,
&wsp[nquad2*nmodes1*numElmt+n*nquad2],1);
// add in (1+c)/2 (1+b)/2 (1+a)/2 component
output[1+n*totmodes] += Blas::Ddot(nquad2,
base2.get()+nquad2,1,
&wsp[nquad2*(nmodes1+1)*numElmt+n*nquad2],1);
}
}
}
/**
*
*/
......
......@@ -82,6 +82,19 @@ void PrismIProduct(bool sortTopVert, int numElmt,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp);
void PyrIProduct(bool sortTopVert, int numElmt,
int nquad0, int nquad1, int nquad2,
int nmodes0, int nmodes1, int nmodes2,
const Array<OneD, const NekDouble> &base0,
const Array<OneD, const NekDouble> &base1,
const Array<OneD, const NekDouble> &base2,
const Array<OneD, const NekDouble> &jac,
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &wsp);
void TetIProduct(bool sortTopEdge, int numElmt,
int nquad0, int nquad1, int nquad2,
int nmodes0, int nmodes1, int nmodes2,
......
......@@ -856,5 +856,99 @@ OperatorKey IProductWRTBase_SumFac_Prism::m_type = GetOperatorFactory().
OperatorKey(ePrism, eIProductWRTBase, eSumFac,false),
IProductWRTBase_SumFac_Prism::create, "IProductWRTBase_SumFac_Prism");
/**
* @brief Backward transform operator using sum-factorisation (Pyr)
*/
class IProductWRTBase_SumFac_Pyr : public Operator
{
public:
OPERATOR_CREATE(IProductWRTBase_SumFac_Pyr)
virtual ~IProductWRTBase_SumFac_Pyr()
{
}
virtual void operator()(
const Array<OneD, const NekDouble> &input,
Array<OneD, NekDouble> &output,
Array<OneD, NekDouble> &output1,
Array<OneD, NekDouble> &output2,
Array<OneD, NekDouble> &wsp)
{
ASSERTL1(wsp.num_elements() == m_wspSize,
"Incorrect workspace size");
PyrIProduct(m_sortTopVertex, m_numElmt,
m_nquad0, m_nquad1, m_nquad2,
m_nmodes0, m_nmodes1, m_nmodes2,
m_base0, m_base1, m_base2,
m_jac,input,output,wsp);
}
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;
bool m_sortTopVertex;
private:
IProductWRTBase_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_jac = pGeomData->GetJacWithStdWeights(pCollExp);
m_wspSize = m_numElmt * m_nquad2
*(max(m_nquad0*m_nquad1,m_nmodes0*m_nmodes1))
+ m_nquad1*m_nquad2*m_numElmt*m_nmodes0;
if(m_stdExp->GetBasis(0)->GetBasisType()
== LibUtilities::eModified_A)
{
m_sortTopVertex = true;
}
else
{
m_sortTopVertex = false;
}
}
};
/// Factory initialisation for the IProductWRTBase_SumFac_Pyr operator
OperatorKey IProductWRTBase_SumFac_Pyr::m_type = GetOperatorFactory().
RegisterCreatorFunction(
OperatorKey(ePyramid, eIProductWRTBase, eSumFac,false),
IProductWRTBase_SumFac_Pyr::create, "IProductWRTBase_SumFac_Pyr");
}
}
......@@ -493,6 +493,10 @@ namespace Nektar
int order1 = m_base[1]->GetNumModes();
int order2 = m_base[2]->GetNumModes();
ASSERTL1(wsp.num_elements() >= nquad1*nquad2*order0 +
nquad2*order0*order1,
"Insufficient workspace size");
Array<OneD, NekDouble > tmp1 = wsp;
Array<OneD, NekDouble > tmp2 = wsp + nquad1*nquad2*order0;
......
......@@ -935,7 +935,6 @@ namespace Nektar
BOOST_CHECK_CLOSE(diff1[i],diff2[i], epsilon);
}
}
#if 0
BOOST_AUTO_TEST_CASE(TestPyrIProductWRTBase_IterPerExp_UniformP_MultiElmt)
{
......@@ -944,9 +943,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);
......@@ -958,7 +956,7 @@ namespace Nektar
Nektar::LibUtilities::BasisType basisTypeDir2 = Nektar::LibUtilities::eModified_A;
const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2,4,PointsKeyDir2);
Nektar::LibUtilities::PointsType PointsTypeDir3 = Nektar::LibUtilities::eGaussLobattoLegendre;
Nektar::LibUtilities::PointsType PointsTypeDir3 = Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
const Nektar::LibUtilities::PointsKey PointsKeyDir3(4, PointsTypeDir3);
Nektar::LibUtilities::BasisType basisTypeDir3 = Nektar::LibUtilities::eModifiedPyr_C;
const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3,4,PointsKeyDir3);
......@@ -1024,9 +1022,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);
......@@ -1038,7 +1035,7 @@ namespace Nektar
Nektar::LibUtilities::BasisType basisTypeDir2 = Nektar::LibUtilities::eModified_A;
const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2,4,PointsKeyDir2);
Nektar::LibUtilities::PointsType PointsTypeDir3 = Nektar::LibUtilities::eGaussLobattoLegendre;
Nektar::LibUtilities::PointsType PointsTypeDir3 = Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
const Nektar::LibUtilities::PointsKey PointsKeyDir3(4, PointsTypeDir3);
Nektar::LibUtilities::BasisType basisTypeDir3 = Nektar::LibUtilities::eModifiedPyr_C;
const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3,4,PointsKeyDir3);
......@@ -1105,9 +1102,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);
......@@ -1119,7 +1115,7 @@ namespace Nektar
Nektar::LibUtilities::BasisType basisTypeDir2 = Nektar::LibUtilities::eModified_A;
const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2,4,PointsKeyDir2);
Nektar::LibUtilities::PointsType PointsTypeDir3 = Nektar::LibUtilities::eGaussLobattoLegendre;
Nektar::LibUtilities::PointsType PointsTypeDir3 = Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
const Nektar::LibUtilities::PointsKey PointsKeyDir3(4, PointsTypeDir3);
Nektar::LibUtilities::BasisType basisTypeDir3 = Nektar::LibUtilities::eModifiedPyr_C;
const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3,4,PointsKeyDir3);
......@@ -1142,7 +1138,9 @@ namespace Nektar
LibUtilities::SessionReaderSharedPtr dummySession;
Collections::CollectionOptimisation colOpt(dummySession, Collections::eSumFac);
Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
//Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
Collections::OperatorImpMap impTypes;
impTypes[Collections::eIProductWRTBase] = Collections::eSumFac;
Collections::Collection c(CollExp, impTypes);
const int nq = Exp->GetTotPoints();
......@@ -1186,9 +1184,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);
......@@ -1200,7 +1197,7 @@ namespace Nektar
Nektar::LibUtilities::BasisType basisTypeDir2 = Nektar::LibUtilities::eModified_A;
const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2,6,PointsKeyDir2);
Nektar::LibUtilities::PointsType PointsTypeDir3 = Nektar::LibUtilities::eGaussLobattoLegendre;
Nektar::LibUtilities::PointsType PointsTypeDir3 = Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
const Nektar::LibUtilities::PointsKey PointsKeyDir3(8, PointsTypeDir3);
Nektar::LibUtilities::BasisType basisTypeDir3 = Nektar::LibUtilities::eModifiedPyr_C;
const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3,8,PointsKeyDir3);
......@@ -1267,9 +1264,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);
......@@ -1281,7 +1277,7 @@ namespace Nektar
Nektar::LibUtilities::BasisType basisTypeDir2 = Nektar::LibUtilities::eModified_A;
const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2,6,PointsKeyDir2);
Nektar::LibUtilities::PointsType PointsTypeDir3 = Nektar::LibUtilities::eGaussLobattoLegendre;
Nektar::LibUtilities::PointsType PointsTypeDir3 = Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
const Nektar::LibUtilities::PointsKey PointsKeyDir3(8, PointsTypeDir3);
Nektar::LibUtilities::BasisType basisTypeDir3 = Nektar::LibUtilities::eModifiedPyr_C;
const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3,8,PointsKeyDir3);
......@@ -1348,9 +1344,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);
......@@ -1362,7 +1357,7 @@ namespace Nektar
Nektar::LibUtilities::BasisType basisTypeDir2 = Nektar::LibUtilities::eModified_A;
const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2,6,PointsKeyDir2);
Nektar::LibUtilities::PointsType PointsTypeDir3 = Nektar::LibUtilities::eGaussLobattoLegendre;
Nektar::LibUtilities::PointsType PointsTypeDir3 = Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
const Nektar::LibUtilities::PointsKey PointsKeyDir3(8, PointsTypeDir3);
Nektar::LibUtilities::BasisType basisTypeDir3 = Nektar::LibUtilities::eModifiedPyr_C;
const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3,8,PointsKeyDir3);
......@@ -1385,7 +1380,9 @@ namespace Nektar
LibUtilities::SessionReaderSharedPtr dummySession;
Collections::CollectionOptimisation colOpt(dummySession, Collections::eSumFac);
Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
//Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
Collections::OperatorImpMap impTypes;
impTypes[Collections::eIProductWRTBase] = Collections::eSumFac;
Collections::Collection c(CollExp, impTypes);
const int nq = Exp->GetTotPoints();
......@@ -1421,6 +1418,7 @@ namespace Nektar
}
}
#if 0
BOOST_AUTO_TEST_CASE(TestPyrIProductWRTDerivBase_IterPerExp_UniformP_MultiElmt)
{
......
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