Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Classes | Functions
CheckXmlFile.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <MultiRegions/ExpList.h>
Include dependency graph for CheckXmlFile.cpp:

Go to the source code of this file.

Classes

class  Ord

Functions

bool CheckTetRotation (Array< OneD, NekDouble > &xc, Array< OneD, NekDouble > &yc, Array< OneD, NekDouble > &xz, std::map< int, SpatialDomains::TetGeomSharedPtr >::iterator &tetIter, int id)
int main (int argc, char *argv[])

Function Documentation

bool CheckTetRotation ( Array< OneD, NekDouble > &  xc,
Array< OneD, NekDouble > &  yc,
Array< OneD, NekDouble > &  xz,
std::map< int, SpatialDomains::TetGeomSharedPtr >::iterator &  tetIter,
int  id 
)

Definition at line 216 of file CheckXmlFile.cpp.

References Ord::x, Ord::y, and Ord::z.

Referenced by main().

{
bool RotationOK = true;
Ord v[4];
NekDouble abx,aby,abz;
v[0].x = xc[(tetIter->second)->GetVid(0)];
v[0].y = yc[(tetIter->second)->GetVid(0)];
v[0].z = zc[(tetIter->second)->GetVid(0)];
v[1].x = xc[(tetIter->second)->GetVid(1)];
v[1].y = yc[(tetIter->second)->GetVid(1)];
v[1].z = zc[(tetIter->second)->GetVid(1)];
v[2].x = xc[(tetIter->second)->GetVid(2)];
v[2].y = yc[(tetIter->second)->GetVid(2)];
v[2].z = zc[(tetIter->second)->GetVid(2)];
v[3].x = xc[(tetIter->second)->GetVid(3)];
v[3].y = yc[(tetIter->second)->GetVid(3)];
v[3].z = zc[(tetIter->second)->GetVid(3)];
// cross product of edge 0 and 2
abx = (v[1].y-v[0].y)*(v[2].z-v[0].z) -
(v[1].z-v[0].z)*(v[2].y-v[0].y);
aby = (v[1].z-v[0].z)*(v[2].x-v[0].x) -
(v[1].x-v[0].x)*(v[2].z-v[0].z);
abz = (v[1].x-v[0].x)*(v[2].y-v[0].y) -
(v[1].y-v[0].y)*(v[2].x-v[0].x);
// inner product of cross product with respect to edge 3 should be positive
if(((v[3].x-v[0].x)*abx + (v[3].y-v[0].y)*aby +
(v[3].z-v[0].z)*abz)<0.0)
{
cerr << "ERROR: Element " << id + 1 << "is NOT counter-clockwise\n" << endl;
RotationOK = false;
}
return RotationOK;
}
int main ( int  argc,
char *  argv[] 
)

Definition at line 13 of file CheckXmlFile.cpp.

References ASSERTL0, CheckTetRotation(), Nektar::LibUtilities::SessionReader::CreateInstance(), Nektar::iterator, and Nektar::SpatialDomains::MeshGraph::Read().

