Commit 75c0ad6e authored by Spencer Sherwin's avatar Spencer Sherwin
Browse files

Added Modified_C basis that can be used with pyramids and also works with a...

Added Modified_C basis that can be used with pyramids and also works with a modified and orothoogonal expansion. Currently orthogonal expansion produced diagonal but with non-unit values
parent bf40c447
......@@ -66,6 +66,9 @@ int main(int argc, char *argv[]){
fprintf(stderr,"\t Modified_A = 4\n");
fprintf(stderr,"\t Modified_B = 5\n");
fprintf(stderr,"\t Modified_C = 6\n");
fprintf(stderr,"\t Ortho_C = 7\n");
fprintf(stderr,"\t Modified_C = 8\n");
fprintf(stderr,"\t Fourier = 7\n");
fprintf(stderr,"\t Lagrange = 8\n");
fprintf(stderr,"\t Gauss Lagrange = 9\n");
......@@ -162,23 +165,23 @@ int main(int argc, char *argv[]){
}
break;
case LibUtilities::ePyramid:
if((btype1 == eOrtho_B) || (btype1 == eOrtho_C)
|| (btype1 == eModified_B) || (btype1 == eModified_C))
if((btype1 == eOrtho_B) || (btype1 == eOrtho_C) || (btype1 == eOrthoPyr_C)
|| (btype1 == eModified_B) || (btype1 == eModified_C) || (btype1 == eModifiedPyr_C))
{
NEKERROR(ErrorUtil::efatal, "Basis 1 cannot be of type Ortho_B, "
"Ortho_C, Modified_B or Modified_C");
}
if((btype2 == eOrtho_B) || (btype2 == eOrtho_C)
|| (btype2 == eModified_B) || (btype2 == eModified_C))
if((btype2 == eOrtho_B) || (btype2 == eOrtho_C) || (btype1 == eOrthoPyr_C)
|| (btype2 == eModified_B) || (btype2 == eModified_C) || (btype1 == eModifiedPyr_C))
{
NEKERROR(ErrorUtil::efatal, "Basis 2 cannot be of type Ortho_B, "
"Ortho_C, Modified_B or Modified_C");
}
if((btype3 == eOrtho_A) || (btype3 == eOrtho_B)
|| (btype3 == eModified_A) || (btype3 == eModified_B))
if((btype3 == eOrtho_A) || (btype3 == eOrtho_B) || (btype3 == eOrtho_C)
|| (btype3 == eModified_A) || (btype3 == eModified_B) || (btype3 == eModified_C))
{
NEKERROR(ErrorUtil::efatal, "Basis 3 cannot be of type Ortho_A, "
"Ortho_B, Modified_A or Modified_B");
"Ortho_B, Ortho_C, Modified_A, Modified_B or Modified_C");
}
break;
case LibUtilities::ePrism:
......
......@@ -240,26 +240,19 @@ namespace Nektar
"than order in 'c' direction.");
// Count number of coefficients explicitly.
const int Pi = Na - 2, Qi = Nb - 2, Ri = Nc - 2;
int nCoeff =
5 + // vertices
Pi * 2 + Qi * 2 + Ri * 4 + // base edges
Pi * Qi + // base quad
Pi * (2*Ri - Pi - 1) + // p-r triangles;
Qi * (2*Ri - Qi - 1); // q-r triangles;
int nCoeff = 0;
// Count number of interior tet modes
for (int a = 0; a < Pi - 1; ++a)
for (int a = 0; a < Na; ++a)
{
for (int b = 0; b < Qi - a - 1; ++b)
for (int b = 0; b < Nb; ++b)
{
for (int c = 0; c < Ri - a - b -1; ++c)
for (int c = 0; c < Nc - std::max(a,b); ++c)
{
++nCoeff;
}
}
}
return nCoeff;
}
......
......@@ -328,6 +328,47 @@ namespace Nektar
}
break;
/** \brief Orthogonal basis C for Pyramid expansion (which is richer than tets)
\f$\tilde \psi_{pqr}^c = \left ( {1 - \eta_3} \over 2 \right)^{p+q} P_r^{2p+2q+2, 0}(\eta_3)\f$ \\
*/
// This is tilde psi_pqr in Spen's book, page 105
// The 4-dimensional array is laid out in memory such that
// 1) Eta_z values are the changing the fastest, then r, q, and finally p.
// 2) r index increases by the stride of numPoints.
case eOrthoPyr_C:
{
int P = numModes - 1, Q = numModes - 1, R = numModes - 1;
NekDouble *mode = m_bdata.data();
for( int p = 0; p <= P; ++p )
{
for( int q = 0; q <= Q; ++q )
{
for( int r = 0; r <= R - max(p,q); ++r, mode += numPoints )
{
Polylib::jacobfd(numPoints, z.data(), mode, NULL, r, 2*p + 2*q + 2.0, 0.0);
for( int k = 0; k < numPoints; ++k )
{
// Note factor of 0.5 is part of normalisation
mode[k] *= pow(0.5*(1.0 - z[k]), p+q);
// finish normalisation
mode[k] *= sqrt(r+p+q+1.5);
}
}
}
}
// Define derivative basis
Blas::Dgemm('n','n',numPoints,numModes*(numModes+1)*
(numModes+2)/6,numPoints,1.0, D, numPoints,
m_bdata.data(),numPoints,0.0,m_dbdata.data(),numPoints);
}
break;
case eModified_A:
// Note the following packing deviates from the
......@@ -452,13 +493,13 @@ namespace Nektar
case eModified_C:
{
// Note the following packing deviates from the
// definition in the Book by Karniadakis in two
// ways. 1) We put the vertex degrees of freedom
// at the lower index range to follow a more
// hierarchic structure. 2) We do not duplicate
// the singular vertex definition (or the
// duplicated face information in the book ) so
// that only a tetrahedral number
// definition in the Book by Karniadakis &
// Sherwin in two ways. 1) We put the vertex
// degrees of freedom at the lower index range to
// follow a more hierarchic structure. 2) We do
// not duplicate the singular vertex definition
// (or the duplicated face information in the book
// ) so that only a tetrahedral number
// (i.e. (modes)*(modes+1)*(modes+2)/6) of modes
// are required consistent with the orthogonal
// basis.
......@@ -505,6 +546,110 @@ namespace Nektar
}
break;
case eModifiedPyr_C:
{
// Note the following packing deviates from the
// definition in the Book by Karniadakis &
// Sherwin in two ways. 1) We put the vertex
// degrees of freedom at the lower index range to
// follow a more hierarchic structure. 2) We do
// not duplicate the singular vertex definition
// so that only a pyramidic number
// (i.e. (modes)*(modes+1)*(2*modes+1)/6) of modes
// are required consistent with the orthogonal
// basis.
// In the current structure the r index runs
// fastest rollowed by q and than the p index so
// that the matrix has a more compact structure
// Generate Modified_B basis;
BasisKey ModBKey(eModified_B,m_basisKey.GetNumModes(),
m_basisKey.GetPointsKey());
BasisSharedPtr ModB = BasisManager()[ModBKey];
Array<OneD, const NekDouble> ModB_data = ModB->GetBdata();
// Copy Modified_B basis into first
// (numModes*(numModes+1)/2)*numPoints entires of
// bdata.
int N;
int B_offset = 0;
int offset = 0;
// Vertex 0,3,4, edges 3,4,7, face 4
N = numPoints*(numModes)*(numModes+1)/2;
Vmath::Vcopy(N, &ModB_data[0],1,&m_bdata[0],1);
offset += N;
B_offset += numPoints*(numModes);
// Vertex 1 edges 5
N = numPoints*(numModes-1);
Vmath::Vcopy(N, &ModB_data[0]+B_offset,1,&m_bdata[0]+offset,1);
offset += N;
// Vertex 2 edges 1,6, face 2
N = numPoints*(numModes-1)*(numModes)/2;
Vmath::Vcopy(N, &ModB_data[0]+B_offset,1,&m_bdata[0]+offset,1);
offset += N;
B_offset += numPoints*(numModes-1);
NekDouble *one_m_z_pow, *one_p_z;
NekDouble *mode;
mode = m_bdata.data() + offset;
for(p = 2; p < numModes; ++p)
{
// edges 0 2, faces 1 3
N = numPoints*(numModes-p);
Vmath::Vcopy(N, &ModB_data[0]+B_offset,1,mode,1);
mode += N;
Vmath::Vcopy(N, &ModB_data[0]+B_offset,1,mode,1);
mode += N;
B_offset += N;
one_m_z_pow = m_bdata.data();
one_p_z = m_bdata.data()+numPoints;
for(q = 2; q < numModes; ++q)
{
// face 0
for(i = 0; i < numPoints; ++i)
{
mode[i] = pow(one_m_z_pow[i],p+q); // [(1-z)/2]^{p+q}
}
one_m_z_pow = mode;
mode += numPoints;
// interior
for(int r = 1; r < numModes-max(p,q); ++r, mode+=numPoints)
{
Polylib::jacobfd(numPoints,z.data(),mode,NULL,r-1,2*p+2*q-1,1.0);
for(i = 0; i < numPoints; ++i)
{
mode[i] *= one_m_z_pow[i]*one_p_z[i];
}
mode += numPoints;
}
}
}
// set up derivative of basis.
Blas::Dgemm('n','n',numPoints,
numModes*(numModes+1)*(2*numModes+1)/6,
numPoints,1.0,D,numPoints,
m_bdata.data(),numPoints,0.0,
m_dbdata.data(),numPoints);
}
break;
case eGLL_Lagrange:
{
mode = m_bdata.data();
......@@ -523,7 +668,7 @@ namespace Nektar
Blas::Dgemm('n', 'n', numPoints, numModes, numPoints, 1.0,
D, numPoints, m_bdata.data(), numPoints, 0.0,
m_dbdata.data(), numPoints);
}//end scope
break;
case eGauss_Lagrange:
......
......@@ -102,7 +102,10 @@ namespace Nektar
case eOrtho_C:
value = m_nummodes*(m_nummodes+1)*(m_nummodes+2)/6;
break;
case eModifiedPyr_C:
case eOrthoPyr_C:
value = m_nummodes*(m_nummodes+1)*(2*m_nummodes+1)/6;
break;
case eOrtho_A:
case eModified_A:
case eFourier:
......
......@@ -49,6 +49,8 @@ namespace Nektar
eModified_A, //!< Principle Modified Functions \f$ \phi^a_p(z_i) \f$
eModified_B, //!< Principle Modified Functions \f$ \phi^b_{pq}(z_i) \f$
eModified_C, //!< Principle Modified Functions \f$ \phi^c_{pqr}(z_i) \f$
eOrthoPyr_C, //!< Principle Orthogonal Functions \f$\widetilde{\psi}^c_{pqr}(z_i) for Pyramids\f$
eModifiedPyr_C, //!< Principle Modified Functions \f$ \phi^c_{pqr}(z_i) for Pyramids\f$
eFourier, //!< Fourier Expansion \f$ \exp(i p\pi z_i)\f$
eGLL_Lagrange, //!< Lagrange for SEM basis \f$ h_p(z_i) \f$
eGauss_Lagrange, //!< Lagrange Polynomials using the Gauss points \f$ h_p(z_i) \f$
......
......@@ -53,6 +53,8 @@ namespace Nektar
"Modified_A",
"Modified_B",
"Modified_C",
"OrthoPyr_C",
"ModifiedPyr_C",
"Fourier",
"GLL_Lagrange",
"Gauss_Lagrange",
......
......@@ -61,17 +61,13 @@ namespace Nektar
Bc.GetNumModes()),
Ba, Bb, Bc)
{
if (Ba.GetNumModes() > Bc.GetNumModes())
{
ASSERTL0(false, "order in 'a' direction is higher "
"than order in 'c' direction");
}
if (Bb.GetNumModes() > Bc.GetNumModes())
{
ASSERTL0(false, "order in 'b' direction is higher "
"than order in 'c' direction");
}
ASSERTL0(Ba.GetNumModes() <= Bc.GetNumModes(), "order in 'a' direction is higher "
"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;
......@@ -79,6 +75,7 @@ namespace Nektar
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++;
......@@ -218,6 +215,7 @@ namespace Nektar
m_idxMap[p][q][r] = pair<int, int>(it->second, rp);
}
}
#endif
}
StdPyrExp::StdPyrExp(const StdPyrExp &T)
......@@ -391,7 +389,14 @@ namespace Nektar
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble>& outarray)
{
Array<OneD, NekDouble> wsp;
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);
v_BwdTrans_SumFacKernel(m_base[0]->GetBdata(),
m_base[1]->GetBdata(),
m_base[2]->GetBdata(),
......@@ -410,6 +415,7 @@ 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();
......@@ -483,6 +489,80 @@ namespace Nektar
}
}
}
#else
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();
int order2 = m_base[2]->GetNumModes();
Array<OneD, NekDouble > tmp = wsp;
Array<OneD, NekDouble > tmp1 = tmp + nquad2*order0*order1;
int i, j, mode, mode1, cnt;
// Perform summation over '2' direction
mode = mode1 = cnt = 0;
for(i = 0; i < order0; ++i)
{
for(j = 0; j < order1; ++j, ++cnt)
{
int ijmax = max(i,j);
Blas::Dgemv('N', nquad2, order2-ijmax,
1.0, base2.get()+mode*nquad2, nquad2,
inarray.get()+mode1, 1,
0.0, tmp.get()+cnt*nquad2, 1);
mode += order2-ijmax;
mode1 += order2-ijmax;
}
//increment mode in case order1!=order2
for(j = order1; j < order2-i; ++j)
{
int ijmax = max(i,j);
mode += order2-ijmax;
}
}
// fix for modified basis by adding split of top singular
// vertex mode - currently (1+c)/2 x (1-b)/2 x (1-a)/2
// component is evaluated
if(m_base[0]->GetBasisType() == LibUtilities::eModified_A)
{
// Not sure why we could not use basis as 1.0
// top singular vertex - (1+c)/2 x (1+b)/2 x (1-a)/2 component
Blas::Daxpy(nquad2,inarray[1],base2.get()+nquad2,1,
&tmp[0]+nquad2,1);
// top singular vertex - (1+c)/2 x (1-b)/2 x (1+a)/2 component
Blas::Daxpy(nquad2,inarray[1],base2.get()+nquad2,1,
&tmp[0]+order1*nquad2,1);
// top singular vertex - (1+c)/2 x (1+b)/2 x (1+a)/2 component
Blas::Daxpy(nquad2,inarray[1],base2.get()+nquad2,1,
&tmp[0]+order1*nquad2+nquad2,1);
}
// Perform summation over '1' direction
mode = 0;
for(i = 0; i < order0; ++i)
{
Blas::Dgemm('N', 'T', nquad1, nquad2, order1,
1.0, base1.get(), nquad1,
tmp.get()+mode*nquad2, nquad2,
0.0, tmp1.get()+i*nquad1*nquad2, nquad1);
mode += order1;
}
// Perform summation over '0' direction
Blas::Dgemm('N', 'T', nquad0, nquad1*nquad2, order0,
1.0, base0.get(), nquad0,
tmp1.get(), nquad1*nquad2,
0.0, outarray.get(), nquad0);
#endif
}
/** \brief Forward transform from physical quadrature space
......@@ -504,14 +584,20 @@ namespace Nektar
v_IProductWRTBase(inarray,outarray);
// get Mass matrix inverse
StdMatrixKey masskey(eInvMass,DetShapeType(),*this);
StdMatrixKey masskey(eMass,DetShapeType(),*this);
DNekMatSharedPtr matsys = GetStdMatrix(masskey);
cout << *matsys << endl;
// get Mass matrix inverse
StdMatrixKey imasskey(eInvMass,DetShapeType(),*this);
DNekMatSharedPtr imatsys = GetStdMatrix(imasskey);
// copy inarray in case inarray == outarray
DNekVec in (m_ncoeffs, outarray);
DNekVec out(m_ncoeffs, outarray, eWrapper);
out = (*matsys)*in;
out = (*imatsys)*in;
}
......@@ -554,7 +640,14 @@ namespace Nektar
Array<OneD, NekDouble>& outarray,
bool multiplybyweights)
{
Array<OneD, NekDouble> wsp;
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(nquad1*nquad2*order0 +
nquad2*order0*order1);
if(multiplybyweights)
{
......@@ -589,6 +682,7 @@ 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();
......@@ -665,6 +759,78 @@ namespace Nektar
}
}
}
#else
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();
int order2 = m_base[2]->GetNumModes();
Array<OneD, NekDouble > tmp1 = wsp;
Array<OneD, NekDouble > tmp2 = wsp + nquad1*nquad2*order0;
int i,j, mode,mode1, cnt;
// Inner product with respect to the '0' direction
Blas::Dgemm('T', 'N', nquad1*nquad2, order0, nquad0,
1.0, inarray.get(), nquad0,
base0.get(), nquad0,
0.0, tmp1.get(), nquad1*nquad2);
// Inner product with respect to the '1' direction
for(mode=i=0; i < order0; ++i)
{
Blas::Dgemm('T', 'N', nquad2, order1, nquad1,
1.0, tmp1.get()+i*nquad1*nquad2, nquad1,
base1.get(), nquad1,
0.0, tmp2.get()+mode*nquad2, nquad2);
mode += order1;
}
// Inner product with respect to the '2' direction
mode = mode1 = cnt = 0;
for(i = 0; i < order0; ++i)
{
for(j = 0; j < order1; ++j, ++cnt)
{
int ijmax = max(i,j);
Blas::Dgemv('T', nquad2, order2-ijmax,
1.0, base2.get()+mode*nquad2, nquad2,
tmp2.get()+cnt*nquad2, 1,
0.0, outarray.get()+mode1, 1);
mode += order2-ijmax;
mode1 += order2-ijmax;
}
//increment mode in case order1!=order2
for(j = order1; j < order2; ++j)
{
int ijmax = max(i,j);
mode += order2-ijmax;
}
}
// fix for modified basis for top singular vertex component
// Already have evaluated (1+c)/2 (1-b)/2 (1-a)/2
if(m_base[0]->GetBasisType() == LibUtilities::eModified_A)
{
// add in (1+c)/2 (1+b)/2 (1-a)/2 component
outarray[1] += Blas::Ddot(nquad2,base2.get()+nquad2,1,
&tmp2[nquad2],1);
// add in (1+c)/2 (1-b)/2 (1+a)/2 component
outarray[1] += Blas::Ddot(nquad2,base2.get()+nquad2,1,
&tmp2[nquad2*order1],1);
// add in (1+c)/2 (1+b)/2 (1+a)/2 component
outarray[1] += Blas::Ddot(nquad2,base2.get()+nquad2,1,
&tmp2[nquad2*order1+nquad2],1);
}
#endif
}
void StdPyrExp::v_IProductWRTDerivBase(
......@@ -940,7 +1106,7 @@ namespace Nektar
ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
GetBasisType(1) == LibUtilities::eGLL_Lagrange,
"BasisType is not a boundary interior form");
ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
ASSERTL1(GetBasisType(2) == LibUtilities::eModifiedPyr_C ||
GetBasisType(2) == LibUtilities::eGLL_Lagrange,
"BasisType is not a boundary interior form");
......@@ -1119,9 +1285,9 @@ namespace Nektar
"Method only implemented if BasisType is identical"
"in x and y directions");
ASSERTL1(GetEdgeBasisType(0) == LibUtilities::eModified_A &&
GetEdgeBasisType(4) == LibUtilities::eModified_C,
GetEdgeBasisType(4) == LibUtilities::eModifiedPyr_C,
"Method only implemented for Modified_A BasisType"
"(x and y direction) and Modified_C BasisType (z "
"(x and y direction) and ModifiedPyr_C BasisType (z "
"direction)");
int i, j, k, p, q, r, nFaceCoeffs;
......
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