Commit 6d54669c authored by Dave Moxey's avatar Dave Moxey
Browse files

Remove the extra basis generation (seems to not be required)

parent 9b6e4c1b
...@@ -72,19 +72,6 @@ namespace Nektar ...@@ -72,19 +72,6 @@ namespace Nektar
"than order in 'c' direction"); "than order in 'c' direction");
} }
// Generate extra modes for interior
LibUtilities::BasisKey Bamod(
Ba.GetBasisType(),
Ba.GetNumModes() + 1,
Ba.GetPointsKey());
LibUtilities::BasisKey Bcmod(
Bc.GetBasisType(),
Bc.GetNumModes() + 1,
Bc.GetPointsKey());
m_base_A = LibUtilities::BasisManager()[Bamod];
m_base_C = LibUtilities::BasisManager()[Bcmod];
// Set up mode mapping which takes 0\leq i\leq N_coeffs -> (p,q,r) // Set up mode mapping which takes 0\leq i\leq N_coeffs -> (p,q,r)
// of the 3D tensor product // of the 3D tensor product
const int P = Ba.GetNumModes() - 1; const int P = Ba.GetNumModes() - 1;
...@@ -93,70 +80,70 @@ namespace Nektar ...@@ -93,70 +80,70 @@ namespace Nektar
int cnt = 0; int cnt = 0;
// Vertices // Vertices
m_map [cnt ] = triple(0, 0, 0, false); m_map [cnt ] = triple(0, 0, 0);
m_rmap[cnt++] = 0; m_rmap[cnt++] = 0;
m_map [cnt ] = triple(1, 0, 0, false); m_map [cnt ] = triple(1, 0, 0);
m_rmap[cnt++] = 0; m_rmap[cnt++] = 0;
m_map [cnt ] = triple(1, 1, 0, false); m_map [cnt ] = triple(1, 1, 0);
m_rmap[cnt++] = 0; m_rmap[cnt++] = 0;
m_map [cnt ] = triple(0, 1, 0, false); m_map [cnt ] = triple(0, 1, 0);
m_rmap[cnt++] = 0; m_rmap[cnt++] = 0;
m_map [cnt ] = triple(0, 0, 1, false); m_map [cnt ] = triple(0, 0, 1);
m_rmap[cnt++] = 1; m_rmap[cnt++] = 1;
// Edge 0 // Edge 0
for (int i = 2; i <= P; ++i) for (int i = 2; i <= P; ++i)
{ {
m_map [cnt ] = triple (i, 0, 0, false); m_map [cnt ] = triple (i, 0, 0);
m_rmap[cnt++] = GetTetMode(i, 0, 0); m_rmap[cnt++] = GetTetMode(i, 0, 0);
} }
// Edge 1 // Edge 1
for (int i = 2; i <= Q; ++i) for (int i = 2; i <= Q; ++i)
{ {
m_map [cnt ] = triple (1, i, 0, false); m_map [cnt ] = triple (1, i, 0);
m_rmap[cnt++] = GetTetMode(0, i, 0); m_rmap[cnt++] = GetTetMode(0, i, 0);
} }
// Edge 2 // Edge 2
for (int i = 2; i <= P; ++i) for (int i = 2; i <= P; ++i)
{ {
m_map [cnt ] = triple (i, 1, 0, false); m_map [cnt ] = triple (i, 1, 0);
m_rmap[cnt++] = GetTetMode(i, 0, 0); m_rmap[cnt++] = GetTetMode(i, 0, 0);
} }
// Edge 3 // Edge 3
for (int i = 2; i <= Q; ++i) for (int i = 2; i <= Q; ++i)
{ {
m_map [cnt ] = triple (0, i, 0, false); m_map [cnt ] = triple (0, i, 0);
m_rmap[cnt++] = GetTetMode(0, i, 0); m_rmap[cnt++] = GetTetMode(0, i, 0);
} }
// Edge 4 // Edge 4
for (int i = 2; i <= R; ++i) for (int i = 2; i <= R; ++i)
{ {
m_map [cnt ] = triple(0, 0, i, false); m_map [cnt ] = triple(0, 0, i);
m_rmap[cnt++] = i; m_rmap[cnt++] = i;
} }
// Edge 5 // Edge 5
for (int i = 2; i <= R; ++i) for (int i = 2; i <= R; ++i)
{ {
m_map [cnt ] = triple(1, 0, i, false); m_map [cnt ] = triple(1, 0, i);
m_rmap[cnt++] = i; m_rmap[cnt++] = i;
} }
// Edge 6 // Edge 6
for (int i = 2; i <= R; ++i) for (int i = 2; i <= R; ++i)
{ {
m_map [cnt ] = triple(1, 1, i, false); m_map [cnt ] = triple(1, 1, i);
m_rmap[cnt++] = i; m_rmap[cnt++] = i;
} }
// Edge 7 // Edge 7
for (int i = 2; i <= R; ++i) for (int i = 2; i <= R; ++i)
{ {
m_map [cnt ] = triple(0, 1, i, false); m_map [cnt ] = triple(0, 1, i);
m_rmap[cnt++] = i; m_rmap[cnt++] = i;
} }
...@@ -165,7 +152,7 @@ namespace Nektar ...@@ -165,7 +152,7 @@ namespace Nektar
{ {
for (int i = 2; i <= P; ++i) for (int i = 2; i <= P; ++i)
{ {
m_map [cnt ] = triple (i, j, 0, false); m_map [cnt ] = triple (i, j, 0);
m_rmap[cnt++] = GetTetMode(2, 0, 0); m_rmap[cnt++] = GetTetMode(2, 0, 0);
} }
} }
...@@ -175,7 +162,7 @@ namespace Nektar ...@@ -175,7 +162,7 @@ namespace Nektar
{ {
for (int j = 1; j <= R-i; ++j) for (int j = 1; j <= R-i; ++j)
{ {
m_map [cnt ] = triple (i, 0, j, false); m_map [cnt ] = triple (i, 0, j);
m_rmap[cnt++] = GetTetMode(i, 0, j); m_rmap[cnt++] = GetTetMode(i, 0, j);
} }
} }
...@@ -185,7 +172,7 @@ namespace Nektar ...@@ -185,7 +172,7 @@ namespace Nektar
{ {
for (int j = 1; j <= R-i; ++j) for (int j = 1; j <= R-i; ++j)
{ {
m_map [cnt ] = triple (1, i, j, false); m_map [cnt ] = triple (1, i, j);
m_rmap[cnt++] = GetTetMode(0, i, j); m_rmap[cnt++] = GetTetMode(0, i, j);
} }
} }
...@@ -195,7 +182,7 @@ namespace Nektar ...@@ -195,7 +182,7 @@ namespace Nektar
{ {
for (int j = 1; j <= R-i; ++j) for (int j = 1; j <= R-i; ++j)
{ {
m_map [cnt ] = triple (i, 1, j, false); m_map [cnt ] = triple (i, 1, j);
m_rmap[cnt++] = GetTetMode(i, 0, j); m_rmap[cnt++] = GetTetMode(i, 0, j);
} }
} }
...@@ -205,7 +192,7 @@ namespace Nektar ...@@ -205,7 +192,7 @@ namespace Nektar
{ {
for (int j = 1; j <= R-i; ++j) for (int j = 1; j <= R-i; ++j)
{ {
m_map [cnt ] = triple (0, i, j, false); m_map [cnt ] = triple (0, i, j);
m_rmap[cnt++] = GetTetMode(0, i, j); m_rmap[cnt++] = GetTetMode(0, i, j);
} }
} }
...@@ -219,8 +206,8 @@ namespace Nektar ...@@ -219,8 +206,8 @@ namespace Nektar
{ {
// need to go to j+1-th mode in the 'b' direction to // need to go to j+1-th mode in the 'b' direction to
// select correct modified_a mode // select correct modified_a mode
m_map [cnt ] = triple (i, j+1, k, true); m_map [cnt ] = triple (i, j+1, k);
m_rmap[cnt++] = GetModTetMode(i-1, j, k); m_rmap[cnt++] = GetTetMode(i-1, j, k);
} }
} }
} }
...@@ -376,15 +363,14 @@ namespace Nektar ...@@ -376,15 +363,14 @@ namespace Nektar
m_base[2]->GetBasisType() != LibUtilities::eModified_C, m_base[2]->GetBasisType() != LibUtilities::eModified_C,
"Basis[2] is not a general tensor type"); "Basis[2] is not a general tensor type");
int Qx = m_base[0]->GetNumPoints(); const int Qx = m_base[0]->GetNumPoints();
int Qy = m_base[1]->GetNumPoints(); const int Qy = m_base[1]->GetNumPoints();
int Qz = m_base[2]->GetNumPoints(); const int Qz = m_base[2]->GetNumPoints();
int s = 0;
Array<OneD, const NekDouble> bx = m_base[0]->GetBdata(); const Array<OneD, const NekDouble> &bx = m_base[0]->GetBdata();
Array<OneD, const NekDouble> by = m_base[1]->GetBdata(); const Array<OneD, const NekDouble> &by = m_base[1]->GetBdata();
Array<OneD, const NekDouble> bz = m_base[2]->GetBdata(); const Array<OneD, const NekDouble> &bz = m_base[2]->GetBdata();
Array<OneD, const NekDouble> bxc = m_base_A->GetBdata();
Array<OneD, const NekDouble> bzc = m_base_C->GetBdata();
int Q = m_base[2]->GetNumModes()-1; int Q = m_base[2]->GetNumModes()-1;
...@@ -392,33 +378,22 @@ namespace Nektar ...@@ -392,33 +378,22 @@ namespace Nektar
{ {
for (int j = 0; j < Qy; ++j) for (int j = 0; j < Qy; ++j)
{ {
for (int i = 0; i < Qx; ++i) for (int i = 0; i < Qx; ++i, ++s)
{ {
NekDouble sum = 0.0; NekDouble sum = 0.0;
for (int cnt = 0; cnt < m_ncoeffs; ++cnt) for (int cnt = 0; cnt < m_ncoeffs; ++cnt)
{ {
triple &idx = m_map[cnt]; triple &idx = m_map[cnt];
const int p = boost::get<0>(idx); const int p = boost::get<0>(idx);
const int q = boost::get<1>(idx); const int q = boost::get<1>(idx);
const int r = boost::get<2>(idx); const int r = boost::get<2>(idx);
const bool in = boost::get<3>(idx);
NekDouble tmp; NekDouble tmp;
if (in) tmp = inarray[cnt]*
{ bx[i + Qx*p]*
tmp = inarray[cnt]* by[j + Qy*q]*
bxc[i + Qx*p]* bz[k + Qz*m_rmap[cnt]];
bxc[j + Qy*q]*
bzc[k + Qz*m_rmap[cnt]];
}
else
{
tmp = inarray[cnt]*
bx[i + Qx*p]*
by[j + Qy*q]*
bz[k + Qz*m_rmap[cnt]];
}
if (r == 0 && p >= 2 && q >= 2) if (r == 0 && p >= 2 && q >= 2)
{ {
...@@ -428,11 +403,6 @@ namespace Nektar ...@@ -428,11 +403,6 @@ namespace Nektar
{ {
tmp *= bz[k + Qz*GetTetMode(blah-2,0,0)]; tmp *= bz[k + Qz*GetTetMode(blah-2,0,0)];
} }
//int blah = min(p+q-2, Q)-p-1;
//if (blah >= 0)
//{
// tmp *= bz[k + Qz*GetTetMode(q-3,0,0)];
//}
} }
sum += tmp; sum += tmp;
} }
...@@ -443,7 +413,7 @@ namespace Nektar ...@@ -443,7 +413,7 @@ namespace Nektar
sum += inarray[m]*bz[k+Qz]*(bx[i+Qx]*by[j+Qy]+ sum += inarray[m]*bz[k+Qz]*(bx[i+Qx]*by[j+Qy]+
bx[i ]*by[j+Qy]+ bx[i ]*by[j+Qy]+
bx[i+Qx]*by[j ]); bx[i+Qx]*by[j ]);
outarray[i + Qx*(j + Qy*k)] = sum; outarray[s] = sum;
} }
} }
} }
...@@ -537,12 +507,10 @@ namespace Nektar ...@@ -537,12 +507,10 @@ namespace Nektar
MultiplyByStdQuadratureMetric(inarray, tmp); MultiplyByStdQuadratureMetric(inarray, tmp);
const Array<OneD, const NekDouble> &bx = m_base[0]->GetBdata(); const Array<OneD, const NekDouble> &bx = m_base[0]->GetBdata();
const Array<OneD, const NekDouble> &by = m_base[1]->GetBdata(); const Array<OneD, const NekDouble> &by = m_base[1]->GetBdata();
const Array<OneD, const NekDouble> &bz = m_base[2]->GetBdata(); const Array<OneD, const NekDouble> &bz = m_base[2]->GetBdata();
const Array<OneD, const NekDouble> &bxc = m_base_A ->GetBdata();
const Array<OneD, const NekDouble> &bzc = m_base_C ->GetBdata();
int Q = m_base[1]->GetNumModes()-1; int Q = m_base[1]->GetNumModes()-1;
// Initial pyramid implementation. We need to iterate over vertices, // Initial pyramid implementation. We need to iterate over vertices,
...@@ -554,7 +522,6 @@ namespace Nektar ...@@ -554,7 +522,6 @@ namespace Nektar
const int p = boost::get<0>(idx); const int p = boost::get<0>(idx);
const int q = boost::get<1>(idx); const int q = boost::get<1>(idx);
const int r = boost::get<2>(idx); const int r = boost::get<2>(idx);
const bool in = boost::get<3>(idx);
NekDouble tmp = 0.0; NekDouble tmp = 0.0;
int s = 0; int s = 0;
...@@ -565,20 +532,10 @@ namespace Nektar ...@@ -565,20 +532,10 @@ namespace Nektar
for (int i = 0; i < Qx; ++i, ++s) for (int i = 0; i < Qx; ++i, ++s)
{ {
NekDouble tmp2; NekDouble tmp2;
if (in) tmp2 = inarray[s]*
{ bx[i + Qx*p]*
tmp2 = inarray[s]* by[j + Qy*q]*
bxc[i + Qx*p]* bz[k + Qz*m_rmap[cnt]];
bxc[j + Qy*q]*
bzc[k + Qz*m_rmap[cnt]];
}
else
{
tmp2 = inarray[s]*
bx[i + Qx*p]*
by[j + Qy*q]*
bz[k + Qz*m_rmap[cnt]];
}
if (r == 0 && p >= 2 && q >= 2) if (r == 0 && p >= 2 && q >= 2)
{ {
...@@ -1117,25 +1074,6 @@ namespace Nektar ...@@ -1117,25 +1074,6 @@ namespace Nektar
return cnt + K; return cnt + K;
} }
int StdPyrExp::GetModTetMode(const int I, const int J, const int K)
{
const int R = m_base[2]->GetNumModes()+1;
int i, j, cnt = 0;
for (i = 0; i < I; ++i)
{
cnt += (R-i)*(R-i+1)/2;
}
i = R-I;
for (j = 0; j < J; ++j)
{
cnt += i;
i--;
}
return cnt + K;
}
void StdPyrExp::v_MultiplyByStdQuadratureMetric( void StdPyrExp::v_MultiplyByStdQuadratureMetric(
const Array<OneD, const NekDouble>& inarray, const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble>& outarray) Array<OneD, NekDouble>& outarray)
......
...@@ -48,7 +48,8 @@ namespace Nektar ...@@ -48,7 +48,8 @@ namespace Nektar
class StdPyrExp : virtual public StdExpansion3D class StdPyrExp : virtual public StdExpansion3D
{ {
public: public:
typedef boost::tuple<unsigned char, unsigned char, unsigned char, bool> triple; typedef boost::tuple<
unsigned char, unsigned char, unsigned char> triple;
STD_REGIONS_EXPORT StdPyrExp(); STD_REGIONS_EXPORT StdPyrExp();
...@@ -76,7 +77,6 @@ namespace Nektar ...@@ -76,7 +77,6 @@ namespace Nektar
return m_rmap; return m_rmap;
} }
STD_REGIONS_EXPORT int GetTetMode(int I, int J, int K); STD_REGIONS_EXPORT int GetTetMode(int I, int J, int K);
STD_REGIONS_EXPORT int GetModTetMode(int I, int J, int K);
protected: protected:
//--------------------------------------- //---------------------------------------
...@@ -188,13 +188,14 @@ namespace Nektar ...@@ -188,13 +188,14 @@ namespace Nektar
//--------------------------------------- //---------------------------------------
STD_REGIONS_EXPORT virtual void v_GetFaceToElementMap( STD_REGIONS_EXPORT virtual void v_GetFaceToElementMap(
const int fid, const int fid,
const Orientation faceOrient, const Orientation faceOrient,
Array<OneD, unsigned int> &maparray, Array<OneD, unsigned int> &maparray,
Array<OneD, int> &signarray, Array<OneD, int> &signarray,
int nummodesA=-1, int nummodesA=-1,
int nummodesB=-1); int nummodesB=-1);
STD_REGIONS_EXPORT virtual int v_GetVertexMap(int localVertexId, STD_REGIONS_EXPORT virtual int v_GetVertexMap(
bool useCoeffPacking = false); int localVertexId,
bool useCoeffPacking = false);
/* /*
STD_REGIONS_EXPORT virtual void v_GetEdgeInteriorMap( STD_REGIONS_EXPORT virtual void v_GetEdgeInteriorMap(
const int eid, const int eid,
...@@ -224,8 +225,6 @@ namespace Nektar ...@@ -224,8 +225,6 @@ namespace Nektar
//--------------------------------------- //---------------------------------------
vector<triple> m_map; vector<triple> m_map;
vector<int> m_rmap; vector<int> m_rmap;
LibUtilities::BasisSharedPtr m_base_A;
LibUtilities::BasisSharedPtr m_base_C;
}; };
typedef boost::shared_ptr<StdPyrExp> StdPyrExpSharedPtr; typedef boost::shared_ptr<StdPyrExp> StdPyrExpSharedPtr;
} //end of namespace } //end of namespace
......
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