Nektar++
Classes | Functions
CheckXmlFile.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <MultiRegions/ExpList.h>
#include <SpatialDomains/MeshGraphIO.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 276 of file CheckXmlFile.cpp.

280{
281 bool RotationOK = true;
282 Ord v[4];
283 NekDouble abx, aby, abz;
284
285 v[0].x = xc[(tetIter->second)->GetVid(0)];
286 v[0].y = yc[(tetIter->second)->GetVid(0)];
287 v[0].z = zc[(tetIter->second)->GetVid(0)];
288
289 v[1].x = xc[(tetIter->second)->GetVid(1)];
290 v[1].y = yc[(tetIter->second)->GetVid(1)];
291 v[1].z = zc[(tetIter->second)->GetVid(1)];
292
293 v[2].x = xc[(tetIter->second)->GetVid(2)];
294 v[2].y = yc[(tetIter->second)->GetVid(2)];
295 v[2].z = zc[(tetIter->second)->GetVid(2)];
296
297 v[3].x = xc[(tetIter->second)->GetVid(3)];
298 v[3].y = yc[(tetIter->second)->GetVid(3)];
299 v[3].z = zc[(tetIter->second)->GetVid(3)];
300
301 // cross product of edge 0 and 2
302 abx = (v[1].y - v[0].y) * (v[2].z - v[0].z) -
303 (v[1].z - v[0].z) * (v[2].y - v[0].y);
304 aby = (v[1].z - v[0].z) * (v[2].x - v[0].x) -
305 (v[1].x - v[0].x) * (v[2].z - v[0].z);
306 abz = (v[1].x - v[0].x) * (v[2].y - v[0].y) -
307 (v[1].y - v[0].y) * (v[2].x - v[0].x);
308
309 // inner product of cross product with respect to edge 3 should be positive
310 if (((v[3].x - v[0].x) * abx + (v[3].y - v[0].y) * aby +
311 (v[3].z - v[0].z) * abz) < 0.0)
312 {
313 cerr << "ERROR: Element " << id + 1 << "is NOT counter-clockwise\n"
314 << endl;
315 RotationOK = false;
316 }
317 return RotationOK;
318}
Ord::x
double x
Definition: CheckXmlFile.cpp:271
Ord::y
double y
Definition: CheckXmlFile.cpp:272
Ord::z
double z
Definition: CheckXmlFile.cpp:273
Nektar::UnitTests::z
std::vector< double > z(NPUPPER)
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43

References Ord::x, Ord::y, Nektar::UnitTests::z(), and Ord::z.

Referenced by main().

◆ main()

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

Definition at line 49 of file CheckXmlFile.cpp.

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

References ASSERTL0, CheckTetRotation(), and Nektar::UnitTests::z().

Generated on Wed Oct 30 2024 16:07:35 for Nektar++ by doxygen 1.9.4