Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CheckXmlFile.cpp
Go to the documentation of this file.
1 #include <cstdio>
2 #include <cstdlib>
3 
4 #include <MultiRegions/ExpList.h>
5 
6 using namespace Nektar;
7 
8 bool CheckTetRotation(Array<OneD, NekDouble> &xc, Array<OneD, NekDouble> &yc,
9  Array<OneD, NekDouble> &xz,
10  std::map<int,SpatialDomains::TetGeomSharedPtr>::iterator &tetIter,
11  int id);
12 
13 int main(int argc, char *argv[])
14 {
15  Array<OneD,NekDouble> fce;
16  Array<OneD,NekDouble> xc0,xc1,xc2;
17 
18  if(argc != 2)
19  {
20  fprintf(stderr,"Usage: CheckXmlFile meshfile.xml\n");
21  exit(1);
22  }
23 
24  LibUtilities::SessionReaderSharedPtr vSession
25  = LibUtilities::SessionReader::CreateInstance(argc, argv);
26 
27  //----------------------------------------------
28  // Read in mesh from input file
29  string meshfile(argv[argc-1]);
30  SpatialDomains::MeshGraphSharedPtr mesh = SpatialDomains::MeshGraph::Read(vSession);
31  //----------------------------------------------
32 
33  //----------------------------------------------
34  // Define Expansion
35  int expdim = mesh->GetMeshDimension();
36 
37  switch(expdim)
38  {
39  case 1:
40  ASSERTL0(false,"1D not set up");
41  break;
42  case 2:
43  {
44  NekDouble x,y,z;
45  string outname(strtok(argv[argc-1],"."));
46  outname += ".dat";
47  FILE *fp = fopen(outname.c_str(),"w");
48 
49  SpatialDomains::TriGeomMap trigeom = mesh->GetAllTriGeoms();
50  SpatialDomains::QuadGeomMap quadgeom = mesh->GetAllQuadGeoms();
51 
52  int nverts = mesh->GetNvertices();
53 
54  Array<OneD, NekDouble> xc(nverts),yc(nverts),zc(nverts);
55 
56  for(int i = 0; i < nverts; ++i)
57  {
58  mesh->GetVertex(i)->GetCoords(x,y,z);
59  xc[i] = x;
60  yc[i] = y;
61  zc[i] = z;
62  }
63 
64  std::map<int,SpatialDomains::TriGeomSharedPtr>::iterator triIter;
65  for(triIter = trigeom.begin(); triIter != trigeom.end(); ++triIter)
66  {
67  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);
68  }
69 
70  std::map<int,SpatialDomains::QuadGeomSharedPtr>::iterator quadIter;
71  for(quadIter = quadgeom.begin(); quadIter != quadgeom.end(); ++quadIter)
72  {
73  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);
74  }
75  }
76  break;
77  case 3:
78  {
79  NekDouble x,y,z;
80  string outname(strtok(argv[argc-1],"."));
81  outname += ".dat";
82 
83  SpatialDomains::TetGeomMap tetgeom = mesh->GetAllTetGeoms();
84  SpatialDomains::PyrGeomMap pyrgeom = mesh->GetAllPyrGeoms();
85  SpatialDomains::PrismGeomMap prismgeom = mesh->GetAllPrismGeoms();
86  SpatialDomains::HexGeomMap hexgeom = mesh->GetAllHexGeoms();
87 
88  int nverts = mesh->GetNvertices();
89 
90  Array<OneD, NekDouble> xc(nverts),yc(nverts),zc(nverts);
91 
92  for(int i = 0; i < nverts; ++i)
93  {
94  // Check element roation
95  mesh->GetVertex(i)->GetCoords(x,y,z);
96  xc[i] = x;
97  yc[i] = y;
98  zc[i] = z;
99 
100  }
101 
102  for (int i = 0; i < nverts; ++i)
103  {
104  for (int j = i+1; j < nverts; ++j)
105  {
106  if ((xc[i]-xc[j])*(xc[i]-xc[j]) +
107  (yc[i]-yc[j])*(yc[i]-yc[j]) +
108  (zc[i]-zc[j])*(zc[i]-zc[j]) < 1e-10)
109  {
110  cout << "Duplicate vertices: " << i << " " << j << endl;
111  }
112  }
113  }
114 
115  std::map<int,SpatialDomains::TetGeomSharedPtr>::iterator tetIter;
116  int cnt = 0;
117  bool NoRotateIssues = true;
118  bool NoOrientationIssues = true;
119  for(tetIter = tetgeom.begin(); tetIter != tetgeom.end(); ++tetIter)
120  {
121  // check rotation and dump
122  NoOrientationIssues = CheckTetRotation(xc,yc,zc,tetIter,cnt++);
123 
124  // Check face rotation
125  if((tetIter->second)->GetFace(0)->GetVid(2) != (tetIter->second)->GetVid(2))
126  {
127  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;
128  NoRotateIssues = false;
129  }
130 
131  for(int i = 1; i < 4; ++i)
132 
133  {
134  if((tetIter->second)->GetFace(i)->GetVid(2) != (tetIter->second)->GetVid(3))
135  {
136  cout << "ERROR: Face " << tetIter->second->GetFid(i) << " is not aligned with top Vertex of Tet " << (tetIter->second)->GetGlobalID() << endl;
137  NoRotateIssues = false;
138  }
139  }
140 
141  }
142  if(NoOrientationIssues)
143  {
144  cout << "All Tet have correct ordering for anticlockwise rotation" << endl;
145  }
146 
147  if(NoRotateIssues)
148  {
149  cout << "All Tet faces are correctly aligned" << endl;
150  }
151 
152 
153  std::map<int,SpatialDomains::PyrGeomSharedPtr>::iterator pyrIter;
154  for(pyrIter = pyrgeom.begin(); pyrIter != pyrgeom.end(); ++pyrIter)
155  {
156  // Put pyramid checks in here
157  }
158 
159 
160  // Put prism checks in here
161  std::map<int,SpatialDomains::PrismGeomSharedPtr>::iterator Iter;
162  NoRotateIssues = true;
163  NoOrientationIssues = true;
164  for(Iter = prismgeom.begin(); Iter != prismgeom.end(); ++Iter)
165  {
166 
167  // Check face rotation
168  if((Iter->second)->GetFace(1)->GetVid(2) != (Iter->second)->GetVid(4))
169  {
170  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;
171  NoRotateIssues = false;
172  }
173 
174 
175  // Check face rotation
176  if((Iter->second)->GetFace(3)->GetVid(2) != (Iter->second)->GetVid(5))
177  {
178  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;
179  NoRotateIssues = false;
180  }
181 
182  }
183 
184  if(NoRotateIssues)
185  {
186  cout << "All Prism Tri faces are correctly aligned" << endl;
187  }
188 
189  std::map<int,SpatialDomains::HexGeomSharedPtr>::iterator hexIter;
190  for(hexIter = hexgeom.begin(); hexIter != hexgeom.end(); ++hexIter)
191  {
192  // PUt Hex checks in here
193  }
194 
195  }
196 
197  break;
198  default:
199  ASSERTL0(false,"Expansion dimension not recognised");
200  break;
201  }
202 
203  //-----------------------------------------------
204 
205  return 0;
206 }
207 
208 class Ord
209 {
210 public:
211  double x;
212  double y;
213  double z;
214 };
215 
216 bool CheckTetRotation(Array<OneD, NekDouble> &xc, Array<OneD, NekDouble> &yc, Array<OneD, NekDouble> &zc, std::map<int,SpatialDomains::TetGeomSharedPtr>::iterator &tetIter, int id)
217 {
218  bool RotationOK = true;
219  Ord v[4];
220  NekDouble abx,aby,abz;
221 
222  v[0].x = xc[(tetIter->second)->GetVid(0)];
223  v[0].y = yc[(tetIter->second)->GetVid(0)];
224  v[0].z = zc[(tetIter->second)->GetVid(0)];
225 
226  v[1].x = xc[(tetIter->second)->GetVid(1)];
227  v[1].y = yc[(tetIter->second)->GetVid(1)];
228  v[1].z = zc[(tetIter->second)->GetVid(1)];
229 
230  v[2].x = xc[(tetIter->second)->GetVid(2)];
231  v[2].y = yc[(tetIter->second)->GetVid(2)];
232  v[2].z = zc[(tetIter->second)->GetVid(2)];
233 
234  v[3].x = xc[(tetIter->second)->GetVid(3)];
235  v[3].y = yc[(tetIter->second)->GetVid(3)];
236  v[3].z = zc[(tetIter->second)->GetVid(3)];
237 
238  // cross product of edge 0 and 2
239  abx = (v[1].y-v[0].y)*(v[2].z-v[0].z) -
240  (v[1].z-v[0].z)*(v[2].y-v[0].y);
241  aby = (v[1].z-v[0].z)*(v[2].x-v[0].x) -
242  (v[1].x-v[0].x)*(v[2].z-v[0].z);
243  abz = (v[1].x-v[0].x)*(v[2].y-v[0].y) -
244  (v[1].y-v[0].y)*(v[2].x-v[0].x);
245 
246  // inner product of cross product with respect to edge 3 should be positive
247  if(((v[3].x-v[0].x)*abx + (v[3].y-v[0].y)*aby +
248  (v[3].z-v[0].z)*abz)<0.0)
249  {
250  cerr << "ERROR: Element " << id + 1 << "is NOT counter-clockwise\n" << endl;
251  RotationOK = false;
252  }
253  return RotationOK;
254 }
Generated on Sun Mar 15 2015 00:11:30 for Nektar++ by   doxygen 1.8.1.2