Commit 24f44507 authored by Dave Moxey's avatar Dave Moxey
Browse files

Add LocalRegions pyramid regression tests

parent 30ce2b35
......@@ -66,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)
......@@ -81,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)
......
......@@ -3,6 +3,7 @@
#include <StdRegions/StdExpansion3D.h>
#include <LocalRegions/HexExp.h>
#include <LocalRegions/PyrExp.h>
#include <LocalRegions/TetExp.h>
#include <LocalRegions/PrismExp.h>
......@@ -54,6 +55,9 @@ SpatialDomains::HexGeomSharedPtr CreateHexGeom(int argc, char *argv[]);
/// Creates a Prism 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[]);
/// Creates a Tet geometry based on the command-line parameters.
SpatialDomains::TetGeomSharedPtr CreateTetGeom(int argc, char *argv[]);
......@@ -86,6 +90,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 Prism = 6\n");
fprintf(stderr,"\t Prism = 7\n");
fprintf(stderr,"\t Hexahedron = 8\n");
......@@ -110,6 +115,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)
{
......@@ -150,6 +156,29 @@ int main(int argc, char *argv[]){
"Ortho_B, Modified_A 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,
"Basis 3 cannot be of type Ortho_A, "
"Ortho_B, Modified_A or Modified_B");
}
break;
case LibUtilities::ePrism:
if((btype1 == eOrtho_B) || (btype1 == eOrtho_C)
|| (btype1 == eModified_B) || (btype1 == eModified_C))
......@@ -234,7 +263,8 @@ int main(int argc, char *argv[]){
if(btype3 != LibUtilities::eFourier)
{
if (regionshape == LibUtilities::eTetrahedron)
if (regionshape == LibUtilities::eTetrahedron ||
regionshape == LibUtilities::ePyramid)
{
Qtype3 = LibUtilities::eGaussRadauMAlpha2Beta0;
}
......@@ -283,6 +313,29 @@ int main(int argc, char *argv[]){
//----------------------------------------------
}
break;
case LibUtilities::ePyramid:
{
const LibUtilities::PointsKey Pkey1(nq1,Qtype1);
const LibUtilities::PointsKey Pkey2(nq2,Qtype2);
const LibUtilities::PointsKey Pkey3(nq3,Qtype3);
const LibUtilities::BasisKey Bkey1(btype1,order1,Pkey1);
const LibUtilities::BasisKey Bkey2(btype2,order2,Pkey2);
const LibUtilities::BasisKey Bkey3(btype3,order3,Pkey3);
SpatialDomains::PyrGeomSharedPtr geom = CreatePyrGeom(argc, argv);
E = new LocalRegions::PyrExp(Bkey1, Bkey2, Bkey3, geom);
E->GetCoords(x,y,z);
//----------------------------------------------
// Define solution to be projected
for(i = 0; i < nq1*nq2*nq3; ++i)
{
sol[i] = Tet_sol(x[i],y[i],z[i],order1,order2,order3);
}
//----------------------------------------------
}
break;
case LibUtilities::ePrism:
{
const LibUtilities::PointsKey Pkey1(nq1,Qtype1);
......@@ -359,6 +412,7 @@ int main(int argc, char *argv[]){
switch(regionshape)
{
case LibUtilities::eTetrahedron:
case LibUtilities::ePyramid:
{
for(i = 0; i < nq1*nq2*nq3; ++i)
{
......@@ -757,6 +811,111 @@ SpatialDomains::PrismGeomSharedPtr CreatePrismGeom(int argc, char *argv[])
return geom;
}
SpatialDomains::PyrGeomSharedPtr CreatePyrGeom(int argc, char *argv[])
{
if (argc != 26)
{
cout << "Insufficient points for a pyramid!" << endl;
}
// /////////////////////////////////////////////////////////////////////
// Set up Pyramid vertex coordinates
// PointGeom (const int coordim, const int vid, double x,
// double y, double z)
const int nVerts = 6;
const double point[][3] = {
{atof(argv[11]),atof(argv[12]),atof(argv[13])},
{atof(argv[14]),atof(argv[15]),atof(argv[16])},
{atof(argv[17]),atof(argv[18]),atof(argv[19])},
{atof(argv[20]),atof(argv[21]),atof(argv[22])},
{atof(argv[23]),atof(argv[24]),atof(argv[25])}
};
// Populate the list of verts
PointGeomSharedPtr verts[nVerts];
const int three = 3;
for( int i = 0; i < nVerts; ++i ) {
verts[i] = MemoryManager<PointGeom>
::AllocateSharedPtr( three, i, point[i][0],
point[i][1], point[i][2] );
}
// /////////////////////////////////////////////////////////////////////
// Set up Pyramid Edges
// SegGeom (int id, const int coordim), EdgeComponent(id, coordim)
const int nEdges = 9;
const int vertexConnectivity[][2] = {
{0,1}, {1,2}, {3,2}, {0,3},
{0,4}, {1,4}, {2,4}, {3,4}
};
// Populate the list of edges
SegGeomSharedPtr edges[nEdges];
for( int i = 0; i < nEdges; ++i ) {
PointGeomSharedPtr vertsArray[2];
for( int j = 0; j < 2; ++j ) {
vertsArray[j] = verts[vertexConnectivity[i][j]];
}
edges[i] = MemoryManager<SegGeom>
::AllocateSharedPtr( i, three, vertsArray);
}
// ////////////////////////////////////////////////////////////////////
// Set up Prism faces
const int nFaces = 5;
const int edgeConnectivity[][4] = {
{0,1,2,3},
{0,5,4,-1}, // Triangular face
{1,6,5,-1}, // Triangular face
{2,6,7,-1}, // Triangular face
{3,7,4,-1} // Triangular face
};
const bool isEdgeFlipped[][4] = {
{0,0,1,1},
{0,0,1,0},
{0,0,1,0},
{0,0,1,0},
{0,0,1,0}
};
// Populate the list of faces
Geometry2DSharedPtr faces[5];
for (int i = 0; i < nFaces; ++i) {
if (i > 0)
{
SegGeomSharedPtr edgeArray[3];
Orientation eorientArray[3];
for (int j = 0; j < 3; ++j) {
edgeArray[j] = edges[edgeConnectivity[i][j]];
eorientArray[j] = isEdgeFlipped[i][j] ? eBackwards : eForwards;
}
faces[i] = MemoryManager<TriGeom>
::AllocateSharedPtr( i, edgeArray, eorientArray);
}
else
{
SegGeomSharedPtr edgeArray[4];
Orientation eorientArray[4];
for (int j = 0; j < 4; ++j) {
edgeArray[j] = edges[edgeConnectivity[i][j]];
eorientArray[j] = isEdgeFlipped[i][j] ? eBackwards : eForwards;
}
faces[i] = MemoryManager<QuadGeom>
::AllocateSharedPtr( i, edgeArray, eorientArray);
}
}
SpatialDomains::PyrGeomSharedPtr geom =
MemoryManager<SpatialDomains::PyrGeom>::AllocateSharedPtr(faces);
geom->SetOwnData();
return geom;
}
SpatialDomains::TetGeomSharedPtr CreateTetGeom(int argc, char *argv[])
{
......
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Project3D Pyramid Modified basis P=6 Q=7</description>
<executable>LocProject3D</executable>
<parameters>6 4 4 6 6 6 6 7 7 6 0 0 0 1 0 0 1 1 0 0 1 0 0.5 0.5 0.866</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-8">1.80317e-12</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-8">1.20881e-11</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>LocProject_Diff3D Reg. Prism Modified Basis, P=6, Q=7</description>
<executable>LocProject_Diff3D</executable>
<parameters>6 4 4 6 6 6 6 7 7 6 0 0 0 1 0 0 1 1 0 0 1 0 0.5 0 1 0.5 0.5 0.866</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-9">5.9107e-12</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-9">3.81348e-11</value>
</metric>
</metrics>
</test>
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