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

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

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

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

Generated on Wed Jan 31 2024 08:26:51 for Nektar++ by doxygen 1.9.4