TexGen
Mesh.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 "Mesh.h"
22#include "Plane.h"
23#include "MeshOctreeClasses.h"
24#include "Textile.h"
25#include "TexGen.h"
26
27using namespace TexGen;
28extern "C"
29{
30#include "../Triangle/triangle.h"
31}
32
33CMesh::CMesh(void)
34{
35}
36
37CMesh::~CMesh(void)
38{
39}
40
41CMesh::CMesh(TiXmlElement &Element)
42{
43 FOR_EACH_TIXMLELEMENT(pNode, Element, "Node")
44 {
45 m_Nodes.push_back(valueify<XYZ>(pNode->Attribute("value")));
46 }
47 FOR_EACH_TIXMLELEMENT(pType, Element, "Element")
48 {
49 ELEMENT_TYPE iType = (ELEMENT_TYPE)valueify<int>(pType->Attribute("type"));
50 FOR_EACH_TIXMLELEMENT(pIndex, *pType, "Index")
51 {
52 m_Indices[iType].push_back(valueify<int>(pIndex->Attribute("value")));
53 }
54 }
55}
56
57void CMesh::PopulateTiXmlElement(TiXmlElement &Element, OUTPUT_TYPE OutputType) const
58{
59 vector<XYZ>::const_iterator itNode;
60 for (itNode = m_Nodes.begin(); itNode != m_Nodes.end(); ++itNode)
61 {
62 TiXmlElement Node("Node");
63 Node.SetAttribute("value", stringify(*itNode));
64 Element.InsertEndChild(Node);
65 }
66 int i;
67 list<int>::const_iterator itIndex;
68 for (i=0; i<NUM_ELEMENT_TYPES; ++i)
69 {
70 TiXmlElement Indices("Element");
71 Indices.SetAttribute("type", i);
72 for (itIndex = m_Indices[i].begin(); itIndex != m_Indices[i].end(); ++itIndex)
73 {
74 TiXmlElement Index("Index");
75 Index.SetAttribute("value", *itIndex);
76 Indices.InsertEndChild(Index);
77 }
78 Element.InsertEndChild(Indices);
79 }
80}
81
82int CMesh::InsertNodes(const CMesh &Mesh, XYZ Offset)
83{
84 int iOffset = (int)m_Nodes.size();
85 vector<XYZ>::const_iterator itNode;
86 for (itNode = Mesh.m_Nodes.begin(); itNode != Mesh.m_Nodes.end(); ++itNode)
87 {
88 m_Nodes.push_back(*itNode+Offset);
89 }
90 return iOffset;
91}
92
93void CMesh::InsertMesh(const CMesh &Mesh, XYZ Offset)
94{
95 int iOffset = InsertNodes(Mesh, Offset);
96 int i;
97 list<int>::const_iterator itIndex;
98 for (i=0; i<NUM_ELEMENT_TYPES; ++i)
99 {
100 for (itIndex = Mesh.m_Indices[i].begin(); itIndex != Mesh.m_Indices[i].end(); ++itIndex)
101 {
102 m_Indices[i].push_back(*itIndex + iOffset);
103 }
104 }
105}
106
107void CMesh::ChangeNodeIndices(int iChangeTo, int iChangeFrom)
108{
109 int i;
110 list<int>::iterator itIndex;
111 for (i=0; i<NUM_ELEMENT_TYPES; ++i)
112 {
113 for (itIndex = m_Indices[i].begin(); itIndex != m_Indices[i].end(); ++itIndex)
114 {
115 if (*itIndex == iChangeFrom)
116 {
117 *itIndex = iChangeTo;
118 }
119 }
120 }
121}
122
123void CMesh::ChangeNodeIndices(int iChangeTo, int iChangeFrom, vector<vector<int*> > &References)
124{
125 vector<int*>::iterator itIntPointer;
126 vector<int*> &ChangeFromPointers = References[iChangeFrom];
127 vector<int*> &ChangeToPointers = References[iChangeTo];
128 for (itIntPointer = ChangeFromPointers.begin(); itIntPointer != ChangeFromPointers.end(); ++itIntPointer)
129 {
130 **itIntPointer = iChangeTo;
131 ChangeToPointers.push_back(*itIntPointer);
132 }
133 ChangeFromPointers.clear();
134}
135
136void CMesh::RemoveOpposingTriangles()
137{
138 list<int> &TriangleIndices = m_Indices[TRI];
139 int i1[3];
140 int i2[3];
141 list<int>::iterator itIter1, itIter2;
142 list<int>::iterator itTriStart1, itTriStart2;
143// int iCommonNodes;
144// int i, j;
145 for (itIter1 = TriangleIndices.begin(); itIter1 != TriangleIndices.end(); )
146 {
147 itTriStart1 = itIter1;
148 i1[0] = *(itIter1++);
149 i1[1] = *(itIter1++);
150 i1[2] = *(itIter1++);
151 for (itIter2 = itIter1; itIter2 != TriangleIndices.end(); )
152 {
153 itTriStart2 = itIter2;
154 i2[0] = *(itIter2++);
155 i2[1] = *(itIter2++);
156 i2[2] = *(itIter2++);
157
158 if (i1[0] == i2[0] && i1[1] == i2[2] && i1[2] == i2[1] ||
159 i1[0] == i2[1] && i1[1] == i2[0] && i1[2] == i2[2] ||
160 i1[0] == i2[2] && i1[1] == i2[1] && i1[2] == i2[0])
161 {
162 TriangleIndices.erase(itTriStart1, itIter1);
163 if (itIter1 == itTriStart2)
164 itIter1 = TriangleIndices.erase(itTriStart2, itIter2);
165 else
166 TriangleIndices.erase(itTriStart2, itIter2);
167 //continue;
168 break;
169 }
170 }
171 }
172}
173
174void CMesh::RemoveOpposingQuads()
175{
176 list<int> &QuadIndices = m_Indices[QUAD];
177 int i1[4];
178 int i2[4];
179 list<int>::iterator itIter1, itIter2;
180 list<int>::iterator itQuadStart1, itQuadStart2;
181// int iCommonNodes;
182// int i, j;
183 for (itIter1 = QuadIndices.begin(); itIter1 != QuadIndices.end(); )
184 {
185 itQuadStart1 = itIter1;
186 i1[0] = *(itIter1++);
187 i1[1] = *(itIter1++);
188 i1[2] = *(itIter1++);
189 i1[3] = *(itIter1++);
190 for (itIter2 = itIter1; itIter2 != QuadIndices.end(); )
191 {
192 itQuadStart2 = itIter2;
193 i2[0] = *(itIter2++);
194 i2[1] = *(itIter2++);
195 i2[2] = *(itIter2++);
196 i2[3] = *(itIter2++);
197
198 if (i1[0] == i2[0] && i1[1] == i2[3] && i1[2] == i2[2] && i1[3] == i2[1] ||
199 i1[0] == i2[1] && i1[1] == i2[0] && i1[2] == i2[3] && i1[3] == i2[2] ||
200 i1[0] == i2[2] && i1[1] == i2[1] && i1[2] == i2[0] && i1[3] == i2[3] ||
201 i1[0] == i2[3] && i1[1] == i2[2] && i1[2] == i2[1] && i1[3] == i2[0])
202 {
203 QuadIndices.erase(itQuadStart1, itIter1);
204 if (itIter1 == itQuadStart2)
205 itIter1 = QuadIndices.erase(itQuadStart2, itIter2);
206 else
207 QuadIndices.erase(itQuadStart2, itIter2);
208 continue;
209 }
210 }
211 }
212}
213
214void CMesh::RemoveDegenerateTriangles()
215{
216 list<int> &TriangleIndices = m_Indices[TRI];
217 list<int>::iterator itIter, itTriStart;
218 int i1;
219 int i2;
220 int i3;
221 for (itIter = TriangleIndices.begin(); itIter != TriangleIndices.end(); )
222 {
223 itTriStart = itIter;
224 i1 = *(itIter++);
225 i2 = *(itIter++);
226 i3 = *(itIter++);
227 if (i1 == i2 || i2 == i3 || i3 == i1)
228 {
229 itIter = TriangleIndices.erase(itTriStart, itIter);
230 }
231 }
232}
233
234void CMesh::RemoveDuplicateElements(CMesh::ELEMENT_TYPE ElementType)
235{
236 list<int> &ElementIndices = m_Indices[ElementType];
237 vector<int> Index1;
238 vector<int> Index2;
239 list<int>::iterator itIter1, itIter2;
240 list<int>::iterator itStart1, itStart2;
241
242 for (itIter1 = ElementIndices.begin(); itIter1 != ElementIndices.end(); )
243 {
244 itStart1 = itIter1;
245 Index1.clear();
246 for ( int i=0; i < GetNumNodes(ElementType); ++i )
247 {
248 Index1.push_back( *(itIter1++));
249 }
250 sort( Index1.begin(), Index1.end() );
251 for (itIter2 = itIter1; itIter2 != ElementIndices.end(); )
252 {
253 itStart2 = itIter2;
254 Index2.clear();
255 for ( int j=0; j < GetNumNodes(ElementType); ++j )
256 {
257 Index2.push_back( *(itIter2++));
258 }
259 sort( Index2.begin(), Index2.end() );
260
261 if ( equal( Index1.begin(), Index1.end(), Index2.begin() ) )
262 {
263 if (itIter1 == itStart2)
264 {
265 itIter1 = ElementIndices.erase(itStart2, itIter2);
266 itIter2 = itIter1;
267 }
268 else
269 itIter2 = ElementIndices.erase(itStart2, itIter2);
270 //continue;
271 }
272 }
273 }
274}
275
276void CMesh::RemoveDuplicateTriangles()
277{
278 list<int> &TriangleIndices = m_Indices[TRI];
279 int i1[3];
280 int i2[3];
281 list<int>::iterator itIter1, itIter2;
282 list<int>::iterator itTriStart1, itTriStart2;
283// int iCommonNodes;
284// int i, j;
285 for (itIter1 = TriangleIndices.begin(); itIter1 != TriangleIndices.end(); )
286 {
287 itTriStart1 = itIter1;
288 i1[0] = *(itIter1++);
289 i1[1] = *(itIter1++);
290 i1[2] = *(itIter1++);
291 for (itIter2 = itIter1; itIter2 != TriangleIndices.end(); )
292 {
293 itTriStart2 = itIter2;
294 i2[0] = *(itIter2++);
295 i2[1] = *(itIter2++);
296 i2[2] = *(itIter2++);
297
298 if ( i1[0] == i2[0] && i1[1] == i2[1] && i1[2] == i2[2] )
299 {
300 //TriangleIndices.erase(itTriStart1, itIter1);
301 if (itIter1 == itTriStart2)
302 {
303 itIter1 = TriangleIndices.erase(itTriStart2, itIter2);
304 itIter2 = itIter1;
305 }
306 else
307 itIter2 = TriangleIndices.erase(itTriStart2, itIter2);
308 //continue;
309 }
310 }
311 }
312}
313
314void CMesh::RemoveDuplicateSegments()
315{
316 list<int> &LineIndices = m_Indices[LINE];
317 list<int>::iterator itIter, itLineStart;
318 list<int>::iterator itCompare;
319 int i1, i2;
320 int j1, j2;
321 for (itIter = LineIndices.begin(); itIter != LineIndices.end(); )
322 {
323 itLineStart = itIter;
324 i1 = *(itIter++);
325 i2 = *(itIter++);
326 for (itCompare = itIter; itCompare != LineIndices.end(); )
327 {
328 j1 = *(itCompare++);
329 j2 = *(itCompare++);
330 if ((i1 == j1 && i2 == j2) ||
331 (i1 == j2 && i2 == j1))
332 {
333 itIter = LineIndices.erase(itLineStart, itIter);
334 break;
335 }
336 }
337 }
338}
339
340pair<XYZ, XYZ> CMesh::GetAABB(double dGrowDistance) const
341{
342 pair<XYZ, XYZ> AABB;
343
344 vector<XYZ>::const_iterator itNode;
345 int iNode;
346 for (itNode = m_Nodes.begin(), iNode=0; itNode != m_Nodes.end(); ++itNode, ++iNode)
347 {
348 if (iNode == 0)
349 {
350 AABB.first = AABB.second = *itNode;
351 }
352 else
353 {
354 AABB.first = Min(AABB.first, *itNode);
355 AABB.second = Max(AABB.second, *itNode);
356/* if (itNode->x < AABB.first.x)
357 AABB.first.x = itNode->x;
358 if (itNode->y < AABB.first.y)
359 AABB.first.y = itNode->y;
360 if (itNode->z < AABB.first.z)
361 AABB.first.z = itNode->z;
362 if (itNode->x > AABB.second.x)
363 AABB.second.x = itNode->x;
364 if (itNode->y > AABB.second.y)
365 AABB.second.y = itNode->y;
366 if (itNode->z > AABB.second.z)
367 AABB.second.z = itNode->z;*/
368 }
369 }
370
371 AABB.first.x -= dGrowDistance;
372 AABB.first.y -= dGrowDistance;
373 AABB.first.z -= dGrowDistance;
374 AABB.second.x += dGrowDistance;
375 AABB.second.y += dGrowDistance;
376 AABB.second.z += dGrowDistance;
377
378 return AABB;
379}
380
381
382int CMesh::MergeNodes(const double TOL)
383{
384 // Merging of nodes is optimised using an octree, that is the space is partitioned to group
385 // nodes that are in proximity to each other. Which reduces the number of checks needed to be performed
386 // when comparing node positions.
387
388 // Get the dimensions of the octree
389
390 pair<XYZ, XYZ> AABB = GetAABB(TOL);
391
392 double dSize = AABB.second.x - AABB.first.x;
393 if (dSize < AABB.second.y - AABB.first.y)
394 dSize = AABB.second.y - AABB.first.y;
395 if (dSize < AABB.second.z - AABB.first.z)
396 dSize = AABB.second.z - AABB.first.z;
397
398 // Prevent problems with creating an octree too small
399 dSize = max(dSize, 1e-3);
400
401 // Create the octree
402
403 XYZ Avg = 0.5*(AABB.second + AABB.first);
404 dSize *= 1.1;
405 Octree<pair<int, XYZ> > Octree(Vector3f(float(Avg.x-0.5*dSize), float(Avg.y-0.5*dSize), float(Avg.z-0.5*dSize)), (float)dSize, 10, 10);
406 COctreeAgentNode Agent(TOL);
407
408 // Add the nodes to the octree, the indices of the nodes need to be stored with the nodes because
409 // the elements refer to the nodes by their index.
410
411 vector<pair<int, XYZ> > NodesWithIndices;
412 NodesWithIndices.resize(m_Nodes.size());
413
414 vector<XYZ>::iterator itNode;
415 int iNode;
416 for (itNode = m_Nodes.begin(), iNode=0; itNode != m_Nodes.end(); ++itNode, ++iNode)
417 {
418 NodesWithIndices[iNode] = pair<int, XYZ>(iNode, *itNode);
419 Octree.insertItem(NodesWithIndices[iNode], Agent);
420 }
421
422 // Visit the octree with the visitor class, this actually does the node merging
423
424 COctreeVisitorMergeNodes Visitor(*this, TOL);
425
426 Octree.visit(Visitor);
427
428 return Visitor.GetNumMerged();
429
430
431 // UNOPTIMISED VERSION (worth keeping for debuging purposes)
432
433 /* //vector<bool> DeletedNodes;
434 set<int> DeletedNodes;
435 //DeletedNodes.resize(m_Nodes.size(), false);
436
437 map<ELEMENT_TYPE, list<int> >::iterator itType;
438 vector<XYZ>::iterator itNode1;
439 vector<XYZ>::iterator itNode2;
440 int iNode1, iNode2;
441 int iNumMerged = 0;
442 for (itNode1 = m_Nodes.begin(), iNode1=0; itNode1 != m_Nodes.end(); ++itNode1, ++iNode1)
443 {
444 for (itNode2 = itNode1+1, iNode2 = iNode1+1; itNode2 != m_Nodes.end(); ++itNode2, ++iNode2)
445 {
446 if (GetLengthSquared(*itNode2, *itNode1) < TOL*TOL)
447 {
448 ++iNumMerged;
449 ChangeNodeIndices(iNode1, iNode2);
450 //DeletedNodes[iNode2] = true;
451 DeletedNodes.insert(iNode2);
452 }
453 }
454 }
455
456 DeleteNodes(DeletedNodes);
457
458 cout << "Num Merged Nodes: " << iNumMerged << endl;
459
460 return iNumMerged;*/
461}
462
463vector<pair<int, int> > CMesh::GetNodePairs(XYZ Vector, const double TOL) const
464{
465 vector<pair<int, int> > NodePairs;
466 vector<XYZ>::const_iterator itNode1;
467 vector<XYZ>::const_iterator itNode2;
468 int iNode1, iNode2;
469 double dTolSquared = TOL*TOL;
470 for (itNode1 = m_Nodes.begin(), iNode1=0; itNode1 != m_Nodes.end(); ++itNode1, ++iNode1)
471 {
472 for (itNode2 = m_Nodes.begin(), iNode2=0; itNode2 != m_Nodes.end(); ++itNode2, ++iNode2)
473 {
474 if (GetLengthSquared(*itNode1, *itNode2-Vector) < dTolSquared)
475 {
476 NodePairs.push_back(make_pair(iNode1, iNode2));
477 break;
478 }
479 }
480 }
481 return NodePairs;
482}
483
484 void CMesh::GetNodePairs(XYZ Vector, vector<pair<int, int> > &NodePairs, const double TOL ) const
485{
486 vector<XYZ>::const_iterator itNode1;
487 vector<XYZ>::const_iterator itNode2;
488
489 vector<XYZ> OffsetNodes;
490 for ( itNode2 = m_Nodes.begin(); itNode2 != m_Nodes.end(); ++itNode2 )
491 {
492 OffsetNodes.push_back( *itNode2 - Vector );
493 }
494
495 bool bProgress = m_Nodes.size() > 10000;
496 int iProgressInc = m_Nodes.size()/10;
497
498 int iNode1, iNode2;
499 double dTolSquared = TOL*TOL;
500 bool bFound;
501 for (itNode1 = m_Nodes.begin(), iNode1=0; itNode1 != m_Nodes.end(); ++itNode1, ++iNode1)
502 {
503 if ( bProgress && !(iNode1 % iProgressInc) )
504 {
505 TGLOG( "GetNodePairs " << iNode1/iProgressInc*10 << "% complete");
506 }
507 bFound = false;
508 for ( iNode2 = iNode1; iNode2 < OffsetNodes.size(); ++ iNode2 ) // Matching node normally later in mesh than itNode1 so start checking here
509 { // to improve speed
510 if (GetLengthSquared(*itNode1, OffsetNodes[iNode2]) < dTolSquared)
511 {
512 NodePairs.push_back(make_pair(iNode1, iNode2));
513 bFound = true;
514 break;
515 }
516 }
517 if ( !bFound )
518 {
519 for ( iNode2 = 0; iNode2 < iNode1; ++iNode2 ) // Check from beginning if node earlier in mesh than one checking against
520 {
521 if (GetLengthSquared(*itNode1, OffsetNodes[iNode2]) < dTolSquared)
522 {
523 NodePairs.push_back(make_pair(iNode1, iNode2));
524 break;
525 }
526 }
527 }
528 }
529 //return NodePairs;
530}
531
532int CMesh::GetClosestNode(XYZ Position) const
533{
534 vector<XYZ>::const_iterator itNode;
535 double dClosestDistSqrd = -1;
536 double dDistSqrd;
537 int iClosestNode = -1, i;
538 for (itNode = m_Nodes.begin(), i=0; itNode != m_Nodes.end(); ++itNode, ++i)
539 {
540 dDistSqrd = GetLengthSquared(Position, *itNode);
541 if ( dDistSqrd < dClosestDistSqrd || dClosestDistSqrd == -1 )
542 {
543 dClosestDistSqrd = dDistSqrd;
544 iClosestNode = i;
545 }
546 }
547 return iClosestNode;
548}
549
550int CMesh::GetClosestNodeDistance(XYZ Position, double dTol ) const
551{
552 vector<XYZ>::const_iterator itNode;
553 double dClosestDistSqrd = -1;
554 double dDistSqrd;
555 int iClosestNode = -1, i;
556 for (itNode = m_Nodes.begin(), i=0; itNode != m_Nodes.end(); ++itNode, ++i)
557 {
558 dDistSqrd = GetLengthSquared(Position, *itNode);
559 if ( dDistSqrd < dClosestDistSqrd || dClosestDistSqrd == -1 )
560 {
561 dClosestDistSqrd = dDistSqrd;
562 if ( dClosestDistSqrd <= dTol )
563 iClosestNode = i;
564 }
565 }
566 return iClosestNode;
567}
568
569void CMesh::ConvertToSurfaceMesh()
570{
571 Convert3Dto2D();
572 RemoveOpposingTriangles();
573 RemoveOpposingQuads();
574}
575
576void CMesh::Convert3Dto2D()
577{
578 ConvertHextoQuad();
579 ConvertWedgeto2D();
580 ConvertTettoTriangle();
581 ConvertPyramidto2D();
582}
583
584void CMesh::ConvertToTetMesh()
585{
586 ConvertWedgetoTet();
587 ConvertHextoTet();
588 ConvertPyramidtoTet();
589 RemoveAllElementsExcept(TET);
590}
591
592void CMesh::ConvertToTriangleMesh()
593{
594 Convert3Dto2D();
595 ConvertQuadstoTriangles();
596 RemoveAllElementsExcept(TRI);
597}
598
599void CMesh::ConvertToSegmentMesh()
600{
601 ConvertToTriangleMesh();
602 ConvertTrianglestoSegments();
603 RemoveDuplicateSegments();
604}
605
606void CMesh::ConvertTriToQuad( double Tolerance )
607{
608 list<int> &TriangleIndices = m_Indices[TRI];
609 vector<int> i1;
610 int i2[3];
611 list<int>::iterator itIter1, itIter2;
612 list<int>::iterator itTriStart1, itTriStart2;
613 vector<int>::iterator iti1;
614
615 int iNumNodes = GetNumNodes(TRI);
616 for (itIter1 = TriangleIndices.begin(); itIter1 != TriangleIndices.end(); )
617 {
618 itTriStart1 = itIter1;
619 i1.clear();
620 i1.push_back(*(itIter1++));
621 i1.push_back(*(itIter1++));
622 i1.push_back(*(itIter1++));
623 for ( itIter2 = itIter1; itIter2 != TriangleIndices.end(); )
624 {
625 itTriStart2 = itIter2;
626 list<int> RemInd;
627 i2[0] = *(itIter2++);
628 RemInd.push_back(i2[0]);
629 i2[1] = *(itIter2++);
630 RemInd.push_back(i2[1]);
631 i2[2] = *(itIter2++);
632 RemInd.push_back(i2[2]);
633 vector<int> CommonInd;
634
635 int i;
636 for ( iti1 = i1.begin(), i = 0; iti1 != i1.end(); iti1++, ++i )
637 {
638
639 for ( int j = 0; j < iNumNodes; ++j )
640 {
641 if ( *iti1 == i2[j] )
642 {
643 CommonInd.push_back(i);
644 RemInd.remove(i2[j]);
645 break;
646 }
647 }
648 }
649 if ( CommonInd.size() == 2 )
650 {
651 // Check normals
652 XYZ N1 = m_Nodes[i1[0]];
653 XYZ N2 = m_Nodes[i1[1]];
654 XYZ N3 = m_Nodes[i1[2]];
655 XYZ Normal1 = CrossProduct(N2-N1, N3-N1);
656 Normal1 /= GetLength(Normal1); // Normalise
657
658 N1 = m_Nodes[i2[0]];
659 N2 = m_Nodes[i2[1]];
660 N3 = m_Nodes[i2[2]];
661 XYZ Normal2 = CrossProduct(N2-N1, N3-N1);
662 Normal2 /= GetLength(Normal2); // Normalise
663
664 if ( Normal1 == Normal2 || Normal1 == -Normal2 )
665 {
666 if ( CommonInd[0] == 0 && CommonInd[1] == 2 )
667 i1.push_back( *(RemInd.begin()) );
668 else
669 i1.insert( i1.begin()+CommonInd[1],*(RemInd.begin()));
670 AddElement( QUAD, i1 );
671 if ( itIter1 == itTriStart2 )
672 itIter1 = TriangleIndices.erase( itTriStart1, itIter2 );
673 else
674 {
675 TriangleIndices.erase( itTriStart2, itIter2 );
676 itIter1 = TriangleIndices.erase( itTriStart1, itIter1 );
677 }
678
679 break;
680 }
681
682 }
683 }
684 }
685}
686
687int CMesh::RemoveUnreferencedNodes()
688{
689/* vector<bool> UnreferencedNodes;
690 UnreferencedNodes.resize(m_Nodes.size(), true);*/
691 set<int> UnreferencedNodes;
692 int i;
693 for (i=0; i<(int)m_Nodes.size(); ++i)
694 {
695 UnreferencedNodes.insert(i);
696 }
697
698 list<int>::iterator itIndex;
699 for (i = 0; i < NUM_ELEMENT_TYPES; ++i)
700 {
701 for (itIndex = m_Indices[i].begin(); itIndex != m_Indices[i].end(); ++itIndex)
702 {
703// UnreferencedNodes[*itIndex] = false;
704 UnreferencedNodes.erase(*itIndex); // Get rid of the nodes which are referenced..
705 }
706 }
707
708 return DeleteNodes(UnreferencedNodes); // and delete the unused ones remaining
709}
710
711void CMesh::RemoveAllElementsExcept(ELEMENT_TYPE Type)
712{
713 int i;
714 for (i = 0; i < NUM_ELEMENT_TYPES; ++i)
715 {
716 if (i == Type)
717 continue;
718 else
719 m_Indices[i].clear();
720 }
721}
722
723void CMesh::RemoveElementType( ELEMENT_TYPE Type )
724{
725 m_Indices[Type].clear();
726}
727
728int CMesh::DeleteNodes(const set<int> &Nodes)
729{
730 vector<int> NewIndices;
731 NewIndices.resize(m_Nodes.size(), 0);
732 int i, j;
733 int iNumNodesDeleted = 0;
734 for (i=0, j=0; i<(int)NewIndices.size(); ++i)
735 {
736 NewIndices[i] = j;
737 if (Nodes.count(i))
738 ++iNumNodesDeleted;
739 else
740 ++j;
741 }
742
743 list<int>::iterator itIndex;
744 for (i = 0; i < NUM_ELEMENT_TYPES; ++i)
745 {
746 for (itIndex = m_Indices[i].begin(); itIndex != m_Indices[i].end(); ++itIndex)
747 {
748 *itIndex = NewIndices[*itIndex];
749 }
750 }
751
752 vector<XYZ>::iterator itNode;
753 vector<XYZ> NewNodes;
754 int iNodeIndex;
755 for (itNode = m_Nodes.begin(), iNodeIndex=0; itNode != m_Nodes.end(); ++itNode, ++iNodeIndex)
756 {
757 if (!Nodes.count(iNodeIndex))
758 NewNodes.push_back(*itNode);
759 }
760 m_Nodes = NewNodes;
761
762 return iNumNodesDeleted;
763}
764
765void CMesh::ConvertHextoQuad()
766{
767 list<int>::iterator itIter;
768 int i1, i2, i3, i4, i5, i6, i7, i8;
769 for (itIter = m_Indices[HEX].begin(); itIter != m_Indices[HEX].end(); )
770 {
771 i1 = *(itIter++);
772 i2 = *(itIter++);
773 i3 = *(itIter++);
774 i4 = *(itIter++);
775 i5 = *(itIter++);
776 i6 = *(itIter++);
777 i7 = *(itIter++);
778 i8 = *(itIter++);
779
780 m_Indices[QUAD].push_back(i4);
781 m_Indices[QUAD].push_back(i3);
782 m_Indices[QUAD].push_back(i2);
783 m_Indices[QUAD].push_back(i1);
784
785 m_Indices[QUAD].push_back(i1);
786 m_Indices[QUAD].push_back(i2);
787 m_Indices[QUAD].push_back(i6);
788 m_Indices[QUAD].push_back(i5);
789
790 m_Indices[QUAD].push_back(i2);
791 m_Indices[QUAD].push_back(i3);
792 m_Indices[QUAD].push_back(i7);
793 m_Indices[QUAD].push_back(i6);
794
795 m_Indices[QUAD].push_back(i3);
796 m_Indices[QUAD].push_back(i4);
797 m_Indices[QUAD].push_back(i8);
798 m_Indices[QUAD].push_back(i7);
799
800 m_Indices[QUAD].push_back(i4);
801 m_Indices[QUAD].push_back(i1);
802 m_Indices[QUAD].push_back(i5);
803 m_Indices[QUAD].push_back(i8);
804
805 m_Indices[QUAD].push_back(i5);
806 m_Indices[QUAD].push_back(i6);
807 m_Indices[QUAD].push_back(i7);
808 m_Indices[QUAD].push_back(i8);
809 }
810 m_Indices[HEX].clear();
811}
812
813void CMesh::ConvertWedgeto2D()
814{
815 list<int>::iterator itIter;
816 int i1, i2, i3, i4, i5, i6;
817 for (itIter = m_Indices[WEDGE].begin(); itIter != m_Indices[WEDGE].end(); )
818 {
819 i1 = *(itIter++);
820 i2 = *(itIter++);
821 i3 = *(itIter++);
822 i4 = *(itIter++);
823 i5 = *(itIter++);
824 i6 = *(itIter++);
825
826 m_Indices[TRI].push_back(i1);
827 m_Indices[TRI].push_back(i2);
828 m_Indices[TRI].push_back(i3);
829
830 m_Indices[QUAD].push_back(i2);
831 m_Indices[QUAD].push_back(i1);
832 m_Indices[QUAD].push_back(i4);
833 m_Indices[QUAD].push_back(i5);
834
835 m_Indices[QUAD].push_back(i3);
836 m_Indices[QUAD].push_back(i2);
837 m_Indices[QUAD].push_back(i5);
838 m_Indices[QUAD].push_back(i6);
839
840 m_Indices[QUAD].push_back(i1);
841 m_Indices[QUAD].push_back(i3);
842 m_Indices[QUAD].push_back(i6);
843 m_Indices[QUAD].push_back(i4);
844
845 m_Indices[TRI].push_back(i6);
846 m_Indices[TRI].push_back(i5);
847 m_Indices[TRI].push_back(i4);
848 }
849 m_Indices[WEDGE].clear();
850}
851
852void CMesh::ConvertTettoTriangle()
853{
854 list<int>::iterator itIter;
855 int i1, i2, i3, i4;
856 for (itIter = m_Indices[TET].begin(); itIter != m_Indices[TET].end(); )
857 {
858 i1 = *(itIter++);
859 i2 = *(itIter++);
860 i3 = *(itIter++);
861 i4 = *(itIter++);
862
863 m_Indices[TRI].push_back(i3);
864 m_Indices[TRI].push_back(i2);
865 m_Indices[TRI].push_back(i1);
866
867 m_Indices[TRI].push_back(i1);
868 m_Indices[TRI].push_back(i2);
869 m_Indices[TRI].push_back(i4);
870
871 m_Indices[TRI].push_back(i2);
872 m_Indices[TRI].push_back(i3);
873 m_Indices[TRI].push_back(i4);
874
875 m_Indices[TRI].push_back(i3);
876 m_Indices[TRI].push_back(i1);
877 m_Indices[TRI].push_back(i4);
878 }
879 m_Indices[TET].clear();
880}
881
882void CMesh::ConvertHextoWedge(bool bQuality)
883{
884 list<int>::iterator itIter;
885 int i0, i1, i2, i3, i4, i5, i6, i7;
886 for (itIter = m_Indices[HEX].begin(); itIter != m_Indices[HEX].end(); )
887 {
888 i0 = *(itIter++);
889 i1 = *(itIter++);
890 i2 = *(itIter++);
891 i3 = *(itIter++);
892 i4 = *(itIter++);
893 i5 = *(itIter++);
894 i6 = *(itIter++);
895 i7 = *(itIter++);
896
897 // Split it up randomly without regard to quality
898 m_Indices[WEDGE].push_back(i4);
899 m_Indices[WEDGE].push_back(i5);
900 m_Indices[WEDGE].push_back(i6);
901 m_Indices[WEDGE].push_back(i0);
902 m_Indices[WEDGE].push_back(i1);
903 m_Indices[WEDGE].push_back(i2);
904
905 m_Indices[WEDGE].push_back(i6);
906 m_Indices[WEDGE].push_back(i7);
907 m_Indices[WEDGE].push_back(i4);
908 m_Indices[WEDGE].push_back(i2);
909 m_Indices[WEDGE].push_back(i3);
910 m_Indices[WEDGE].push_back(i0);
911
912 // Todo: If the bQuality flag is set then use a more intelligent
913 // way to decide how to split the element up
914 }
915 m_Indices[HEX].clear();
916}
917
918void CMesh::ConvertWedgetoTetandPyramid(bool bQuality)
919{
920 list<int>::iterator itIter;
921 int i0, i1, i2, i3, i4, i5;
922 for (itIter = m_Indices[WEDGE].begin(); itIter != m_Indices[WEDGE].end(); )
923 {
924 i0 = *(itIter++);
925 i1 = *(itIter++);
926 i2 = *(itIter++);
927 i3 = *(itIter++);
928 i4 = *(itIter++);
929 i5 = *(itIter++);
930
931 // Split it up randomly without regard to quality
932 m_Indices[PYRAMID].push_back(i0);
933 m_Indices[PYRAMID].push_back(i1);
934 m_Indices[PYRAMID].push_back(i4);
935 m_Indices[PYRAMID].push_back(i3);
936 m_Indices[PYRAMID].push_back(i2);
937
938 m_Indices[TET].push_back(i3);
939 m_Indices[TET].push_back(i4);
940 m_Indices[TET].push_back(i5);
941 m_Indices[TET].push_back(i2);
942
943 // Todo: If the bQuality flag is set then use a more intelligent
944 // way to decide how to split the element up
945 }
946 m_Indices[WEDGE].clear();
947}
948
949void CMesh::ConvertPyramidtoTet(bool bQuality)
950{
951 list<int>::iterator itIter;
952 int i0, i1, i2, i3, i4;
953 for (itIter = m_Indices[PYRAMID].begin(); itIter != m_Indices[PYRAMID].end(); )
954 {
955 i0 = *(itIter++);
956 i1 = *(itIter++);
957 i2 = *(itIter++);
958 i3 = *(itIter++);
959 i4 = *(itIter++);
960
961 double d1 = 0, d2 = 0;
962 if (bQuality)
963 {
964 // If we don't care about the quality then just skip this calculation
965 d1 = GetLengthSquared(m_Nodes[i0], m_Nodes[i2]);
966 d2 = GetLengthSquared(m_Nodes[i1], m_Nodes[i3]);
967 }
968 if (d1<d2)
969 {
970 m_Indices[TET].push_back(i0);
971 m_Indices[TET].push_back(i1);
972 m_Indices[TET].push_back(i2);
973 m_Indices[TET].push_back(i4);
974
975 m_Indices[TET].push_back(i0);
976 m_Indices[TET].push_back(i2);
977 m_Indices[TET].push_back(i3);
978 m_Indices[TET].push_back(i4);
979 }
980 else
981 {
982 m_Indices[TET].push_back(i1);
983 m_Indices[TET].push_back(i2);
984 m_Indices[TET].push_back(i3);
985 m_Indices[TET].push_back(i4);
986
987 m_Indices[TET].push_back(i3);
988 m_Indices[TET].push_back(i0);
989 m_Indices[TET].push_back(i1);
990 m_Indices[TET].push_back(i4);
991 }
992 }
993 m_Indices[PYRAMID].clear();
994}
995
996void CMesh::ConvertWedgetoTet(bool bQuality)
997{
998 ConvertWedgetoTetandPyramid(bQuality);
999 ConvertPyramidtoTet(bQuality);
1000}
1001
1002void CMesh::ConvertHextoTet(bool bQuality)
1003{
1004 ConvertHextoWedge(bQuality);
1005 ConvertWedgetoTet(bQuality);
1006}
1007
1008void CMesh::ConvertPyramidto2D()
1009{
1010 list<int>::iterator itIter;
1011 int i0, i1, i2, i3, i4;
1012 for (itIter = m_Indices[PYRAMID].begin(); itIter != m_Indices[PYRAMID].end(); )
1013 {
1014 i0 = *(itIter++);
1015 i1 = *(itIter++);
1016 i2 = *(itIter++);
1017 i3 = *(itIter++);
1018 i4 = *(itIter++);
1019
1020 m_Indices[QUAD].push_back(i3);
1021 m_Indices[QUAD].push_back(i2);
1022 m_Indices[QUAD].push_back(i1);
1023 m_Indices[QUAD].push_back(i0);
1024
1025 m_Indices[TRI].push_back(i0);
1026 m_Indices[TRI].push_back(i1);
1027 m_Indices[TRI].push_back(i4);
1028
1029 m_Indices[TRI].push_back(i1);
1030 m_Indices[TRI].push_back(i2);
1031 m_Indices[TRI].push_back(i4);
1032
1033 m_Indices[TRI].push_back(i2);
1034 m_Indices[TRI].push_back(i3);
1035 m_Indices[TRI].push_back(i4);
1036
1037 m_Indices[TRI].push_back(i3);
1038 m_Indices[TRI].push_back(i0);
1039 m_Indices[TRI].push_back(i4);
1040 }
1041 m_Indices[PYRAMID].clear();
1042}
1043
1044list<int>::iterator CMesh::ConvertQuadtoTriangles(list<int>::iterator itQuad)
1045{
1046 int i0, i1, i2, i3;
1047
1048 list<int>::iterator itEnd = itQuad;
1049
1050 i0 = *(itEnd++);
1051 i1 = *(itEnd++);
1052 i2 = *(itEnd++);
1053 i3 = *(itEnd++);
1054
1055 m_Indices[TRI].push_back(i0);
1056 m_Indices[TRI].push_back(i1);
1057 m_Indices[TRI].push_back(i2);
1058
1059 m_Indices[TRI].push_back(i0);
1060 m_Indices[TRI].push_back(i2);
1061 m_Indices[TRI].push_back(i3);
1062
1063 return m_Indices[QUAD].erase(itQuad, itEnd);
1064}
1065
1066void CMesh::ConvertTrianglestoSegments()
1067{
1068 list<int>::iterator itIter;
1069 int i0, i1, i2;
1070 for (itIter = m_Indices[TRI].begin(); itIter != m_Indices[TRI].end(); )
1071 {
1072 i0 = *(itIter++);
1073 i1 = *(itIter++);
1074 i2 = *(itIter++);
1075
1076 m_Indices[LINE].push_back(i0);
1077 m_Indices[LINE].push_back(i1);
1078
1079 m_Indices[LINE].push_back(i1);
1080 m_Indices[LINE].push_back(i2);
1081
1082 m_Indices[LINE].push_back(i2);
1083 m_Indices[LINE].push_back(i0);
1084 }
1085 m_Indices[TRI].clear();
1086}
1087
1088void CMesh::ConvertQuadstoTriangles(bool bQuality)
1089{
1090 list<int>::iterator itIter;
1091 int i0, i1, i2, i3;
1092 for (itIter = m_Indices[QUAD].begin(); itIter != m_Indices[QUAD].end(); )
1093 {
1094 i0 = *(itIter++);
1095 i1 = *(itIter++);
1096 i2 = *(itIter++);
1097 i3 = *(itIter++);
1098
1099 if (!bQuality)
1100 {
1101 m_Indices[TRI].push_back(i0);
1102 m_Indices[TRI].push_back(i1);
1103 m_Indices[TRI].push_back(i2);
1104
1105 m_Indices[TRI].push_back(i0);
1106 m_Indices[TRI].push_back(i2);
1107 m_Indices[TRI].push_back(i3);
1108 }
1109 else
1110 {
1111 double d1, d2;
1112 d1 = GetLengthSquared(m_Nodes[i0], m_Nodes[i2]);
1113 d2 = GetLengthSquared(m_Nodes[i1], m_Nodes[i3]);
1114 if (d1<d2)
1115 {
1116 m_Indices[TRI].push_back(i0);
1117 m_Indices[TRI].push_back(i1);
1118 m_Indices[TRI].push_back(i2);
1119
1120 m_Indices[TRI].push_back(i0);
1121 m_Indices[TRI].push_back(i2);
1122 m_Indices[TRI].push_back(i3);
1123 }
1124 else
1125 {
1126 m_Indices[TRI].push_back(i1);
1127 m_Indices[TRI].push_back(i2);
1128 m_Indices[TRI].push_back(i3);
1129
1130 m_Indices[TRI].push_back(i3);
1131 m_Indices[TRI].push_back(i0);
1132 m_Indices[TRI].push_back(i1);
1133 }
1134 }
1135 }
1136 m_Indices[QUAD].clear();
1137}
1138
1139void CMesh::Translate(XYZ Vector)
1140{
1141 vector<XYZ>::iterator itNode;
1142 for (itNode = m_Nodes.begin(); itNode != m_Nodes.end(); ++itNode)
1143 {
1144 *itNode += Vector;
1145 }
1146}
1147
1148void CMesh::Rotate(WXYZ Rotation, XYZ Origin)
1149{
1150 vector<XYZ>::iterator itNode;
1151 for (itNode = m_Nodes.begin(); itNode != m_Nodes.end(); ++itNode)
1152 {
1153 (*itNode) = Rotation * (*itNode-Origin) + Origin;
1154 }
1155}
1156/*
1157void CMesh::Rotate(XYZ X, XYZ Y, XYZ Z)
1158{
1159 vector<XYZ>::iterator itNode;
1160 for (itNode = m_Nodes.begin(); itNode != m_Nodes.end(); ++itNode)
1161 {
1162 *itNode = X * itNode->x + Y * itNode->y + Z * itNode->z;
1163 }
1164}
1165*/
1166void CMesh::FlipNormals()
1167{
1168 list<int>::iterator itIter;
1169 list<int>::iterator it0, it1, it2, it3;
1170 int i0, i1, i2, i3;
1171 for (itIter = m_Indices[QUAD].begin(); itIter != m_Indices[QUAD].end(); )
1172 {
1173 it0 = itIter++;
1174 it1 = itIter++;
1175 it2 = itIter++;
1176 it3 = itIter++;
1177 i0 = *it0;
1178 i1 = *it1;
1179 i2 = *it2;
1180 i3 = *it3;
1181 *it0 = i3;
1182 *it1 = i2;
1183 *it2 = i1;
1184 *it3 = i0;
1185 }
1186 for (itIter = m_Indices[TRI].begin(); itIter != m_Indices[TRI].end(); )
1187 {
1188 it0 = itIter++;
1189 it1 = itIter++;
1190 it2 = itIter++;
1191 i0 = *it0;
1192 i1 = *it1;
1193 i2 = *it2;
1194 *it0 = i2;
1195 *it1 = i1;
1196 *it2 = i0;
1197 }
1198}
1199
1200void CMesh::MeshClosedLoop(const XYZ &Normal, const vector<int> &ClosedLoopVector, bool bQuality)
1201{
1202 // A more efficient algorithm is described in http://www.cs.umd.edu/~mount/754/Lects/754lects.pdf
1203 // worth implementing if this function needs optimising. Although will need extending to 3d.
1204 const double TOL = 1e-10;
1205
1206 list<int> ClosedLoop(ClosedLoopVector.begin(), ClosedLoopVector.end());
1207
1208 list<int> &TriIndices = m_Indices[CMesh::TRI];
1209
1210 list<int>::iterator itPrev;
1211 list<int>::iterator itCurrent;
1212 list<int>::iterator itNext;
1213 list<int>::iterator itCompare;
1214 list<int>::iterator itBest;
1215
1216 XYZ V1, V2, V3;
1217 XYZ P1, P2, P3, P;
1218
1219 PLANE PL1, PL2, PL3;
1220 double l1, l2;
1221
1222 double dBestAngle;
1223 double dAngle;
1224
1225 bool bFound;
1226
1227 // Keep going while there are at least 3 indices in the loop
1228 while (ClosedLoop.size() >= 3)
1229 {
1230 // Three iterators are used to keep track of 3 consecutive points
1231 bFound = false;
1232 itPrev = ClosedLoop.end(); --itPrev;
1233 itCurrent = ClosedLoop.begin();
1234 itNext = ClosedLoop.begin(); ++itNext;
1235 // Go round increment the three points until all combination of triangles have been looked at
1236 while (itCurrent != ClosedLoop.end())
1237 {
1238 // Get the point coordinates
1239 P1 = m_Nodes[*itPrev];
1240 P2 = m_Nodes[*itCurrent];
1241 P3 = m_Nodes[*itNext];
1242
1243 // Form vectors from the edges of the triangle
1244 V1 = P2 - P1;
1245 V2 = P3 - P2;
1246 V3 = P1 - P3;
1247
1248 // Check that the triangle is facing the right way, this is important
1249 // for concave loops
1250 if (DotProduct(Normal, CrossProduct(V1, V3)) > -TOL)
1251 {
1252 PL1.Normal = CrossProduct(V1, Normal);
1253 Normalise(PL1.Normal);
1254 PL1.d = DotProduct(PL1.Normal, P1);
1255
1256 PL2.Normal = CrossProduct(V2, Normal);
1257 Normalise(PL2.Normal);
1258 PL2.d = DotProduct(PL2.Normal, P2);
1259
1260 PL3.Normal = CrossProduct(V3, Normal);
1261 Normalise(PL3.Normal);
1262 PL3.d = DotProduct(PL3.Normal, P3);
1263
1264 // Make sure that none of the other points lie within the triangle, again important
1265 // for concave loops
1266 for (itCompare = ClosedLoop.begin(); itCompare != ClosedLoop.end(); ++itCompare)
1267 {
1268 if (itCompare == itPrev || itCompare == itCurrent || itCompare == itNext)
1269 continue;
1270
1271 // Get the coordinate of the point
1272 P = m_Nodes[*itCompare];
1273
1274 // If true then the point lies within the triangle, we must not form a triangle with this
1275 // set of points
1276 if (DotProduct(P, PL1.Normal) > PL1.d && DotProduct(P, PL2.Normal) > PL2.d && DotProduct(P, PL3.Normal) > PL3.d)
1277 {
1278 break;
1279 }
1280 }
1281
1282 // If this is true then no other point lies within the proposed triangle
1283 if (itCompare == ClosedLoop.end())
1284 {
1285 l1 = GetLength(V1);
1286 l2 = GetLength(V2);
1287
1288 // Find the angle between V1 and V2
1289 dAngle = acos(DotProduct(V1, V2)/(l1*l2));
1290
1291 // If this angle is larger than the rest of the valid triangles
1292 // then consider this as our best canditate for creating a triangle
1293 if (!bFound || dAngle > dBestAngle)
1294 {
1295 dBestAngle = dAngle;
1296 itBest = itCurrent;
1297 bFound = true;
1298 if (!bQuality)
1299 break;
1300 }
1301 }
1302 }
1303
1304 // Go to the next set of 3 points
1305 ++itPrev;
1306 ++itCurrent;
1307 ++itNext;
1308 if (itPrev == ClosedLoop.end())
1309 itPrev = ClosedLoop.begin();
1310 if (itNext == ClosedLoop.end())
1311 itNext = ClosedLoop.begin();
1312 }
1313
1314 // Didn't find a best triangle, in which case just choose one at random.
1315 // this should never happen however.
1316 if (!bFound) // Does happen if slice through yarn at edge of domain
1317 { // but don't want the whole program to keel over so taken assert out!
1318 // Just doesn't render that end section
1319 TGERROR( "Couldn't mesh closed loop on domain edge" );
1320 return;
1321 //itBest = ClosedLoop.begin();
1322 //assert(false);
1323 }
1324
1325 // The best 3 points will be used to form a triangle, the middle point will
1326 // then be removed and the remaining polygon will continue to be meshed.
1327
1328 if (itBest == ClosedLoop.begin())
1329 {itPrev = ClosedLoop.end(); --itPrev;}
1330 else
1331 {itPrev = itBest; --itPrev;}
1332 itCurrent = itBest;
1333 itNext = itBest; ++itNext;
1334 if (itNext == ClosedLoop.end())
1335 itNext = ClosedLoop.begin();
1336 TriIndices.push_back(*itPrev);
1337 TriIndices.push_back(*itCurrent);
1338 TriIndices.push_back(*itNext);
1339 ClosedLoop.erase(itCurrent);
1340 }
1341}
1342
1343
1344/*void CMesh::MeshClosedLoop(const XYZ &Normal, const vector<int> &ClosedLoop)
1345{
1346 // NOTE: Triangle is not very robust, if anything fucks up, the program is likely to just crash
1347 // thus the code has been replaced with my own meshing algorithms. This is kept here in case
1348 // it needs to be re-implemented.
1349
1350 // Triangle will be used to triangulate the gaps
1351
1352 if (ClosedLoop.size() < 3) // Triangle will quit if no vertices are given to it
1353 return;
1354
1355 bool bReverse = false;
1356
1357 char szSwitches[] = "pzQPB";
1358
1359 triangulateio TriangleInput;
1360 triangulateio TriangleOutput;
1361 memset(&TriangleInput, 0, sizeof(TriangleInput));
1362 memset(&TriangleOutput, 0, sizeof(TriangleOutput));
1363
1364 TriangleInput.pointlist = new REAL [ClosedLoop.size()*2];
1365 TriangleInput.numberofpoints = (int)ClosedLoop.size();
1366
1367 // Ignore one of the coordinates, the one with the highest value for the plane normal
1368 // This is a very basic projection of the 3D coordinates to a 2d plane. One of the three planes
1369 // XY, YZ or XZ is chosen. The plane is chosen to give minimum distortion. There will still be some
1370 // distortion, but this is not an issue since a quality mesh is not needed.
1371 int iIgnoreCoordinate;
1372 if (abs(Normal.x) >= abs(Normal.y) && abs(Normal.x) >= abs(Normal.z))
1373 iIgnoreCoordinate = 0;
1374 else if (abs(Normal.y) >= abs(Normal.z))
1375 iIgnoreCoordinate = 1;
1376 else
1377 iIgnoreCoordinate = 2;
1378
1379 if (Normal[iIgnoreCoordinate] > 0)
1380 bReverse = true;
1381
1382 // Fill the nodes into the TriangleInput structure
1383 int i;
1384 for (i=0; i<(int)ClosedLoop.size(); ++i)
1385 {
1386 // Ignore one of the coordinates, this is decided further up
1387 TriangleInput.pointlist[i*2] = m_Nodes[ClosedLoop[i]][(iIgnoreCoordinate+1)%3];
1388 TriangleInput.pointlist[i*2+1] = m_Nodes[ClosedLoop[i]][(iIgnoreCoordinate+2)%3];
1389 }
1390
1391 TriangleInput.segmentlist = new int [ClosedLoop.size()*2];
1392 TriangleInput.numberofsegments = (int)ClosedLoop.size();
1393
1394 // Pass in the segments with the new node indices
1395 for (i=0; i<(int)ClosedLoop.size(); ++i)
1396 {
1397 TriangleInput.segmentlist[i*2] = i;
1398 TriangleInput.segmentlist[i*2+1] = (i+1)%ClosedLoop.size();
1399 }
1400
1401 // Triangulate
1402 triangulate(szSwitches, &TriangleInput, &TriangleOutput, NULL);
1403
1404 delete [] TriangleInput.pointlist;
1405 delete [] TriangleInput.segmentlist;
1406
1407 list<int> &TriIndices = m_Indices[CMesh::TRI];
1408
1409 int i1, i2, i3;
1410 for (i=0; i<TriangleOutput.numberoftriangles; ++i)
1411 {
1412 i1 = TriangleOutput.trianglelist[i*3];
1413 i2 = TriangleOutput.trianglelist[i*3+1];
1414 i3 = TriangleOutput.trianglelist[i*3+2];
1415 // Convert the node from those passed to triangle to the Mesh indices
1416 if (!bReverse)
1417 {
1418 TriIndices.push_back(ClosedLoop[i1]);
1419 TriIndices.push_back(ClosedLoop[i2]);
1420 TriIndices.push_back(ClosedLoop[i3]);
1421 }
1422 else
1423 {
1424 TriIndices.push_back(ClosedLoop[i1]);
1425 TriIndices.push_back(ClosedLoop[i3]);
1426 TriIndices.push_back(ClosedLoop[i2]);
1427 }
1428 }
1429
1430 // Free up the memory allocated by triangle
1431 trifree(TriangleOutput.pointlist);
1432 trifree(TriangleOutput.trianglelist);
1433
1434 return;
1435}*/
1436
1437void CMesh::CopySelfToRange(XYZ Vector, int iLowerLimit, int iUpperLimit)
1438{
1439 CMesh Original = *this;
1440 Clear();
1441 int i;
1442 for (i=iLowerLimit; i<=iUpperLimit; ++i)
1443 {
1444 InsertMesh(Original, i*Vector);
1445 }
1446}
1447
1448void CMesh::Clear()
1449{
1450 m_Nodes.clear();
1451 int i;
1452 for (i = 0; i < NUM_ELEMENT_TYPES; ++i)
1453 {
1454 m_Indices[i].clear();
1455 }
1456/* map<ELEMENT_TYPE, list<int> >::iterator itType;
1457 for (itType = m_Indices.begin(); itType != m_Indices.end(); ++itType)
1458 {
1459 itType->second.clear();
1460 }*/
1461}
1462
1463int CMesh::OutputNodes(ostream &Output, int iStartIndex, string Seperator, bool bSCIRun) const
1464{
1465 int iNodeIndex;
1466 vector<XYZ>::const_iterator itNode;
1467
1468 if ( bSCIRun )
1469 Output << m_Nodes.size() << endl;
1470
1471 for (itNode = m_Nodes.begin(), iNodeIndex=0; itNode != m_Nodes.end(); ++itNode, ++iNodeIndex)
1472 {
1473 if (iStartIndex >= 0 && !bSCIRun)
1474 Output << iNodeIndex+iStartIndex << Seperator;
1475 Output << *itNode << endl;
1476 }
1477 return iNodeIndex+iStartIndex;
1478}
1479
1480int CMesh::OutputElements(ostream &Output, CMesh::ELEMENT_TYPE ElementType, int iStartIndex, int iIndexOffset, string Seperator, bool bSCIRun) const
1481{
1482 assert(ElementType >= 0 && ElementType < NUM_ELEMENT_TYPES);
1483 int i;
1484 int iElementNumber = 0;
1485 list<int>::const_iterator itIndex;
1486
1487 if ( bSCIRun )
1488 Output << m_Indices[ElementType].size()/GetNumNodes(ElementType) << endl;
1489 for (itIndex = m_Indices[ElementType].begin(); itIndex != m_Indices[ElementType].end(); ++iElementNumber)
1490 {
1491 if (iStartIndex >= 0 && !bSCIRun)
1492 Output << iElementNumber+iStartIndex << Seperator;
1493 for (i=0; i<GetNumNodes(ElementType); ++i, ++itIndex)
1494 {
1495 if (i>0)
1496 Output << Seperator;
1497 Output << (*itIndex)+iIndexOffset;
1498 }
1499 Output << endl;
1500 }
1501 return iElementNumber+iStartIndex;
1502}
1503
1504void CMesh::GetNodeElementReferences(vector<vector<int*> > &References)
1505{
1506 References.clear();
1507 References.resize(m_Nodes.size());
1508 int i;
1509 list<int>::iterator itIndex;
1510 for (i = 0; i < NUM_ELEMENT_TYPES; ++i)
1511 {
1512 for (itIndex = m_Indices[i].begin(); itIndex != m_Indices[i].end(); ++itIndex)
1513 {
1514 References[*itIndex].push_back(&(*itIndex));
1515 }
1516 }
1517}
1518
1519double CMesh::CalculateVolume() const
1520{
1521 int i;
1522 for (i = 0; i < NUM_ELEMENT_TYPES; ++i)
1523 {
1524 if (i != TRI && !m_Indices[i].empty())
1525 {
1526 CMesh TriMesh(*this);
1527 TriMesh.ConvertToTriangleMesh(); // Converts whole mesh to triangle not just current element type
1528 return TriMesh.CalculateVolume(); // Which is why can return after CalculateVolume
1529 }
1530 }
1531
1532 double dVolume = 0;
1533 // Code from http://www.gamedev.net/reference/articles/article2247.asp
1534 const list<int> &TriangleIndices = m_Indices[TRI];
1535 XYZ P1, P2, P3;
1536 list<int>::const_iterator itIter;
1537 for (itIter = TriangleIndices.begin(); itIter != TriangleIndices.end(); )
1538 {
1539 P1 = m_Nodes[*(itIter++)];
1540 P2 = m_Nodes[*(itIter++)];
1541 P3 = m_Nodes[*(itIter++)];
1542
1543 dVolume += ((P2.y-P1.y)*(P3.z-P1.z)-(P2.z-P1.z)*(P3.y-P1.y)) * (P1.x+P2.x+P3.x);
1544 }
1545
1546 return dVolume / 6;
1547}
1548
1549vector<XYZ> CMesh::GetElementCenters() const
1550{
1551 TGLOG("Get Element Centres");
1552 vector<XYZ> ElementCenters;
1553 int i, j, iNumNodes;
1554 list<int>::const_iterator itIndex;
1555 for (i = 0; i < NUM_ELEMENT_TYPES; ++i)
1556 {
1557 iNumNodes = GetNumNodes((ELEMENT_TYPE)i);
1558 if ( iNumNodes != -1 )
1559 {
1560 for (itIndex = m_Indices[i].begin(); itIndex != m_Indices[i].end(); )
1561 {
1562 XYZ Center;
1563 for (j=0; j<iNumNodes; ++j)
1564 {
1565 Center += m_Nodes[*itIndex];
1566 ++itIndex;
1567 }
1568 Center /= iNumNodes;
1569 ElementCenters.push_back(Center);
1570 }
1571 }
1572 }
1573 return ElementCenters;
1574}
1575
1576vector<XYZ> CMesh::GetElementCenters( ELEMENT_TYPE type )
1577{
1578 vector<XYZ> ElementCentres;
1579 int iNumNodes = GetNumNodes(type);
1580 list<int>::const_iterator itIndex;
1581 for ( itIndex = m_Indices[type].begin(); itIndex != m_Indices[type].end(); )
1582 {
1583 XYZ Centre;
1584 for (int i=0; i<iNumNodes; ++i)
1585 {
1586 Centre += m_Nodes[*itIndex];
1587 ++itIndex;
1588 }
1589 Centre /= iNumNodes;
1590 ElementCentres.push_back(Centre);
1591 }
1592 return ElementCentres;
1593}
1594
1595int CMesh::CountInvertedElements() const
1596{
1597 int i;
1598 int iNumInvertedElements = 0;
1599 for (i = 0; i < NUM_ELEMENT_TYPES; ++i)
1600 {
1601 iNumInvertedElements += CountInvertedElements((ELEMENT_TYPE)i);
1602 }
1603 return iNumInvertedElements;
1604}
1605
1606int CMesh::CountInvertedElements(ELEMENT_TYPE ElementType) const
1607{
1608 if (ElementType != TET && ElementType != PYRAMID && ElementType != WEDGE)
1609 {
1610 return 0;
1611 }
1612 list<int>::const_iterator itIter;
1613 int i, iNumInvertedElements = 0;
1614 CMesh CopiedMesh;
1615 CopiedMesh.m_Nodes = m_Nodes;
1616 const list<int> &Indices = m_Indices[ElementType];
1617 int iNumNodes = GetNumNodes(ElementType);
1618 for (itIter = Indices.begin(); itIter != Indices.end(); )
1619 {
1620 CopiedMesh.m_Indices[ElementType].clear();
1621 for (i=0; i<iNumNodes; ++i)
1622 {
1623 CopiedMesh.m_Indices[ElementType].push_back(*(itIter++));
1624 }
1625 if (CopiedMesh.CalculateVolume() < 0)
1626 {
1627 ++iNumInvertedElements;
1628 }
1629 }
1630 return iNumInvertedElements;
1631}
1632
1633int CMesh::IntersectLine(const XYZ &P1, const XYZ &P2, vector<pair<double, XYZ> > &IntersectionPoints, pair<bool, bool> TrimResults, bool bForceFind/*, const CElementsOctree *pOctree*/) const
1634{
1635 int i;
1636 for (i = 0; i < NUM_ELEMENT_TYPES; ++i)
1637 {
1638 if (i != TRI && !m_Indices[i].empty())
1639 {
1640 TGERROR("Warning: IntersectLine only works for triangle mesh, mesh contains other element types");
1641 break;
1642 }
1643 }
1644
1645 IntersectionPoints.clear();
1646
1647 XYZ T1, T2, T3;
1648 XYZ Normal, Intersection;
1649 double dU;
1650 const double dTolerance = 1e-9;
1651/* if (pOctree && pOctree->GetOctree())
1652 {
1653 list<vector<int> > Elements;
1654 list<vector<int> >::iterator itIter;
1655
1656 COctreeVisitorElementNearLine Visitor(P1, P2, Elements, TrimResults);
1657
1658 pOctree->GetOctree()->visit(Visitor);
1659
1660 for (itIter = Elements.begin(); itIter != Elements.end(); ++itIter)
1661 {
1662 T1 = m_Nodes[(*itIter)[0]];
1663 T2 = m_Nodes[(*itIter)[1]];
1664 T3 = m_Nodes[(*itIter)[2]];
1665
1666 Normal = CrossProduct(T2-T1, T3-T1);
1667 Normalise(Normal);
1668
1669 if (GetIntersectionLinePlane(P1, P2, T1, Normal, Intersection, &dU))
1670 {
1671 if (TrimResults.first && dU < 0)
1672 ; // Do nothing
1673 else if(TrimResults.second && dU > 1)
1674 ; // Do nothing
1675 else if (PointInsideTriangle(T1, T2, T3, Intersection, Normal))
1676 {
1677 IntersectionPoints.push_back(pair<double, XYZ>(dU, Normal));
1678 }
1679 }
1680 }
1681 }
1682 else
1683 {*/
1684 bool bFirst = true;
1685 double dMin;
1686 double dAccuracy;
1687 double dBestU;
1688 XYZ BestNormal;
1689
1690 const list<int> &TriangleIndices = m_Indices[TRI];
1691 list<int>::const_iterator itIter;
1692 for (itIter = TriangleIndices.begin(); itIter != TriangleIndices.end(); )
1693 {
1694 T1 = m_Nodes[*(itIter++)];
1695 T2 = m_Nodes[*(itIter++)];
1696 T3 = m_Nodes[*(itIter++)];
1697
1698 Normal = CrossProduct(T2-T1, T3-T1);
1699 if (Normal)
1700 {
1701 Normalise(Normal);
1702
1703 if (GetIntersectionLinePlane(P1, P2, T1, Normal, Intersection, &dU))
1704 {
1705 if (TrimResults.first && dU < 0)
1706 continue; // Do nothing
1707 if(TrimResults.second && dU > 1)
1708 continue; // Do nothing
1709 dMin = PointInsideTriangleAccuracy(T1, T2, T3, Intersection, Normal);
1710 if (dMin >= -dTolerance)
1711 {
1712 IntersectionPoints.push_back(pair<double, XYZ>(dU, Normal));
1713 }
1714 if (bFirst || dMin > dAccuracy)
1715 {
1716 dAccuracy = dMin;
1717 bFirst = false;
1718 dBestU = dU;
1719 BestNormal = Normal;
1720 }
1721/* if (PointInsideTriangle(T1, T2, T3, Intersection, Normal))
1722 {
1723 IntersectionPoints.push_back(pair<double, XYZ>(dU, Normal));
1724 }*/
1725 }
1726 }
1727 }
1728// }
1729
1730 if (bForceFind && IntersectionPoints.empty())
1731 {
1732 IntersectionPoints.push_back(pair<double, XYZ>(dBestU, BestNormal));
1733 }
1734
1735 sort(IntersectionPoints.begin(), IntersectionPoints.end(), LessPairDoubleXYZ());
1736
1737 return IntersectionPoints.size();
1738}
1739
1740void CMesh::AddOrCancel(list<pair<int, int> > &EdgeStack, pair<int, int> Edge)
1741{
1742 list<pair<int, int> >::iterator FindResult;
1743 FindResult = find(EdgeStack.begin(), EdgeStack.end(), Edge);
1744 if (FindResult != EdgeStack.end())
1745 {
1746 EdgeStack.erase(FindResult);
1747 return;
1748 }
1749 FindResult = find(EdgeStack.begin(), EdgeStack.end(), pair<int, int>(Edge.second, Edge.first));
1750 if (FindResult != EdgeStack.end())
1751 {
1752 EdgeStack.erase(FindResult);
1753 return;
1754 }
1755 EdgeStack.push_back(Edge);
1756}
1757
1758void CMesh::MeshConvexHull()
1759{
1760 // Algorithm and code taken from: (the incremental method is implemented since it seems to be
1761 // the easiest to implement and quite robust, I tried to implement the gift wrap method however
1762 // it doesn't seem as robust, problems occur when points are collinear and/or coplanar)
1763 // http://www.cse.unsw.edu.au/~lambert/java/3d/hull.html
1764 // http://www.cse.unsw.edu.au/~lambert/java/3d/incremental.html
1765 // http://www.cse.unsw.edu.au/~lambert/java/3d/source/Incremental.java
1766 const double TOL = 1e-7;
1767
1768 list<pair<int, int> > EdgeStack;
1769 list<pair<int, int> >::iterator itEdge;
1770
1771 MergeNodes();
1772
1773 int i;
1774 for (i = 0; i < NUM_ELEMENT_TYPES; ++i)
1775 {
1776 m_Indices[i].clear();
1777 }
1778
1779 if (m_Nodes.size() < 3)
1780 return;
1781
1782 list<int> &TriangleIndices = m_Indices[TRI];
1783 // Generate a list of triangle normals
1784 list<PLANE> TrianglePlanes;
1785 PLANE Plane;
1786 XYZ P1, P2, P3, P;
1787
1788 // Find the first set of 3 points that are not collinear to create a triangle
1789 for (i=0; i+2<(int)m_Nodes.size(); ++i)
1790 {
1791 P1 = m_Nodes[i+0];
1792 P2 = m_Nodes[i+1];
1793 P3 = m_Nodes[i+2];
1794 Plane.Normal = CrossProduct(P2 - P1, P3 - P1);
1795 if (GetLengthSquared(Plane.Normal) > TOL*TOL)
1796 {
1797 break;
1798 }
1799 }
1800 if (i+2 == (int)m_Nodes.size())
1801 {
1802 TGERROR("Unable to create convex hull, all points are collinear");
1803 assert(false);
1804 return;
1805 }
1806
1807 TriangleIndices.push_back(i+0);
1808 TriangleIndices.push_back(i+1);
1809 TriangleIndices.push_back(i+2);
1810 Normalise(Plane.Normal);
1811 Plane.d = DotProduct(Plane.Normal, P1);
1812 TrianglePlanes.push_back(Plane);
1813
1814 TriangleIndices.push_back(i+0);
1815 TriangleIndices.push_back(i+2);
1816 TriangleIndices.push_back(i+1);
1817 Plane.Normal = -Plane.Normal;
1818 Plane.d = -Plane.d;
1819 TrianglePlanes.push_back(Plane);
1820
1821 list<int>::iterator itIter, itTriStart;
1822 list<PLANE>::iterator itPlane;
1823 int i1, i2, i3;
1824 bool bFacesAdded;
1825 // This loop will be iterated until no more faces are added, this is necessary when
1826 // we have more than 3 nodes that are coplanar. If not some nodes will be missed.
1827 do
1828 {
1829 bFacesAdded = false;
1830 for (i=0; i<(int)m_Nodes.size(); ++i)
1831 {
1832 P = m_Nodes[i];
1833 // Delete faces that this vertex can see
1834 for (itIter = TriangleIndices.begin(), itPlane = TrianglePlanes.begin(); itIter != TriangleIndices.end(); )
1835 {
1836 itTriStart = itIter;
1837 i1 = *(itIter++);
1838 i2 = *(itIter++);
1839 i3 = *(itIter++);
1840 // If the vertex can see the plane (with a tolerance) we don't want to remove
1841 // triangles that are in the same place as the vertex or the algorithm will fail.
1842 if (DotProduct(itPlane->Normal, P) > itPlane->d+TOL)
1843 {
1844 // If the edge already exist in the edge stack then it should cancel with it
1845 // (remove the edge rather than adding it) otherwise add the edge.
1846 AddOrCancel(EdgeStack, pair<int, int>(i1, i2));
1847 AddOrCancel(EdgeStack, pair<int, int>(i2, i3));
1848 AddOrCancel(EdgeStack, pair<int, int>(i3, i1));
1849 // Delete the triangle
1850 itIter = TriangleIndices.erase(itTriStart, itIter);
1851 itPlane = TrianglePlanes.erase(itPlane);
1852 }
1853 else
1854 ++itPlane;
1855 }
1856 // Create new triangles and calculate the planes of the new triangle.
1857 // Not only is it an optimisation to calculate the triangle planes only once
1858 // it is also for consitency.
1859 for (itEdge = EdgeStack.begin(); itEdge != EdgeStack.end(); ++itEdge)
1860 {
1861 i1 = itEdge->first;
1862 i2 = itEdge->second;
1863 i3 = i;
1864 TriangleIndices.push_back(i1);
1865 TriangleIndices.push_back(i2);
1866 TriangleIndices.push_back(i3);
1867 P1 = m_Nodes[i1];
1868 P2 = m_Nodes[i2];
1869 P3 = m_Nodes[i3];
1870 Plane.Normal = CrossProduct(P2 - P1, P3 - P1);
1871 Normalise(Plane.Normal);
1872 Plane.d = DotProduct(Plane.Normal, P1);
1873 TrianglePlanes.push_back(Plane);
1874 bFacesAdded = true;
1875 }
1876 EdgeStack.clear();
1877 }
1878 } while (bFacesAdded);
1879}
1880
1881void CMesh::BuildGrid(XYZ Min, XYZ Max, int iNumX, int iNumY, int iNumZ)
1882{
1883 XYZ P;
1884 int i, j, k;
1885 for (i=0; i<iNumX; ++i)
1886 {
1887 for (j=0; j<iNumY; ++j)
1888 {
1889 for (k=0; k<iNumZ; ++k)
1890 {
1891 P.x = i/double(iNumX-1);
1892 P.y = j/double(iNumY-1);
1893 P.z = k/double(iNumZ-1);
1894 P = Min + (Max-Min)*P;
1895 m_Nodes.push_back(P);
1896 // Create hex elements out of the grid
1897 if (i < iNumX-1 && j < iNumY-1 && k < iNumZ-1)
1898 {
1899 m_Indices[HEX].push_back((k+0) + (j+0)*iNumZ + (i+0)*iNumZ*iNumY);
1900 m_Indices[HEX].push_back((k+0) + (j+0)*iNumZ + (i+1)*iNumZ*iNumY);
1901 m_Indices[HEX].push_back((k+0) + (j+1)*iNumZ + (i+1)*iNumZ*iNumY);
1902 m_Indices[HEX].push_back((k+0) + (j+1)*iNumZ + (i+0)*iNumZ*iNumY);
1903 m_Indices[HEX].push_back((k+1) + (j+0)*iNumZ + (i+0)*iNumZ*iNumY);
1904 m_Indices[HEX].push_back((k+1) + (j+0)*iNumZ + (i+1)*iNumZ*iNumY);
1905 m_Indices[HEX].push_back((k+1) + (j+1)*iNumZ + (i+1)*iNumZ*iNumY);
1906 m_Indices[HEX].push_back((k+1) + (j+1)*iNumZ + (i+0)*iNumZ*iNumY);
1907 }
1908 }
1909 }
1910 }
1911}
1912
1913void CMesh::BuildGrid(XYZ Min, XYZ Max, double dPointsPerUnit)
1914{
1915 int iNumX, iNumY, iNumZ;
1916 iNumX = Round((Max.x-Min.x)*dPointsPerUnit);
1917 iNumY = Round((Max.y-Min.y)*dPointsPerUnit);
1918 iNumZ = Round((Max.z-Min.z)*dPointsPerUnit);
1919 BuildGrid(Min, Max, iNumX, iNumY, iNumZ);
1920}
1921
1922void CMesh::WriteBinaryXYZ(ostream &Output, XYZ Vector)
1923{
1924 float val;
1925 int i;
1926 for (i=0; i<3; ++i)
1927 {
1928 val = (float)Vector[i];
1929 Output.write((char*)&val, 4);
1930 }
1931}
1932
1933bool CMesh::SaveToSTL(string Filename, bool bBinary) const
1934{
1935 AddExtensionIfMissing(Filename, ".stl");
1936
1937 int i;
1938 for (i = 0; i < NUM_ELEMENT_TYPES; ++i)
1939 {
1940 if (i != TRI && !m_Indices[i].empty())
1941 {
1942 CMesh TriMesh(*this);
1943 TriMesh.ConvertToTriangleMesh();
1944 return TriMesh.SaveToSTL(Filename, bBinary);
1945 }
1946 }
1947
1948 ofstream Output;
1949 if (bBinary)
1950 Output.open(Filename.c_str(), ios::out|ios::binary);
1951 else
1952 Output.open(Filename.c_str(), ios::out);
1953 if (!Output)
1954 return false;
1955
1956 const list<int> &TriangleIndices = m_Indices[TRI];
1957
1958 if (bBinary)
1959 {
1960 char szHeader[80];
1961 int iNumFacets = TriangleIndices.size()/3;
1962 strncpy(szHeader, Filename.c_str(), 80);
1963 Output.write(szHeader, 80);
1964 Output.write((char*)&iNumFacets, 4);
1965 }
1966 else
1967 Output << "solid " << Filename << endl;
1968
1969 XYZ T1, T2, T3, Normal;
1970 list<int>::const_iterator itIndex;
1971 short int Padding = 0;
1972 for (itIndex = TriangleIndices.begin(); itIndex != TriangleIndices.end(); )
1973 {
1974 T1 = m_Nodes[*(itIndex++)];
1975 T2 = m_Nodes[*(itIndex++)];
1976 T3 = m_Nodes[*(itIndex++)];
1977
1978 Normal = CrossProduct(T2-T1, T3-T1);
1979 Normalise(Normal);
1980
1981 if (bBinary)
1982 {
1983 WriteBinaryXYZ(Output, Normal);
1984 WriteBinaryXYZ(Output, T1);
1985 WriteBinaryXYZ(Output, T2);
1986 WriteBinaryXYZ(Output, T3);
1987 Output.write((char*)&Padding, 2);
1988 }
1989 else
1990 {
1991 Output << " facet normal " << Normal.x << " " << Normal.y << " " << Normal.z << endl;
1992 Output << " outer loop" << endl;
1993 Output << " vertex " << T1.x << " " << T1.y << " " << T1.z << endl;
1994 Output << " vertex " << T2.x << " " << T2.y << " " << T2.z << endl;
1995 Output << " vertex " << T3.x << " " << T3.y << " " << T3.z << endl;
1996 Output << " endloop" << endl;
1997 Output << " endfacet" << endl;
1998 }
1999 }
2000
2001 if (!bBinary)
2002 Output << "endsolid " << Filename << endl;
2003
2004 Output.close();
2005
2006 return true;
2007}
2008
2009bool CMesh::SaveToSMESH(string Filename) const
2010{
2011 AddExtensionIfMissing(Filename, ".smesh");
2012
2013 ofstream Output(Filename.c_str(), ios::out);
2014 if (!Output)
2015 return false;
2016 Output << "# node count, 3 dim, no attribute, no boundary marker" << endl;
2017 Output << m_Nodes.size() << " 3 0 0" << endl;
2018
2019 Output << "# node index, node coordinates" << endl;
2020 int iNodeIndex;
2021 vector<XYZ>::const_iterator itNode;
2022 for (itNode = m_Nodes.begin(), iNodeIndex=0; itNode != m_Nodes.end(); ++itNode, ++iNodeIndex)
2023 {
2024 Output << iNodeIndex << " " << itNode->x << " " << itNode->y << " " << itNode->z << endl;
2025 }
2026
2027 int iNumTriangles = m_Indices[TRI].size()/3;
2028 int iNumQuads = m_Indices[QUAD].size()/4;
2029
2030 Output << "# facet count, no boundary marker" << endl;
2031 Output << iNumTriangles+iNumQuads << " 0" << endl;
2032 Output << "# facets" << endl;
2033 list<int>::const_iterator itIndex;
2034 int iNumNodesPerElement;
2035 int i, j;
2036 ELEMENT_TYPE ElemType;
2037 for (j=0; j<2; ++j)
2038 {
2039 if (j==0)
2040 ElemType = QUAD;
2041 else
2042 ElemType = TRI;
2043 iNumNodesPerElement = GetNumNodes(ElemType);
2044 for (itIndex = m_Indices[ElemType].begin(); itIndex != m_Indices[ElemType].end(); )
2045 {
2046 Output << iNumNodesPerElement << " ";
2047 for (i=0; i<iNumNodesPerElement; ++i, ++itIndex)
2048 {
2049 if (i>0)
2050 Output << " ";
2051 Output << (*itIndex);
2052 }
2053 Output << endl;
2054 }
2055 }
2056 Output << "# Part 3 - hole list" << endl;
2057 Output << "0 # no hole" << endl;
2058 Output << "# Part 4 - region list" << endl;
2059 Output << "0 # no region" << endl;
2060
2061 return true;
2062}
2063/*
2064bool CMesh::ReadFromTetGen(string NodesFile, string ElementsFile, string FacesFile)
2065{
2066 // Development paused, implemented in python more easily
2067 ifstream Nodes(NodesFile.c_str());
2068 ifstream Elements(ElementsFile.c_str());
2069 ifstream Faces(FacesFile.c_str());
2070
2071 if (Nodes)
2072 {
2073 int iNumNodes;
2074 }
2075 else
2076 return false; // Not much point in reading the other files if we don't have nodal information
2077 if (Elements)
2078 {
2079
2080 }
2081 if (Faces)
2082 {
2083
2084 }
2085 return true;
2086}
2087*/
2088
2089bool CMesh::SaveToABAQUS(string Filename, const vector<POINT_INFO> *pElementInfo, bool bCreateStep, bool bCreateMaterial, int iElementType)
2090{
2091 AddExtensionIfMissing(Filename, ".inp");
2092
2093 ofstream Output(Filename.c_str());
2094
2095 if (!Output)
2096 {
2097 TGERROR("Unable to output mesh to ABAQUS file format, could not open file: " << Filename);
2098 return false;
2099 }
2100
2101 TGLOG("Saving mesh data to " << Filename);
2102
2103 Output << "*Heading" << endl;
2104 Output << "File generated by TexGen v" << TEXGEN.GetVersion() << endl;
2105
2106 Output << "************" << endl;
2107 Output << "*** MESH ***" << endl;
2108 Output << "************" << endl;
2109 Output << "*Node" << endl;
2110 OutputNodes(Output, 1);
2111
2112 int iStartIndex = 1;
2113
2114 // Linear elements
2115 if (!m_Indices[TET].empty())
2116 {
2117 Output << "*Element, Type=C3D4" << endl;
2118 iStartIndex = OutputElements(Output, TET, iStartIndex, 1);
2119 }
2120
2121 if (!m_Indices[WEDGE].empty())
2122 {
2123 Output << "*Element, Type=C3D6" << endl;
2124 // NOTE: The element ordering is different between ABAQUS and VTK...
2125 // and since the index ordering in this class is based on VTK we need
2126 // to reorder the elements :(. So we won't use the OutputElements function
2127 // in this case
2128// iStartIndex = OutputElements(Output, WEDGE, iStartIndex, 1);
2129 int i1, i2, i3, i4, i5, i6;
2130 int iElementNumber = 0;
2131 list<int>::const_iterator itIndex;
2132 for (itIndex = m_Indices[WEDGE].begin(); itIndex != m_Indices[WEDGE].end(); ++iElementNumber)
2133 {
2134 i1 = *(itIndex++)+1;
2135 i2 = *(itIndex++)+1;
2136 i3 = *(itIndex++)+1;
2137 i4 = *(itIndex++)+1;
2138 i5 = *(itIndex++)+1;
2139 i6 = *(itIndex++)+1;
2140 Output << iElementNumber+iStartIndex
2141 << ", " << i4 << ", " << i5 << ", " << i6
2142 << ", " << i1 << ", " << i2 << ", " << i3
2143 << endl;
2144 }
2145 iStartIndex += iElementNumber;
2146 }
2147
2148 if (!m_Indices[HEX].empty())
2149 {
2150 if ( !iElementType )
2151 {
2152 Output << "*Element, Type=C3D8R" << endl;
2153 }
2154 else
2155 {
2156 Output << "*Element, Type=C3D8" << endl;
2157 }
2158 iStartIndex = OutputElements(Output, HEX, iStartIndex, 1);
2159 }
2160
2161 // Quadratic elements
2162 if (!m_Indices[QUADRATIC_TET].empty())
2163 {
2164 Output << "*Element, Type=C3D10" << endl;
2165 iStartIndex = OutputElements(Output, QUADRATIC_TET, iStartIndex, 1);
2166 }
2167
2168 // Shell elements
2169 if (!m_Indices[QUAD].empty())
2170 {
2171 Output << "*Element, Type=S4R" << endl;
2172 iStartIndex = OutputElements(Output, QUAD, iStartIndex, 1);
2173 }
2174
2175 if (pElementInfo)
2176 {
2177 string OrientationsFilename = Filename;
2178 OrientationsFilename.replace(OrientationsFilename.end()-4, OrientationsFilename.end(), ".ori");
2179 ofstream OriOutput(OrientationsFilename.c_str());
2180 string ElementDataFilename = Filename;
2181 ElementDataFilename.replace(ElementDataFilename.end()-4, ElementDataFilename.end(), ".eld");
2182 ofstream DataOutput(ElementDataFilename.c_str());
2183
2184 if (!OriOutput)
2185 {
2186 TGERROR("Unable to output orientations, could not open file: " << OrientationsFilename);
2187 return false;
2188 }
2189 if (!DataOutput)
2190 {
2191 TGERROR("Unable to output additional element data, could not open file: " << ElementDataFilename);
2192 return false;
2193 }
2194
2195 TGLOG("Saving element orientations data to " << OrientationsFilename);
2196 TGLOG("Saving additional element data to " << ElementDataFilename);
2197
2198 WriteOrientationsHeader( Output );
2199 Output << "*Distribution Table, Name=TexGenOrientationVectors" << endl;
2200 Output << "COORD3D,COORD3D" << endl;
2201 Output << "*Distribution, Location=Element, Table=TexGenOrientationVectors, Name=TexGenOrientationVectors, Input=" << StripPath(OrientationsFilename) << endl;
2202 Output << "*Orientation, Name=TexGenOrientations, Definition=coordinates" << endl;
2203 Output << "TexGenOrientationVectors" << endl;
2204 Output << "1, 0" << endl;
2205
2206 // Default orientation
2207 WriteOrientationsHeader( OriOutput );
2208 OriOutput << ", 1.0, 0.0, 0.0, 0.0, 1.0, 0.0" << endl;
2209
2210 int i;
2211
2212 WriteElementsHeader( DataOutput );
2213
2214 map<int, vector<int> > ElementSets;
2215 vector<POINT_INFO>::const_iterator itData;
2216 for (itData = pElementInfo->begin(), i=1; itData != pElementInfo->end(); ++itData, ++i)
2217 {
2218 if (itData->iYarnIndex != -1)
2219 {
2220 if ( GetLength(itData->Up ) )
2221 {
2222 XYZ Up = itData->Up;
2223 XYZ Dir = itData->Orientation;
2224
2225 XYZ Perp = CrossProduct(Dir, Up);
2226 Normalise(Perp);
2227 OriOutput << i << ", " << Dir << ", " << Perp << endl;
2228 }
2229 else
2230 {
2231 // Default orientation
2232 OriOutput << i << ", 1.0, 0.0, 0.0, 0.0, 1.0, 0.0" << endl;
2233 }
2234 }
2235 else
2236 {
2237 // Default orientation
2238 OriOutput << i << ", 1.0, 0.0, 0.0, 0.0, 1.0, 0.0" << endl;
2239 }
2240 DataOutput << i;
2241 DataOutput << ", " << itData->iYarnIndex;
2242 DataOutput << ", " << itData->Location; // This counts as 2 DepVars
2243 DataOutput << ", " << itData->dVolumeFraction;
2244 DataOutput << ", " << itData->dSurfaceDistance;
2245 DataOutput << endl;
2246 ElementSets[itData->iYarnIndex].push_back(i);
2247 }
2248
2249 // Output element sets
2250 Output << "********************" << endl;
2251 Output << "*** ELEMENT SETS ***" << endl;
2252 Output << "********************" << endl;
2253 Output << "** TexGen generates a number of element sets:" << endl;
2254 Output << "** All - Contains all elements" << endl;
2255 Output << "** Matrix - Contains all elements belonging to the matrix" << endl;
2256 Output << "** YarnX - Where X represents the yarn index" << endl;
2257 Output << "*ElSet, ElSet=All, Generate" << endl;
2258 Output << "1, " << GetNumElements() << ", 1" << endl;
2259 vector<int>::iterator itIndices;
2260 map<int, vector<int> >::iterator itElementSet;
2261 for (itElementSet = ElementSets.begin(); itElementSet != ElementSets.end(); ++itElementSet)
2262 {
2263 if (itElementSet->first == -1)
2264 Output << "*ElSet, ElSet=Matrix" << endl;
2265 else
2266 Output << "*ElSet, ElSet=Yarn" << itElementSet->first << endl;
2267 int iLinePos = 0;
2268 for (itIndices = itElementSet->second.begin(); itIndices != itElementSet->second.end(); ++itIndices)
2269 {
2270 if (iLinePos == 0)
2271 {
2272 // Do nothing...
2273 }
2274 else if (iLinePos < 16)
2275 {
2276 Output << ", ";
2277 }
2278 else
2279 {
2280 Output << endl;
2281 iLinePos = 0;
2282 }
2283 Output << *itIndices;
2284 ++iLinePos;
2285 }
2286 Output << endl;
2287 }
2288
2289 // Output materials, this is only here because it is necessary in order to associate
2290 // orientations and other things to the yarns. It just creates a linear elastic isotropic
2291 // material with young's modulus of 1.
2292 if (bCreateMaterial)
2293 {
2294 Output << "*****************" << endl;
2295 Output << "*** MATERIALS ***" << endl;
2296 Output << "*****************" << endl;
2297 Output << "*Material, Name=TexGenGenerated" << endl;
2298 Output << "*Elastic" << endl;
2299 Output << "1.0, 0.0" << endl;
2300 Output << "*DepVar" << endl;
2301 Output << "5" << endl;
2302 // Output << "*Solid Section, ElSet=All, Material=TexGenGenerated, Orientation=TexGenOrientations" << endl;
2303 // Output << "1.0," << endl;
2304 for (itElementSet = ElementSets.begin(); itElementSet != ElementSets.end(); ++itElementSet)
2305 {
2306 // Let's create 1 section definition for the matrix and each yarn so that it is easy
2307 // to apply different material properties to different yarns in ABAQUS CAE
2308 ostringstream SetName;
2309 if (itElementSet->first == -1)
2310 {
2311 SetName << "Matrix";
2312 }
2313 else
2314 {
2315 SetName << "Yarn" << itElementSet->first;
2316 }
2317 Output << "*Solid Section, ElSet=" << SetName.str() << ", Material=TexGenGenerated, Orientation=TexGenOrientations" << endl;
2318 Output << "1.0," << endl;
2319 }
2320 Output << "** Note: Additional element data are stored as a depvars:" << endl;
2321 Output << "** 1 - Yarn Index (-1 for matrix, first yarn starting at 0)" << endl;
2322 Output << "** 2/3 - Location (x and y cross-section coordinates of element relative to yarn centerline)" << endl;
2323 Output << "** 4 - Volume fraction" << endl;
2324 Output << "** 5 - Distance of element from the surface of the yarn (for yarn elements only, distance is negative)" << endl;
2325 Output << "*Initial Conditions, Type=Solution, Input=" << StripPath(ElementDataFilename) << endl;
2326 }
2327
2328 if (bCreateStep)
2329 {
2330 Output << "************" << endl;
2331 Output << "*** STEP ***" << endl;
2332 Output << "************" << endl;
2333 Output << "*Step, Name=TexGenGenerated" << endl;
2334 Output << "*Static" << endl;
2335 Output << "*Output, Field, Variable=PRESELECT" << endl;
2336 // Output << "*Output, History, Variable=PRESELECT" << endl;
2337 Output << "*Element Output" << endl;
2338 Output << "SDV" << endl;
2339 Output << "*Node Print" << endl;
2340 Output << "U" << endl;
2341 Output << "*End Step" << endl;
2342 }
2343 }
2344
2345 return true;
2346}
2347
2348bool CMesh::SaveToVTK(string Filename, const vector<CMeshDataBase*> *pMeshData) const
2349{
2350 AddExtensionIfMissing(Filename, ".vtu");
2351
2352 TiXmlDocument doc;
2353 TiXmlDeclaration decl("1.0", "", "");
2354 doc.InsertEndChild(decl);
2355 TiXmlElement VTKFile("VTKFile");
2356 {
2357 VTKFile.SetAttribute("type", "UnstructuredGrid");
2358 VTKFile.SetAttribute("version", "0.1");
2359 VTKFile.SetAttribute("byte_order", "LittleEndian");
2360 TiXmlElement UnstructuredGrid("UnstructuredGrid");
2361 {
2362 TiXmlElement Piece("Piece");
2363 {
2364 TiXmlElement Points("Points");
2365 int iNumPoints = FillVTKPointData(Points);
2366 Piece.InsertEndChild(Points);
2367
2368 TiXmlElement Cells("Cells");
2369 int iNumCells = FillVTKCellData(Cells);
2370 Piece.InsertEndChild(Cells);
2371
2372 if (pMeshData)
2373 {
2374 TiXmlElement PointData("PointData");
2375 TiXmlElement CellData("CellData");
2376 vector<CMeshDataBase*>::const_iterator itData;
2377 for (itData = pMeshData->begin(); itData != pMeshData->end(); ++itData)
2378 {
2379 if ((*itData)->m_DataType == CMeshDataBase::NODE)
2380 (*itData)->InsertVTKData(PointData);
2381 else
2382 (*itData)->InsertVTKData(CellData);
2383 }
2384 TGLOG("Finished iterator");
2385 Piece.InsertEndChild(PointData);
2386 Piece.InsertEndChild(CellData);
2387 TGLOG("End InsertVTKData");
2388 }
2389
2390 Piece.SetAttribute("NumberOfPoints", stringify(iNumPoints));
2391 Piece.SetAttribute("NumberOfCells", stringify(iNumCells));
2392 }
2393 UnstructuredGrid.InsertEndChild(Piece);
2394 }
2395 VTKFile.InsertEndChild(UnstructuredGrid);
2396 }
2397
2398 doc.InsertEndChild(VTKFile);
2399 TGLOG("Return call SaveFile");
2400 return doc.SaveFile(Filename);
2401}
2402/*
2403bool CMesh::SaveToVTK(string Filename, const vector<TiXmlElement> *pAdditionalData) const
2404{
2405 AddExtensionIfMissing(Filename, ".vtu");
2406
2407 TiXmlDocument doc;
2408 TiXmlDeclaration decl("1.0", "", "");
2409 doc.InsertEndChild(decl);
2410 TiXmlElement VTKFile("VTKFile");
2411 {
2412 VTKFile.SetAttribute("type", "UnstructuredGrid");
2413 VTKFile.SetAttribute("version", "0.1");
2414 VTKFile.SetAttribute("byte_order", "LittleEndian");
2415 TiXmlElement UnstructuredGrid("UnstructuredGrid");
2416 {
2417 TiXmlElement Piece("Piece");
2418 {
2419 TiXmlElement Points("Points");
2420 int iNumPoints = FillVTKPointData(Points);
2421 Piece.InsertEndChild(Points);
2422
2423 TiXmlElement Cells("Cells");
2424 int iNumCells = FillVTKCellData(Cells);
2425 Piece.InsertEndChild(Cells);
2426
2427 if (pAdditionalData)
2428 {
2429 vector<TiXmlElement>::const_iterator itElement;
2430 for (itElement = pAdditionalData->begin(); itElement != pAdditionalData->end(); ++itElement)
2431 Piece.InsertEndChild(*itElement);
2432 }
2433
2434 Piece.SetAttribute("NumberOfPoints", stringify(iNumPoints));
2435 Piece.SetAttribute("NumberOfCells", stringify(iNumCells));
2436 }
2437 UnstructuredGrid.InsertEndChild(Piece);
2438 }
2439 VTKFile.InsertEndChild(UnstructuredGrid);
2440 }
2441
2442 doc.InsertEndChild(VTKFile);
2443
2444 return doc.SaveFile(Filename);
2445}
2446*/
2447int CMesh::FillVTKPointData(TiXmlElement &Points) const
2448{
2449 TGLOG("FillVTKPointData");
2450 ostringstream PointsData;
2451 vector<XYZ>::const_iterator itNode;
2452 for (itNode = m_Nodes.begin(); itNode != m_Nodes.end(); ++itNode)
2453 {
2454 PointsData << itNode->x << " " << itNode->y << " " << itNode->z << " ";
2455 }
2456
2457 TiXmlElement DataArray("DataArray");
2458 {
2459 DataArray.SetAttribute("type", "Float32");
2460 DataArray.SetAttribute("NumberOfComponents", "3");
2461 DataArray.SetAttribute("format", "ascii");
2462 DataArray.InsertEndChild(TiXmlText(PointsData.str()));
2463 }
2464 Points.InsertEndChild(DataArray);
2465 TGLOG("return");
2466 return m_Nodes.size();
2467}
2468
2469int CMesh::FillVTKCellData(TiXmlElement &Cells) const
2470{
2471 TGLOG("FillVTKCellData");
2472 ostringstream ConnectivityData;
2473 ostringstream OffsetsData;
2474 ostringstream TypesData;
2475
2476 int i, iOffset = 0;
2477 int iTotalNumCells = 0;
2478 list<int>::const_iterator itIndex;
2479 int j;
2480 for (j = 0; j < NUM_ELEMENT_TYPES; ++j)
2481 {
2482 int iVTKType = 0;
2483 switch (j)
2484 {
2485 case TRI:
2486 iVTKType = 5;
2487 break;
2488 case QUAD:
2489 iVTKType = 9;
2490 break;
2491 case TET:
2492 iVTKType = 10;
2493 break;
2494 case WEDGE:
2495 iVTKType = 13;
2496 break;
2497 case PYRAMID:
2498 iVTKType = 14;
2499 break;
2500 case HEX:
2501 iVTKType = 12;
2502 break;
2503 case LINE:
2504 iVTKType = 3;
2505 break;
2506 case POLYLINE:
2507 iVTKType = 4;
2508 break;
2509 case QUADRATIC_TET:
2510 iVTKType = 24;
2511 break;
2512 }
2513 if (iVTKType != 0)
2514 {
2515 int iNumNodesPerElement = GetNumNodes((ELEMENT_TYPE)j);
2516 int iNumElements = m_Indices[j].size()/iNumNodesPerElement;
2517 for (itIndex = m_Indices[j].begin(), i=0; itIndex != m_Indices[j].end(); ++itIndex, ++i)
2518 {
2519 ConnectivityData << *itIndex << " ";
2520 }
2521 for (i=0; i<iNumElements; ++i)
2522 {
2523 TypesData << iVTKType << " ";
2524 }
2525 for (i=0; i<iNumElements; ++i)
2526 {
2527 iOffset += iNumNodesPerElement;
2528 OffsetsData << iOffset << " ";
2529 }
2530 iTotalNumCells += iNumElements;
2531 }
2532 }
2533
2534
2535 TiXmlElement Connectivity("DataArray");
2536 {
2537 Connectivity.SetAttribute("type", "Int32");
2538 Connectivity.SetAttribute("Name", "connectivity");
2539 Connectivity.SetAttribute("format", "ascii");
2540 Connectivity.InsertEndChild(TiXmlText(ConnectivityData.str()));
2541 }
2542 Cells.InsertEndChild(Connectivity);
2543
2544 TiXmlElement Offsets("DataArray");
2545 {
2546 Offsets.SetAttribute("type", "Int32");
2547 Offsets.SetAttribute("Name", "offsets");
2548 Offsets.SetAttribute("format", "ascii");
2549 Offsets.InsertEndChild(TiXmlText(OffsetsData.str()));
2550 }
2551 Cells.InsertEndChild(Offsets);
2552
2553 TiXmlElement Types("DataArray");
2554 {
2555 Types.SetAttribute("type", "Int32");
2556 Types.SetAttribute("Name", "types");
2557 Types.SetAttribute("format", "ascii");
2558 Types.InsertEndChild(TiXmlText(TypesData.str()));
2559 }
2560 Cells.InsertEndChild(Types);
2561 TGLOG("return");
2562 return iTotalNumCells;
2563}
2564
2565bool CMesh::AddElement(ELEMENT_TYPE Type, const vector<int> &Indices)
2566{
2567 if (Type < 0 || Type >= NUM_ELEMENT_TYPES)
2568 {
2569 TGERROR("Tried to add element of invalid type: " << Type);
2570 return false;
2571 }
2572 if ((int)Indices.size() != GetNumNodes(Type) && Type != CMesh::POLYGON )
2573 {
2574 TGERROR("Tried to add element of type " << Type << " with invalid number of indices: " << Indices.size());
2575 return false;
2576 }
2577 if ( Type == CMesh::POLYGON && *(Indices.begin()) != *(Indices.end()-1) )
2578 {
2579 TGERROR("Tried to add unclosed POLYGON to mesh" );
2580 return false;
2581 }
2582 m_Indices[Type].insert(m_Indices[Type].end(), Indices.begin(), Indices.end());
2583 return true;
2584}
2585
2586int CMesh::GetNumElements(ELEMENT_TYPE Type) const
2587{
2588 assert(Type >= 0 && Type < NUM_ELEMENT_TYPES);
2589 return m_Indices[Type].size()/GetNumNodes(Type);
2590}
2591
2592int CMesh::GetNumElements() const
2593{
2594 int iNumElements = 0;
2595 int i;
2596 for (i = 0; i < NUM_ELEMENT_TYPES; ++i)
2597 {
2598 if ( (ELEMENT_TYPE)i != CMesh::POLYGON )
2599 iNumElements += GetNumElements((ELEMENT_TYPE)i);
2600 }
2601 return iNumElements;
2602}
2603
2604vector<XYZ>::const_iterator CMesh::NodesBegin() const
2605{
2606 return m_Nodes.begin();
2607}
2608
2609vector<XYZ>::const_iterator CMesh::NodesEnd() const
2610{
2611 return m_Nodes.end();
2612}
2613
2614vector<XYZ>::iterator CMesh::NodesBegin()
2615{
2616 return m_Nodes.begin();
2617}
2618
2619vector<XYZ>::iterator CMesh::NodesEnd()
2620{
2621 return m_Nodes.end();
2622}
2623
2624const int CMesh::AddNode(XYZ Node)
2625{
2626 m_Nodes.push_back(Node);
2627 return (int)m_Nodes.size()-1;
2628}
2629
2630void CMesh::SetNode(int iIndex, XYZ Node)
2631{
2632 assert(iIndex >= 0 && iIndex<(int)m_Nodes.size());
2633 m_Nodes[iIndex] = Node;
2634}
2635
2636const XYZ& CMesh::GetNode(int iIndex) const
2637{
2638 return m_Nodes[iIndex];
2639}
2640
2641vector<XYZ>::iterator CMesh::DeleteNode(vector<XYZ>::iterator it)
2642{
2643 return m_Nodes.erase(it);
2644}
2645
2646int CMesh::GetNumNodes() const
2647{
2648 return (int)m_Nodes.size();
2649}
2650
2651const bool CMesh::NodesEmpty() const
2652{
2653 return m_Nodes.empty();
2654}
2655
2656const vector<XYZ>& CMesh::GetNodes() const
2657{
2658 return m_Nodes;
2659}
2660
2661vector<XYZ>& CMesh::GetNodes()
2662{
2663 return m_Nodes;
2664}
2665
2666void CMesh::SetNumNodes(int NumNodes)
2667{
2668 m_Nodes.resize(NumNodes);
2669}
2670
2671const list<int>& CMesh::GetIndices(ELEMENT_TYPE ElemType) const
2672{
2673 assert(ElemType >= 0 && ElemType < NUM_ELEMENT_TYPES);
2674 return m_Indices[ElemType];
2675}
2676
2677list<int>& CMesh::GetIndices(ELEMENT_TYPE ElemType)
2678{
2679 assert(ElemType >= 0 && ElemType < NUM_ELEMENT_TYPES);
2680 return m_Indices[ElemType];
2681}
2682
2683void CMesh::ConvertElementListToVector( ELEMENT_TYPE ElementType, vector<int> &Indices )
2684{
2685 assert(ElementType >= 0 && ElementType < NUM_ELEMENT_TYPES);
2686
2687 Indices.insert( Indices.begin(), m_Indices[ElementType].begin(), m_Indices[ElementType].end() );
2688}
2689
2690bool CMesh::SaveToSCIRun(string Filename)
2691{
2692 Filename = RemoveExtension( Filename, ".pts" );
2693 AddExtensionIfMissing(Filename, ".tri.pts");
2694
2695 ofstream Output(Filename.c_str());
2696
2697 if (!Output)
2698 {
2699 TGERROR("Unable to output mesh to SCIRun .pts file format, could not open file: " << Filename);
2700 return false;
2701 }
2702
2703 TGLOG("Saving mesh data to " << Filename);
2704
2705
2706 OutputNodes(Output, 1, ",", true );
2707
2708 int iStartIndex = 1;
2709
2710 Filename = RemoveExtension( Filename, ".pts" );
2711 AddExtensionIfMissing( Filename, ".fac" );
2712
2713 ofstream ElementOutput( Filename.c_str());
2714
2715 if (!ElementOutput)
2716 {
2717 TGERROR("Unable to output mesh to SCIRun .fac file format, could not open file: " << Filename);
2718 return false;
2719 }
2720
2721 TGLOG("Saving mesh data to " << Filename);
2722
2723 ConvertQuadstoTriangles();
2724
2725 // Tri elements
2726 if (!m_Indices[TRI].empty())
2727 {
2728 iStartIndex = OutputElements(ElementOutput, TRI, iStartIndex, 1, ", ", true);
2729 }
2730 return true;
2731}
2732
2733
TOL
#define TOL
Definition: AdjustMeshInterference.cpp:27
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
FOR_EACH_TIXMLELEMENT
#define FOR_EACH_TIXMLELEMENT(CHILDELEMENT, PARENTELEMENT, ELEMENTNAME)
Macro to enable looping over tinyxml easier.
Definition: Misc.h:45
TEXGEN
#define TEXGEN
Helper macro to get the texgen instance.
Definition: TexGen.h:76
TexGen::CMesh
Defines the nodes and elements of a surface or volume mesh.
Definition: Mesh.h:58
TexGen::CMesh::RemoveOpposingQuads
void RemoveOpposingQuads()
Remove quads that have the same indices but opposite normals.
Definition: Mesh.cpp:174
TexGen::CMesh::FillVTKPointData
int FillVTKPointData(TiXmlElement &Points) const
Definition: Mesh.cpp:2447
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::ConvertHextoWedge
void ConvertHextoWedge(bool bQuality=true)
Convert Hexahedral elements to Wedge elements.
Definition: Mesh.cpp:882
TexGen::CMesh::GetNumNodes
int GetNumNodes() const
Return the number of nodes.
Definition: Mesh.cpp:2646
TexGen::CMesh::GetNumElements
int GetNumElements() const
Get the total number of elements in the mesh.
Definition: Mesh.cpp:2592
TexGen::CMesh::RemoveAllElementsExcept
void RemoveAllElementsExcept(ELEMENT_TYPE Type)
Remove all elements except those of given type.
Definition: Mesh.cpp:711
TexGen::CMesh::ChangeNodeIndices
void ChangeNodeIndices(int iChangeTo, int iChangeFrom)
Change all the indices from one number to another.
Definition: Mesh.cpp:107
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::GetNodeElementReferences
void GetNodeElementReferences(vector< vector< int * > > &References)
Get a list of elements which reference each node.
Definition: Mesh.cpp:1504
TexGen::CMesh::RemoveDuplicateSegments
void RemoveDuplicateSegments()
Remove segments which have the same indices.
Definition: Mesh.cpp:314
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::ConvertQuadstoTriangles
void ConvertQuadstoTriangles(bool bQuality=true)
Convert the quad elements to triangles.
Definition: Mesh.cpp:1088
TexGen::CMesh::RemoveDegenerateTriangles
void RemoveDegenerateTriangles()
Remove triangles which have two equal corner indices.
Definition: Mesh.cpp:214
TexGen::CMesh::GetNodePairs
vector< pair< int, int > > GetNodePairs(XYZ Vector, const double Tolerance=1e-6) const
Return a list of node pairs (A, B) where A == B + Vector.
Definition: Mesh.cpp:463
TexGen::CMesh::ConvertHextoQuad
void ConvertHextoQuad()
Convert Hex elements to Quad elements representing their surface.
Definition: Mesh.cpp:765
TexGen::CMesh::GetNodes
const vector< XYZ > & GetNodes() const
Get a const reference to the nodes.
Definition: Mesh.cpp:2656
TexGen::CMesh::RemoveDuplicateTriangles
void RemoveDuplicateTriangles()
Definition: Mesh.cpp:276
TexGen::CMesh::SaveToSMESH
bool SaveToSMESH(string Filename) const
Save the mesh to a .smesh file to be used in conjunction with tetgen.
Definition: Mesh.cpp:2009
TexGen::CMesh::Convert3Dto2D
void Convert3Dto2D()
Convert all 3D elements to 2D elements representing their surface.
Definition: Mesh.cpp:576
TexGen::CMesh::SaveToABAQUS
bool SaveToABAQUS(string Filename, const vector< POINT_INFO > *pElementInfo=NULL, bool bCreateStep=true, bool bCreateMaterial=true, int iElementType=0)
Save the mesh to ABAQUS input file format with information such as yarn tangents.
Definition: Mesh.cpp:2089
TexGen::CMesh::WriteBinaryXYZ
static void WriteBinaryXYZ(ostream &Output, XYZ Vector)
Definition: Mesh.cpp:1922
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::PopulateTiXmlElement
void PopulateTiXmlElement(TiXmlElement &Element, OUTPUT_TYPE OutputType) const
Used for saving data to XML.
Definition: Mesh.cpp:57
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::GetElementCenters
vector< XYZ > GetElementCenters() const
Get a vector of element centers, one entry for each element.
Definition: Mesh.cpp:1549
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::ConvertPyramidto2D
void ConvertPyramidto2D()
Convert Pyramid elements to Quad and Triangle elements representing their surface.
Definition: Mesh.cpp:1008
TexGen::CMesh::ConvertTriToQuad
void ConvertTriToQuad(double Tolerance=1e-6)
Convert Triangel elements to Quads.
Definition: Mesh.cpp:606
TexGen::CMesh::ConvertToTetMesh
void ConvertToTetMesh()
Definition: Mesh.cpp:584
TexGen::CMesh::OutputElements
int OutputElements(ostream &Output, ELEMENT_TYPE ElementType, int iStartIndex=1, int iIndexOffset=1, string Seperator=", ", bool bSCIRun=false) const
Output the element indices in the mesh each on a new line, indices seperated by commas.
Definition: Mesh.cpp:1480
TexGen::CMesh::m_Indices
list< int > m_Indices[NUM_ELEMENT_TYPES]
Map of indices into the nodes.
Definition: Mesh.h:515
TexGen::CMesh::MeshConvexHull
void MeshConvexHull()
Create a triangular convex hull of the nodes contained within the mesh.
Definition: Mesh.cpp:1758
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::ConvertTrianglestoSegments
void ConvertTrianglestoSegments()
Convert triangle elements to segments.
Definition: Mesh.cpp:1066
TexGen::CMesh::Translate
void Translate(XYZ Vector)
Translate whole mesh by given vector.
Definition: Mesh.cpp:1139
TexGen::CMesh::SaveToSCIRun
bool SaveToSCIRun(string Filename)
Save the mesh to SCIRun format.
Definition: Mesh.cpp:2690
TexGen::CMesh::~CMesh
~CMesh(void)
Definition: Mesh.cpp:37
TexGen::CMesh::SaveToSTL
bool SaveToSTL(string Filename, bool bBinary=true) const
Save the mesh to STL file.
Definition: Mesh.cpp:1933
TexGen::CMesh::ConvertHextoTet
void ConvertHextoTet(bool bQuality=true)
Convert Hexahedral elements to Tetrahedral elements.
Definition: Mesh.cpp:1002
TexGen::CMesh::GetClosestNode
int GetClosestNode(XYZ Position) const
Get the index of the node closest to the given position space.
Definition: Mesh.cpp:532
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::RemoveOpposingTriangles
void RemoveOpposingTriangles()
Remove triangles that have the same indices but opposite normals.
Definition: Mesh.cpp:136
TexGen::CMesh::DeleteNodes
int DeleteNodes(const set< int > &Nodes)
Delete nodes and adjust node indices.
Definition: Mesh.cpp:728
TexGen::CMesh::ELEMENT_TYPE
ELEMENT_TYPE
Each element type is represented by a unique integer value.
Definition: Mesh.h:66
TexGen::CMesh::NUM_ELEMENT_TYPES
@ NUM_ELEMENT_TYPES
Definition: Mesh.h:77
TexGen::CMesh::PYRAMID
@ PYRAMID
Definition: Mesh.h:70
TexGen::CMesh::POLYLINE
@ POLYLINE
Definition: Mesh.h:74
TexGen::CMesh::QUADRATIC_TET
@ QUADRATIC_TET
Definition: Mesh.h:75
TexGen::CMesh::POLYGON
@ POLYGON
Definition: Mesh.h:76
TexGen::CMesh::ConvertPyramidtoTet
void ConvertPyramidtoTet(bool bQuality=true)
Convert Pyramid elements to Tetrahedral elements.
Definition: Mesh.cpp:949
TexGen::CMesh::OutputNodes
int OutputNodes(ostream &Output, int iStartIndex=1, string Seperator=", ", bool bSCIRun=false) const
Output the node coordinates in the mesh each on a new line, components seperated by commas.
Definition: Mesh.cpp:1463
TexGen::CMesh::SetNode
void SetNode(int iIndex, XYZ Node)
Set the node at given index.
Definition: Mesh.cpp:2630
TexGen::CMesh::MeshClosedLoop
void MeshClosedLoop(const XYZ &Normal, const vector< int > &ClosedLoopVector, bool bQuality=false)
Definition: Mesh.cpp:1200
TexGen::CMesh::CopySelfToRange
void CopySelfToRange(XYZ Vector, int iLowerLimit, int iUpperLimit)
Copy the mesh to the range given by the lower limit and upper limit with the given repeat vector.
Definition: Mesh.cpp:1437
TexGen::CMesh::ConvertWedgetoTetandPyramid
void ConvertWedgetoTetandPyramid(bool bQuality=true)
Convert Wedge elements to a combination of Tetrahedral and Pyramid elements.
Definition: Mesh.cpp:918
TexGen::CMesh::ConvertToTriangleMesh
void ConvertToTriangleMesh()
Definition: Mesh.cpp:592
TexGen::CMesh::ConvertWedgeto2D
void ConvertWedgeto2D()
Convert Wedge elements to Quad and Triangle elements representing their surface.
Definition: Mesh.cpp:813
TexGen::CMesh::m_Nodes
vector< XYZ > m_Nodes
List of nodes.
Definition: Mesh.h:513
TexGen::CMesh::BuildGrid
void BuildGrid(XYZ Min, XYZ Max, int iNumX, int iNumY, int iNumZ)
Build grid of points.
Definition: Mesh.cpp:1881
TexGen::CMesh::ConvertWedgetoTet
void ConvertWedgetoTet(bool bQuality=true)
Convert Wedge elements to Tetrahedral elements.
Definition: Mesh.cpp:996
TexGen::CMesh::SetNumNodes
void SetNumNodes(int NumNodes)
Resize the vector size.
Definition: Mesh.cpp:2666
TexGen::CMesh::FlipNormals
void FlipNormals()
Flip the normals of the mesh for triangles and quads.
Definition: Mesh.cpp:1166
TexGen::CMesh::ConvertElementListToVector
void ConvertElementListToVector(ELEMENT_TYPE ElementType, vector< int > &Indices)
Convert element index list to vector so can access by index.
Definition: Mesh.cpp:2683
TexGen::CMesh::CalculateVolume
double CalculateVolume() const
Calculate the volume of the mesh.
Definition: Mesh.cpp:1519
TexGen::CMesh::GetClosestNodeDistance
int GetClosestNodeDistance(XYZ Position, double dTol) const
Get the index of the node within a given tolerance distance to the given position space.
Definition: Mesh.cpp:550
TexGen::CMesh::CMesh
CMesh(void)
Definition: Mesh.cpp:33
TexGen::CMesh::GetAABB
pair< XYZ, XYZ > GetAABB(double dGrowDistance=0) const
Get an axis aligned bounding box for the mesh.
Definition: Mesh.cpp:340
TexGen::CMesh::Clear
void Clear()
Empty mesh nodes and indices.
Definition: Mesh.cpp:1448
TexGen::CMesh::Rotate
void Rotate(WXYZ Rotation, XYZ Origin=XYZ(0, 0, 0))
Rotate the whole mesh by given quaternion.
Definition: Mesh.cpp:1148
TexGen::CMesh::AddOrCancel
static void AddOrCancel(list< pair< int, int > > &EdgeStack, pair< int, int > Edge)
Add an edge to the stack of edges if it doesn't already exist otherwise if it exists delete it.
Definition: Mesh.cpp:1740
TexGen::CMesh::FillVTKCellData
int FillVTKCellData(TiXmlElement &Cells) const
Definition: Mesh.cpp:2469
TexGen::CMesh::ConvertToSegmentMesh
void ConvertToSegmentMesh()
Convert all elements to segments.
Definition: Mesh.cpp:599
TexGen::CMesh::NodesBegin
vector< XYZ >::const_iterator NodesBegin() const
Definition: Mesh.cpp:2604
TexGen::CMesh::NodesEmpty
const bool NodesEmpty() const
Returns true if the nodes array is empty.
Definition: Mesh.cpp:2651
TexGen::CMesh::DeleteNode
vector< XYZ >::iterator DeleteNode(vector< XYZ >::iterator it)
Delete a node given iterator.
Definition: Mesh.cpp:2641
TexGen::CMesh::CountInvertedElements
int CountInvertedElements() const
Check and count if any of the elements are inverted.
Definition: Mesh.cpp:1595
TexGen::CMesh::InsertNodes
int InsertNodes(const CMesh &Mesh, XYZ Offset=XYZ(0, 0, 0))
Add the nodes of Mesh to this mesh.
Definition: Mesh.cpp:82
TexGen::CMesh::ConvertToSurfaceMesh
void ConvertToSurfaceMesh()
Convert a volume mesh into a surface mesh (interior surfaces are elliminated)
Definition: Mesh.cpp:569
TexGen::CMesh::NodesEnd
vector< XYZ >::const_iterator NodesEnd() const
Definition: Mesh.cpp:2609
TexGen::CMesh::RemoveDuplicateElements
void RemoveDuplicateElements(CMesh::ELEMENT_TYPE ElementType)
Remove duplicate elements which have the same indices (leaves one copy of element)
Definition: Mesh.cpp:234
TexGen::CMesh::IntersectLine
int IntersectLine(const XYZ &P1, const XYZ &P2, vector< pair< double, XYZ > > &IntersectionPoints, pair< bool, bool > TrimResults=make_pair(false, false), bool bForceFind=false) const
Find the points where a line intersects the mesh.
Definition: Mesh.cpp:1633
TexGen::CMesh::ConvertTettoTriangle
void ConvertTettoTriangle()
Convert Tet elements to Triangle elements representing their surface.
Definition: Mesh.cpp:852
TexGen::COctreeAgentNode
Octree agent used to add indexed nodes to an octree.
Definition: MeshOctreeClasses.h:31
TexGen::COctreeVisitorMergeNodes
Octree visitor used to merge nodes together within a given tolerance.
Definition: MeshOctreeClasses.h:45
TexGen
Namespace containing a series of customised math operations not found in the standard c++ library.
Definition: AdjustMeshInterference.h:27
TexGen::PointInsideTriangleAccuracy
double PointInsideTriangleAccuracy(const XYZ &P1, const XYZ &P2, const XYZ &P3, const XYZ &P, const XYZ &N)
P1, P2, P3 are the three triangle corners, P is the points to be tested, N is the normal to the trian...
Definition: mymath.h:1285
TexGen::RemoveExtension
std::string RemoveExtension(std::string &Filename, std::string Extension)
Strip the extension from the filename.
Definition: Misc.cpp:122
TexGen::StripPath
std::string StripPath(std::string Filename)
Strip the path from the filename (e.g. "c:\folder\file.ext -> file.ext")
Definition: Misc.cpp:107
TexGen::GetLengthSquared
double GetLengthSquared(const XYZ &Point1, const XYZ &Point2)
Get the length squared between two points.
Definition: mymath.h:548
TexGen::WriteElementsHeader
void WriteElementsHeader(std::ostream &Output)
Write elements header for ABAQUS .eld files.
Definition: Misc.cpp:167
TexGen::WriteOrientationsHeader
void WriteOrientationsHeader(std::ostream &Output)
Write orientations header for ABAQUS .ori files.
Definition: Misc.cpp:157
TexGen::GetIntersectionLinePlane
bool GetIntersectionLinePlane(const XYZ &p1, const XYZ &p2, const XYZ &p3, const XYZ &N, XYZ &Intersection, double *pdU=NULL)
P1 and P2 are two points on the line, P3 is a point on the plane, N is the plane normal,...
Definition: mymath.h:1025
TexGen::OUTPUT_TYPE
OUTPUT_TYPE
Definition: Misc.h:105
TexGen::Round
int Round(double dValue)
Definition: mymath.h:1169
TexGen::Max
double Max(XYZ &Vector)
Get maximum element of vector and return it.
Definition: mymath.h:642
TexGen::GetLength
double GetLength(const XYZ &Point1, const XYZ &Point2)
Get the length between two points.
Definition: mymath.h:540
TexGen::stringify
std::string stringify(const T &x, int iPrecision=12, bool bScientific=true)
Function to convert a value (e.g. int, double, etc...) to a string.
Definition: Misc.h:50
TexGen::Min
XYZ Min(const XYZ &P1, const XYZ &P2)
Given two points, return a new point who's coordinates are the smaller of the two.
Definition: mymath.h:1142
TexGen::Normalise
WXYZ & Normalise(WXYZ &Quaternion)
Normalise the quaternion and return it.
Definition: mymath.h:609
TexGen::AddExtensionIfMissing
void AddExtensionIfMissing(std::string &Filename, std::string Extension)
Adds an extension to the filename if it is missing otherwise do nothing (e.g. picture -> picture....
Definition: Misc.cpp:116
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::LessPairDoubleXYZ
Used to sort double-XYZ pairs.
Definition: Misc.h:92
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::WXYZ
Struct for representing a quaternion.
Definition: mymath.h:38
TexGen::XYZ
Struct for representing points in 3D space.
Definition: mymath.h:56
TexGen::XYZ::z
double z
Definition: mymath.h:57
TexGen::XYZ::x
double x
Definition: mymath.h:57
TexGen::XYZ::y
double y
Definition: mymath.h:57
Generated by doxygen 1.9.2