Commit 31d2b68a authored by Spencer Sherwin's avatar Spencer Sherwin

Merge branch 'master' into feature/AvgFilterMultiOutput

parents bdbf0976 2c7169fd
......@@ -14,7 +14,6 @@ ADD_NEKTAR_EXECUTABLE(LocProject1D demos LocProject1DSource)
TARGET_LINK_LIBRARIES(LocProject1D ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocProject1D)
SET(LocProject2DSource LocProject2D.cpp)
ADD_NEKTAR_EXECUTABLE(LocProject2D demos LocProject2DSource)
TARGET_LINK_LIBRARIES(LocProject2D ${LinkLibraries})
......@@ -25,7 +24,6 @@ ADD_NEKTAR_EXECUTABLE(LocProject3D demos LocProject3DSource)
TARGET_LINK_LIBRARIES(LocProject3D ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocProject3D)
SET(LocProject_Diff1DSource LocProject_Diff1D.cpp)
ADD_NEKTAR_EXECUTABLE(LocProject_Diff1D demos LocProject_Diff1DSource)
TARGET_LINK_LIBRARIES(LocProject_Diff1D ${LinkLibraries})
......@@ -41,46 +39,6 @@ ADD_NEKTAR_EXECUTABLE(LocProject_Diff3D demos LocProject_Diff3DSource)
TARGET_LINK_LIBRARIES(LocProject_Diff3D ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocProject_Diff3D)
SET(LocTetExpDemoSource LocTetExpDemo.cpp)
ADD_NEKTAR_EXECUTABLE(LocTetExpDemo demos LocTetExpDemoSource)
TARGET_LINK_LIBRARIES(LocTetExpDemo ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocTetExpDemo)
SET(LocPyrExpDemoSource LocPyrExpDemo.cpp)
ADD_NEKTAR_EXECUTABLE(LocPyrExpDemo demos LocPyrExpDemoSource)
TARGET_LINK_LIBRARIES(LocPyrExpDemo ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocPyrExpDemo)
SET(LocPrismExpDemoSource LocPrismExpDemo.cpp)
ADD_NEKTAR_EXECUTABLE(LocPrismExpDemo demos LocPrismExpDemoSource)
TARGET_LINK_LIBRARIES(LocPrismExpDemo ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocPrismExpDemo)
SET(LocHexExpDemoSource LocHexExpDemo.cpp)
ADD_NEKTAR_EXECUTABLE(LocHexExpDemo demos LocHexExpDemoSource)
TARGET_LINK_LIBRARIES(LocHexExpDemo ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocHexExpDemo)
SET(LocHexExpDemoSource LocHexExpDemo_2.cpp)
ADD_NEKTAR_EXECUTABLE(LocHexExpDemo_2 demos LocHexExpDemoSource)
TARGET_LINK_LIBRARIES(LocHexExpDemo_2 ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocHexExpDemo_2)
SET(LocPrismExpDemoSource LocPrismExpDemo_2.cpp)
ADD_NEKTAR_EXECUTABLE(LocPrismExpDemo_2 demos LocPrismExpDemoSource)
TARGET_LINK_LIBRARIES(LocPrismExpDemo_2 ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocPrismExpDemo_2)
SET(LocPyrExpDemoSource LocPyrExpDemo_2.cpp)
ADD_NEKTAR_EXECUTABLE(LocPyrExpDemo_2 demos LocPyrExpDemoSource)
TARGET_LINK_LIBRARIES(LocPyrExpDemo_2 ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocPyrExpDemo_2)
SET(LocTetExpDemoSource LocTetExpDemo_2.cpp)
ADD_NEKTAR_EXECUTABLE(LocTetExpDemo_2 demos LocTetExpDemoSource)
TARGET_LINK_LIBRARIES(LocTetExpDemo_2 ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocTetExpDemo_2)
# Generate list of available subdirectories
FILE(GLOB dir_list "*")
FOREACH(dir ${dir_list})
......@@ -108,6 +66,7 @@ ADD_NEKTAR_TEST(LocProject3D_Tet_Ortho_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject3D_Tet_Mod_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject3D_Prism_Ortho_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject3D_Prism_Mod_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject3D_Pyr_Mod_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject3D_Hex_Ortho_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject3D_Hex_Mod_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject3D_Hex_Lagrange_Basis_P6_Q7)
......@@ -123,6 +82,7 @@ ADD_NEKTAR_TEST(LocProject_Diff3D_Reg_Hex_Mod_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject_Diff3D_Reg_Hex_Ortho_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject_Diff3D_Reg_Prism_Mod_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject_Diff3D_Reg_Prism_Ortho_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject_Diff3D_Reg_Pyr_Mod_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject_Diff3D_Reg_Tet_Mod_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject_Diff3D_Reg_Tet_Ortho_Basis_P6_Q7)
......
This diff is collapsed.
#include <cstdlib>
#include <LocalRegions/HexExp.h>
#include <SpatialDomains/MeshGraph3D.h>
#include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
using namespace Nektar;
/// modification to deal with exact solution. Return 1 if integer < 0
static double pow_loc(const double val, const int i)
{
return (i < 0)? 1.0: pow(val,i);
}
/// Solution to excite all modes in Hex expansion.
NekDouble Hex_sol( NekDouble x, NekDouble y, NekDouble z,
int order1, int order2, int order3,
LibUtilities::BasisType btype1,
LibUtilities::BasisType btype2,
LibUtilities::BasisType btype3);
/// Derivative of solution in Hex expansion.
NekDouble Hex_Diff_Sol( NekDouble x, NekDouble y, NekDouble z,
int P, int Q, int R,
LibUtilities::BasisType btype1,
LibUtilities::BasisType btype2,
LibUtilities::BasisType btype3, int dir);
int main(int argc, char *argv[])
{
if(argc != 8)
{
cerr << "usage: LocHexExpDemo MeshFile nummodes0 nummodes1 nummodes2 "
"nx ny nz" << endl;
exit(1);
}
int i;
string in(argv[1]);
int nummodes0 = atoi(argv[2]);
int nummodes1 = atoi(argv[3]);
int nummodes2 = atoi(argv[4]);
int nquad0 = atoi(argv[5]);
int nquad1 = atoi(argv[6]);
int nquad2 = atoi(argv[7]);
int P = nummodes0 - 1;
int Q = nummodes1 - 1;
int R = nummodes2 - 1;
int ntotquad = nquad0*nquad1*nquad2;
SpatialDomains::MeshGraph3D graph3D;
graph3D.ReadGeometry(in);
LibUtilities::BasisType bType = LibUtilities::eModified_A;
LibUtilities::PointsType qtype = LibUtilities::eGaussLobattoLegendre;
const LibUtilities::PointsKey pointsKey0( nquad0, qtype );
const LibUtilities::PointsKey pointsKey1( nquad1, qtype );
const LibUtilities::PointsKey pointsKey2( nquad2, qtype );
const LibUtilities::BasisKey basisKey0( bType, nummodes0, pointsKey0 );
const LibUtilities::BasisKey basisKey1( bType, nummodes1, pointsKey1 );
const LibUtilities::BasisKey basisKey2( bType, nummodes2, pointsKey2 );
SpatialDomains::HexGeomSharedPtr geom;
if(!(geom = boost::dynamic_pointer_cast<SpatialDomains::HexGeom>(
graph3D.GetCompositeItem(0,0))))
{
cerr << "Could not find HexGeom in input file" << endl;
exit(1);
}
LocalRegions::HexExpSharedPtr E = MemoryManager<LocalRegions::HexExp>
::AllocateSharedPtr(basisKey0,basisKey1,basisKey2,geom);
Array<OneD, NekDouble> sol( ntotquad );
Array<OneD, NekDouble> dx( ntotquad, 0.0 );
Array<OneD, NekDouble> dy( ntotquad, 0.0 );
Array<OneD, NekDouble> dz( ntotquad, 0.0 );
Array<OneD, NekDouble> diff_solution_x( ntotquad, 0.0 );
Array<OneD, NekDouble> diff_solution_y( ntotquad, 0.0 );
Array<OneD, NekDouble> diff_solution_z( ntotquad, 0.0 );
Array<OneD, NekDouble> derivatives( ntotquad, 0.0 );
Array<OneD, NekDouble> x( ntotquad );
Array<OneD, NekDouble> y( ntotquad );
Array<OneD, NekDouble> z( ntotquad );
E->GetCoords(x,y,z);
//----------------------------------------------
// Define solution to be projected
for(i = 0; i < ntotquad; i++)
{
sol[i] = Hex_sol(x[i],y[i],z[i],P,Q,R, bType, bType, bType);
}
//----------------------------------------------
//----------------------------------------------
// Define the derivative solution
for(int n = 0; n < ntotquad; ++n){
diff_solution_x[n] = Hex_Diff_Sol( x[n], y[n], z[n], P,Q,R,
bType, bType, bType, 1);
diff_solution_y[n] = Hex_Diff_Sol( x[n], y[n], z[n], P,Q,R,
bType, bType, bType, 2);
diff_solution_z[n] = Hex_Diff_Sol( x[n], y[n], z[n], P,Q,R,
bType, bType, bType, 3);
}
//---------------------------------------------
//---------------------------------------------
// Project onto Expansion
Array<OneD, NekDouble> coeffs(E->GetNcoeffs());
Array<OneD, NekDouble> phys (ntotquad);
E->FwdTrans( sol, coeffs );
//---------------------------------------------
//-------------------------------------------
// Backward Transform Solution to get projected values
E->BwdTrans( coeffs, phys );
//-------------------------------------------
//--------------------------------------------
// Calculate L_p error
cout << "L infinity error: " << E->Linf(phys, sol) << endl;
cout << "L 2 error: " << E->L2 (phys, sol) << endl;
//--------------------------------------------
//-------------------------------------------
// Evaulate solution at x = y = z = 0 and print error
Array<OneD, NekDouble> t = Array<OneD, NekDouble>(3);
t[0] = -0.39;
t[1] = -0.25;
t[2] = 0.5;
NekDouble exact_sol = 0.0;
exact_sol = Hex_sol( t[0], t[1], t[2], nummodes0, nummodes1, nummodes2,
bType, bType, bType);
NekDouble numericSolution = E->PhysEvaluate(t, phys);
cout << "Numeric solution = " << numericSolution << endl;
cout << "Exact solution = " << exact_sol << endl;
cout << "Difference at x = ( "
<< t[0] << ", " << t[1] << ", " << t[2] << " ): "
<< numericSolution - exact_sol << endl;
//-------------------------------------------
E->PhysDeriv( sol, dx, dy, dz);
double error_x = 0, error_y=0, error_z=0;
for( int n = 0; n < ntotquad; ++n ) {
error_x += fabs(diff_solution_x[n] - dx[n]);
error_y += fabs(diff_solution_y[n] - dy[n]);
error_z += fabs(diff_solution_z[n] - dz[n]);
}
cout << "L 1 error of derivatives X = " << error_x << endl;
cout << "L 1 error of derivatives Y = " << error_y << endl;
cout << "L 1 error of derivatives Z = " << error_z << endl;
return 0;
}
NekDouble Hex_sol(NekDouble x, NekDouble y, NekDouble z,
int order1, int order2, int order3,
LibUtilities::BasisType btype1,
LibUtilities::BasisType btype2,
LibUtilities::BasisType btype3)
{
int i,j,k;
NekDouble sol = 0.0;
int Nx = (btype1 == LibUtilities::eFourier ? order1/2 : order1);
int Ny = (btype2 == LibUtilities::eFourier ? order2/2 : order2);
int Nz = (btype3 == LibUtilities::eFourier ? order3/2 : order3);
bool Fx = (btype1 == LibUtilities::eFourier);
bool Fy = (btype2 == LibUtilities::eFourier);
bool Fz = (btype3 == LibUtilities::eFourier);
NekDouble a;
for (i = 0; i < Nx; ++i)
{
for (j = 0; j < Ny; ++j)
{
for (k = 0; k < Nz; ++k)
{
a = (Fx ? sin(M_PI*i*x) + cos(M_PI*i*x) : pow_loc(x,i));
a *= (Fy ? sin(M_PI*j*y) + cos(M_PI*j*y) : pow_loc(y,j));
a *= (Fz ? sin(M_PI*k*z) + cos(M_PI*k*z) : pow_loc(z,k));
sol += a;
}
}
}
return sol;
}
NekDouble Hex_Diff_Sol(NekDouble x, NekDouble y, NekDouble z,
int order1, int order2, int order3,
LibUtilities::BasisType btype1,
LibUtilities::BasisType btype2,
LibUtilities::BasisType btype3, int dir){
int i,j,k;
NekDouble sol = 0.0;
int Nx = (btype1 == LibUtilities::eFourier ? order1/2 : order1);
int Ny = (btype2 == LibUtilities::eFourier ? order2/2 : order2);
int Nz = (btype3 == LibUtilities::eFourier ? order3/2 : order3);
bool Fx = (btype1 == LibUtilities::eFourier);
bool Fy = (btype2 == LibUtilities::eFourier);
bool Fz = (btype3 == LibUtilities::eFourier);
NekDouble a;
for (i = 0; i < Nx; ++i)
{
for (j = 0; j < Ny; ++j)
{
for (k = 0; k < Nz; ++k)
{
switch (dir)
{
case 1:
a = i*(Fx ? M_PI*(cos(M_PI*i*x) - sin(M_PI*i*x))
: pow_loc(x,i-1));
a *= (Fy ? (sin(M_PI*j*y) + cos(M_PI*j*y)) : pow_loc(y,j));
a *= (Fz ? (sin(M_PI*k*z) + cos(M_PI*k*z)) : pow_loc(z,k));
sol += a;
break;
case 2:
a = (Fx ? (sin(M_PI*i*x) + cos(M_PI*i*x)) : pow_loc(x,i));
a *= j*(Fy ? M_PI*(cos(M_PI*j*y) - sin(M_PI*j*y))
: pow_loc(y,j-1));
a *= (Fz ? (sin(M_PI*k*z) + cos(M_PI*k*z)) : pow_loc(z,k));
sol += a;
break;
case 3:
a = (Fx ? (sin(M_PI*i*x) + cos(M_PI*i*x)) : pow_loc(x,i));
a *= (Fy ? (sin(M_PI*j*y) + cos(M_PI*j*y)) : pow_loc(y,j));
a *= k*(Fz ? M_PI*(cos(M_PI*k*z) - sin(M_PI*k*z))
: pow_loc(z,k-1));
sol += a;
break;
}
}
}
}
return sol;
}
This diff is collapsed.
#include <cstdlib>
#include <LocalRegions/PrismExp.h>
#include <SpatialDomains/MeshGraph3D.h>
#include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
using namespace Nektar;
NekDouble Prism_sol(NekDouble x, NekDouble y, NekDouble z, int order1, int order2, int order3);
NekDouble Prism_sol2(NekDouble x, NekDouble y, NekDouble z, int P, int Q, int R,
LibUtilities::BasisType bType_x,
LibUtilities::BasisType bType_y,
LibUtilities::BasisType bType_z );
int main(int argc, char *argv[])
{
if(argc != 5)
{
cerr << "usage: LocPrismExpDemo MeshFile nummodes0 nummodes1 nummodes2" << endl;
exit(1);
}
int i;
string in(argv[1]);
int nummodes0 = atoi(argv[2]);
int nummodes1 = atoi(argv[3]);
int nummodes2 = atoi(argv[4]);
int nquad0 = nummodes0 + 1;
int nquad1 = nummodes1 + 1;
int nquad2 = nummodes2 + 1;
int ntotquad = nquad0*nquad1*nquad2;
SpatialDomains::MeshGraph3D graph3D;
graph3D.ReadGeometry(in);
LibUtilities::BasisType bTypeX = LibUtilities::eModified_A;
LibUtilities::BasisType bTypeY = LibUtilities::eModified_A;
LibUtilities::BasisType bTypeZ = LibUtilities::eModified_B;
LibUtilities::PointsType qtypeX = LibUtilities::eGaussLobattoLegendre;
LibUtilities::PointsType qtypeY = LibUtilities::eGaussLobattoLegendre;
LibUtilities::PointsType qtypeZ = LibUtilities::eGaussRadauMAlpha1Beta0;
const LibUtilities::PointsKey pointsKey0( nquad0, qtypeX );
const LibUtilities::PointsKey pointsKey1( nquad1, qtypeY );
const LibUtilities::PointsKey pointsKey2( nquad2, qtypeZ );
const LibUtilities::BasisKey basisKey0( bTypeX, nummodes0, pointsKey0 );
const LibUtilities::BasisKey basisKey1( bTypeY, nummodes1, pointsKey1 );
const LibUtilities::BasisKey basisKey2( bTypeZ, nummodes2, pointsKey2 );
SpatialDomains::PrismGeomSharedPtr geom;
if(!(geom = boost::dynamic_pointer_cast<SpatialDomains::PrismGeom>(graph3D.GetCompositeItem(0,0))))
{
cerr << "Could not find PrismGeom in input file" << endl;
exit(1);
}
LocalRegions::PrismExpSharedPtr E = MemoryManager<LocalRegions::PrismExp>::
AllocateSharedPtr(basisKey0,basisKey1,basisKey2,geom);
Array<OneD,NekDouble> x(ntotquad);
Array<OneD,NekDouble> y(ntotquad);
Array<OneD,NekDouble> z(ntotquad);
E->GetCoords(x,y,z);
Array<OneD,NekDouble> sol(ntotquad);
//----------------------------------------------
// Define solution to be projected
for(i = 0; i < ntotquad; i++)
{
sol[i] = Prism_sol(x[i],y[i],z[i],nummodes0,nummodes1,nummodes2);
}
//----------------------------------------------
//---------------------------------------------
// Project onto Expansion
Array<OneD, NekDouble> coeffs(E->GetNcoeffs());
Array<OneD, NekDouble> phys (ntotquad);
E->FwdTrans( sol, coeffs );
//---------------------------------------------
//-------------------------------------------
// Backward Transform Solution to get projected values
E->BwdTrans( coeffs, phys );
//-------------------------------------------
//--------------------------------------------
// Calculate L_p error
cout << "L infinity error: " << E->Linf(phys, sol) << endl;
cout << "L 2 error: " << E->L2 (phys, sol) << endl;
//--------------------------------------------
return 0;
}
NekDouble Prism_sol(NekDouble x, NekDouble y, NekDouble z, int nummodes1, int nummodes2, int nummodes3)
{
NekDouble sol = 0;
for(int i = 0; i < nummodes1; ++i)
{
for(int j = 0; j < nummodes2; ++j)
{
for(int k = 0; k < nummodes3; ++k)
{
sol += pow(x,i) * pow(y,j) * pow(z,k);
}
}
}
return sol;
}
NekDouble Prism_sol2(NekDouble x, NekDouble y, NekDouble z, int P, int Q, int R,
LibUtilities::BasisType bType_x,
LibUtilities::BasisType bType_y,
LibUtilities::BasisType bType_z )
{
NekDouble sol = 0;
// case 1 -- Common case
if( (bType_x != LibUtilities::eFourier) && (bType_y != LibUtilities::eFourier) && (bType_z != LibUtilities::eFourier) )
{
for(int p = 0; p <= P; ++p) {
for(int q = 0; q <= Q; ++q) {
for(int r = 0; r <= R - p; ++r) {
sol += pow(x,p) * pow(y,q) * pow(z,r);
}
}
}
} else // case 2
if((bType_x != LibUtilities::eFourier) && (bType_y != LibUtilities::eFourier) && (bType_z == LibUtilities::eFourier))
{
for(int i = 0; i <= P; ++i) {
for(int j = 0; j <= Q; ++j) {
for(int k = 0; k <= R/2; ++k) {
sol += pow(x,i) * pow(y,j) * sin(M_PI*k*z) + pow(x,i) * pow(y,j) * cos(M_PI*k*z);
}
}
}
}else // case 3
if((bType_x == LibUtilities::eFourier) && (bType_y != LibUtilities::eFourier) && (bType_z != LibUtilities::eFourier))
{
for(int i = 0; i <= P/2; ++i) {
for(int j = 0; j <= Q; ++j) {
for(int k = 0; k <= R - i; ++k) {
sol += sin(M_PI*i*x)* pow(y,j) * pow(z,k) + cos(M_PI*i*x)* pow(y,j) * pow(z,k);
}
}
}
}
else // case 4
if((bType_x != LibUtilities::eFourier) && (bType_y == LibUtilities::eFourier) && (bType_z != LibUtilities::eFourier))
{
for(int i = 0; i <= P; ++i) {
for(int j = 0; j <= Q/2; ++j) {
for(int k = 0; k <= R - i; ++k) {
sol += pow(x,i)*sin(M_PI*j*y)*pow(z,k) + pow(x,i)*cos(M_PI*j*y)*pow(z,k);
}
}
}
}
return sol;
}
......@@ -4,6 +4,7 @@
#include <StdRegions/StdExpansion3D.h>
#include <LocalRegions/HexExp.h>
#include <LocalRegions/PrismExp.h>
#include <LocalRegions/PyrExp.h>
#include <LocalRegions/TetExp.h>
#include <LibUtilities/Foundations/Foundations.hpp>
......@@ -39,6 +40,9 @@ SpatialDomains::HexGeomSharedPtr CreateHexGeom(int argc, char *argv[]);
/// Generates a Hex geometry based on the command-line parameters.
SpatialDomains::PrismGeomSharedPtr CreatePrismGeom(int argc, char *argv[]);
/// Generates a Pyramid geometry based on the command-line parameters.
SpatialDomains::PyrGeomSharedPtr CreatePyrGeom(int argc, char *argv[]);
/// Generates a Tet geometry based on the command-line parameters.
SpatialDomains::TetGeomSharedPtr CreateTetGeom(int argc, char *argv[]);
......@@ -48,7 +52,6 @@ static double pow_loc(const double val, const int i)
return (i < 0)? 1.0: pow(val,i);
}
/// This routine projects a polynomial or trigonmetric functions which
/// has energy in all modes of the expansions and reports and error
int main(int argc, char *argv[]){
......@@ -70,6 +73,7 @@ int main(int argc, char *argv[]){
fprintf(stderr,"Where RegionShape is an integer value which "
"dictates the region shape:\n");
fprintf(stderr,"\t Tetrahedron = 5\n");
fprintf(stderr,"\t Pyramid = 6\n");
fprintf(stderr,"\t Prism = 7\n");
fprintf(stderr,"\t Hexahedron = 8\n");
......@@ -95,6 +99,7 @@ int main(int argc, char *argv[]){
// Check to see if 3D region
if (regionshape != LibUtilities::eTetrahedron &&
regionshape != LibUtilities::ePyramid &&
regionshape != LibUtilities::ePrism &&
regionshape != LibUtilities::eHexahedron)
{
......@@ -134,6 +139,29 @@ int main(int argc, char *argv[]){
"or Modified_B");
}
break;
case LibUtilities::ePyramid:
if((btype1 == eOrtho_B) || (btype1 == eOrtho_C)
|| (btype1 == eModified_B) || (btype1 == eModified_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))
{
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))
{
NEKERROR(ErrorUtil::efatal,