{
Array<OneD,NekDouble> fce;
Array<OneD,NekDouble> xc0,xc1,xc2;
if(argc != 2)
{
fprintf(stderr,"Usage: CheckXmlFile meshfile.xml\n");
exit(1);
}
= LibUtilities::SessionReader::CreateInstance(argc, argv);
//----------------------------------------------
// Read in mesh from input file
string meshfile(argv[argc-1]);
SpatialDomains::MeshGraphSharedPtr mesh = SpatialDomains::MeshGraph::Read(vSession);
//----------------------------------------------
//----------------------------------------------
// Define Expansion
int expdim = mesh->GetMeshDimension();
switch(expdim)
{
case 1:
ASSERTL0(false,"1D not set up");
break;
case 2:
{
NekDouble x,y,z;
string outname(strtok(argv[argc-1],"."));
outname += ".dat";
FILE *fp = fopen(outname.c_str(),"w");
SpatialDomains::TriGeomMap trigeom = mesh->GetAllTriGeoms();
SpatialDomains::QuadGeomMap quadgeom = mesh->GetAllQuadGeoms();
int nverts = mesh->GetNvertices();
Array<OneD, NekDouble> xc(nverts),yc(nverts),zc(nverts);
for(int i = 0; i < nverts; ++i)
{
mesh->GetVertex(i)->GetCoords(x,y,z);
xc[i] = x;
yc[i] = y;
zc[i] = z;
}
for(triIter = trigeom.begin(); triIter != trigeom.end(); ++triIter)
{
fprintf(fp,"%d %d %d %d\n",(triIter->second)->GetVid(0)+1,(triIter->second)->GetVid(1)+1,(triIter->second)->GetVid(2)+1,(triIter->second)->GetVid(2)+1);
}
for(quadIter = quadgeom.begin(); quadIter != quadgeom.end(); ++quadIter)
{
fprintf(fp,"%d %d %d %d\n",(quadIter->second)->GetVid(0)+1,(quadIter->second)->GetVid(1)+1,(quadIter->second)->GetVid(2)+1,(quadIter->second)->GetVid(3)+1);
}
}
break;
case 3:
{
NekDouble x,y,z;
string outname(strtok(argv[argc-1],"."));
outname += ".dat";
SpatialDomains::TetGeomMap tetgeom = mesh->GetAllTetGeoms();
SpatialDomains::PyrGeomMap pyrgeom = mesh->GetAllPyrGeoms();
SpatialDomains::PrismGeomMap prismgeom = mesh->GetAllPrismGeoms();
SpatialDomains::HexGeomMap hexgeom = mesh->GetAllHexGeoms();
int nverts = mesh->GetNvertices();
Array<OneD, NekDouble> xc(nverts),yc(nverts),zc(nverts);
for(int i = 0; i < nverts; ++i)
{
// Check element roation
mesh->GetVertex(i)->GetCoords(x,y,z);
xc[i] = x;
yc[i] = y;
zc[i] = z;
}
for (int i = 0; i < nverts; ++i)
{
for (int j = i+1; j < nverts; ++j)
{
if ((xc[i]-xc[j])*(xc[i]-xc[j]) +
(yc[i]-yc[j])*(yc[i]-yc[j]) +
(zc[i]-zc[j])*(zc[i]-zc[j]) < 1e-10)
{
cout << "Duplicate vertices: " << i << " " << j << endl;
}
}
}
int cnt = 0;
bool NoRotateIssues = true;
bool NoOrientationIssues = true;
for(tetIter = tetgeom.begin(); tetIter != tetgeom.end(); ++tetIter)
{
// check rotation and dump
NoOrientationIssues = CheckTetRotation(xc,yc,zc,tetIter,cnt++);
// Check face rotation
if((tetIter->second)->GetFace(0)->GetVid(2) != (tetIter->second)->GetVid(2))
{
cout << "ERROR: Face " << tetIter->second->GetFid(0) << " (vert "<< (tetIter->second)->GetFace(0)->GetVid(2) << ") is not aligned with base vertex of Tet " << (tetIter->second)->GetGlobalID() << " (vert "<< (tetIter->second)->GetVid(2) <<")" << endl;
NoRotateIssues = false;
}
for(int i = 1; i < 4; ++i)
{
if((tetIter->second)->GetFace(i)->GetVid(2) != (tetIter->second)->GetVid(3))
{
cout << "ERROR: Face " << tetIter->second->GetFid(i) << " is not aligned with top Vertex of Tet " << (tetIter->second)->GetGlobalID() << endl;
NoRotateIssues = false;
}
}
}
if(NoOrientationIssues)
{
cout << "All Tet have correct ordering for anticlockwise rotation" << endl;
}
if(NoRotateIssues)
{
cout << "All Tet faces are correctly aligned" << endl;
}
for(pyrIter = pyrgeom.begin(); pyrIter != pyrgeom.end(); ++pyrIter)
{
// Put pyramid checks in here
}
// Put prism checks in here
NoRotateIssues = true;
NoOrientationIssues = true;
for(Iter = prismgeom.begin(); Iter != prismgeom.end(); ++Iter)
{
// Check face rotation
if((Iter->second)->GetFace(1)->GetVid(2) != (Iter->second)->GetVid(4))
{
cout << "ERROR: Face " << Iter->second->GetFid(1) << " (vert "<< (Iter->second)->GetFace(1)->GetVid(2) << ") not aligned to face 1 singular vert of Prism " << (Iter->second)->GetGlobalID() << " (vert "<< (Iter->second)->GetVid(4) <<")" << endl;
NoRotateIssues = false;
}
// Check face rotation
if((Iter->second)->GetFace(3)->GetVid(2) != (Iter->second)->GetVid(5))
{
cout << "ERROR: Face " << Iter->second->GetFid(3) << " (vert "<< (Iter->second)->GetFace(3)->GetVid(2) << ") not aligned to face 3 singular vert of Prism " << (Iter->second)->GetGlobalID() << " (vert "<< (Iter->second)->GetVid(5) <<")" << endl;
NoRotateIssues = false;
}
}
if(NoRotateIssues)
{
cout << "All Prism Tri faces are correctly aligned" << endl;
}
for(hexIter = hexgeom.begin(); hexIter != hexgeom.end(); ++hexIter)
{
// PUt Hex checks in here
}
}
break;
default:
ASSERTL0(false,"Expansion dimension not recognised");
break;
}
//-----------------------------------------------
return 0;
}