Commit f0e78ed8 authored by Chris Cantwell's avatar Chris Cantwell

Generalised num coefficient calculations.

parent 50d99e2f
......@@ -513,25 +513,32 @@ namespace Nektar
switch (m_meshComposites[cId].type)
{
case 'A':
weight = StdTetData::getNumberOfCoefficients(na, nb, nc);
bndWeight = 2 * (na-2) * (na-1) + 6 * (na-1) + 4;
weight = StdTetData::getNumberOfCoefficients(na, nb, nc);
bndWeight = StdTetData::getNumberOfBndCoefficients(na, nb, nc);
break;
case 'R':
weight = StdPrismData::getNumberOfCoefficients(na, nb, nc);
bndWeight = (na-2)*(na-1) + 3*(na-1)*(na-1) + 9*(na-1) + 6;
weight = StdPrismData::getNumberOfCoefficients(na, nb, nc);
bndWeight = StdPrismData::getNumberOfBndCoefficients(na, nb, nc);
break;
case 'H':
weight = StdHexData::getNumberOfCoefficients(na, nb, nc);
bndWeight = 6*(na-1)*(na-1) + 12*(na-1) + 8;
weight = StdHexData::getNumberOfCoefficients(na, nb, nc);
bndWeight = StdHexData::getNumberOfBndCoefficients(na, nb, nc);
break;
case 'P':
weight = StdPyrData::getNumberOfCoefficients(na, nb, nc);
weight = StdPyrData::getNumberOfCoefficients(na, nb, nc);
bndWeight = StdPyrData::getNumberOfBndCoefficients(na, nb, nc);
break;
case 'Q':
weight = StdQuadData::getNumberOfCoefficients(na, nb);
weight = StdQuadData::getNumberOfCoefficients(na, nb);
bndWeight = StdQuadData::getNumberOfBndCoefficients(na, nb);
break;
case 'T':
weight = StdTriData::getNumberOfCoefficients(na, nb);
weight = StdTriData::getNumberOfCoefficients(na, nb);
bndWeight = StdTriData::getNumberOfBndCoefficients(na, nb);
break;
case 'S':
weight = StdSegData::getNumberOfCoefficients(na);
bndWeight = StdSegData::getNumberOfBndCoefficients(na);
break;
default:
break;
......@@ -654,8 +661,12 @@ namespace Nektar
" does not correspond to mesh dimension");
int na = it->second[0];
int nb = it->second[1];
int nb = 0;
int nc = 0;
if (m_dim >= 2)
{
nb = it->second[1];
}
if (m_dim == 3)
{
nc = it->second[2];
......@@ -666,25 +677,32 @@ namespace Nektar
switch (m_meshComposites[cId].type)
{
case 'A':
weight = StdTetData::getNumberOfCoefficients(na, nb, nc);
bndWeight = 2 * (na-2) * (na-1) + 6 * (na-1) + 4;
weight = StdTetData::getNumberOfCoefficients(na, nb, nc);
bndWeight = StdTetData::getNumberOfBndCoefficients(na, nb, nc);
break;
case 'R':
weight = StdPrismData::getNumberOfCoefficients(na, nb, nc);
bndWeight = (na-2)*(na-1) + 3*(na-1)*(na-1) + 9*(na-1) + 6;
weight = StdPrismData::getNumberOfCoefficients(na, nb, nc);
bndWeight = StdPrismData::getNumberOfBndCoefficients(na, nb, nc);
break;
case 'H':
weight = StdHexData::getNumberOfCoefficients(na, nb, nc);
bndWeight = 6*(na-1)*(na-1) + 12*(na-1) + 8;
weight = StdHexData::getNumberOfCoefficients(na, nb, nc);
bndWeight = StdHexData::getNumberOfBndCoefficients(na, nb, nc);
break;
case 'P':
weight = StdPyrData::getNumberOfCoefficients(na, nb, nc);
weight = StdPyrData::getNumberOfCoefficients(na, nb, nc);
bndWeight = StdPyrData::getNumberOfBndCoefficients(na, nb, nc);
break;
case 'Q':
weight = StdQuadData::getNumberOfCoefficients(na, nb);
weight = StdQuadData::getNumberOfCoefficients(na, nb);
bndWeight = StdQuadData::getNumberOfBndCoefficients(na, nb);
break;
case 'T':
weight = StdTriData::getNumberOfCoefficients(na, nb);
weight = StdTriData::getNumberOfCoefficients(na, nb);
bndWeight = StdTriData::getNumberOfBndCoefficients(na, nb);
break;
case 'S':
weight = StdSegData::getNumberOfCoefficients(na);
bndWeight = StdSegData::getNumberOfBndCoefficients(na);
break;
default:
break;
......
......@@ -92,23 +92,56 @@ namespace Nektar
};
namespace StdSegData
{
inline int getNumberOfCoefficients(int Na)
{
return Na;
}
inline int getNumberOfBndCoefficients(int Na)
{
return 2;
}
}
// Dimensions of coefficients for each space
namespace StdTriData
{
inline int getNumberOfCoefficients(int Na, int Nb)
{
ASSERTL0(Na <= Nb, "order in 'a' direction is higher "
ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL1(Na <= Nb, "order in 'a' direction is higher "
"than order in 'b' direction");
return Na*(Na+1)/2 + Na*(Nb-Na);
}
inline int getNumberOfBndCoefficients(int Na, int Nb)
{
ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL1(Na <= Nb, "order in 'a' direction is higher "
"than order in 'b' direction");
return (Na-1) + 2*(Nb-1);
}
}
namespace StdQuadData
{
inline int getNumberOfCoefficients(int Na, int Nb)
{
ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
return Na*Nb;
}
inline int getNumberOfBndCoefficients(int Na, int Nb)
{
ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
return 2*(Na-1) + 2*(Nb-1);
}
}
......@@ -117,8 +150,20 @@ namespace Nektar
{
inline int getNumberOfCoefficients( int Na, int Nb, int Nc )
{
ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL2(Nc > 1, "Order in 'c' direction must be > 1.");
return Na*Nb*Nc;
}
inline int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
{
ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL2(Nc > 1, "Order in 'c' direction must be > 1.");
return 2*Na*Nb + 2*Na*Nc + 2*Nb*Nc
- 4*(Na + Nb + Nc) + 8;
}
}
namespace StdTetData
......@@ -140,6 +185,13 @@ namespace Nektar
*/
inline int getNumberOfCoefficients(int Na, int Nb, int Nc)
{
ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL2(Nc > 1, "Order in 'c' direction must be > 1.");
ASSERTL1(Na <= Nb, "order in 'a' direction is higher "
"than order in 'b' direction");
ASSERTL1(Nb <= Nc, "order in 'b' direction is higher "
"than order in 'c' direction");
int nCoef = 0;
for (int a = 0; a < Na; ++a)
{
......@@ -153,6 +205,25 @@ namespace Nektar
}
return nCoef;
}
inline int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
{
ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL2(Nc > 1, "Order in 'c' direction must be > 1.");
ASSERTL1(Na <= Nb, "order in 'a' direction is higher "
"than order in 'b' direction");
ASSERTL1(Nb <= Nc, "order in 'b' direction is higher "
"than order in 'c' direction");
int nCoef = Na*(Na+1)/2 + (Nb-Na)*Na // base
+ Na*(Na+1)/2 + (Nc-Na)*Na // front
+ 2*(Nb*(Nb+1)/2 + (Nc-Nb)*Nb)// 2 other sides
- Na - 2*Nb - 3*Nc // less edges
+ 4; // plus vertices
return nCoef;
}
}
......@@ -160,6 +231,13 @@ namespace Nektar
{
inline int getNumberOfCoefficients(int Na, int Nb, int Nc)
{
ASSERTL1(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL1(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL1(Nc > 1, "Order in 'c' direction must be > 1.");
ASSERTL1(Na <= Nc, "Order in 'a' direction is higher "
"than order in 'c' direction.");
ASSERTL1(Nb <= Nc, "Order in 'b' direction is higher "
"than order in 'c' direction.");
int nCoef = 0;
for (int c = 0; c < Nc; ++c)
{
......@@ -171,31 +249,54 @@ namespace Nektar
}
}
}
/*
for (int a = 0; a < Na; ++a)
{
for (int b = 0; b < Nb; ++b)
{
for (int c = 0; c < Nc - a - b; ++c)
{
++nCoef;
}
}
}
*/
//std::cout << "Na = " << Na << " Nb = " << Nb << " Nc = " << Nc << " nCoef = " << nCoef << std::endl;
return nCoef;
}
inline int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
{
ASSERTL1(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL1(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL1(Nc > 1, "Order in 'c' direction must be > 1.");
ASSERTL1(Na <= Nc, "Order in 'a' direction is higher "
"than order in 'c' direction.");
ASSERTL1(Nb <= Nc, "Order in 'b' direction is higher "
"than order in 'c' direction.");
return Na*Nb // base
+ 2*(Na*(Na+1)/2 + (Nc-Na)*Na) // front and back
+ 2*(Nb*(Nb+1)/2 + (Nc-Nb)*Nb) // sides
- 2*Na - 2*Nb - 4*Nc // less edges
+ 5; // plus vertices
}
}
namespace StdPrismData
{
inline int getNumberOfCoefficients( int Na, int Nb, int Nc )
{
ASSERTL1(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL1(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL1(Nc > 1, "Order in 'c' direction must be > 1.");
ASSERTL1(Na <= Nc, "Order in 'a' direction is higher "
"than order in 'c' direction.");
return Nb*StdTriData::getNumberOfCoefficients(Na,Nc);
}
}
inline int getNumberOfBndCoefficients( int Na, int Nb, int Nc)
{
ASSERTL1(Na > 1, "Order in 'a' direction must be > 1.");
ASSERTL1(Nb > 1, "Order in 'b' direction must be > 1.");
ASSERTL1(Nc > 1, "Order in 'c' direction must be > 1.");
ASSERTL1(Na <= Nc, "Order in 'a' direction is higher "
"than order in 'c' direction.");
return Na*Nb + 2*Nb*Nc // rect faces
+ 2*( Na*(Na+1)/2 + (Nc-Na)*Na ) // tri faces
- 2*Na - 3*Nb - 4*Nc // less edges
+ 6; // plus vertices
}
}
inline int GetNumberOfCoefficients(ShapeType shape, std::vector<unsigned int> &modes, int offset)
{
......
......@@ -34,6 +34,7 @@
///////////////////////////////////////////////////////////////////////////////
#include <StdRegions/StdPrismExp.h>
#include <LibUtilities/BasicUtils/ShapeType.hpp>
namespace Nektar
{
......@@ -944,14 +945,12 @@ namespace Nektar
GetBasisType(2) == LibUtilities::eGLL_Lagrange,
"BasisType is not a boundary interior form");
int P = m_base[0]->GetNumModes() - 1;
int Q = m_base[1]->GetNumModes() - 1;
int R = m_base[2]->GetNumModes() - 1;
int P = m_base[0]->GetNumModes();
int Q = m_base[1]->GetNumModes();
int R = m_base[2]->GetNumModes();
return (P+1)*(Q+1) + 2*(Q+1)*(R+1) // 3 rect. faces in p-q and q-r planes
+ 2*(R+1) + P*(1 + 2*R - P) // 2 tri. faces in p-r plane
- 3*(Q+1) - 2*(P+1) - 4*(R+1) // subtract double counted points on edge
+ 6; // add vertices
return LibUtilities::StdPrismData::
getNumberOfBndCoefficients(P,Q,R);
}
int StdPrismExp::v_NumDGBndryCoeffs() const
......
......@@ -733,15 +733,12 @@ namespace Nektar
GetBasisType(2) == LibUtilities::eGLL_Lagrange,
"BasisType is not a boundary interior form");
int P = m_base[0]->GetNumModes() - 1;
int Q = m_base[1]->GetNumModes() - 1;
int R = m_base[2]->GetNumModes() - 1;
int P = m_base[0]->GetNumModes();
int Q = m_base[1]->GetNumModes();
int R = m_base[2]->GetNumModes();
return (P+1)*(Q+1) // 1 rect. face in p-q plane
+ 2*(R+1) + P*(1+2*R-P) // 2 tri. faces in p-r plane
+ 2*(R+1) + Q*(1+2*R-Q) // 2 tri. faces in q-r plane
- 2*(P+1)-2*(Q+1)-4*(R+1) // subtract double counted edges
+ 5; // add vertices
return LibUtilities::StdPyrData::
getNumberOfBndCoefficients(P, Q, R);
}
int StdPyrExp::v_GetEdgeNcoeffs(const int i) const
......
......@@ -1160,20 +1160,8 @@ namespace Nektar
int Q = m_base[1]->GetNumModes();
int R = m_base[2]->GetNumModes();
int p_hat, k;
// All modes in the first layer are boundary modes
int tot = P*(P+1)/2 + (Q-P)*P;
// Loop over each plane in the stack
for (int i = 1; i < R - 1; ++i)
{
p_hat = min(P, R-i);
k = min(Q-P, max(0, Q-i-1));
// First two columns and bottom row are boundary modes
tot += (p_hat + k) + (p_hat + k - 1) + p_hat - 2;
}
// Add on top vertex mode
return tot + 1;
return LibUtilities::StdTetData::
getNumberOfBndCoefficients(P, Q, R);
}
int StdTetExp::v_NumDGBndryCoeffs() const
......
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