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 217 of file CheckXmlFile.cpp.

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

Referenced by main().

218 {
219  bool RotationOK = true;
220  Ord v[4];
221  NekDouble abx,aby,abz;
222 
223  v[0].x = xc[(tetIter->second)->GetVid(0)];
224  v[0].y = yc[(tetIter->second)->GetVid(0)];
225  v[0].z = zc[(tetIter->second)->GetVid(0)];
226 
227  v[1].x = xc[(tetIter->second)->GetVid(1)];
228  v[1].y = yc[(tetIter->second)->GetVid(1)];
229  v[1].z = zc[(tetIter->second)->GetVid(1)];
230 
231  v[2].x = xc[(tetIter->second)->GetVid(2)];
232  v[2].y = yc[(tetIter->second)->GetVid(2)];
233  v[2].z = zc[(tetIter->second)->GetVid(2)];
234 
235  v[3].x = xc[(tetIter->second)->GetVid(3)];
236  v[3].y = yc[(tetIter->second)->GetVid(3)];
237  v[3].z = zc[(tetIter->second)->GetVid(3)];
238 
239  // cross product of edge 0 and 2
240  abx = (v[1].y-v[0].y)*(v[2].z-v[0].z) -
241  (v[1].z-v[0].z)*(v[2].y-v[0].y);
242  aby = (v[1].z-v[0].z)*(v[2].x-v[0].x) -
243  (v[1].x-v[0].x)*(v[2].z-v[0].z);
244  abz = (v[1].x-v[0].x)*(v[2].y-v[0].y) -
245  (v[1].y-v[0].y)*(v[2].x-v[0].x);
246 
247  // inner product of cross product with respect to edge 3 should be positive
248  if(((v[3].x-v[0].x)*abx + (v[3].y-v[0].y)*aby +
249  (v[3].z-v[0].z)*abz)<0.0)
250  {
251  cerr << "ERROR: Element " << id + 1 << "is NOT counter-clockwise\n" << endl;
252  RotationOK = false;
253  }
254  return RotationOK;
255 }
double x
double y
double z
double NekDouble
int main ( int  argc,
char *  argv[] 
)

Definition at line 14 of file CheckXmlFile.cpp.

References ASSERTL0, CheckTetRotation(), and Nektar::iterator.

