TexGen
DomainPrism.cpp
Go to the documentation of this file.
1/*=============================================================================
2TexGen: Geometric textile modeller.
3Copyright (C) 2006 Martin Sherburn
4
5This program is free software; you can redistribute it and/or
6modify it under the terms of the GNU General Public License
7as published by the Free Software Foundation; either version 2
8of the License, or (at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13GNU General Public License for more details.
14
15You should have received a copy of the GNU General Public License
16along with this program; if not, write to the Free Software
17Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18=============================================================================*/
19
20#include "PrecompiledHeaders.h"
21#include "DomainPrism.h"
22#include "Yarn.h"
23#include "SectionPolygon.h"
24using namespace TexGen;
25
26CDomainPrism::CDomainPrism(const vector<XY> &Points, XYZ &start, XYZ &end)
27 :m_Points(Points)
28{
29 // Prism is described by a yarn with two nodes (to give the orientation )
30 // and a cross-section defined by a polygon section created using the points input
31 CNode Node = CNode(start);
32 m_Yarn.AddNode(Node);
33 Node = CNode(end);
34 m_Yarn.AddNode(Node);
35 // Last parameter in CSectionPolygon forces input points only to be used to generate the polygon
36 m_Yarn.AssignSection(CYarnSectionConstant(CSectionPolygon(Points, false, true)));
37 m_Yarn.SetResolution(2, (int)Points.size()); // Need better definition of resolution
38 BuildMesh();
39 GeneratePlanes();
40}
41
42CDomainPrism::~CDomainPrism(void)
43{
44}
45
46CDomainPrism::CDomainPrism(TiXmlElement &Element)
47 : CDomain(Element)
48{
49 TiXmlElement* pYarn = Element.FirstChildElement("Yarn");
50 m_Yarn = CYarn(*pYarn);
51 // Build the domain mesh if this is empty
52 if (m_Mesh.GetNumNodes() == 0)
53 {
54 BuildMesh();
55 }
56 GeneratePlanes();
57}
58
59void CDomainPrism::PopulateTiXmlElement(TiXmlElement &Element, OUTPUT_TYPE OutputType) const
60{
61 CDomain::PopulateTiXmlElement(Element, OutputType);
62 TiXmlElement Yarn("Yarn");
63 m_Yarn.PopulateTiXmlElement(Yarn, OutputType);
64 Element.InsertEndChild(Yarn);
65}
66
67void CDomainPrism::BuildMesh()
68{
69 // Build a surface mesh for the yarn representing the domain
70 m_Mesh.Clear();
71 m_Yarn.AddSurfaceToMesh( m_Mesh );
72 m_Mesh.MergeNodes();
73 // m_Mesh.SaveToVTK("DomainMesh"); // Useful for debugging to check mesh created
74}
75
76void CDomainPrism::GeneratePlanes()
77{
78 if (m_Mesh.GetNumNodes() == 0)
79 {
80 BuildMesh();
81 }
82
83 // Generate set of planes to form the domain, one for each element of the surface mesh of the domain 'yarn'
84
85 list<int>::iterator itInt;
86 list<int>::iterator itStart;
87 XYZ points[4];
88 m_ElementPlanes.clear();
89
90 list<int> &QuadIndices = m_Mesh.GetIndices(CMesh::QUAD);
91 for (itInt = QuadIndices.begin(); itInt != QuadIndices.end(); )
92 {
93 itStart = itInt;
94
95 points[0] = m_Mesh.GetNode(*(itInt++));
96 points[1] = m_Mesh.GetNode(*(itInt++));
97 points[2] = m_Mesh.GetNode(*(itInt++));
98 points[3] = m_Mesh.GetNode(*(itInt++));
99
100 PLANE ElementPlane;
101 if (GetPlane( points, ElementPlane))
102 m_ElementPlanes.push_back(ElementPlane);
103 }
104
105 list<int> &TriIndices = m_Mesh.GetIndices(CMesh::TRI);
106 for (itInt = TriIndices.begin(); itInt != TriIndices.end(); )
107 {
108 itStart = itInt;
109
110 points[0] = m_Mesh.GetNode(*(itInt++));
111 points[1] = m_Mesh.GetNode(*(itInt++));
112 points[2] = m_Mesh.GetNode(*(itInt++));
113
114 PLANE ElementPlane;
115 if (GetPlane(points, ElementPlane))
116 m_ElementPlanes.push_back(ElementPlane);
117 }
118
119 list<int> &Indices = m_Mesh.GetIndices(CMesh::POLYGON);
120 list<int>::iterator itIndices;
121 int StartIndex;
122 list<int>::iterator itStartIndex;
123
124 for (itIndices = Indices.begin(); itIndices != Indices.end(); )
125 {
126 StartIndex = (*(itIndices));
127 itStartIndex = itIndices;
128 vector<int> PolygonIndex;
129 XYZ p;
130
131 int i = 0, iPoints = 0;
132
133 do {
134 if (i < 3)
135 {
136 points[i] = m_Mesh.GetNode(*(itIndices));
137 }
138 ++i;
139 ++itIndices;
140 if (i == 3)
141 {
142 PLANE ElementPlane;
143 if (GetPlane(points, ElementPlane))
144 m_ElementPlanes.push_back(ElementPlane);
145 }
146 } while ((*itIndices) != StartIndex);
147 ++itIndices;
148 }
149 RemoveDuplicatePlanes();
150}
151
152bool CDomainPrism::GetPlane(XYZ *p, PLANE &plane)
153{
154 plane.Normal = CrossProduct(p[2] - p[0], p[1] - p[0]);
155
156 if (GetLength(plane.Normal))
157 {
158 Normalise(plane.Normal);
159 plane.d = DotProduct(plane.Normal, p[0] - XYZ(0, 0, 0));
160 return(true);
161 }
162 return(false);
163}
164
165void CDomainPrism::RemoveDuplicatePlanes()
166{
167 vector<PLANE>::iterator itPlanes1, itPlanes2;
168
169 for (itPlanes1 = m_ElementPlanes.begin(); itPlanes1 != m_ElementPlanes.end(); ++itPlanes1)
170 {
171 for (itPlanes2 = itPlanes1 + 1; itPlanes2 != m_ElementPlanes.end(); )
172 {
173 if (PlaneEqual(*itPlanes1, *itPlanes2))
174 {
175 itPlanes2 = m_ElementPlanes.erase(itPlanes2);
176 }
177 else
178 ++itPlanes2;
179 }
180 }
181}
182
183bool CDomainPrism::PlaneEqual(PLANE Plane1, PLANE Plane2)
184{
185 if (fabs(GetLength(Plane1.Normal, Plane2.Normal)) > 1e-9)
186 return false;
187 if (Plane1.d != Plane2.d)
188 return false;
189 return true;
190}
191
192void CDomainPrism::ClipIntersectMeshToDomain(CMesh &Mesh, bool bFillGaps) const
193{
194 // Clip mesh elements which are known to intersect with the domain
195 // Checks each mesh element against each domain plane
196 const double TOL = 1e-9;
197
198 TGLOGINDENT("Clipping mesh to domain");
199
200 vector<XYZ> NewTriangles;
201 vector<bool> NewTrianglesFlipped;
202 vector<XYZ>::iterator itXYZ;
203 vector<XYZ> NewQuads;
204
205 vector<PLANE>::const_iterator itPlane;
206 list<int>::iterator itStart;
207 list<int>::iterator itInt;
208 const XYZ *p1, *p2, *p3, *p4;
209 double d1, d2, d3, d4; // d represents the distance of the point to the plane (i.e. +ve inside, -ve outside, 0 on top)
210 double dTemp;
211 const XYZ *pTemp;
212 double u;
213 int i;
214 int iLastNodeIndex;
215 bool bFlipped;
216
217 // Deal with surface elements
218 for (itPlane = m_ElementPlanes.begin(); itPlane != m_ElementPlanes.end(); ++itPlane)
219 {
220 list<int> &QuadIndices = Mesh.GetIndices(CMesh::QUAD);
221 for (itInt = QuadIndices.begin(); itInt != QuadIndices.end(); )
222 {
223 itStart = itInt;
224
225 p1 = &Mesh.GetNode(*(itInt++));
226 p2 = &Mesh.GetNode(*(itInt++));
227 p3 = &Mesh.GetNode(*(itInt++));
228 p4 = &Mesh.GetNode(*(itInt++));
229
230 d1 = DotProduct(itPlane->Normal, *p1) - itPlane->d;
231 d2 = DotProduct(itPlane->Normal, *p2) - itPlane->d;
232 d3 = DotProduct(itPlane->Normal, *p3) - itPlane->d;
233 d4 = DotProduct(itPlane->Normal, *p4) - itPlane->d;
234
235 if (d1 <= TOL && d2 <= TOL && d3 <= TOL && d4 <= TOL) // The quad lies completely on or outside the plane
236 {
237 if ( fabs(d1) < TOL || fabs(d2) < TOL || fabs(d3) < TOL || fabs(d4) < TOL)
238 itInt = QuadIndices.erase(itStart, itInt); // Delete the quad
239 }
240 else if (d1 >= -TOL && d2 >= -TOL && d3 >= -TOL && d4 >= -TOL) // The quad lies completely inside the plane
241 {
242 // Do nothing
243 }
244 else if (d1 < TOL && d2 < TOL && d3 > TOL && d4 > TOL) // Points 1 & 2 outside plane
245 {
246 itInt = QuadIndices.erase(itStart, itInt); // Delete the quad
247 u = d4 / (d4 - d1);
248 NewQuads.push_back(*p4 + (*p1 - *p4) * u);
249 u = d3 / (d3 - d2);
250 NewQuads.push_back(*p3 + (*p2 - *p3) * u);
251 NewQuads.push_back(*p3);
252 NewQuads.push_back(*p4);
253 }
254 else if (d2 < TOL && d3 < TOL && d4 > TOL && d1 > TOL) // Points 2 & 3 outside plane
255 {
256 itInt = QuadIndices.erase(itStart, itInt); // Delete the quad
257 NewQuads.push_back(*p1);
258 u = d1 / (d1 - d2);
259 NewQuads.push_back(*p1 + (*p2 - *p1) * u);
260 u = d4 / (d4 - d3);
261 NewQuads.push_back(*p4 + (*p3 - *p4) * u);
262 NewQuads.push_back(*p4);
263 }
264 else if (d3 < TOL && d4 < TOL && d1 > TOL && d2 > TOL) // Points 3 & 4 outside plane
265 {
266 itInt = QuadIndices.erase(itStart, itInt); // Delete the quad
267 NewQuads.push_back(*p1);
268 NewQuads.push_back(*p2);
269 u = d2 / (d2 - d3);
270 NewQuads.push_back(*p2 + (*p3 - *p2) * u);
271 u = d1 / (d1 - d4);
272 NewQuads.push_back(*p1 + (*p4 - *p1) * u);
273 }
274 else if (d4 < TOL && d1 < TOL && d2 > TOL && d3 > TOL) // Points 4 & 1 outside plane
275 {
276 itInt = QuadIndices.erase(itStart, itInt); // Delete the quad
277 u = d2 / (d2 - d1);
278 NewQuads.push_back(*p2 + (*p1 - *p2) * u);
279 NewQuads.push_back(*p2);
280 NewQuads.push_back(*p3);
281 u = d3 / (d3 - d4);
282 NewQuads.push_back(*p3 + (*p4 - *p3) * u);
283
284 }
285 else // Convert the quad to a triangle for trimming if 1 or 3 points inside plane
286 {
287 itInt = Mesh.ConvertQuadtoTriangles(itStart);
288 }
289 }
290
291 // Add the new quads to the mesh, and clear the new quad list
292 iLastNodeIndex = int(Mesh.GetNumNodes());
293 for (i = 0; i < int(NewQuads.size() / 4); ++i)
294 {
295
296 QuadIndices.push_back(4 * i + iLastNodeIndex);
297 QuadIndices.push_back(4 * i + 1 + iLastNodeIndex);
298 QuadIndices.push_back(4 * i + 2 + iLastNodeIndex);
299 QuadIndices.push_back(4 * i + 3 + iLastNodeIndex);
300 }
301 for (itXYZ = NewQuads.begin(); itXYZ != NewQuads.end(); ++itXYZ)
302 {
303 Mesh.AddNode(*itXYZ);
304 }
305 NewQuads.clear();
306
307 list<int> &TriIndices = Mesh.GetIndices(CMesh::TRI);
308 for (itInt = TriIndices.begin(); itInt != TriIndices.end(); )
309 {
310 itStart = itInt;
311
312 p1 = &Mesh.GetNode(*(itInt++));
313 p2 = &Mesh.GetNode(*(itInt++));
314 p3 = &Mesh.GetNode(*(itInt++));
315
316 d1 = DotProduct(itPlane->Normal, *p1) - itPlane->d;
317 d2 = DotProduct(itPlane->Normal, *p2) - itPlane->d;
318 d3 = DotProduct(itPlane->Normal, *p3) - itPlane->d;
319
320 if (d1 <= TOL && d2 <= TOL && d3 <= TOL) // The triangle lies completely outside or on the plane
321 {
322 itInt = TriIndices.erase(itStart, itInt); // Delete the triangle
323 }
324 else if (d1 >= -TOL && d2 >= -TOL && d3 >= -TOL) // The triangle lies completely inside the plane
325 {
326 // Do nothing
327 }
328 else
329 {
330 itInt = TriIndices.erase(itStart, itInt); // Delete the triangle, will need to be seperated into smaller ones
331 // Order points such that d1 >= d2 >= d3
332 bFlipped = false; // Keep track of whether the triangle is flipped or not
333 if (d2 > d1)
334 {
335 dTemp = d2; d2 = d1; d1 = dTemp;
336 pTemp = p2; p2 = p1; p1 = pTemp;
337 bFlipped = !bFlipped;
338 }
339 if (d3 > d2)
340 {
341 dTemp = d3; d3 = d2; d2 = dTemp;
342 pTemp = p3; p3 = p2; p2 = pTemp;
343 bFlipped = !bFlipped;
344 if (d2 > d1)
345 {
346 dTemp = d2; d2 = d1; d1 = dTemp;
347 pTemp = p2; p2 = p1; p1 = pTemp;
348 bFlipped = !bFlipped;
349 }
350 }
351
352 if (d2 <= TOL) // One point of the triangle lies inside the plane, the other two are outside or on the plane
353 {
354 NewTriangles.push_back(*p1);
355 u = d1 / (d1 - d2);
356 NewTriangles.push_back(*p1 + (*p2 - *p1) * u);
357 u = d1 / (d1 - d3);
358 NewTriangles.push_back(*p1 + (*p3 - *p1) * u);
359 NewTrianglesFlipped.push_back(bFlipped);
360 }
361 else if (d3 <= TOL) // Two points of the triangle lie inside the plane, the other is outside or on the plane
362 {
363 NewTriangles.push_back(*p1);
364 NewTriangles.push_back(*p2);
365 u = d2 / (d2 - d3);
366 NewTriangles.push_back(*p2 + (*p3 - *p2) * u);
367 NewTrianglesFlipped.push_back(bFlipped);
368
369 NewTriangles.push_back(*p2 + (*p3 - *p2) * u);
370 u = d1 / (d1 - d3);
371 NewTriangles.push_back(*p1 + (*p3 - *p1) * u);
372 NewTriangles.push_back(*p1);
373 NewTrianglesFlipped.push_back(bFlipped);
374 }
375 }
376 }
377
378 // Add the new triangles to the mesh, and clear the new triangles list
379 iLastNodeIndex = int(Mesh.GetNumNodes());
380 for (i = 0; i < int(NewTriangles.size() / 3); ++i)
381 {
382 // The order of the vertices determines the direction of the normal,
383 // the normal should always point outwards from the yarn. Thus if the
384 // normal was flipped we must make sure the indices are swapped so
385 // the normal is flipped back to its original state.
386 if (!NewTrianglesFlipped[i])
387 {
388 TriIndices.push_back(3 * i + iLastNodeIndex);
389 TriIndices.push_back(3 * i + 1 + iLastNodeIndex);
390 TriIndices.push_back(3 * i + 2 + iLastNodeIndex);
391 }
392 else
393 {
394 TriIndices.push_back(3 * i + iLastNodeIndex);
395 TriIndices.push_back(3 * i + 2 + iLastNodeIndex);
396 TriIndices.push_back(3 * i + 1 + iLastNodeIndex);
397 }
398 }
399 for (itXYZ = NewTriangles.begin(); itXYZ != NewTriangles.end(); ++itXYZ)
400 {
401 Mesh.AddNode(*itXYZ);
402 }
403 NewTriangles.clear();
404 NewTrianglesFlipped.clear();
405
406 vector<int> ClosedLoop;
407 if (bFillGaps)
408 FillGaps(Mesh, *itPlane, ClosedLoop);
409 }
410
411 // TODO Reinstate volume elements if want to clip through elements rather than using
412 // centre point to determine whether inside domain and retaining whole element
413
414 // Polygon mesh not currently sent to intersect mesh clip - reinstate if added to intersection mesh
415
416 Mesh.RemoveUnreferencedNodes();
417}
418
419bool CDomainPrism::ClipIntersectMeshToDomain(CMesh &Mesh, vector<CMesh> &DomainMeshes, bool bFillGaps) const
420{
421 // Clip mesh elements which are known to intersect with the domain
422 // Checks each mesh element against each domain plane
423 const double TOL = 1e-9;
424
425 TGLOGINDENT("Clipping mesh to domain");
426
427 vector<XYZ> NewTriangles;
428 vector<bool> NewTrianglesFlipped;
429 vector<XYZ>::iterator itXYZ;
430 vector<XYZ> NewQuads;
431
432 vector<PLANE>::const_iterator itPlane;
433 list<int>::iterator itStart;
434 list<int>::iterator itInt;
435 const XYZ *p1, *p2, *p3, *p4;
436 double d1, d2, d3, d4; // d represents the distance of the point to the plane (i.e. +ve inside, -ve outside, 0 on top)
437 double dTemp;
438 const XYZ *pTemp;
439 double u;
440 int i;
441 int iLastNodeIndex;
442 bool bFlipped;
443
444 vector<CMesh>::iterator itDomainMeshes;
445 bool bSaveDomainMeshes = false;
446 if (DomainMeshes.size() != 0)
447 bSaveDomainMeshes = true;
448
449 // Deal with surface elements
450 for (itPlane = m_ElementPlanes.begin(); itPlane != m_ElementPlanes.end(); ++itPlane)
451 {
452 list<int> &QuadIndices = Mesh.GetIndices(CMesh::QUAD);
453 for (itInt = QuadIndices.begin(); itInt != QuadIndices.end(); )
454 {
455 itStart = itInt;
456
457 p1 = &Mesh.GetNode(*(itInt++));
458 p2 = &Mesh.GetNode(*(itInt++));
459 p3 = &Mesh.GetNode(*(itInt++));
460 p4 = &Mesh.GetNode(*(itInt++));
461
462 d1 = DotProduct(itPlane->Normal, *p1) - itPlane->d;
463 d2 = DotProduct(itPlane->Normal, *p2) - itPlane->d;
464 d3 = DotProduct(itPlane->Normal, *p3) - itPlane->d;
465 d4 = DotProduct(itPlane->Normal, *p4) - itPlane->d;
466
467 if (d1 <= TOL && d2 <= TOL && d3 <= TOL && d4 <= TOL) // The quad lies completely on or outside the plane
468 {
469 if (fabs(d1) < TOL || fabs(d2) < TOL || fabs(d3) < TOL || fabs(d4) < TOL)
470 itInt = QuadIndices.erase(itStart, itInt); // Delete the quad
471 }
472 else if (d1 >= -TOL && d2 >= -TOL && d3 >= -TOL && d4 >= -TOL) // The quad lies completely inside the plane
473 {
474 // Do nothing
475 }
476 else if (d1 < TOL && d2 < TOL && d3 > TOL && d4 > TOL) // Points 1 & 2 outside plane
477 {
478 itInt = QuadIndices.erase(itStart, itInt); // Delete the quad
479 u = d4 / (d4 - d1);
480 NewQuads.push_back(*p4 + (*p1 - *p4) * u);
481 u = d3 / (d3 - d2);
482 NewQuads.push_back(*p3 + (*p2 - *p3) * u);
483 NewQuads.push_back(*p3);
484 NewQuads.push_back(*p4);
485 }
486 else if (d2 < TOL && d3 < TOL && d4 > TOL && d1 > TOL) // Points 2 & 3 outside plane
487 {
488 itInt = QuadIndices.erase(itStart, itInt); // Delete the quad
489 NewQuads.push_back(*p1);
490 u = d1 / (d1 - d2);
491 NewQuads.push_back(*p1 + (*p2 - *p1) * u);
492 u = d4 / (d4 - d3);
493 NewQuads.push_back(*p4 + (*p3 - *p4) * u);
494 NewQuads.push_back(*p4);
495 }
496 else if (d3 < TOL && d4 < TOL && d1 > TOL && d2 > TOL) // Points 3 & 4 outside plane
497 {
498 itInt = QuadIndices.erase(itStart, itInt); // Delete the quad
499 NewQuads.push_back(*p1);
500 NewQuads.push_back(*p2);
501 u = d2 / (d2 - d3);
502 NewQuads.push_back(*p2 + (*p3 - *p2) * u);
503 u = d1 / (d1 - d4);
504 NewQuads.push_back(*p1 + (*p4 - *p1) * u);
505 }
506 else if (d4 < TOL && d1 < TOL && d2 > TOL && d3 > TOL) // Points 4 & 1 outside plane
507 {
508 itInt = QuadIndices.erase(itStart, itInt); // Delete the quad
509 u = d2 / (d2 - d1);
510 NewQuads.push_back(*p2 + (*p1 - *p2) * u);
511 NewQuads.push_back(*p2);
512 NewQuads.push_back(*p3);
513 u = d3 / (d3 - d4);
514 NewQuads.push_back(*p3 + (*p4 - *p3) * u);
515
516 }
517 else // Convert the quad to a triangle for trimming if 1 or 3 points inside plane
518 {
519 itInt = Mesh.ConvertQuadtoTriangles(itStart);
520 }
521 }
522
523 // Add the new quads to the mesh, and clear the new quad list
524 iLastNodeIndex = int(Mesh.GetNumNodes());
525 for (i = 0; i < int(NewQuads.size() / 4); ++i)
526 {
527
528 QuadIndices.push_back(4 * i + iLastNodeIndex);
529 QuadIndices.push_back(4 * i + 1 + iLastNodeIndex);
530 QuadIndices.push_back(4 * i + 2 + iLastNodeIndex);
531 QuadIndices.push_back(4 * i + 3 + iLastNodeIndex);
532 }
533 for (itXYZ = NewQuads.begin(); itXYZ != NewQuads.end(); ++itXYZ)
534 {
535 Mesh.AddNode(*itXYZ);
536 }
537 NewQuads.clear();
538
539 list<int> &TriIndices = Mesh.GetIndices(CMesh::TRI);
540 for (itInt = TriIndices.begin(); itInt != TriIndices.end(); )
541 {
542 itStart = itInt;
543
544 p1 = &Mesh.GetNode(*(itInt++));
545 p2 = &Mesh.GetNode(*(itInt++));
546 p3 = &Mesh.GetNode(*(itInt++));
547
548 d1 = DotProduct(itPlane->Normal, *p1) - itPlane->d;
549 d2 = DotProduct(itPlane->Normal, *p2) - itPlane->d;
550 d3 = DotProduct(itPlane->Normal, *p3) - itPlane->d;
551
552 if (d1 <= TOL && d2 <= TOL && d3 <= TOL) // The triangle lies completely outside or on the plane
553 {
554 itInt = TriIndices.erase(itStart, itInt); // Delete the triangle
555 }
556 else if (d1 >= -TOL && d2 >= -TOL && d3 >= -TOL) // The triangle lies completely inside the plane
557 {
558 // Do nothing
559 }
560 else
561 {
562 itInt = TriIndices.erase(itStart, itInt); // Delete the triangle, will need to be seperated into smaller ones
563 // Order points such that d1 >= d2 >= d3
564 bFlipped = false; // Keep track of whether the triangle is flipped or not
565 if (d2 > d1)
566 {
567 dTemp = d2; d2 = d1; d1 = dTemp;
568 pTemp = p2; p2 = p1; p1 = pTemp;
569 bFlipped = !bFlipped;
570 }
571 if (d3 > d2)
572 {
573 dTemp = d3; d3 = d2; d2 = dTemp;
574 pTemp = p3; p3 = p2; p2 = pTemp;
575 bFlipped = !bFlipped;
576 if (d2 > d1)
577 {
578 dTemp = d2; d2 = d1; d1 = dTemp;
579 pTemp = p2; p2 = p1; p1 = pTemp;
580 bFlipped = !bFlipped;
581 }
582 }
583
584 if (d2 <= TOL) // One point of the triangle lies inside the plane, the other two are outside or on the plane
585 {
586 NewTriangles.push_back(*p1);
587 u = d1 / (d1 - d2);
588 NewTriangles.push_back(*p1 + (*p2 - *p1) * u);
589 u = d1 / (d1 - d3);
590 NewTriangles.push_back(*p1 + (*p3 - *p1) * u);
591 NewTrianglesFlipped.push_back(bFlipped);
592 }
593 else if (d3 <= TOL) // Two points of the triangle lie inside the plane, the other is outside or on the plane
594 {
595 NewTriangles.push_back(*p1);
596 NewTriangles.push_back(*p2);
597 u = d2 / (d2 - d3);
598 NewTriangles.push_back(*p2 + (*p3 - *p2) * u);
599 NewTrianglesFlipped.push_back(bFlipped);
600
601 NewTriangles.push_back(*p2 + (*p3 - *p2) * u);
602 u = d1 / (d1 - d3);
603 NewTriangles.push_back(*p1 + (*p3 - *p1) * u);
604 NewTriangles.push_back(*p1);
605 NewTrianglesFlipped.push_back(bFlipped);
606 }
607 }
608 }
609
610 // Add the new triangles to the mesh, and clear the new triangles list
611 iLastNodeIndex = int(Mesh.GetNumNodes());
612 for (i = 0; i < int(NewTriangles.size() / 3); ++i)
613 {
614 // The order of the vertices determines the direction of the normal,
615 // the normal should always point outwards from the yarn. Thus if the
616 // normal was flipped we must make sure the indices are swapped so
617 // the normal is flipped back to its original state.
618 if (!NewTrianglesFlipped[i])
619 {
620 TriIndices.push_back(3 * i + iLastNodeIndex);
621 TriIndices.push_back(3 * i + 1 + iLastNodeIndex);
622 TriIndices.push_back(3 * i + 2 + iLastNodeIndex);
623 }
624 else
625 {
626 TriIndices.push_back(3 * i + iLastNodeIndex);
627 TriIndices.push_back(3 * i + 2 + iLastNodeIndex);
628 TriIndices.push_back(3 * i + 1 + iLastNodeIndex);
629 }
630 }
631 for (itXYZ = NewTriangles.begin(); itXYZ != NewTriangles.end(); ++itXYZ)
632 {
633 Mesh.AddNode(*itXYZ);
634 }
635 NewTriangles.clear();
636 NewTrianglesFlipped.clear();
637
638 vector<int> ClosedLoop;
639 if (bFillGaps)
640 {
641 if (!FillGaps(Mesh, *itPlane, ClosedLoop, false))
642 return false;
643 }
644 if (bSaveDomainMeshes)
645 {
646 if (ClosedLoop.size() > 0)
647 {
648 for (itDomainMeshes = DomainMeshes.begin(); itDomainMeshes != DomainMeshes.end(); itDomainMeshes++)
649 {
650 XYZ Points[3];
651 XYZ Normal;
652 double dPlane;
653 for (int i = 0; i < 3; i++)
654 {
655 Points[i] = (*itDomainMeshes).GetNode(i);
656 }
657 Normal = CrossProduct((Points[0] - Points[1]), Points[2] - Points[1]);
658 dPlane = DotProduct(Points[2], Normal);
659 if (fabs(DotProduct(Mesh.GetNode(ClosedLoop[0]), Normal) - dPlane) < 0.000001)
660 {
661 int iIndex = (*itDomainMeshes).GetNumNodes();
662 int iPolyStart = iIndex;
663 vector<int>::iterator itClosedLoop;
664 vector<int> Indices;
665 int iStartInd = ClosedLoop[0];
666 for (itClosedLoop = ClosedLoop.begin(); itClosedLoop != ClosedLoop.end(); itClosedLoop++)
667 {
668 if (*itClosedLoop == iStartInd && iIndex > iPolyStart) // Back to start of loop
669 {
670 Indices.push_back(iPolyStart);
671 }
672 else
673 {
674 Indices.push_back(iIndex++);
675 (*itDomainMeshes).AddNode(Mesh.GetNode(*itClosedLoop));
676 }
677 }
678 if (ClosedLoop[ClosedLoop.size() - 1] != iStartInd) // Close loop if not already closed
679 Indices.push_back(iPolyStart);
680 (*itDomainMeshes).AddElement(CMesh::POLYGON, Indices);
681 break;
682 }
683 }
684 }
685 }
686 }
687
688 // TODO Reinstate volume elements if want to clip through elements rather than using
689 // centre point to determine whether inside domain and retaining whole element
690
691 // Polygon mesh not currently sent to intersect mesh clip - reinstate if added to intersection mesh
692
693 Mesh.RemoveUnreferencedNodes();
694 return true;
695}
696
697void CDomainPrism::ClipMeshToDomain(CMesh &Mesh, bool bFillGaps) const
698{
699 const double TOL = 1e-9;
700
701 TGLOGINDENT("Clipping mesh to domain");
702
703 vector<XYZ> NewTriangles;
704
705 vector<XYZ>::iterator itXYZ;
706 vector<XYZ> NewQuads;
707
708 vector<PLANE>::const_iterator itPlane;
709 list<int>::iterator itStart;
710 list<int>::iterator itInt;
711 const XYZ *p1, *p2, *p3, *p4;
712
713 int i;
714 int iLastNodeIndex;
715 bool b1, b2, b3, b4;
716
717 CMesh IntersectMesh; // Mesh used to save Mesh elements (from yarn) which intersect with the domain planes
718 vector<XYZ> IntersectQuads;
719
720 // Deal with surface elements
721
722 // Quad elements
723 list<int> &QuadIndices = Mesh.GetIndices(CMesh::QUAD);
724 for (itInt = QuadIndices.begin(); itInt != QuadIndices.end(); )
725 {
726 itStart = itInt;
727
728 double d1= 0, d2=0, d3=0, d4 = 0;
729 p1 = &Mesh.GetNode(*(itInt++));
730 p2 = &Mesh.GetNode(*(itInt++));
731 p3 = &Mesh.GetNode(*(itInt++));
732 p4 = &Mesh.GetNode(*(itInt++));
733
734 b1 = m_Yarn.PointInsideYarn(*p1, NULL, NULL, NULL, &d1);
735 b2 = m_Yarn.PointInsideYarn(*p2, NULL, NULL, NULL, &d2);
736 b3 = m_Yarn.PointInsideYarn(*p3, NULL, NULL, NULL, &d3);
737 b4 = m_Yarn.PointInsideYarn(*p4, NULL, NULL, NULL, &d4);
738
739
740 if (!b1 && !b2 && !b3 && !b4) // The quad lies completely on or outside the plane
741 {
742 itInt = QuadIndices.erase(itStart, itInt); // Delete the quad
743 }
744 else if (b1 && b2 && b3 && b4) // The quad lies completely inside the plane
745 {
746 // Do nothing
747 }
748 else // quad to intersections
749 {
750 // For concave shapes PointInsideYarn may return incorrect value (point may be on 'incorrect' side of plane for one part of shape even though actually inside)
751 // So, save points for questionable elements to IntersectQuads vector ( to allow later checking of quad against each domain plane )
752 itInt = QuadIndices.erase(itStart, itInt); // Delete the quad
753 IntersectQuads.push_back(*p1);
754 IntersectQuads.push_back(*p2);
755 IntersectQuads.push_back(*p3);
756 IntersectQuads.push_back(*p4);
757 }
758 }
759
760 // Add points for intersecting quads to intersection mesh
761 list<int> &IntersectQuadIndices = IntersectMesh.GetIndices(CMesh::QUAD);
762 for (i = 0; i < int(IntersectQuads.size() / 4); ++i)
763 {
764 IntersectQuadIndices.push_back(4 * i);
765 IntersectQuadIndices.push_back(4 * i + 1);
766 IntersectQuadIndices.push_back(4 * i + 2);
767 IntersectQuadIndices.push_back(4 * i + 3);
768 }
769 for (itXYZ = IntersectQuads.begin(); itXYZ != IntersectQuads.end(); ++itXYZ)
770 {
771 IntersectMesh.AddNode(*itXYZ);
772 }
773 IntersectQuads.clear();
774
775
776 // Triangle elements
777 vector<XYZ> IntersectTriangles;
778
779 list<int> &TriIndices = Mesh.GetIndices(CMesh::TRI);
780 for (itInt = TriIndices.begin(); itInt != TriIndices.end(); )
781 {
782 itStart = itInt;
783
784 double d1 = 0, d2 = 0, d3 = 0, d4 = 0;
785 p1 = &Mesh.GetNode(*(itInt++));
786 p2 = &Mesh.GetNode(*(itInt++));
787 p3 = &Mesh.GetNode(*(itInt++));
788
789 b1 = m_Yarn.PointInsideYarn(*p1, NULL, NULL, NULL, &d1);
790 b2 = m_Yarn.PointInsideYarn(*p2, NULL, NULL, NULL, &d2);
791 b3 = m_Yarn.PointInsideYarn(*p3, NULL, NULL, NULL, &d3);
792
793 if (!b1 && !b2 && !b3) // The triangle lies completely outside or on the plane
794 {
795 itInt = TriIndices.erase(itStart, itInt); // Delete the triangle
796 }
797 else if (b1 && b2 && b3) // The triangle lies completely inside the plane
798 {
799 // Do nothing
800 }
801 else
802 {
803 itInt = TriIndices.erase(itStart, itInt);
804 IntersectTriangles.push_back(*p1);
805 IntersectTriangles.push_back(*p2);
806 IntersectTriangles.push_back(*p3);
807 }
808 }
809
810 // Add points for intersecting triangular elements to intersection mesh
811 list<int> &IntersectTriangleIndices = IntersectMesh.GetIndices(CMesh::TRI);
812 iLastNodeIndex = IntersectMesh.GetNumNodes();
813 for (i = 0; i < int(IntersectQuads.size() / 4); ++i)
814 {
815 IntersectTriangleIndices.push_back(4 * i + iLastNodeIndex);
816 IntersectTriangleIndices.push_back(4 * i + 1 + iLastNodeIndex);
817 IntersectTriangleIndices.push_back(4 * i + 2 + iLastNodeIndex);
818 IntersectTriangleIndices.push_back(4 * i + 3 + iLastNodeIndex);
819 }
820 for (itXYZ = IntersectTriangles.begin(); itXYZ != IntersectTriangles.end(); ++itXYZ)
821 {
822 IntersectMesh.AddNode(*itXYZ);
823 }
824 IntersectTriangles.clear();
825
826 // Deal with volume elements
827 int iNumNodes;
828 double d;
829 vector<CMesh::ELEMENT_TYPE> VolumeElements;
830 vector<CMesh::ELEMENT_TYPE>::iterator itElementType;
831 VolumeElements.push_back(CMesh::TET);
832 VolumeElements.push_back(CMesh::PYRAMID);
833 VolumeElements.push_back(CMesh::WEDGE);
834 VolumeElements.push_back(CMesh::HEX);
835 VolumeElements.push_back(CMesh::QUADRATIC_TET);
836 //for (itPlane = m_Planes.begin(); itPlane != m_Planes.end(); ++itPlane)
837 //{
838 for (itElementType = VolumeElements.begin(); itElementType != VolumeElements.end(); ++itElementType)
839 {
840 list<int> &Indices = Mesh.GetIndices(*itElementType);
841 for (itInt = Indices.begin(); itInt != Indices.end(); )
842 {
843 iNumNodes = CMesh::GetNumNodes(*itElementType);
844 itStart = itInt;
845
846 XYZ Center;
847 for (i = 0; i < iNumNodes; ++i)
848 {
849 Center += Mesh.GetNode(*(itInt++));
850 }
851 Center /= iNumNodes;
852
853 b1 = m_Yarn.PointInsideYarn(Center, NULL, NULL, NULL, &d);
854 // d = DotProduct(itPlane->Normal, Center) - itPlane->d;
855 // if (d < 0)
856 if (!b1)
857 itInt = Indices.erase(itStart, itInt); // Delete the volume element
858 }
859 }
860 //}
861
862
863 // Does it need to deal with polygon elements?
864
865 IntersectMesh.SaveToVTK("IntersectionMesh"); // Debug check for viewing of intersection mesh
866
867 // Clip the elements of the yarn mesh which intersect with the domain and add back into the yarn mesh
868 ClipIntersectMeshToDomain(IntersectMesh);
869 Mesh.InsertMesh(IntersectMesh);
870
871}
872
873bool CDomainPrism::ClipMeshToDomain(CMesh &Mesh, vector<CMesh> &DomainMeshes, bool bFillGaps) const
874{
875 const double TOL = 1e-9;
876
877 TGLOGINDENT("Clipping mesh to domain");
878
879 vector<XYZ> NewTriangles;
880
881 vector<XYZ>::iterator itXYZ;
882 vector<XYZ> NewQuads;
883
884 vector<PLANE>::const_iterator itPlane;
885 list<int>::iterator itStart;
886 list<int>::iterator itInt;
887 const XYZ *p1, *p2, *p3, *p4;
888
889 int i;
890 int iLastNodeIndex;
891 bool b1, b2, b3, b4;
892
893
894 CMesh IntersectMesh; // Mesh used to save Mesh elements (from yarn) which intersect with the domain planes
895 vector<XYZ> IntersectQuads;
896
897 // Deal with surface elements
898
899 // Quad elements
900 list<int> &QuadIndices = Mesh.GetIndices(CMesh::QUAD);
901 for (itInt = QuadIndices.begin(); itInt != QuadIndices.end(); )
902 {
903 itStart = itInt;
904
905 double d1 = 0, d2 = 0, d3 = 0, d4 = 0;
906 p1 = &Mesh.GetNode(*(itInt++));
907 p2 = &Mesh.GetNode(*(itInt++));
908 p3 = &Mesh.GetNode(*(itInt++));
909 p4 = &Mesh.GetNode(*(itInt++));
910
911 b1 = m_Yarn.PointInsideYarn(*p1, NULL, NULL, NULL, &d1);
912 b2 = m_Yarn.PointInsideYarn(*p2, NULL, NULL, NULL, &d2);
913 b3 = m_Yarn.PointInsideYarn(*p3, NULL, NULL, NULL, &d3);
914 b4 = m_Yarn.PointInsideYarn(*p4, NULL, NULL, NULL, &d4);
915
916
917 if (!b1 && !b2 && !b3 && !b4) // The quad lies completely on or outside the plane
918 {
919 itInt = QuadIndices.erase(itStart, itInt); // Delete the quad
920 }
921 else if (b1 && b2 && b3 && b4) // The quad lies completely inside the plane
922 {
923 // Do nothing
924 }
925 else // quad to intersections
926 {
927 // For concave shapes PointInsideYarn may return incorrect value (point may be on 'incorrect' side of plane for one part of shape even though actually inside)
928 // So, save points for questionable elements to IntersectQuads vector ( to allow later checking of quad against each domain plane )
929 itInt = QuadIndices.erase(itStart, itInt); // Delete the quad
930 IntersectQuads.push_back(*p1);
931 IntersectQuads.push_back(*p2);
932 IntersectQuads.push_back(*p3);
933 IntersectQuads.push_back(*p4);
934 }
935 }
936
937 // Add points for intersecting quads to intersection mesh
938 list<int> &IntersectQuadIndices = IntersectMesh.GetIndices(CMesh::QUAD);
939 for (i = 0; i < int(IntersectQuads.size() / 4); ++i)
940 {
941 IntersectQuadIndices.push_back(4 * i);
942 IntersectQuadIndices.push_back(4 * i + 1);
943 IntersectQuadIndices.push_back(4 * i + 2);
944 IntersectQuadIndices.push_back(4 * i + 3);
945 }
946 for (itXYZ = IntersectQuads.begin(); itXYZ != IntersectQuads.end(); ++itXYZ)
947 {
948 IntersectMesh.AddNode(*itXYZ);
949 }
950 IntersectQuads.clear();
951
952
953 // Triangle elements
954 vector<XYZ> IntersectTriangles;
955
956 list<int> &TriIndices = Mesh.GetIndices(CMesh::TRI);
957 for (itInt = TriIndices.begin(); itInt != TriIndices.end(); )
958 {
959 itStart = itInt;
960
961 double d1 = 0, d2 = 0, d3 = 0, d4 = 0;
962 p1 = &Mesh.GetNode(*(itInt++));
963 p2 = &Mesh.GetNode(*(itInt++));
964 p3 = &Mesh.GetNode(*(itInt++));
965
966 b1 = m_Yarn.PointInsideYarn(*p1, NULL, NULL, NULL, &d1);
967 b2 = m_Yarn.PointInsideYarn(*p2, NULL, NULL, NULL, &d2);
968 b3 = m_Yarn.PointInsideYarn(*p3, NULL, NULL, NULL, &d3);
969
970 if (!b1 && !b2 && !b3) // The triangle lies completely outside or on the plane
971 {
972 itInt = TriIndices.erase(itStart, itInt); // Delete the triangle
973 }
974 else if (b1 && b2 && b3) // The triangle lies completely inside the plane
975 {
976 // Do nothing
977 }
978 else
979 {
980 itInt = TriIndices.erase(itStart, itInt);
981 IntersectTriangles.push_back(*p1);
982 IntersectTriangles.push_back(*p2);
983 IntersectTriangles.push_back(*p3);
984 }
985 }
986
987 // Add points for intersecting triangular elements to intersection mesh
988 list<int> &IntersectTriangleIndices = IntersectMesh.GetIndices(CMesh::TRI);
989 iLastNodeIndex = IntersectMesh.GetNumNodes();
990 for (i = 0; i < int(IntersectQuads.size() / 4); ++i)
991 {
992 IntersectTriangleIndices.push_back(4 * i + iLastNodeIndex);
993 IntersectTriangleIndices.push_back(4 * i + 1 + iLastNodeIndex);
994 IntersectTriangleIndices.push_back(4 * i + 2 + iLastNodeIndex);
995 IntersectTriangleIndices.push_back(4 * i + 3 + iLastNodeIndex);
996 }
997 for (itXYZ = IntersectTriangles.begin(); itXYZ != IntersectTriangles.end(); ++itXYZ)
998 {
999 IntersectMesh.AddNode(*itXYZ);
1000 }
1001 IntersectTriangles.clear();
1002
1003 // Deal with volume elements
1004 int iNumNodes;
1005 double d;
1006 vector<CMesh::ELEMENT_TYPE> VolumeElements;
1007 vector<CMesh::ELEMENT_TYPE>::iterator itElementType;
1008 VolumeElements.push_back(CMesh::TET);
1009 VolumeElements.push_back(CMesh::PYRAMID);
1010 VolumeElements.push_back(CMesh::WEDGE);
1011 VolumeElements.push_back(CMesh::HEX);
1012 VolumeElements.push_back(CMesh::QUADRATIC_TET);
1013 //for (itPlane = m_Planes.begin(); itPlane != m_Planes.end(); ++itPlane)
1014 //{
1015 for (itElementType = VolumeElements.begin(); itElementType != VolumeElements.end(); ++itElementType)
1016 {
1017 list<int> &Indices = Mesh.GetIndices(*itElementType);
1018 for (itInt = Indices.begin(); itInt != Indices.end(); )
1019 {
1020 iNumNodes = CMesh::GetNumNodes(*itElementType);
1021 itStart = itInt;
1022
1023 XYZ Center;
1024 for (i = 0; i < iNumNodes; ++i)
1025 {
1026 Center += Mesh.GetNode(*(itInt++));
1027 }
1028 Center /= iNumNodes;
1029
1030 b1 = m_Yarn.PointInsideYarn(Center, NULL, NULL, NULL, &d);
1031 // d = DotProduct(itPlane->Normal, Center) - itPlane->d;
1032 // if (d < 0)
1033 if (!b1)
1034 itInt = Indices.erase(itStart, itInt); // Delete the volume element
1035 }
1036 }
1037 //}
1038
1039
1040 // Does it need to deal with polygon elements?
1041
1042 IntersectMesh.SaveToVTK("IntersectionMesh"); // Debug check for viewing of intersection mesh
1043
1044 // Clip the elements of the yarn mesh which intersect with the domain and add back into the yarn mesh
1045 if (!ClipIntersectMeshToDomain(IntersectMesh, DomainMeshes, bFillGaps))
1046 return false;
1047 Mesh.InsertMesh(IntersectMesh);
1048 return true;
1049}
1050
1051bool CDomainPrism::FillGaps(CMesh &Mesh, const PLANE &Plane, vector<int> &Polygon, bool bMeshGaps)
1052{
1053 const double TOL = 1e-9;
1054
1055 int i1, i2, i3, i4;
1056 const XYZ *p1, *p2, *p3, *p4;
1057 double d1, d2, d3, d4; // d represents the distance of the point to the plane (i.e. +ve inside, -ve outside, 0 on top)
1058
1059 vector<pair<int, int> > Segments;
1060
1061 // Merge the nodes together and remove degenerate triangles
1062 Mesh.MergeNodes();
1063 Mesh.RemoveDegenerateTriangles();
1064
1065 // Build a list of segments which lie on the plane
1066 list<int>::iterator itInt;
1067 // Check each quad to see if any of the edges lie on the plane
1068 list<int> &QuadIndices = Mesh.GetIndices(CMesh::QUAD);
1069 for (itInt = QuadIndices.begin(); itInt != QuadIndices.end(); )
1070 {
1071 i1 = *(itInt++);
1072 i2 = *(itInt++);
1073 i3 = *(itInt++);
1074 i4 = *(itInt++);
1075
1076 p1 = &Mesh.GetNode(i1);
1077 p2 = &Mesh.GetNode(i2);
1078 p3 = &Mesh.GetNode(i3);
1079 p4 = &Mesh.GetNode(i4);
1080
1081 d1 = DotProduct(Plane.Normal, *p1) - Plane.d;
1082 d2 = DotProduct(Plane.Normal, *p2) - Plane.d;
1083 d3 = DotProduct(Plane.Normal, *p3) - Plane.d;
1084 d4 = DotProduct(Plane.Normal, *p4) - Plane.d;
1085
1086 d1 = abs(d1);
1087 d2 = abs(d2);
1088 d3 = abs(d3);
1089 d4 = abs(d4);
1090
1091 // Add the segments which lie on the plane
1092 // The order of the segment indices is important, it tells us
1093 // which side is inside and which side is outside
1094 // Here the segments are added clockwise when viewed along the plane normal
1095 if (d1 <= TOL && d2 <= TOL)
1096 {
1097 Segments.push_back(pair<int, int>(i1, i2));
1098 }
1099 if (d2 <= TOL && d3 <= TOL)
1100 {
1101 Segments.push_back(pair<int, int>(i2, i3));
1102 }
1103 if (d3 <= TOL && d4 <= TOL)
1104 {
1105 Segments.push_back(pair<int, int>(i3, i4));
1106 }
1107 if (d4 <= TOL && d1 <= TOL)
1108 {
1109 Segments.push_back(pair<int, int>(i4, i1));
1110 }
1111 }
1112 // Check each triangle to find edges that lie on the plane
1113 list<int> &TriIndices = Mesh.GetIndices(CMesh::TRI);
1114 for (itInt = TriIndices.begin(); itInt != TriIndices.end(); )
1115 {
1116 i1 = *(itInt++);
1117 i2 = *(itInt++);
1118 i3 = *(itInt++);
1119
1120 p1 = &Mesh.GetNode(i1);
1121 p2 = &Mesh.GetNode(i2);
1122 p3 = &Mesh.GetNode(i3);
1123
1124 d1 = DotProduct(Plane.Normal, *p1) - Plane.d;
1125 d2 = DotProduct(Plane.Normal, *p2) - Plane.d;
1126 d3 = DotProduct(Plane.Normal, *p3) - Plane.d;
1127
1128 d1 = abs(d1);
1129 d2 = abs(d2);
1130 d3 = abs(d3);
1131
1132 // Add the segments which lie on the plane
1133 // The order of the segment indices is important, it tells us
1134 // which side is inside and which side is outside
1135 // Here the segments are added clockwise when viewed along the plane normal
1136 if (d1 <= TOL && d2 <= TOL)
1137 {
1138 Segments.push_back(pair<int, int>(i1, i2));
1139 }
1140 if (d2 <= TOL && d3 <= TOL)
1141 {
1142 Segments.push_back(pair<int, int>(i2, i3));
1143 }
1144 if (d3 <= TOL && d1 <= TOL)
1145 {
1146 Segments.push_back(pair<int, int>(i3, i1));
1147 }
1148 }
1149
1150 vector<pair<int, int> >::iterator itSegment;
1151
1152 // Find closed loops
1153 int iIndex, iFirstIndex;
1154 bool bFound;
1155
1156#ifdef _DEBUG
1157 // This code will help to find out where the source of the problem comes from.
1158 // It basically checks how many times each node has been referenced by the segments.
1159 // Each node should be referenced twice, otherwise there will be a problem.
1160
1161 map<int, int> Indices;
1162 map<int, int>::iterator itIndex;
1163
1164 for (itSegment = Segments.begin(); itSegment != Segments.end(); ++itSegment)
1165 {
1166 Indices[itSegment->first]++;
1167 Indices[itSegment->second]++;
1168 }
1169#endif
1170
1171 while (Segments.size())
1172 {
1173 vector<int> ClosedLoop;
1174 // Start at a random index and go round counterclockwise until a full circle is done
1175 // Indices are added to the ClosedLoop list with each index specified only once.
1176 // As segments are followed they are removed from the segments list.
1177 iFirstIndex = iIndex = Segments.begin()->first;
1178 do
1179 {
1180 bFound = false;
1181 for (itSegment = Segments.begin(); itSegment != Segments.end(); ++itSegment)
1182 {
1183 // Follow segments counter-clockwise only
1184 if (itSegment->second == iIndex)
1185 {
1186 // Adjust the index to go follow the segment
1187 iIndex = itSegment->first;
1188 // Delete the segment, it is no longer needed
1189 itSegment = Segments.erase(itSegment);
1190 // Found the segment, now stop searching
1191 bFound = true;
1192 break;
1193 }
1194 // This part is commented because we only want to go counter-clockwise
1195 /* else if (itSegment->first == iIndex)
1196 {
1197 iIndex = itSegment->second;
1198 itSegment = Segments.erase(itSegment);
1199 bFound = true;
1200 break;
1201 }*/
1202 }
1203 if (bFound)
1204 {
1205 // Add the index to the loop
1206 ClosedLoop.push_back(iIndex);
1207 // If the index is the same as the first index then a full circle
1208 // has been made, its time to exit the loop
1209 if (iFirstIndex == iIndex)
1210 break;
1211 }
1212 } while (bFound);
1213
1214 // If a dead end was reached the bFound will be false. This means that the segments
1215 // do not form a fully closed loop. This can occur if nodes are not merged together
1216 // correctly (two nodes may occupy the same position where each one is only referenced
1217 // once which will cause a break in the loop).
1218 if (!bFound)
1219 {
1220 // Report the error
1221 TGERROR("Unable to fill gaps satisfactorily");
1222
1223#ifdef _DEBUG
1224 // Print out the total number of nodes
1225 TGLOG("Number of nodes: " << Indices.size());
1226
1227 // Print out the number of nodes referenced more than once
1228 // for (itIndex = Indices.begin(); itIndex != Indices.end(); ++itIndex)
1229 // {
1230 // if (itIndex->second > 1)
1231 // cout << "Node " << itIndex->first << " referenced " << itIndex->second << " times. (" << Mesh.m_Nodes[itIndex->first] << ")" << endl;
1232 // }
1233 // Print out the number of nodes referenced only once
1234 for (itIndex = Indices.begin(); itIndex != Indices.end(); ++itIndex)
1235 {
1236 if (itIndex->second != 2)
1237 TGLOG("Node " << itIndex->first << " referenced " << itIndex->second << " times. (" << Mesh.GetNode(itIndex->first) << ")");
1238 }
1239
1240 //assert(false);
1241#endif
1242 return false;
1243 }
1244
1245 // Check for two points next to each other which are the same
1246 // Happens on issue #2 - would be better to find cause
1247 vector<int>::iterator itCurrent = ClosedLoop.begin();
1248 vector<int>::iterator itNext = itCurrent + 1;
1249 while (itNext != ClosedLoop.end())
1250 {
1251 if (*itCurrent == *itNext)
1252 {
1253 itNext = ClosedLoop.erase(itNext);
1254 }
1255 else
1256 {
1257 ++itCurrent;
1258 ++itNext;
1259 }
1260 }
1261
1262 // Check for spike in loop where points two apart are the same
1263 // Would be better to find the root cause of this happening
1264
1265 vector<int>::iterator itPrev = ClosedLoop.begin();
1266 itCurrent = itPrev + 1;
1267 itNext = itCurrent + 1;
1268
1269 while (itNext != ClosedLoop.end())
1270 {
1271 if (*itPrev == *itNext)
1272 {
1273 ClosedLoop.erase(itCurrent, itNext + 1);
1274 if (itPrev != ClosedLoop.begin())
1275 {
1276 itCurrent = itPrev;
1277 --itPrev;
1278 itNext = itCurrent + 1;
1279 }
1280 }
1281 else
1282 {
1283 // Go to the next set of 3 points
1284 ++itPrev;
1285 ++itCurrent;
1286 ++itNext;
1287 if (itPrev == ClosedLoop.end())
1288 itPrev = ClosedLoop.begin();
1289 }
1290 }
1291
1292 Polygon.insert(Polygon.begin(), ClosedLoop.begin(), ClosedLoop.end());
1293 // Mesh the closed loop
1294 if (bMeshGaps) // If just saving closed loop don't want to fill in end
1295 Mesh.MeshClosedLoop(Plane.Normal, ClosedLoop);
1296 };
1297
1298 return true;
1299
1300}
1301
1302void CDomainPrism::GetPolygonLimits( XYZ &StartPoint, XYZ *SizeVecs )
1303{
1304 pair<XY, XY> XYCorners;
1305 XYZ Node0, Node1;
1306
1307 Node0 = m_Yarn.GetNode(0)->GetPosition();
1308 Node1 = m_Yarn.GetNode(1)->GetPosition();
1309 GetMinMaxXY(m_Points, XYCorners.first, XYCorners.second);
1310
1311 XY CornerXDir(XYCorners.second.x, XYCorners.first.y);
1312 XY CornerZDir(XYCorners.first.x, XYCorners.second.y);
1313
1314 const CNode *node = m_Yarn.GetNode(0);
1315
1316 const vector<CSlaveNode> &SlaveNodes = m_Yarn.GetSlaveNodes(CYarn::SURFACE);
1317
1318 XYZ Up = SlaveNodes[0].GetUp();
1319 XYZ Side = SlaveNodes[0].GetSide();
1320
1321 // Rotate the 2d section point to the global 3d coordinate system
1322 StartPoint = Side * XYCorners.first.x;
1323 StartPoint += Up * XYCorners.first.y;
1324 SizeVecs[0] = Side * CornerXDir.x;
1325 SizeVecs[0] += Up * CornerXDir.y;
1326 SizeVecs[2] = Side * CornerZDir.x;
1327 SizeVecs[2] += Up * CornerZDir.y;
1328
1329 // Translate the point to its global position
1330 StartPoint += Node0;
1331 SizeVecs[0] += Node0;
1332 SizeVecs[0] -= StartPoint;
1333 SizeVecs[2] += Node0;
1334 SizeVecs[2] -= StartPoint;
1335
1336 SizeVecs[1] = Node1 - Node0;
1337}
1338
1339void CDomainPrism::GetMeshWithPolygonEnd(CMesh &Mesh)
1340{
1341 vector<int> StartIndices, EndIndices;
1342 int NumSectionPoints = m_Points.size();
1343 // Create vector of indices for polygon ends. Will need to update if change from 2 slave nodes in m_Yarn
1344 for (int i = 0; i < NumSectionPoints; ++i)
1345 {
1346 StartIndices.push_back(i);
1347 EndIndices.push_back(i + NumSectionPoints);
1348 }
1349 StartIndices.push_back(0);
1350 EndIndices.push_back(NumSectionPoints);
1351
1352 if (m_Mesh.GetNumNodes() == 0)
1353 {
1354 BuildMesh();
1355 }
1356 Mesh = m_Mesh;
1357
1358 Mesh.AddElement(CMesh::POLYGON, StartIndices);
1359 Mesh.AddElement(CMesh::POLYGON, EndIndices);
1360 Mesh.RemoveElementType(CMesh::TRI); // Don't need triangles on end mesh
1361}
TOL
#define TOL
Definition: AdjustMeshInterference.cpp:27
TGLOGINDENT
#define TGLOGINDENT(MESSAGE)
Combines the TGLOG macro and TGLOGAUTOINDENT macro in one.
Definition: Logger.h:68
TGERROR
#define TGERROR(MESSAGE)
Macros used to report the file name and line number to the TexGenError and TexGenLog functions.
Definition: Logger.h:29
TGLOG
#define TGLOG(MESSAGE)
Definition: Logger.h:36
NULL
#define NULL
Definition: ShinyConfig.h:50
TexGen::CDomain
Abstract base class representing the domain in which a textile cell may lie.
Definition: Domain.h:34
TexGen::CDomain::PopulateTiXmlElement
virtual void PopulateTiXmlElement(TiXmlElement &Element, OUTPUT_TYPE OutputType=OUTPUT_STANDARD) const
Used for saving data to XML.
Definition: Domain.cpp:43
TexGen::CDomain::m_Mesh
CMesh m_Mesh
A mesh representing the domain as a surface mesh.
Definition: Domain.h:112
TexGen::CDomainPrism::m_ElementPlanes
vector< PLANE > m_ElementPlanes
Planes corresponding to mesh elements.
Definition: DomainPrism.h:127
TexGen::CDomainPrism::RemoveDuplicatePlanes
void RemoveDuplicatePlanes()
Iterate through m_ElementPlanes to remove duplicates.
Definition: DomainPrism.cpp:165
TexGen::CDomainPrism::GetPlane
bool GetPlane(XYZ *points, PLANE &plane)
Generate a plane from a vector of co-planar points.
Definition: DomainPrism.cpp:152
TexGen::CDomainPrism::GetMeshWithPolygonEnd
void GetMeshWithPolygonEnd(CMesh &Mesh)
Definition: DomainPrism.cpp:1339
TexGen::CDomainPrism::BuildMesh
void BuildMesh()
Get the limits for a single given repeat vector and surface mesh.
Definition: DomainPrism.cpp:67
TexGen::CDomainPrism::ClipIntersectMeshToDomain
void ClipIntersectMeshToDomain(CMesh &Mesh, bool bFillGaps=true) const
Definition: DomainPrism.cpp:192
TexGen::CDomainPrism::GetPolygonLimits
void GetPolygonLimits(XYZ &StartPoint, XYZ *SizeVecs)
Definition: DomainPrism.cpp:1302
TexGen::CDomainPrism::GeneratePlanes
void GeneratePlanes()
Generate a set of planes corresponding to the mesh elements.
Definition: DomainPrism.cpp:76
TexGen::CDomainPrism::PlaneEqual
bool PlaneEqual(PLANE Plane1, PLANE Plane2)
Test if both normal and d of two PLANE structures are equal.
Definition: DomainPrism.cpp:183
TexGen::CDomainPrism::m_Yarn
CYarn m_Yarn
Create the domain as yarn with constant polygon section and two nodes.
Definition: DomainPrism.h:124
TexGen::CDomainPrism::FillGaps
static bool FillGaps(CMesh &Mesh, const PLANE &Plane, vector< int > &Polygon, bool bMeshGaps=true)
Definition: DomainPrism.cpp:1051
TexGen::CDomainPrism::CDomainPrism
CDomainPrism(const vector< XY > &Points, XYZ &start, XYZ &end)
Constructor.
Definition: DomainPrism.cpp:26
TexGen::CDomainPrism::ClipMeshToDomain
void ClipMeshToDomain(CMesh &Mesh, bool bFillGaps=true) const
Definition: DomainPrism.cpp:697
TexGen::CDomainPrism::PopulateTiXmlElement
void PopulateTiXmlElement(TiXmlElement &Element, OUTPUT_TYPE OutputType=OUTPUT_STANDARD) const
Used for saving data to XML.
Definition: DomainPrism.cpp:59
TexGen::CDomainPrism::m_Points
vector< XY > m_Points
Prism section points.
Definition: DomainPrism.h:130
TexGen::CMesh
Defines the nodes and elements of a surface or volume mesh.
Definition: Mesh.h:58
TexGen::CMesh::AddElement
bool AddElement(ELEMENT_TYPE Type, const vector< int > &Indices)
Add an element to the mesh of given type with node number checking.
Definition: Mesh.cpp:2565
TexGen::CMesh::GetNumNodes
int GetNumNodes() const
Return the number of nodes.
Definition: Mesh.cpp:2646
TexGen::CMesh::MergeNodes
int MergeNodes(const double Tolerance=1e-8)
If any nodes share the same coordinates merge the nodes together and adjust indices accordingly.
Definition: Mesh.cpp:382
TexGen::CMesh::InsertMesh
void InsertMesh(const CMesh &Mesh, XYZ Offset=XYZ(0, 0, 0))
Add the contents of Mesh to this mesh.
Definition: Mesh.cpp:93
TexGen::CMesh::RemoveDegenerateTriangles
void RemoveDegenerateTriangles()
Remove triangles which have two equal corner indices.
Definition: Mesh.cpp:214
TexGen::CMesh::AddNode
const int AddNode(XYZ Node)
Append a node to the list of nodes, the integer returns the index of the node
Definition: Mesh.cpp:2624
TexGen::CMesh::GetNode
const XYZ & GetNode(int iIndex) const
Get the node with given ID.
Definition: Mesh.cpp:2636
TexGen::CMesh::RemoveUnreferencedNodes
int RemoveUnreferencedNodes()
Remove nodes that are not referenced by any elements.
Definition: Mesh.cpp:687
TexGen::CMesh::ConvertQuadtoTriangles
list< int >::iterator ConvertQuadtoTriangles(list< int >::iterator itQuad)
Convert a specific quad element to two triangles.
Definition: Mesh.cpp:1044
TexGen::CMesh::SaveToVTK
bool SaveToVTK(string Filename, const vector< CMeshDataBase * > *pMeshData=NULL) const
Save the mesh to VTK unstructured grid file format (.vtu)
Definition: Mesh.cpp:2348
TexGen::CMesh::GetNumNodes
static int GetNumNodes(ELEMENT_TYPE Type)
Get the number of nodes a particular element type contains.
Definition: Mesh.h:81
TexGen::CMesh::RemoveElementType
void RemoveElementType(ELEMENT_TYPE Type)
Remove elements of given type.
Definition: Mesh.cpp:723
TexGen::CMesh::GetIndices
const list< int > & GetIndices(ELEMENT_TYPE ElemType) const
Get the element indices of a given element type.
Definition: Mesh.cpp:2671
TexGen::CMesh::PYRAMID
@ PYRAMID
Definition: Mesh.h:70
TexGen::CMesh::QUADRATIC_TET
@ QUADRATIC_TET
Definition: Mesh.h:75
TexGen::CMesh::POLYGON
@ POLYGON
Definition: Mesh.h:76
TexGen::CMesh::MeshClosedLoop
void MeshClosedLoop(const XYZ &Normal, const vector< int > &ClosedLoopVector, bool bQuality=false)
Definition: Mesh.cpp:1200
TexGen::CMesh::Clear
void Clear()
Empty mesh nodes and indices.
Definition: Mesh.cpp:1448
TexGen::CNode
Represents a point on the centreline of a yarn.
Definition: Node.h:28
TexGen::CNode::GetPosition
XYZ GetPosition() const
Definition: Node.h:57
TexGen::CSectionPolygon
Creates a polygonal section, where a list of points are given to form the closed polygon.
Definition: SectionPolygon.h:36
TexGen::CYarn
Represents a yarn consisting of master nodes, section and interpolation function.
Definition: Yarn.h:49
TexGen::CYarn::AddSurfaceToMesh
bool AddSurfaceToMesh(CMesh &Mesh, bool bAddEndCaps=true) const
Create surface mesh for this yarn and add it to the surface mesh object.
Definition: Yarn.cpp:860
TexGen::CYarn::PopulateTiXmlElement
void PopulateTiXmlElement(TiXmlElement &Element, OUTPUT_TYPE OutputType)
Used for saving data to XML.
Definition: Yarn.cpp:106
TexGen::CYarn::GetNode
const CNode * GetNode(int iIndex) const
Get a master node by index.
Definition: Yarn.cpp:299
TexGen::CYarn::AssignSection
void AssignSection(const CYarnSection &YarnSection)
Assign a section to the yarn.
Definition: Yarn.cpp:649
TexGen::CYarn::GetSlaveNodes
const vector< CSlaveNode > & GetSlaveNodes(BUILD_TYPE Usage) const
Get the slave nodes and build them if necessary.
Definition: Yarn.cpp:1749
TexGen::CYarn::SURFACE
@ SURFACE
Definition: Yarn.h:60
TexGen::CYarn::AddNode
void AddNode(const CNode &Node)
Add a node to the end of the list of nodes (note the nodes must be ordered)
Definition: Yarn.cpp:196
TexGen::CYarn::SetResolution
void SetResolution(int iNumSlaveNodes, int iNumSectionPoints)
Set the resolution of the mesh created.
Definition: Yarn.cpp:306
TexGen::CYarn::PointInsideYarn
bool PointInsideYarn(const XYZ &Point, XYZ *pTangent=NULL, XY *pLoc=NULL, double *pVolumeFraction=NULL, double *pDistanceToSurface=NULL, double dTolerance=1e-9, XYZ *pOrientation=NULL, XYZ *pUp=NULL, bool bSurface=false) const
Determine if the given point lies within the yarn (this function doesn't take the repeats into accoun...
Definition: Yarn.cpp:1336
TexGen::CYarnSectionConstant
Creates a section which is constant all along the yarn.
Definition: YarnSectionConstant.h:28
TexGen
Namespace containing a series of customised math operations not found in the standard c++ library.
Definition: AdjustMeshInterference.h:27
TexGen::GetMinMaxXY
void GetMinMaxXY(const std::vector< XY > &Points, XY &Min, XY &Max)
Definition: Misc.h:244
TexGen::OUTPUT_TYPE
OUTPUT_TYPE
Definition: Misc.h:105
TexGen::GetLength
double GetLength(const XYZ &Point1, const XYZ &Point2)
Get the length between two points.
Definition: mymath.h:540
TexGen::Normalise
WXYZ & Normalise(WXYZ &Quaternion)
Normalise the quaternion and return it.
Definition: mymath.h:609
TexGen::DotProduct
double DotProduct(const XYZ &left, const XYZ &right)
Get the dot product of two vectors.
Definition: mymath.h:512
TexGen::CrossProduct
XYZ CrossProduct(const XYZ &left, const XYZ &right)
Get the cross product of two vectors.
Definition: mymath.h:524
TexGen::PLANE
Struct for representing a Plane.
Definition: Plane.h:25
TexGen::PLANE::Normal
XYZ Normal
Definition: Plane.h:30
TexGen::PLANE::d
double d
Definition: Plane.h:31
TexGen::XY
Struct for representing points in 2D space.
Definition: mymath.h:103
TexGen::XY::x
double x
Definition: mymath.h:104
TexGen::XY::y
double y
Definition: mymath.h:104
TexGen::XYZ
Struct for representing points in 3D space.
Definition: mymath.h:56
TexGen::XYZ::x
double x
Definition: mymath.h:57
TexGen::XYZ::y
double y
Definition: mymath.h:57
Generated by doxygen 1.9.2