Commit ada521ac authored by Spencer Sherwin's avatar Spencer Sherwin
Browse files

Updated with mapping array for GetFwdBwdTracePhys and GetTracePhys. Currently...

Updated with mapping array for GetFwdBwdTracePhys and GetTracePhys. Currently have left timing in IncNS and also do not think release mode is executing correctly
parent 267a6f53
......@@ -58,7 +58,7 @@ int main(int argc, char *argv[])
//BwdTrans comparison
cout << "BwdTrans Op: Ntest = " << Ntest << endl;
for(int N = 2; N < 11; ++N)
for(int N = 2; N < 6; ++N)
{
graph3D->SetExpansionsToPolyOrder(N);
......@@ -103,7 +103,7 @@ int main(int argc, char *argv[])
// IProductWRTBase Comparison
cout << "IProductWRTBase Op: Ntest = " << Ntest << endl;
for(int N = 2; N < 11; ++N)
for(int N = 2; N < 6; ++N)
{
graph3D->SetExpansionsToPolyOrder(N);
......@@ -155,7 +155,7 @@ int main(int argc, char *argv[])
// PhysDeriv Comparison
cout << "PhysDeriv Op: Ntest = " << Ntest << endl;
for(int N = 2; N < 11; ++N)
for(int N = 2; N < 6; ++N)
{
graph3D->SetExpansionsToPolyOrder(N);
......@@ -206,7 +206,7 @@ int main(int argc, char *argv[])
// IProductWRTDerivBase Comparison
cout << "IProductWRTDerivBase Op: Ntest = " << Ntest << endl;
for(int N = 2; N < 11; ++N)
for(int N = 2; N < 6; ++N)
{
graph3D->SetExpansionsToPolyOrder(N);
......
......@@ -141,6 +141,14 @@ namespace Nektar
const PointsKey &tpoints1,
NekDouble *to)
{
// default interpolation
if((fpoints0 == tpoints0)&&(fpoints1 == tpoints1))
{
Vmath::Vcopy(tpoints0.GetNumPoints()*tpoints1.GetNumPoints(),
from,1,to,1);
return;
}
DNekMatSharedPtr I0,I1;
Array<OneD, NekDouble> wsp(tpoints1.GetNumPoints()*fpoints0.GetNumPoints()); // fnp0*tnp1
......
This diff is collapsed.
......@@ -97,6 +97,11 @@ namespace Nektar
inline SpatialDomains::Geometry3DSharedPtr GetGeom3D() const;
LOCAL_REGIONS_EXPORT void ReOrientFacePhysMap(const int nvert,
const StdRegions::Orientation orient,
const int nq0,
const int nq1,
Array<OneD, int> &idmap);
protected:
virtual void v_DGDeriv(
const int dir,
......@@ -119,6 +124,13 @@ namespace Nektar
virtual NekDouble v_Integrate(
const Array<OneD, const NekDouble>& inarray);
virtual void v_GetFacePhysVals(
const int face,
const StdRegions::StdExpansionSharedPtr &FaceExp,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
StdRegions::Orientation orient);
//-----------------------------
// Low Energy Basis functions
//-----------------------------
......@@ -139,6 +151,17 @@ namespace Nektar
LOCAL_REGIONS_EXPORT virtual DNekMatSharedPtr v_BuildVertexMatrix(
const DNekScalMatSharedPtr &r_bnd);
LOCAL_REGIONS_EXPORT void ReOrientTriFacePhysMap(const StdRegions::Orientation orient,
const int nq0,
const int nq1,
Array<OneD, int> &idmap);
LOCAL_REGIONS_EXPORT void ReOrientQuadFacePhysMap(const StdRegions::Orientation orient,
const int nq0,
const int nq1,
Array<OneD, int> &idmap);
private:
// Do not add members here since it may lead to conflicts.
// Only use this class for member functions
......
......@@ -689,6 +689,7 @@ namespace Nektar
}
#if 0
///Returns the physical values at the quadrature points of a face
void HexExp::v_GetFacePhysVals(
const int face,
......@@ -1178,8 +1179,128 @@ namespace Nektar
break;
}
}
#endif
void HexExp::v_GetFacePhysMap(const int face,
Array<OneD, int> &outarray)
{
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nquad2 = m_base[2]->GetNumPoints();
int nq0 = 0;
int nq1 = 0;
switch(face)
{
case 0:
nq0 = nquad0;
nq1 = nquad1;
//Directions A and B positive
if(outarray.num_elements()!=nq0*nq1)
{
outarray = Array<OneD, int>(nq0*nq1);
}
for (int i = 0; i < nquad0*nquad1; ++i)
{
outarray[i] = i;
}
break;
case 1:
nq0 = nquad0;
nq1 = nquad2;
//Direction A and B positive
if(outarray.num_elements()!=nq0*nq1)
{
outarray = Array<OneD, int>(nq0*nq1);
}
//Direction A and B positive
for (int k = 0; k < nquad2; k++)
{
for(int i = 0; i < nquad0; ++i)
{
outarray[k*nquad0 + i] = nquad0*nquad1*k + i;
}
}
break;
case 2:
nq0 = nquad1;
nq1 = nquad2;
//Direction A and B positive
if(outarray.num_elements()!=nq0*nq1)
{
outarray = Array<OneD, int>(nq0*nq1);
}
for (int i = 0; i < nquad1*nquad2; i++)
{
outarray[i] = nquad0-1 + i*nquad0;
}
break;
case 3:
nq0 = nquad0;
nq1 = nquad2;
//Direction A and B positive
if(outarray.num_elements()!=nq0*nq1)
{
outarray = Array<OneD, int>(nq0*nq1);
}
for (int k = 0; k < nquad2; k++)
{
for (int i = 0; i < nquad0; i++)
{
outarray[k*nquad0 + i] = (nquad0*(nquad1-1))+(k*nquad0*nquad1) + i;
}
}
break;
case 4:
nq0 = nquad1;
nq1 = nquad2;
//Direction A and B positive
if(outarray.num_elements()!=nq0*nq1)
{
outarray = Array<OneD, int>(nq0*nq1);
}
for (int i = 0; i < nquad1*nquad2; i++)
{
outarray[i] = i*nquad0;
}
break;
case 5:
nq0 = nquad0;
nq1 = nquad1;
//Directions A and B positive
if(outarray.num_elements()!=nq0*nq1)
{
outarray = Array<OneD, int>(nq0*nq1);
}
for (int i = 0; i < nquad0*nquad1; i++)
{
outarray[i] = nquad0*nquad1*(nquad2-1) + i;
}
break;
default:
ASSERTL0(false,"face value (> 5) is out of range");
break;
}
}
void HexExp::v_ComputeFaceNormal(const int face)
{
int i;
......
......@@ -168,13 +168,14 @@ namespace Nektar
LOCAL_REGIONS_EXPORT virtual
bool v_GetFaceDGForwards(const int i) const;
#if 0
LOCAL_REGIONS_EXPORT virtual void v_GetFacePhysVals(
const int face,
const StdRegions::StdExpansionSharedPtr &FaceExp,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
StdRegions::Orientation orient);
#endif
LOCAL_REGIONS_EXPORT virtual void v_GetTracePhysVals(
const int face,
......@@ -183,6 +184,10 @@ namespace Nektar
Array<OneD, NekDouble> &outarray,
StdRegions::Orientation orient);
LOCAL_REGIONS_EXPORT virtual void v_GetFacePhysMap(
const int face,
Array<OneD, int> &outarray);
LOCAL_REGIONS_EXPORT void v_ComputeFaceNormal(const int face);
//---------------------------------------
......
......@@ -588,6 +588,7 @@ namespace Nektar
v_GetFacePhysVals(face,FaceExp,inarray,outarray,orient);
}
#if 0 // have put general method in Expandion3D
void PrismExp::v_GetFacePhysVals(
const int face,
const StdRegions::StdExpansionSharedPtr &FaceExp,
......@@ -878,6 +879,110 @@ namespace Nektar
break;
}
}
#endif
void PrismExp::v_GetFacePhysMap(const int face,
Array<OneD, int> &outarray)
{
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nquad2 = m_base[2]->GetNumPoints();
int nq0 = 0;
int nq1 = 0;
switch(face)
{
case 0:
nq0 = nquad0;
nq1 = nquad1;
if(outarray.num_elements()!=nq0*nq1)
{
outarray = Array<OneD, int>(nq0*nq1);
}
//Directions A and B positive
for(int i = 0; i < nquad0*nquad1; ++i)
{
outarray[i] = i;
}
break;
case 1:
nq0 = nquad0;
nq1 = nquad2;
if(outarray.num_elements()!=nq0*nq1)
{
outarray = Array<OneD, int>(nq0*nq1);
}
//Direction A and B positive
for (int k=0; k<nquad2; k++)
{
for(int i = 0; i < nquad0; ++i)
{
outarray[k*nquad0+i] = (nquad0*nquad1*k)+i;
}
}
break;
case 2:
nq0 = nquad1;
nq1 = nquad2;
if(outarray.num_elements()!=nq0*nq1)
{
outarray = Array<OneD, int>(nq0*nq1);
}
//Directions A and B positive
for(int j = 0; j < nquad1*nquad2; ++j)
{
outarray[j] = nquad0-1 + j*nquad0;
}
break;
case 3:
nq0 = nquad0;
nq1 = nquad2;
if(outarray.num_elements()!=nq0*nq1)
{
outarray = Array<OneD, int>(nq0*nq1);
}
//Direction A and B positive
for (int k=0; k<nquad2; k++)
{
for(int i = 0; i < nquad0; ++i)
{
outarray[k*nquad0+i] = nquad0*(nquad1-1) + (nquad0*nquad1*k)+i;
}
}
break;
case 4:
nq0 = nquad1;
nq1 = nquad2;
if(outarray.num_elements()!=nq0*nq1)
{
outarray = Array<OneD, int>(nq0*nq1);
}
//Directions A and B positive
for(int j = 0; j < nquad1*nquad2; ++j)
{
outarray[j] = j*nquad0;
}
break;
default:
ASSERTL0(false,"face value (> 4) is out of range");
break;
}
}
void PrismExp::v_ComputeFaceNormal(const int face)
{
......
......@@ -140,12 +140,14 @@ namespace Nektar
NekDouble * coeffs);
LOCAL_REGIONS_EXPORT virtual
StdRegions::Orientation v_GetFaceOrient(int face);
#if 0
LOCAL_REGIONS_EXPORT virtual void v_GetFacePhysVals(
const int face,
const StdRegions::StdExpansionSharedPtr &FaceExp,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
StdRegions::Orientation orient);
#endif
LOCAL_REGIONS_EXPORT virtual void v_GetTracePhysVals(
const int face,
......@@ -154,6 +156,10 @@ namespace Nektar
Array<OneD, NekDouble> &outarray,
StdRegions::Orientation orient);
LOCAL_REGIONS_EXPORT virtual void v_GetFacePhysMap(
const int face,
Array<OneD, int> &outarray);
LOCAL_REGIONS_EXPORT void v_ComputeFaceNormal(const int face);
//---------------------------------------
......
......@@ -369,7 +369,7 @@ namespace Nektar
{
v_GetFacePhysVals(face,FaceExp,inarray,outarray,orient);
}
#if 0
void PyrExp::v_GetFacePhysVals(
const int face,
const StdRegions::StdExpansionSharedPtr &FaceExp,
......@@ -516,6 +516,107 @@ namespace Nektar
}
}
}
#endif
void PyrExp::v_GetFacePhysMap(const int face,
Array<OneD, int> &outarray)
{
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nquad2 = m_base[2]->GetNumPoints();
int nq0 = 0;
int nq1 = 0;
switch(face)
{
case 0:
nq0 = nquad0;
nq1 = nquad1;
if(outarray.num_elements()!=nq0*nq1)
{
outarray = Array<OneD, int>(nq0*nq1);
}
//Directions A and B positive
for(int i = 0; i < nquad0*nquad1; ++i)
{
outarray[i] = i;
}
break;
case 1:
nq0 = nquad0;
nq1 = nquad2;
if(outarray.num_elements()!=nq0*nq1)
{
outarray = Array<OneD, int>(nq0*nq1);
}
//Direction A and B positive
for (int k=0; k<nquad2; k++)
{
for(int i = 0; i < nquad0; ++i)
{
outarray[k*nquad0+i] = (nquad0*nquad1*k)+i;
}
}
break;
case 2:
nq0 = nquad1;
nq1 = nquad2;
if(outarray.num_elements()!=nq0*nq1)
{
outarray = Array<OneD, int>(nq0*nq1);
}
//Directions A and B positive
for(int j = 0; j < nquad1*nquad2; ++j)
{
outarray[j] = nquad0-1 + j*nquad0;
}
break;
case 3:
nq0 = nquad0;
nq1 = nquad2;
if(outarray.num_elements()!=nq0*nq1)
{
outarray = Array<OneD, int>(nq0*nq1);
}
//Direction A and B positive
for (int k=0; k<nquad2; k++)
{
for(int i = 0; i < nquad0; ++i)
{
outarray[k*nquad0+i] = nquad0*(nquad1-1) + (nquad0*nquad1*k)+i;
}
}
case 4:
nq0 = nquad1;
nq1 = nquad2;
if(outarray.num_elements()!=nq0*nq1)
{
outarray = Array<OneD, int>(nq0*nq1);
}
//Directions A and B positive
for(int j = 0; j < nquad1*nquad2; ++j)
{
outarray[j] = j*nquad0;
}
break;
default:
ASSERTL0(false,"face value (> 4) is out of range");
break;
}
}
void PyrExp::v_ComputeFaceNormal(const int face)
{
......
......@@ -122,13 +122,15 @@ namespace Nektar
LOCAL_REGIONS_EXPORT virtual
StdRegions::Orientation v_GetFaceOrient(int face);
LOCAL_REGIONS_EXPORT void v_ComputeFaceNormal(const int face);
#if 0
LOCAL_REGIONS_EXPORT virtual void v_GetFacePhysVals(
const int face,
const StdRegions::StdExpansionSharedPtr &FaceExp,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
StdRegions::Orientation orient);
#endif
LOCAL_REGIONS_EXPORT virtual void v_GetTracePhysVals(
const int face,
const StdRegions::StdExpansionSharedPtr &FaceExp,
......@@ -136,6 +138,10 @@ namespace Nektar
Array<OneD, NekDouble> &outarray,
StdRegions::Orientation orient);
LOCAL_REGIONS_EXPORT virtual void v_GetFacePhysMap(
const int face,
Array<OneD, int> &outarray);
//---------------------------------------
// Matrix creation functions
//---------------------------------------
......
......@@ -66,11 +66,11 @@ namespace Nektar
Expansion (geom),
Expansion3D (geom),
m_matrixManager(
boost::bind(&TetExp::CreateMatrix, this, _1),
std::string("TetExpMatrix")),
boost::bind(&TetExp::CreateMatrix, this, _1),
std::string("TetExpMatrix")),
m_staticCondMatrixManager(
boost::bind(&TetExp::CreateStaticCondMatrix, this, _1),
std::string("TetExpStaticCondMatrix"))
boost::bind(&TetExp::CreateStaticCondMatrix, this, _1),
std::string("TetExpStaticCondMatrix"))
{
}
......@@ -112,7 +112,7 @@ namespace Nektar
* point.
*/
NekDouble TetExp::v_Integral(
const Array<OneD, const NekDouble> &inarray)
const Array<OneD, const NekDouble> &inarray)
{
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
......@@ -153,16 +153,16 @@ namespace Nektar
* @param out_d2 Derivative in third coordinate direction.
*/
void TetExp::v_PhysDeriv(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &out_d0,
Array<OneD, NekDouble> &out_d1,
Array<OneD, NekDouble> &out_d2)
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &out_d0,
Array<OneD, NekDouble> &out_d1,
Array<OneD, NekDouble> &out_d2)
{
int TotPts = m_base[0]->GetNumPoints()*m_base[1]->GetNumPoints()*
m_base[2]->GetNumPoints();
Array<TwoD, const NekDouble> df =
m_metricinfo->GetDerivFactors(GetPointsKeys());
m_metricinfo->GetDerivFactors(GetPointsKeys());
Array<OneD,NekDouble> Diff0 = Array<OneD,NekDouble>(3*TotPts);
Array<OneD,NekDouble> Diff1 = Diff0 + TotPts;
Array<OneD,NekDouble> Diff2 = Diff1 + TotPts;
......@@ -231,8 +231,8 @@ namespace Nektar