Nektar++
Classes | Functions
CheckXmlFile.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <MultiRegions/ExpList.h>

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

◆ CheckTetRotation()

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

245 {
246  bool RotationOK = true;
247  Ord v[4];
248  NekDouble abx, aby, abz;
249 
250  v[0].x = xc[(tetIter->second)->GetVid(0)];
251  v[0].y = yc[(tetIter->second)->GetVid(0)];
252  v[0].z = zc[(tetIter->second)->GetVid(0)];
253 
254  v[1].x = xc[(tetIter->second)->GetVid(1)];
255  v[1].y = yc[(tetIter->second)->GetVid(1)];
256  v[1].z = zc[(tetIter->second)->GetVid(1)];
257 
258  v[2].x = xc[(tetIter->second)->GetVid(2)];
259  v[2].y = yc[(tetIter->second)->GetVid(2)];
260  v[2].z = zc[(tetIter->second)->GetVid(2)];
261 
262  v[3].x = xc[(tetIter->second)->GetVid(3)];
263  v[3].y = yc[(tetIter->second)->GetVid(3)];
264  v[3].z = zc[(tetIter->second)->GetVid(3)];
265 
266  // cross product of edge 0 and 2
267  abx = (v[1].y - v[0].y) * (v[2].z - v[0].z) -
268  (v[1].z - v[0].z) * (v[2].y - v[0].y);
269  aby = (v[1].z - v[0].z) * (v[2].x - v[0].x) -
270  (v[1].x - v[0].x) * (v[2].z - v[0].z);
271  abz = (v[1].x - v[0].x) * (v[2].y - v[0].y) -
272  (v[1].y - v[0].y) * (v[2].x - v[0].x);
273 
274  // inner product of cross product with respect to edge 3 should be positive
275  if (((v[3].x - v[0].x) * abx + (v[3].y - v[0].y) * aby +
276  (v[3].z - v[0].z) * abz) < 0.0)
277  {
278  cerr << "ERROR: Element " << id + 1 << "is NOT counter-clockwise\n"
279  << endl;
280  RotationOK = false;
281  }
282  return RotationOK;
283 }
Ord::x
double x
Definition: CheckXmlFile.cpp:236
Ord::y
double y
Definition: CheckXmlFile.cpp:237
Ord::z
double z
Definition: CheckXmlFile.cpp:238
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43

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

Referenced by main().

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 14 of file CheckXmlFile.cpp.

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

References ASSERTL0, and CheckTetRotation().

Generated on Tue Aug 23 2022 07:27:16 for Nektar++ by doxygen 1.9.1