15 {
17  Array<OneD,NekDouble> xc0,xc1,xc2;
18 
19  if(argc != 2)
20  {
21  fprintf(stderr,"Usage: CheckXmlFile meshfile.xml\n");
22  exit(1);
23  }
24 
26  = LibUtilities::SessionReader::CreateInstance(argc, argv);
27 
28  //----------------------------------------------
29  // Read in mesh from input file
30  string meshfile(argv[argc-1]);
31  SpatialDomains::MeshGraphSharedPtr mesh = SpatialDomains::MeshGraph::Read(vSession);
32  //----------------------------------------------
33 
34  //----------------------------------------------
35  // Define Expansion
36  int expdim = mesh->GetMeshDimension();
37 
38  switch(expdim)
39  {
40  case 1:
41  ASSERTL0(false,"1D not set up");
42  break;
43  case 2:
44  {
45  NekDouble x,y,z;
46  string outname(strtok(argv[argc-1],"."));
47  outname += ".dat";
48  FILE *fp = fopen(outname.c_str(),"w");
49 
50  SpatialDomains::TriGeomMap trigeom = mesh->GetAllTriGeoms();
51  SpatialDomains::QuadGeomMap quadgeom = mesh->GetAllQuadGeoms();
52 
53  int nverts = mesh->GetNvertices();
54 
55  Array<OneD, NekDouble> xc(nverts),yc(nverts),zc(nverts);
56 
57  for(int i = 0; i < nverts; ++i)
58  {
59  mesh->GetVertex(i)->GetCoords(x,y,z);
60  xc[i] = x;
61  yc[i] = y;
62  zc[i] = z;
63  }
64 
66  for(triIter = trigeom.begin(); triIter != trigeom.end(); ++triIter)
67  {
68  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);
69  }
70 
72  for(quadIter = quadgeom.begin(); quadIter != quadgeom.end(); ++quadIter)
73  {
74  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);
75  }
76  }
77  break;
78  case 3:
79  {
80  NekDouble x,y,z;
81  string outname(strtok(argv[argc-1],"."));
82  outname += ".dat";
83 
84  SpatialDomains::TetGeomMap tetgeom = mesh->GetAllTetGeoms();
85  SpatialDomains::PyrGeomMap pyrgeom = mesh->GetAllPyrGeoms();
86  SpatialDomains::PrismGeomMap prismgeom = mesh->GetAllPrismGeoms();
87  SpatialDomains::HexGeomMap hexgeom = mesh->GetAllHexGeoms();
88 
89  int nverts = mesh->GetNvertices();
90 
91  Array<OneD, NekDouble> xc(nverts),yc(nverts),zc(nverts);
92 
93  for(int i = 0; i < nverts; ++i)
94  {
95  // Check element roation
96  mesh->GetVertex(i)->GetCoords(x,y,z);
97  xc[i] = x;
98  yc[i] = y;
99  zc[i] = z;
100 
101  }
102 
103  for (int i = 0; i < nverts; ++i)
104  {
105  for (int j = i+1; j < nverts; ++j)
106  {
107  if ((xc[i]-xc[j])*(xc[i]-xc[j]) +
108  (yc[i]-yc[j])*(yc[i]-yc[j]) +
109  (zc[i]-zc[j])*(zc[i]-zc[j]) < 1e-10)
110  {
111  cout << "Duplicate vertices: " << i << " " << j << endl;
112  }
113  }
114  }
115 
117  int cnt = 0;
118  bool NoRotateIssues = true;
119  bool NoOrientationIssues = true;
120  for(tetIter = tetgeom.begin(); tetIter != tetgeom.end(); ++tetIter)
121  {
122  // check rotation and dump
123  NoOrientationIssues = CheckTetRotation(xc,yc,zc,tetIter,cnt++);
124 
125  // Check face rotation
126  if((tetIter->second)->GetFace(0)->GetVid(2) != (tetIter->second)->GetVid(2))
127  {
128  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;
129  NoRotateIssues = false;
130  }
131 
132  for(int i = 1; i < 4; ++i)
133 
134  {
135  if((tetIter->second)->GetFace(i)->GetVid(2) != (tetIter->second)->GetVid(3))
136  {
137  cout << "ERROR: Face " << tetIter->second->GetFid(i) << " is not aligned with top Vertex of Tet " << (tetIter->second)->GetGlobalID() << endl;
138  NoRotateIssues = false;
139  }
140  }
141 
142  }
143  if(NoOrientationIssues)
144  {
145  cout << "All Tet have correct ordering for anticlockwise rotation" << endl;
146  }
147 
148  if(NoRotateIssues)
149  {
150  cout << "All Tet faces are correctly aligned" << endl;
151  }
152 
153 
155  for(pyrIter = pyrgeom.begin(); pyrIter != pyrgeom.end(); ++pyrIter)
156  {
157  // Put pyramid checks in here
158  }
159 
160 
161  // Put prism checks in here
163  NoRotateIssues = true;
164  NoOrientationIssues = true;
165  for(Iter = prismgeom.begin(); Iter != prismgeom.end(); ++Iter)
166  {
167 
168  // Check face rotation
169  if((Iter->second)->GetFace(1)->GetVid(2) != (Iter->second)->GetVid(4))
170  {
171  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;
172  NoRotateIssues = false;
173  }
174 
175 
176  // Check face rotation
177  if((Iter->second)->GetFace(3)->GetVid(2) != (Iter->second)->GetVid(5))
178  {
179  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;
180  NoRotateIssues = false;
181  }
182 
183  }
184 
185  if(NoRotateIssues)
186  {
187  cout << "All Prism Tri faces are correctly aligned" << endl;
188  }
189 
191  for(hexIter = hexgeom.begin(); hexIter != hexgeom.end(); ++hexIter)
192  {
193  // PUt Hex checks in here
194  }
195 
196  }
197 
198  break;
199  default:
200  ASSERTL0(false,"Expansion dimension not recognised");
201  break;
202  }
203 
204  //-----------------------------------------------
205 
206  return 0;
207 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
std::map< int, PrismGeomSharedPtr > PrismGeomMap
Definition: PrismGeom.h:112
std::map< int, PyrGeomSharedPtr > PyrGeomMap
Definition: PyrGeom.h:87
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
bool CheckTetRotation(Array< OneD, NekDouble > &xc, Array< OneD, NekDouble > &yc, Array< OneD, NekDouble > &xz, std::map< int, SpatialDomains::TetGeomSharedPtr >::iterator &tetIter, int id)
std::map< int, TetGeomSharedPtr > TetGeomMap
Definition: TetGeom.h:109
std::map< int, TriGeomSharedPtr > TriGeomMap
Definition: TriGeom.h:62
double NekDouble
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
std::map< int, QuadGeomSharedPtr > QuadGeomMap
Definition: QuadGeom.h:57
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
std::map< int, HexGeomSharedPtr > HexGeomMap
Definition: HexGeom.h:113