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

Version which now passes all the library regression tests locally

parent 4b7010f8
......@@ -8,10 +8,10 @@
</files>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">5.36868e-07</value>
<value tolerance="1e-12">5.22808e-07</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">3.78881e-06</value>
<value tolerance="1e-12">4.16465e-06</value>
</metric>
</metrics>
</test>
......@@ -8,7 +8,7 @@
</files>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">4.61818e-05</value>
<value tolerance="1e-12">4.1907e-05</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">0.000329233</value>
......
......@@ -619,7 +619,11 @@ namespace Nektar
// face 0
for(i = 0; i < numPoints; ++i)
{
mode[i] = pow(m_bdata[i],p+q); // [(1-z)/2]^{p+q}
// [(1-z)/2]^{p+q-2} Note in book it
// seems to suggest p+q-1 but that
// does not seem to give complete
// polynomial space for pyramids
mode[i] = pow(m_bdata[i],p+q-2);
}
one_m_z_pow = mode;
......@@ -628,7 +632,7 @@ namespace Nektar
// interior
for(int r = 1; r < numModes-max(p,q); ++r)
{
Polylib::jacobfd(numPoints,z.data(),mode,NULL,r-1,2*p+2*q-1,1.0);
Polylib::jacobfd(numPoints,z.data(),mode,NULL,r-1,2*p+2*q-3,1.0);
for(i = 0; i < numPoints; ++i)
{
......
......@@ -3306,7 +3306,7 @@ namespace Nektar
returnval.push_back(bkey);
const LibUtilities::PointsKey pkey1(nummodes+quadoffset-1, LibUtilities::eGaussRadauMAlpha2Beta0);
LibUtilities::BasisKey bkey1(LibUtilities::eModified_C, nummodes, pkey1);
LibUtilities::BasisKey bkey1(LibUtilities::eModifiedPyr_C, nummodes, pkey1);
returnval.push_back(bkey1);
}
break;
......
......@@ -600,7 +600,7 @@ namespace Nektar
if (fabs(elementAaxis_length*faceBaxis_length
- fabs(dotproduct1)) > NekConstants::kNekZeroTol)
{
cout << "Warning: Prism axes not parallel" << endl;
cout << "Warning: Pyramid axes not parallel" << endl;
}
// if the result is negative, both axis point in reverse
......@@ -621,7 +621,7 @@ namespace Nektar
if (fabs(elementBaxis_length*faceAaxis_length
- fabs(dotproduct2)) > NekConstants::kNekZeroTol)
{
cout << "Warning: Prism axes not parallel" << endl;
cout << "Warning: Pyramid axes not parallel" << endl;
}
if( dotproduct2 < 0.0 )
......@@ -710,7 +710,7 @@ namespace Nektar
LibUtilities::PointsKey(
order1+1, LibUtilities::eGaussLobattoLegendre));
const LibUtilities::BasisKey C(
LibUtilities::eModified_C, order2,
LibUtilities::eModifiedPyr_C, order2,
LibUtilities::PointsKey(
order2, LibUtilities::eGaussRadauMAlpha2Beta0));
......
......@@ -422,6 +422,7 @@ namespace Nektar
}
case LibUtilities::eModified_B:
case LibUtilities::eModified_C:
case LibUtilities::eModifiedPyr_C:
{
switch (facedir)
{
......@@ -484,6 +485,7 @@ namespace Nektar
case LibUtilities::eOrtho_A:
case LibUtilities::eOrtho_B:
case LibUtilities::eOrtho_C:
case LibUtilities::eOrthoPyr_C:
{
switch (facedir)
{
......
......@@ -29,7 +29,7 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: pyramadic routines built upon StdExpansion3D
// Description: pyramidic routines built upon StdExpansion3D
//
///////////////////////////////////////////////////////////////////////////////
......@@ -66,156 +66,10 @@ namespace Nektar
"than order in 'c' direction");
ASSERTL0(Bb.GetNumModes() <= Bc.GetNumModes(), "order in 'b' direction is higher "
"than order in 'c' direction");
#if 0
// Set up mode mapping which takes 0\leq i\leq N_coeffs -> (p,q,r)
// of the 3D tensor product
const int P = Ba.GetNumModes() - 1;
const int Q = Bb.GetNumModes() - 1;
const int R = Bc.GetNumModes() - 1;
int cnt = 0;
// Vertices
m_map[Mode(0, 0, 0, 0)] = cnt++;
m_map[Mode(1, 0, 0, 0)] = cnt++;
m_map[Mode(1, 1, 0, 0)] = cnt++;
m_map[Mode(0, 1, 0, 0)] = cnt++;
m_map[Mode(0, 0, 1, 1)] = cnt++;
// Edge 0
for (int i = 2; i <= P; ++i)
{
m_map[Mode(i, 0, 0, GetTetMode(i, 0, 0))] = cnt++;
}
// Edge 1
for (int i = 2; i <= Q; ++i)
{
m_map[Mode(1, i, 0, GetTetMode(0, i, 0))] = cnt++;
}
// Edge 2
for (int i = 2; i <= P; ++i)
{
m_map[Mode(i, 1, 0, GetTetMode(i, 0, 0))] = cnt++;
}
// Edge 3
for (int i = 2; i <= Q; ++i)
{
m_map[Mode(0, i, 0, GetTetMode(0, i, 0))] = cnt++;
}
// Edge 4
for (int i = 2; i <= R; ++i)
{
m_map[Mode(0, 0, i, i)] = cnt++;
}
// Edge 5
for (int i = 2; i <= R; ++i)
{
m_map[Mode(1, 0, i, i)] = cnt++;
}
// Edge 6
for (int i = 2; i <= R; ++i)
{
m_map[Mode(1, 1, i, i)] = cnt++;
}
// Edge 7
for (int i = 2; i <= R; ++i)
{
m_map[Mode(0, 1, i, i)] = cnt++;
}
// Face 0 - TODO check this
for (int j = 2; j <= Q; ++j)
{
for (int i = 2; i <= P; ++i)
{
m_map[Mode(i, j, 0, GetTetMode((i-2+j-2) % (Q-1) + 2, 0, 0))] = cnt++;
}
}
// Face 1
for (int i = 2; i <= P; ++i)
{
for (int j = 1; j <= R-i; ++j)
{
m_map[Mode(i, 0, j, GetTetMode(i, 0, j))] = cnt++;
}
}
// Face 2
for (int i = 2; i <= Q; ++i)
{
for (int j = 1; j <= R-i; ++j)
{
m_map[Mode(1, i, j, GetTetMode(0, i, j))] = cnt++;
}
}
// Face 3
for (int i = 2; i <= P; ++i)
{
for (int j = 1; j <= R-i; ++j)
{
m_map[Mode(i, 1, j, GetTetMode(i, 0, j))] = cnt++;
}
}
// Face 4
for (int i = 2; i <= Q; ++i)
{
for (int j = 1; j <= R-i; ++j)
{
m_map[Mode(0, i, j, GetTetMode(0, i, j))] = cnt++;
}
}
// Interior (tetrahedral modes)
for (int i = 2; i <= P+1; ++i)
{
for (int j = 1; j <= Q-i+1; ++j)
{
for (int k = 1; k <= R-i-j+1; ++k)
{
// need to go to j+1-th mode in the 'b' direction to
// select correct modified_a mode
m_map[Mode(i, j+1, k, GetTetMode(i-1, j, k))] = cnt++;
}
}
}
ASSERTL0(m_map.size() == m_ncoeffs,
"Duplicate coefficient entries in map");
map<Mode, unsigned int, cmpop>::iterator it;
for (it = m_map.begin(); it != m_map.end(); ++it)
{
const int p = it->first.get<0>();
const int q = it->first.get<1>();
const int r = it->first.get<2>();
const int rp = it->first.get<3>();
if (m_idxMap.find(p) == m_idxMap.end())
{
m_idxMap[p] = map<int, map<int, pair<int, int> > >();
}
if (m_idxMap[p].find(q) == m_idxMap[p].end())
{
m_idxMap[p][q] = map<int, pair<int, int> >();
}
if (m_idxMap[p][q].find(r) == m_idxMap[p][q].end())
{
m_idxMap[p][q][r] = pair<int, int>(it->second, rp);
}
}
#endif
ASSERTL1(Bc.GetBasisType() == LibUtilities::eModifiedPyr_C ||
Bc.GetBasisType() == LibUtilities::eOrthoPyr_C,
"Expected basis type in 'c' direction to be ModifiedPyr_C or OrthoPyr_C");
}
StdPyrExp::StdPyrExp(const StdPyrExp &T)
......@@ -389,13 +243,14 @@ 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();
Array<OneD, NekDouble> wsp(nquad2*order0*order1+
nquad2*nquad1*order0);
nquad2*nquad1*nquad0);
v_BwdTrans_SumFacKernel(m_base[0]->GetBdata(),
m_base[1]->GetBdata(),
......@@ -415,81 +270,6 @@ namespace Nektar
bool doCheckCollDir1,
bool doCheckCollDir2)
{
#if 0
const int Qx = m_base[0]->GetNumPoints();
const int Qy = m_base[1]->GetNumPoints();
const int Qz = m_base[2]->GetNumPoints();
const NekDouble *bx = base0.get();
const NekDouble *by = base1.get();
const NekDouble *bz = base2.get();
// Need to count coeffs for storage...
map<int, map<int, map<int, pair<int, int> > > >::iterator it_p;
map<int, map<int, pair<int, int> > > ::iterator it_q;
map<int, pair<int, int> > ::iterator it_r;
int pqcnt = 0;
for (it_p = m_idxMap.begin(); it_p != m_idxMap.end(); ++it_p)
{
for (it_q = it_p->second.begin(); it_q != it_p->second.end(); ++it_q)
{
pqcnt++;
}
}
Array<OneD, NekDouble> fpq(pqcnt);
Array<OneD, NekDouble> fp (m_base[0]->GetNumModes());
int i ,j, k, s = 0, cnt = 0, cnt2 = 0;
for (k = 0; k < Qz; ++k)
{
NekDouble bz1 = bz[k+Qz];
cnt = 0;
for (it_p = m_idxMap.begin(); it_p != m_idxMap.end(); ++it_p)
{
for (it_q = it_p->second.begin(); it_q != it_p->second.end(); ++it_q)
{
NekDouble sum = 0.0;
for (it_r = it_q->second.begin(); it_r != it_q->second.end(); ++it_r)
{
sum += inarray[it_r->second.first] * bz[k + Qz*it_r->second.second];
}
fpq[cnt++] = sum;
}
}
for (j = 0; j < Qy; ++j)
{
NekDouble by0 = bz1*by[j];
NekDouble by1 = bz1*by[j+Qy];
cnt = cnt2 = 0;
for (it_p = m_idxMap.begin(); it_p != m_idxMap.end(); ++it_p)
{
NekDouble sum = 0.0;
for (it_q = it_p->second.begin(); it_q != it_p->second.end(); ++it_q)
{
sum += by[j + Qy*it_q->first] * fpq[cnt++];
}
fp[cnt2++] = sum;
}
for (i = 0; i < Qx; ++i, ++s)
{
cnt2 = 0;
NekDouble sum = 0.0;
for (it_p = m_idxMap.begin(); it_p != m_idxMap.end(); ++it_p)
{
sum += bx[i + Qx*it_p->first] * fp[cnt2++];
}
sum += inarray[4]*(by1*(bx[i] + bx[i+Qx]) + by0*bx[i+Qx]);
outarray[s] = sum;
}
}
}
#else
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nquad2 = m_base[2]->GetNumPoints();
......@@ -562,7 +342,6 @@ namespace Nektar
tmp1.get(), nquad1*nquad2,
0.0, outarray.get(), nquad0);
#endif
}
/** \brief Forward transform from physical quadrature space
......@@ -677,84 +456,6 @@ namespace Nektar
bool doCheckCollDir1,
bool doCheckCollDir2)
{
#if 0
int i, j, k, s;
int Qx = m_base[0]->GetNumPoints();
int Qy = m_base[1]->GetNumPoints();
int Qz = m_base[2]->GetNumPoints();
const NekDouble *bx = base0.get();
const NekDouble *by = base1.get();
const NekDouble *bz = base2.get();
map<int, map<int, map<int, pair<int, int> > > >::iterator it_p;
map<int, map<int, pair<int, int> > > ::iterator it_q;
map<int, pair<int, int> > ::iterator it_r;
Array<OneD, NekDouble> f (Qy*Qz);
Array<OneD, NekDouble> fb(Qz);
for (it_p = m_idxMap.begin(); it_p != m_idxMap.end(); ++it_p)
{
const int p = it_p->first;
s = 0;
for (k = 0; k < Qz; ++k)
{
for (j = 0; j < Qy; ++j)
{
NekDouble sum = 0.0;
for (i = 0; i < Qx; ++i, ++s)
{
sum += bx[i + Qx*p]*inarray[s];
}
f[j+Qy*k] = sum;
}
}
for (it_q = it_p->second.begin(); it_q != it_p->second.end(); ++it_q)
{
const int q = it_q->first;
for (k = 0; k < Qz; ++k)
{
NekDouble sum = 0.0;
for (j = 0; j < Qy; ++j)
{
sum += by[j + Qy*q]*f[j+Qy*k];
}
fb[k] = sum;
}
for (it_r = it_q->second.begin(); it_r != it_q->second.end(); ++it_r)
{
const int rpqr = it_r->second.second;
NekDouble sum = 0.0;
for (k = 0; k < Qz; ++k)
{
sum += bz[k + Qz*rpqr]*fb[k];
}
outarray[it_r->second.first] = sum;
}
}
}
// Correct for top mode
s = 0;
for (k = 0; k < Qz; ++k)
{
for (j = 0; j < Qy; ++j)
{
for (i = 0; i < Qx; ++i, ++s)
{
outarray[4] += inarray[s] * bz[k+Qz]*(
bx[i+Qx]*by[j+Qy] +
bx[i+Qx]*by[j ] +
bx[i ]*by[j+Qy]);
}
}
}
#else
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nquad2 = m_base[2]->GetNumPoints();
......@@ -825,7 +526,6 @@ namespace Nektar
outarray[1] += Blas::Ddot(nquad2,base2.get()+nquad2,1,
&tmp2[nquad2*order1+nquad2],1);
}
#endif
}
void StdPyrExp::v_IProductWRTDerivBase(
......@@ -857,7 +557,13 @@ namespace Nektar
Array<OneD, NekDouble> gfac1(nquad1);
Array<OneD, NekDouble> gfac2(nquad2);
Array<OneD, NekDouble> tmp0 (nqtot);
Array<OneD, NekDouble> wsp;
int order0 = m_base[0]->GetNumModes();
int order1 = m_base[1]->GetNumModes();
Array<OneD, NekDouble> wsp(nquad1*nquad2*order0 +
nquad2*order0*order1);
const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
......@@ -1285,7 +991,7 @@ namespace Nektar
"(x and y direction) and ModifiedPyr_C BasisType (z "
"direction)");
int i, j, k, p, q, r, nFaceCoeffs;
int i, j, k, p, q, r, nFaceCoeffs, idx = 0;
int nummodesA=0, nummodesB=0;
int order0 = m_base[0]->GetNumModes();
......@@ -1368,299 +1074,95 @@ namespace Nektar
// Set up ordering inside each 2D face. Also for triangular faces,
// populate signarray.
int cnt = 0, cnt2;
switch (fid)
{
case 0: // Bottom quad
// Fill in vertices
maparray[arrayindx[0]] = 0;
maparray[arrayindx[1]] = 1;
maparray[arrayindx[P+1]] = 2;
maparray[arrayindx[P]] = 3;
// Edge 0
cnt = 5;
for (p = 2; p < P; ++p)
{
maparray[arrayindx[p]] = p-2 + cnt;
}
// Edge 1
cnt += P-2;
for (q = 2; q < Q; ++q)
{
maparray[arrayindx[q*P+1]] = q-2 + cnt;
}
// Edge 2
cnt += Q-2;
for (p = 2; p < P; ++p)
{
maparray[arrayindx[P+p]] = p-2 + cnt;
}
// Edge 3
cnt += P-2;
for (q = 2; q < Q; ++q)
{
maparray[arrayindx[q*P]] = q-2 + cnt;
}
// Interior
cnt += Q-2 + 4*(P-2);
for (q = 2; q < Q; ++q)
{
for (p = 2; p < P; ++p)
{
maparray[arrayindx[q*P+p]] = cnt + (q-2)*P+(p-2);
}
}
break;
case 1: // Left triangle
// Vertices
maparray[0] = 0;
maparray[1] = 4;
maparray[Q] = 1;
// Edge 0 (pyramid edge 0)
cnt = 5;
q = 2*Q-1;
for (p = 2; p < P; q += Q-p, ++p)
{
maparray[q] = cnt++;
if ((int)faceOrient == 7)
{
signarray[q] = p % 2 ? -1 : 1;
}
}
// Edge 1 (pyramid edge 5)
cnt = 5 + 2*(order0-2) + 2*(order1-2) + (order2-2);
for (q = 2; q < Q; ++q)
{
maparray[q] = cnt++;
if ((int)faceOrient == 7)
{
signarray[q] = q % 2 ? -1 : 1;
}
}
// Edge 2 (pyramid edge 4)
cnt = 5 + 2*(order0-2) + 2*(order1-2);
for (q = 2; q < Q; ++q)
{
maparray[Q+q-1] = cnt++;
if ((int)faceOrient == 7)
{
signarray[Q+q-1] = q % 2 ? -1 : 1;
}
}
// Interior
cnt = 5 + 2*(order0-2) + 2*(order1-2) + 4*(order2-2)
+ v_GetFaceIntNcoeffs(0);
cnt2 = 2*Q + 1;
for (p = 2; p < P; ++p)
{
for (r = 2; r < Q-p; ++r)
{
maparray[cnt2] = cnt++;
if ((int)faceOrient == 7 && p > 1)
{
signarray[cnt2++] = p % 2 ? -1 : 1;
}
}
cnt2++;
}
break;
case 2:
// Vertices
maparray[0] = 1;
maparray[1] = 4;
maparray[Q] = 2;
// Edge 0 (pyramid edge 1)
cnt = 5 + (order0-2);
q = 2*Q-1;
for (p = 2; p < P; q += Q-p, ++p)
{
maparray[q] = cnt++;
if ((int)faceOrient == 7)
{
signarray[q] = p % 2 ? -1 : 1;
}
}