TexGen
VoxelMesh.cpp
Go to the documentation of this file.
1/*=============================================================================
2TexGen: Geometric textile modeller.
3Copyright (C) 2010 Louise Brown
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 "VoxelMesh.h"
22#include "TexGen.h"
23#include "PeriodicBoundaries.h"
27#include <iterator>
28//#define SHINY_PROFILER TRUE
29
30using namespace TexGen;
31
32CVoxelMesh::CVoxelMesh(string Type)
33{
34 if ( Type == "CShearedPeriodicBoundaries" )
36 else if ( Type == "CStaggeredPeriodicBoundaries" )
38 else if ( Type == "CRotatedPeriodicBoundaries" )
40 else if ( Type == "CBendingPeriodicBoundaries" )
42 else
44}
45
47{
49}
50
51void CVoxelMesh::SaveVoxelMesh(CTextile &Textile, string OutputFilename, int XVoxNum, int YVoxNum, int ZVoxNum,
52 bool bOutputMatrix, bool bOutputYarns, int iBoundaryConditions, int iElementType, int FileType)
53{
54 //PROFILE_SHARED_DEFINE(ProfileTest)
55 //PROFILE_FUNC()
56
57 const CDomain* pDomain = Textile.GetDomain();
58 if (!pDomain)
59 {
60 TGERROR("Unable to create ABAQUS input file: No domain specified");
61 return;
62 }
63 //PROFILE_SHARED_BEGIN(ProfileTest);
64 m_XVoxels = XVoxNum;
65 m_YVoxels = YVoxNum;
66 m_ZVoxels = ZVoxNum;
67 TGLOG("Calculating voxel sizes");
68 if (!CalculateVoxelSizes(Textile))
69 {
70 TGERROR("Unable to create ABAQUS input file: Error calculating voxel sizes");
71 return;
72 }
73 TGLOG("Replacing spaces in filename with underscore for ABAQUS compatibility");
74 OutputFilename = ReplaceFilenameSpaces(OutputFilename);
75 //GetYarnGridIntersections(Textile);
76 if (FileType == INP_EXPORT)
77 {
78 CTimer timer;
79 timer.start("Timing SaveToAbaqus");
80 SaveToAbaqus(OutputFilename, Textile, bOutputMatrix, bOutputYarns, iBoundaryConditions, iElementType);
81 timer.check("End of SaveToAbaqus");
82 timer.stop();
83 }
84 else
85 SaveVoxelMeshToVTK(OutputFilename, Textile);
86
87 m_ElementsInfo.clear(); // Clear point_info data as otherwise retains memory space until create another voxel mesh or exit program
88
89 // PROFILE_END();
90 // PROFILER_UPDATE();
91 // PROFILER_OUTPUT("ProfileOutput.txt");
92 //SaveToSCIRun( OutputFilename, Textile );
93}
94
95/*void CVoxelMesh::SaveShearedVoxelMesh( CTextile &Textile, string OutputFilename, int XVoxNum, int YVoxNum, int ZVoxNum,
96 bool bOutputMatrix, bool bOutputYarns, int iBoundariesUntiedZ, int iElementType )
97{
98 const CDomain* pDomain = Textile.GetDomain();
99 if (!pDomain)
100 {
101 TGERROR("Unable to create ABAQUS input file: No domain specified");
102 return;
103 }
104 //PROFILE_BEGIN(Begin);
105 m_XVoxels = XVoxNum;
106 m_YVoxels = YVoxNum;
107 m_ZVoxels = ZVoxNum;
108 TGLOG("Calculating voxel sizes");
109
110 CalculateShearedVoxelSizes( Textile );
111 TGLOG("Replacing spaces in filename with underscore for ABAQUS compatibility");
112 OutputFilename = ReplaceFilenameSpaces( OutputFilename );
113 //GetYarnGridIntersections(Textile);
114 CTimer timer;
115 timer.start("Timing SaveToAbaqus");
116 SaveToAbaqus( OutputFilename, Textile, bOutputMatrix, bOutputYarns, iBoundariesUntiedZ, iElementType );
117 timer.check("End of SaveToAbaqus");
118 timer.stop();
119
120 //PROFILE_END();
121 //PROFILER_UPDATE();
122 //PROFILER_OUTPUT("ProfileOutput.txt");
123 //SaveToSCIRun( OutputFilename, Textile );
124}*/
125
126// Leaving this in case want to utilise for rendering mesh at some stage (memory allowing)
127/*void CVoxelMesh::SaveVoxelMesh(CTextile &Textile, string OutputFilename, int XVoxNum, int YVoxNum, int ZVoxNum, bool bOutputMatrix, bool bOutputYarns )
128{
129 vector<POINT_INFO> ElementsInfo;
130 vector<POINT_INFO>::iterator itElementInfo;
131
132 list<int>::iterator itIndex;
133 list<int>::iterator itStartIndex;
134 int i;
135 int iNumNodes = m_Mesh.GetNumNodes(CMesh::HEX);
136 int iNumElements;
137
138 m_XVoxels = XVoxNum;
139 m_YVoxels = YVoxNum;
140 m_ZVoxels = ZVoxNum;
141 TGLOG("Calculating voxel sizes");
142 CalculateVoxelSizes( Textile );
143 TGLOG("Adding nodes");
144
145 AddNodes();
146 TGLOG("Adding elements");
147 AddElements();
148 Textile.GetPointInformation(m_Mesh.GetElementCenters(), ElementsInfo);
149 list<int> &Indices = m_Mesh.GetIndices(CMesh::HEX);
150 TGLOG("Deleting fibre elements");
151 iNumElements = ElementsInfo.size();
152 if ( !bOutputYarns )
153 {
154 for ( i = 0, itIndex = Indices.begin(); i < iNumElements && itIndex != Indices.end(); ++i)
155 {
156 if ( ElementsInfo[i].iYarnIndex == -1 )
157 {
158 advance( itIndex, iNumNodes );
159 }
160 else
161 {
162 itStartIndex = itIndex;
163 advance( itIndex, iNumNodes);
164 itIndex = Indices.erase( itStartIndex, itIndex );
165 }
166 }
167 }
168 else if ( !bOutputMatrix )
169 {
170 for ( i = 0, itIndex = Indices.begin(), itElementInfo = ElementsInfo.begin(); i < iNumElements && itIndex != Indices.end(); ++i)
171 {
172 if ( itElementInfo->iYarnIndex == -1 )
173 {
174 itStartIndex = itIndex;
175 advance( itIndex, iNumNodes);
176 itIndex = Indices.erase( itStartIndex, itIndex );
177 itElementInfo = ElementsInfo.erase( itElementInfo );
178 }
179 else
180 {
181 advance( itIndex, iNumNodes );
182 ++itElementInfo;
183 }
184 }
185 }
186 // else outputting both so don't need to delete anything
187
188 TGLOG("Save to ABAQUS");
189 //SaveVoxelMeshToVTK( "TestVoxelMesh", ElementsInfo );
190
191 m_Mesh.SaveToABAQUS(OutputFilename, bOutputYarns ? &ElementsInfo : NULL, false, false);
192 TGLOG("Finished save voxel mesh");
193}
194
195
196
197
198
199void CVoxelMesh::AddNodes()
200{
201 int x,y,z;
202
203 for ( z = 0; z <= m_ZVoxels; ++z )
204 {
205 for ( y = 0; y <= m_YVoxels; ++y )
206 {
207 for ( x = 0; x <=m_XVoxels; ++x )
208 {
209 XYZ Point;
210 Point.x = m_DomainAABB.first.x + m_VoxSize[0] * x;
211 Point.y = m_DomainAABB.first.y + m_VoxSize[1] * y;
212 Point.z = m_DomainAABB.first.z + m_VoxSize[2] * z;
213 m_Mesh.AddNode( Point );
214 }
215 }
216 }
217}
218
219void CVoxelMesh::AddElements()
220{
221 int numx = m_XVoxels + 1;
222 int numy = m_YVoxels + 1;
223 int x,y,z;
224
225 for ( z = 0; z < m_ZVoxels; ++z )
226 {
227 for ( y = 0; y < m_YVoxels; ++y )
228 {
229 for ( x = 0; x < m_XVoxels; ++x )
230 {
231 vector<int> Indices;
232 Indices.push_back(x + y*numx + z*numx*numy);
233 Indices.push_back((x+1) + y*numx + z*numx*numy);
234 Indices.push_back((x+1) +y*numx + (z+1)*numx*numy);
235 Indices.push_back(x +y*numx + (z+1)*numx*numy);
236 Indices.push_back(x + (y+1)*numx + z*numx*numy);
237 Indices.push_back((x+1) + (y+1)*numx + z*numx*numy);
238 Indices.push_back((x+1) +(y+1)*numx + (z+1)*numx*numy);
239 Indices.push_back(x +(y+1)*numx + (z+1)*numx*numy);
240 m_Mesh.AddElement(CMesh::HEX, Indices );
241 }
242 }
243 }
244}*/
245
246void CVoxelMesh::SaveVoxelMeshToVTK(string Filename, CTextile &Textile )
247{
248 AddExtensionIfMissing(Filename, ".vtu");
249 ofstream Output(Filename.c_str());
250
251 OutputNodes(Output, Textile, VTU_EXPORT);
252
253 OutputHexElements(Output, true, true, VTU_EXPORT);
254
255 CMeshData<int> YarnIndex("YarnIndex", CMeshDataBase::ELEMENT);
256 CMeshData<XYZ> YarnTangent("YarnTangent", CMeshDataBase::ELEMENT);
257 CMeshData<XY> Location("Location", CMeshDataBase::ELEMENT);
258 CMeshData<double> VolumeFraction("VolumeFraction", CMeshDataBase::ELEMENT);
259 CMeshData<double> SurfaceDistance("SurfaceDistance", CMeshDataBase::ELEMENT);
260 CMeshData<XYZ> Orientation("Orientation", CMeshDataBase::ELEMENT);
261
262 vector<POINT_INFO>::iterator itElementInfo;
263
264 for (itElementInfo = m_ElementsInfo.begin(); itElementInfo != m_ElementsInfo.end(); ++itElementInfo)
265 {
266 YarnIndex.m_Data.push_back(itElementInfo->iYarnIndex);
267 YarnTangent.m_Data.push_back(itElementInfo->YarnTangent);
268 Location.m_Data.push_back(itElementInfo->Location);
269 VolumeFraction.m_Data.push_back(itElementInfo->dVolumeFraction);
270 SurfaceDistance.m_Data.push_back(itElementInfo->dSurfaceDistance);
271 Orientation.m_Data.push_back(itElementInfo->Orientation);
272 }
273
274
275 vector<CMeshDataBase*> MeshData;
276
277 MeshData.push_back(&YarnIndex);
278 MeshData.push_back(&YarnTangent);
279 MeshData.push_back(&Location);
280 MeshData.push_back(&VolumeFraction);
281 MeshData.push_back(&SurfaceDistance);
282 MeshData.push_back(&Orientation);
283
284 m_Mesh.SaveToVTK(Filename, &MeshData);
285}
286
287void CVoxelMesh::SaveToAbaqus( string Filename, CTextile &Textile, bool bOutputMatrix, bool bOutputYarn, int iBoundaryConditions, int iElementType )
288{
289 //PROFILE_FUNC();
290 AddExtensionIfMissing(Filename, ".inp");
291
292 ofstream Output(Filename.c_str());
293
294 if (!Output)
295 {
296 TGERROR("Unable to output voxel mesh to ABAQUS file format, could not open file: " << Filename);
297 return;
298 }
299
300 TGLOG("Saving voxel mesh data to " << Filename);
301
302 Output << "*Heading" << "\n";
303 Output << "File generated by TexGen v" << TEXGEN.GetVersion() << "\n";
304
305 Output << "************" << "\n";
306 Output << "*** MESH ***" << "\n";
307 Output << "************" << "\n";
308 Output << "*Node" << "\n";
309 //PROFILE_BEGIN(OutputNodes);
310 OutputNodes(Output, Textile);
311 //PROFILE_END();
312 TGLOG("Outputting hex elements");
313 //Output the voxel HEX elements
314 int iNumHexElements = 0;
315 if ( !iElementType )
316 {
317 Output << "*Element, Type=C3D8R" << "\n";
318 }
319 else
320 {
321 Output << "*Element, Type=C3D8" << "\n";
322 }
323 //PROFILE_BEGIN(OutputHexElements);
324 iNumHexElements = OutputHexElements( Output, bOutputMatrix, bOutputYarn );
325 //PROFILE_END();
326 bool bMatrixOnly = false;
327 if ( bOutputMatrix && !bOutputYarn )
328 bMatrixOnly = true;
329
330 if ( bOutputYarn )
331 {
332 TGLOG("Outputting orientations & element sets");
333 //PROFILE_BEGIN(OutputOrientations);
334 OutputOrientationsAndElementSets( Filename, Output );
335 //PROFILE_END();
336 }
337 else if ( bMatrixOnly )
338 {
339 OutputMatrixElementSet( Filename, Output, iNumHexElements, bMatrixOnly );
340 }
341 //PROFILE_BEGIN(OutputNodeSets);
342 OutputAllNodesSet( Filename, Output );
343
344 // Output material properties
345 m_Materials.SetupMaterials( Textile );
346 m_Materials.OutputMaterials( Output, Textile.GetNumYarns(), bMatrixOnly );
347
348 //PROFILE_END();
349 if ( iBoundaryConditions != NO_BOUNDARY_CONDITIONS )
350 {
351 //PROFILE_BEGIN(OutputPBCs);
352 OutputPeriodicBoundaries( Output, Textile, iBoundaryConditions, bMatrixOnly );
353 //PROFILE_END();
354 }
355 TGLOG("Finished saving to Abaqus");
356}
357
358void CVoxelMesh::SaveToSCIRun( string Filename, CTextile &Textile )
359{
360 Filename = RemoveExtension( Filename, ".inp" );
361 Filename += ".hx.pts";
362
363 ofstream NodesFile(Filename.c_str());
364
365 if (!NodesFile)
366 {
367 TGERROR("Unable to output voxel mesh to SCIRun file format, could not open file: " << Filename);
368 return;
369 }
370
371 TGLOG("Saving voxel mesh data points to " << Filename);
372
373 OutputNodes(NodesFile, Textile, SCIRUN_EXPORT );
374
375 Filename = RemoveExtension( Filename, ".pts" );
376 Filename += ".hex";
377
378 ofstream OutputElements(Filename.c_str());
379
380 if (!OutputElements)
381 {
382 TGERROR("Unable to output voxel mesh to SCIRun file format, could not open file: " << Filename);
383 return;
384 }
385
386 TGLOG("Outputting hex elements");
387 //Output the voxel HEX elements
388 int iNumHexElements = 0;
389
390 iNumHexElements = OutputHexElements( OutputElements, true, true, SCIRUN_EXPORT );
391
392 TGLOG("Finished saving to SCIRun format");
393}
394
395
396
397
398
399int CVoxelMesh::OutputHexElements(ostream &Output, bool bOutputMatrix, bool bOutputYarn, int Filetype )
400{
401 int numx = m_XVoxels + 1;
402 int numy = m_YVoxels + 1;
403 int x,y,z;
404 vector<POINT_INFO>::iterator itElementInfo = m_ElementsInfo.begin();
405 int iElementNumber = 1;
406
407 vector<POINT_INFO> NewElementInfo;
408
409 if ( Filetype == SCIRUN_EXPORT )
410 Output << m_XVoxels*m_YVoxels*m_ZVoxels << "\n";
411
412 for ( z = 0; z < m_ZVoxels; ++z )
413 {
414 for ( y = 0; y < m_YVoxels; ++y )
415 {
416 for ( x = 0; x < m_XVoxels; ++x )
417 {
418 if ( (itElementInfo->iYarnIndex == -1 && bOutputMatrix)
419 || (itElementInfo->iYarnIndex >=0 && bOutputYarn) )
420 {
421 if ( Filetype == INP_EXPORT )
422 {
423 Output << iElementNumber << ", ";
424 Output << (x+1) +y*numx + z*numx*numy + 1 << ", " << (x+1) + (y+1)*numx + z*numx*numy + 1 << ", ";
425 Output << x + (y+1)*numx + z*numx*numy + 1 << ", " << x + y*numx + z*numx*numy + 1 << ", ";
426 Output << (x+1) +y*numx + (z+1)*numx*numy + 1 << ", " << (x+1) +(y+1)*numx + (z+1)*numx*numy + 1 << ", ";
427 Output << x +(y+1)*numx + (z+1)*numx*numy + 1 << ", " << x +y*numx + (z+1)*numx*numy + 1 << "\n";
428 }
429 else if ( Filetype == SCIRUN_EXPORT)
430 {
431 Output << x +y*numx + z*numx*numy + 1 << ", " << (x+1) + y*numx + z*numx*numy + 1 << ", ";
432 Output << x + y*numx + (z+1)*numx*numy + 1 << ", " << (x+1) + y*numx + (z+1)*numx*numy + 1 << ", ";
433 Output << x + (y+1)*numx + z*numx*numy + 1 << ", " << (x+1) +(y+1)*numx + z*numx*numy + 1 << ", ";
434 Output << x +(y+1)*numx + (z+1)*numx*numy + 1 << ", " << (x+1) + (y+1)*numx + (z+1)*numx*numy + 1 << "\n";
435 }
436 else // VTU export
437 {
438 vector<int> Indices;
439 Indices.push_back(x + y*numx + z*numx*numy);
440 Indices.push_back((x + 1) + y*numx + z*numx*numy);
441 Indices.push_back((x + 1) + y*numx + (z + 1)*numx*numy);
442 Indices.push_back(x + y*numx + (z + 1)*numx*numy);
443 Indices.push_back(x + (y + 1)*numx + z*numx*numy);
444 Indices.push_back((x + 1) + (y + 1)*numx + z*numx*numy);
445 Indices.push_back((x + 1) + (y + 1)*numx + (z + 1)*numx*numy);
446 Indices.push_back(x + (y + 1)*numx + (z + 1)*numx*numy);
447 m_Mesh.AddElement(CMesh::HEX, Indices);
448 }
449 ++iElementNumber;
450 if ( bOutputYarn && !bOutputMatrix ) // Just saving yarn so need to make element array with just yarn info
451 {
452 NewElementInfo.push_back( *itElementInfo );
453 }
454
455 }
456 ++itElementInfo;
457 }
458 }
459 }
460
461
462 if ( bOutputYarn && !bOutputMatrix )
463 {
464 m_ElementsInfo.clear();
465 m_ElementsInfo = NewElementInfo;
466 }
467 return ( iElementNumber-1 );
468}
469
471{
472 AddExtensionIfMissing(Filename, ".inp");
473
474 ofstream Output(Filename.c_str(), ofstream::app);
475 OutputOrientationsAndElementSets( Filename, Output );
476}
477
478void CVoxelMesh::OutputOrientationsAndElementSets( string Filename, ostream &Output )
479{
480 string OrientationsFilename = Filename;
481 OrientationsFilename.replace(OrientationsFilename.end()-4, OrientationsFilename.end(), ".ori");
482 ofstream OriOutput(OrientationsFilename.c_str());
483 string ElementDataFilename = Filename;
484 ElementDataFilename.replace(ElementDataFilename.end()-4, ElementDataFilename.end(), ".eld");
485 ofstream DataOutput(ElementDataFilename.c_str());
486
487 if (!OriOutput)
488 {
489 TGERROR("Unable to output orientations, could not open file: " << OrientationsFilename);
490 return;
491 }
492 if (!DataOutput)
493 {
494 TGERROR("Unable to output additional element data, could not open file: " << ElementDataFilename);
495 return;
496 }
497
498 TGLOG("Saving element orientations data to " << OrientationsFilename);
499 TGLOG("Saving additional element data to " << ElementDataFilename);
500
501 WriteOrientationsHeader( Output );
502 Output << "*Distribution Table, Name=TexGenOrientationVectors" << "\n";
503 Output << "COORD3D,COORD3D" << "\n";
504 Output << "*Distribution, Location=Element, Table=TexGenOrientationVectors, Name=TexGenOrientationVectors, Input=" << StripPath(OrientationsFilename) << "\n";
505 Output << "*Orientation, Name=TexGenOrientations, Definition=coordinates" << "\n";
506 Output << "TexGenOrientationVectors" << "\n";
507 Output << "1, 0" << "\n";
508
509 // Default orientation
510 WriteOrientationsHeader( OriOutput );
511 OriOutput << ", 1.0, 0.0, 0.0, 0.0, 1.0, 0.0" << "\n";
512
513 WriteElementsHeader( DataOutput );
514
515 int i;
516
517 map<int, vector<int> > ElementSets;
518 vector<POINT_INFO>::iterator itData;
519 for (itData = m_ElementsInfo.begin(), i=1; itData != m_ElementsInfo.end(); ++itData, ++i)
520 {
521 if (itData->iYarnIndex != -1)
522 {
523 XYZ Up = itData->Up;
524 XYZ Dir = itData->Orientation;
525
526 XYZ Perp = CrossProduct(Dir, Up);
527 Normalise(Perp);
528 OriOutput << i << ", " << Dir << ", " << Perp << "\n";
529 }
530 else
531 {
532 // Default orientation
533 OriOutput << i << ", 1.0, 0.0, 0.0, 0.0, 1.0, 0.0" << "\n";
534 }
535 DataOutput << i;
536 DataOutput << ", " << itData->iYarnIndex;
537 DataOutput << ", " << itData->Location; // This counts as 2 DepVars
538 DataOutput << ", " << itData->dVolumeFraction;
539 DataOutput << ", " << itData->dSurfaceDistance;
540 DataOutput << "\n";
541 ElementSets[itData->iYarnIndex].push_back(i);
542 }
543
544 // Output element sets
545 Output << "********************" << "\n";
546 Output << "*** ELEMENT SETS ***" << "\n";
547 Output << "********************" << "\n";
548 Output << "** TexGen generates a number of element sets:" << "\n";
549 Output << "** All - Contains all elements" << "\n";
550 Output << "** Matrix - Contains all elements belonging to the matrix" << "\n";
551 Output << "** YarnX - Where X represents the yarn index" << "\n";
552 Output << "*ElSet, ElSet=AllElements, Generate" << "\n";
553 Output << "1, " << m_ElementsInfo.size() << ", 1" << "\n";
554 vector<int>::iterator itIndices;
555 map<int, vector<int> >::iterator itElementSet;
556 for (itElementSet = ElementSets.begin(); itElementSet != ElementSets.end(); ++itElementSet)
557 {
558 if (itElementSet->first == -1)
559 Output << "*ElSet, ElSet=Matrix" << "\n";
560 else
561 Output << "*ElSet, ElSet=Yarn" << itElementSet->first << "\n";
562
563 WriteValues(Output, itElementSet->second, 16);
564 }
565}
566
567void CVoxelMesh::OutputMatrixElementSet( string Filename, ostream &Output, int iNumHexElements, bool bMatrixOnly )
568{
569 if ( !bMatrixOnly )
570 return;
571
572 // Output element sets
573 Output << "********************" << "\n";
574 Output << "*** ELEMENT SETS ***" << "\n";
575 Output << "********************" << "\n";
576 Output << "** TexGen generates a number of element sets:" << "\n";
577 Output << "** All - Contains all elements" << "\n";
578 Output << "** Matrix - Contains all elements belonging to the matrix" << "\n";
579 Output << "** YarnX - Where X represents the yarn index" << "\n";
580 Output << "*ElSet, ElSet=Matrix, Generate" << "\n";
581 Output << "1, " << iNumHexElements << ", 1" << "\n";
582}
583
584void CVoxelMesh::OutputAllNodesSet( string Filename, ostream &Output )
585{
586 Output << "*****************" << "\n";
587 Output << "*** NODE SETS ***" << "\n";
588 Output << "*****************" << "\n";
589 Output << "** AllNodes - Node set containing all elements" << "\n";
590 Output << "*NSet, NSet=AllNodes, Generate" << "\n";
591 Output << "1, " << (m_XVoxels+1)*(m_YVoxels+1)*(m_ZVoxels+1) << ", 1" << "\n";
592}
593
594void CVoxelMesh::OutputPeriodicBoundaries(ostream &Output, CTextile& Textile, int iBoundaryConditions, bool bMatrixOnly )
595{
597
598 vector<int> GroupA;
599 vector<int> GroupB;
600 pair< vector<int>, vector<int> > Faces;
601
602 int numx = m_XVoxels + 1;
603 int numy = m_YVoxels + 1;
604 int numz = m_ZVoxels + 1;
605 //int iDummyNodeNum = numx*numy*numz + 1;
606
607 // Create a set of nodes for opposite faces of the domain, then output
608 for ( int z = 1; z < numz-1; ++z )
609 {
610 for ( int y = 1; y < numy-1; ++y )
611 {
612 GroupA.push_back( numx-1 + y*numx + z*numx*numy +1 );
613 GroupB.push_back( y*numx + z*numx*numy + 1);
614
615 }
616 }
617 m_PeriodicBoundaries->SetFaceA( GroupA, GroupB );
618
619 GroupA.clear();
620 GroupB.clear();
621 for ( int z = 1; z < numz-1; ++z )
622 {
623 for ( int x = 1; x < numx-1; ++x )
624 {
625 GroupA.push_back( x + (numy-1)*numx + z*numx*numy +1 );
626 GroupB.push_back( x + z*numx*numy + 1);
627
628 }
629 }
630 m_PeriodicBoundaries->SetFaceB( GroupA, GroupB );
631
632 GroupA.clear();
633 GroupB.clear();
634 for ( int y = 1; y < numy-1; ++y )
635 {
636 for ( int x = 1; x < numx-1; ++x )
637 {
638 GroupA.push_back( x + y*numx + (numz-1)*numx*numy + 1);
639 GroupB.push_back( x + y*numx + 1 );
640
641 }
642 }
643 m_PeriodicBoundaries->SetFaceC( GroupA, GroupB );
644
645 // Create edge node sets
646 vector<int> Edges[4];
647
648 for ( int z = 1; z < numz-1; ++z )
649 {
650 Edges[0].push_back( z*numx*numy + 1 );
651 Edges[3].push_back( numx*(numy-1) + z*numx*numy +1 );
652 Edges[2].push_back( (z+1)*numx*numy );
653 Edges[1].push_back( numx + z*numx*numy );
654 }
655
656 for ( int i = 0; i < 4; ++i )
657 {
658 m_PeriodicBoundaries->SetEdges( Edges[i] );
659 Edges[i].clear();
660 }
661
662 for ( int y = 1; y < numy-1; ++y )
663 {
664 Edges[3].push_back( y*numx + (numz-1)*numx*numy + 1 );
665 Edges[2].push_back( (y+1)*numx +(numz-1)*numx*numy );
666 Edges[1].push_back( (y+1)*numx );
667 Edges[0].push_back( y*numx + 1 );
668 }
669
670 for ( int i = 0; i < 4; ++i )
671 {
672 m_PeriodicBoundaries->SetEdges( Edges[i] );
673 Edges[i].clear();
674 }
675
676 for ( int x = 1; x < numx-1; ++x )
677 {
678 Edges[3].push_back( x + (numz-1)*numx*numy + 1 );
679 Edges[2].push_back( x + numx*(numy-1) + (numz-1)*numx*numy +1 );
680 Edges[1].push_back( x + numx*(numy-1) + 1 );
681 Edges[0].push_back( x + 1 );
682 }
683
684 for ( int i = 0; i < 4; ++i )
685 {
686 m_PeriodicBoundaries->SetEdges( Edges[i] );
687 Edges[i].clear();
688 }
689
690 // Create vertex sets
691 // Vertices rearranged to tie up with Li notation (previous numbering in comments)
694 m_PeriodicBoundaries->SetVertex( numx*numy );//7
695 m_PeriodicBoundaries->SetVertex( (numy-1)*numx + 1 );//6
696 m_PeriodicBoundaries->SetVertex( (numz-1)*numx*numy + 1 );//1
697 m_PeriodicBoundaries->SetVertex( numx + (numz-1)*numx*numy );//4
698 m_PeriodicBoundaries->SetVertex( numx*numy*numz );//3
699 m_PeriodicBoundaries->SetVertex( numx*(numy-1) + (numz-1)*numx*numy + 1 );//2
700
701 m_PeriodicBoundaries->CreatePeriodicBoundaries( Output, numx*numy*numz + 1, Textile, iBoundaryConditions, bMatrixOnly );
702}
703
704void CVoxelMesh::AddElementInfo(vector<POINT_INFO> &RowInfo)
705{
706 m_ElementsInfo.insert(m_ElementsInfo.end(), RowInfo.begin(), RowInfo.end() );
707}
708
710/*template <typename T>
711void CVoxelMesh::WriteValues( ostream &Output, T &Values, int iMaxPerLine)
712{
713 int iLinePos = 0;
714 typename T::const_iterator itValue;
715 for (itValue = Values.begin(); itValue != Values.end(); ++itValue)
716 {
717 if (iLinePos == 0)
718 {
719 // Do nothing...
720 }
721 else if (iLinePos < iMaxPerLine)
722 {
723 Output << ", ";
724 }
725 else
726 {
727 Output << endl;
728 iLinePos = 0;
729 }
730 Output << *itValue;
731 ++iLinePos;
732 }
733 Output << endl;
734}*/
735
736/*void CVoxelMesh::GetYarnGridIntersections( CTextile &Textile )
737{
738 int x,y,z;
739 int iNodeIndex = 1;
740 vector<XYZ> CentrePoints;
741 vector<POINT_INFO> RowInfo;
742
743 vector<pair<XYZ, XYZ> > xyLines;
744 vector<pair<XYZ, XYZ> > xzLines;
745 vector<pair<XYZ, XYZ> > yzLines;
746
747 // Create lines in z direction to test for intersections
748 for ( y = 0; y <= m_YVoxels; ++y )
749 {
750 for ( x = 0; x <=m_XVoxels; ++x )
751 {
752 XYZ Point1, Point2;
753 Point1.x = Point2.x = m_DomainAABB.first.x + m_VoxSize[0] * x;
754 Point1.y = Point2.y = m_DomainAABB.first.y + m_VoxSize[1] * y;
755 Point1.z = 0.0;
756 Point2.z = 1.0;
757 xyLines.push_back( make_pair(Point1,Point2) );
758 }
759 }
760
761 // Create lines in x direction to test for intersections
762 for ( z = 0; z <= m_ZVoxels; ++z )
763 {
764 for ( y = 0; y <=m_YVoxels; ++y )
765 {
766 XYZ Point1, Point2;
767 Point1.x = 0.0;
768 Point2.x = 1.0;
769 Point1.y = Point2.y = m_DomainAABB.first.y + m_VoxSize[1] * y;
770 Point1.z = Point2.z = m_DomainAABB.first.z + m_VoxSize[2] * z;
771 yzLines.push_back( make_pair(Point1,Point2) );
772 }
773 }
774
775 // Create lines in y direction to test for intersections
776 for ( z = 0; z <= m_ZVoxels; ++z )
777 {
778 for ( x = 0; x <=m_XVoxels; ++x )
779 {
780 XYZ Point1, Point2;
781 Point1.x = Point2.x = m_DomainAABB.first.x + m_VoxSize[0] * x;
782 Point1.y = 0.0;
783 Point2.y = 1.0;
784 Point1.z = Point2.z = m_DomainAABB.first.z + m_VoxSize[2] * z;
785 xzLines.push_back( make_pair(Point1,Point2) );
786 }
787 }
788
789 vector<CYarn> &Yarns = Textile.GetYarns();
790 vector<CYarn>::iterator itYarn;
791 const CDomain* pDomain = Textile.GetDomain();
792
793 //vector< pair<int, pair<double,double>>> xyIntersections;
794 vector<vector<pair<int, double> > > xyIntersections;
795 vector<vector<pair<int, double> > > xzIntersections;
796 vector<vector<pair<int, double> > > yzIntersections;
797 xyIntersections.clear();
798 xyIntersections.resize(xyLines.size());
799 xzIntersections.clear();
800 xzIntersections.resize(xzLines.size());
801 yzIntersections.clear();
802 yzIntersections.resize(yzLines.size());
803
804 int i = 0;
805 for ( itYarn = Yarns.begin(); itYarn != Yarns.end(); ++itYarn, ++i)
806 {
807 CMesh Mesh;
808 vector<pair<XYZ, XYZ> >::iterator itLines;
809
810 itYarn->AddSurfaceToMesh(Mesh, *pDomain);
811 Mesh.ConvertQuadstoTriangles();
812
813 CalculateIntersections( Mesh, xyLines, i, xyIntersections );
814 CalculateIntersections( Mesh, xzLines, i, xzIntersections );
815 CalculateIntersections( Mesh, yzLines, i, yzIntersections );
816 }
817}
818
819void CVoxelMesh::CalculateIntersections( CMesh &Mesh, vector<pair<XYZ,XYZ> > &Lines, int YarnIndex, vector<vector<pair<int, double> > > &Intersections )
820{
821 vector< pair<double, XYZ> > IntersectionPoints;
822 vector< pair<double, XYZ> >::iterator itIntersectionPoints;
823 vector<pair<XYZ, XYZ>>::iterator itLines;
824 int j = 0;
825 double dTol = 1e-8;
826
827 for( itLines = Lines.begin(); itLines != Lines.end(); ++itLines, ++j )
828 {
829 int numIntersections;
830 IntersectionPoints.clear();
831 numIntersections = Mesh.IntersectLine( itLines->first, itLines->second, IntersectionPoints );
832
833 pair<int,double> Intersection;
834 Intersection.first = YarnIndex;
835 double PrevDist = 0.0;
836 for( itIntersectionPoints = IntersectionPoints.begin(); itIntersectionPoints != IntersectionPoints.end(); ++itIntersectionPoints )
837 {
838 if ( itIntersectionPoints != IntersectionPoints.begin() && fabs(PrevDist - itIntersectionPoints->first) < dTol )
839 {
840 continue;
841 }
842 Intersection.second = itIntersectionPoints->first;
843 Intersections[j].push_back(Intersection);
844 PrevDist = itIntersectionPoints->first;
845 }
846 }
847}*/
#define TGERROR(MESSAGE)
Macros used to report the file name and line number to the TexGenError and TexGenLog functions.
Definition: Logger.h:29
#define TGLOG(MESSAGE)
Definition: Logger.h:36
#define TEXGEN
Helper macro to get the texgen instance.
Definition: TexGen.h:76
Class used to generate Abaqus input file for periodic boundary conditions for a unit cell in tension ...
Abstract base class representing the domain in which a textile cell may lie.
Definition: Domain.h:34
const CMesh & GetMesh() const
Get the mesh representing the domain as a surface mesh.
Definition: Domain.h:73
vector< T > m_Data
Definition: MeshData.h:87
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
bool SaveToVTK(string Filename, const vector< CMeshDataBase * > *pMeshData=NULL) const
Save the mesh to VTK unstructured grid file format (.vtu)
Definition: Mesh.cpp:2348
Class used to generate Abaqus output for periodic boundary conditions.
void SetFaceB(vector< int > &B1, vector< int > &B2)
void CreatePeriodicBoundaries(ostream &Output, int iDummyNodeNum, CTextile &Textile, int iBoundarConditions, bool bMatrixOnly)
void SetEdges(vector< int > &Edge)
virtual void SetDomainSize(const CMesh &Mesh)
void SetFaceA(vector< int > &A1, vector< int > &A2)
void SetFaceC(vector< int > &C1, vector< int > &C2)
Class used to generate Abaqus input file for periodic boundary conditions for a textile with rotated ...
Class used to generate Abaqus input file for periodic boundary conditions for a textile with sheared ...
Class used to generate Abaqus input file for periodic boundary conditions for a textile with staggere...
Represents a textile cell containing yarns.
Definition: Textile.h:39
const CDomain * GetDomain() const
Definition: Textile.h:287
int GetNumYarns() const
Definition: Textile.cpp:704
void OutputMaterials(ostream &Output, int iNumYarns, bool bMatrixOnly)
Output materials and assign to yarn element sets.
Definition: Materials.cpp:187
void SetupMaterials(CTextile &Textile)
Set up a material for each unique set of material constants.
Definition: Materials.cpp:31
Class used to meaure the amount of time it takes to perform a certain task.
Definition: Timer.h:12
void stop(const char *msg=0)
Stop the timer and print an optional message.
Definition: Timer.h:86
void check(const char *msg=0)
Definition: Timer.h:96
void start(const char *msg=0)
Definition: Timer.h:57
int m_XVoxels
Number of voxels along x,y and z axes.
Definition: VoxelMesh.h:115
void SaveToSCIRun(string Filename, CTextile &Textile)
Save voxel mesh in SCIRun .pts and .hex format without boundary conditions.
Definition: VoxelMesh.cpp:358
virtual void OutputNodes(ostream &Output, CTextile &Textile, int Filetype=INP_EXPORT)=0
Outputs nodes to .inp file and gets element information.
CTextileMaterials m_Materials
Definition: VoxelMesh.h:60
virtual void OutputPeriodicBoundaries(ostream &Output, CTextile &Textile, int iBoundaryConditions, bool bMatrixOnly)
Output periodic boundary conditions to .inp file.
Definition: VoxelMesh.cpp:594
virtual bool CalculateVoxelSizes(CTextile &Textile)=0
Calculate voxel size based on number of voxels on each axis and domain size.
virtual int OutputHexElements(ostream &Output, bool bOutputMatrix, bool bOutputYarn, int Filetype=INP_EXPORT)
Output hex elements to .inp file.
Definition: VoxelMesh.cpp:399
CPeriodicBoundaries * m_PeriodicBoundaries
Definition: VoxelMesh.h:126
void SaveToAbaqus(string Filename, CTextile &Textile, bool bOutputMatrix, bool bOutputYarn, int iBoundaryConditions, int iElementType)
Definition: VoxelMesh.cpp:287
void OutputMatrixElementSet(string Filename, ostream &Output, int iNumHexElements, bool bMatrixOnly)
Outputs all elements when only outputting matrix.
Definition: VoxelMesh.cpp:567
void OutputOrientationsAndElementSets(string Filename)
Outputs yarn orientations and element sets to .ori and .eld files.
Definition: VoxelMesh.cpp:470
void SaveVoxelMeshToVTK(string Filename, CTextile &Textile)
Add nodes for whole domain.
Definition: VoxelMesh.cpp:246
vector< POINT_INFO > m_ElementsInfo
Element information as calculated by GetPointInformation.
Definition: VoxelMesh.h:123
void AddElementInfo(vector< POINT_INFO > &RowInfo)
Add a row of element information.
Definition: VoxelMesh.cpp:704
void OutputAllNodesSet(string Filename, ostream &Output)
Output node set containing all nodes.
Definition: VoxelMesh.cpp:584
CMesh m_Mesh
Find intersections of yarn surfaces with grid of lines from node points in each axis.
Definition: VoxelMesh.h:112
virtual void SaveVoxelMesh(CTextile &Textile, string OutputFilename, int XVoxNum, int YVoxNum, int ZVoxNum, bool bOutputMatrix, bool bOutputYarns, int iBoundaryConditions, int iElementType=0, int FileType=INP_EXPORT)
Definition: VoxelMesh.cpp:51
virtual ~CVoxelMesh(void)
Definition: VoxelMesh.cpp:46
Namespace containing a series of customised math operations not found in the standard c++ library.
std::string RemoveExtension(std::string &Filename, std::string Extension)
Strip the extension from the filename.
Definition: Misc.cpp:122
std::string StripPath(std::string Filename)
Strip the path from the filename (e.g. "c:\folder\file.ext -> file.ext")
Definition: Misc.cpp:107
std::string ReplaceFilenameSpaces(std::string Filename)
Replaces spaces in filename with underscore.
Definition: Misc.cpp:130
void WriteElementsHeader(std::ostream &Output)
Write elements header for ABAQUS .eld files.
Definition: Misc.cpp:167
void WriteOrientationsHeader(std::ostream &Output)
Write orientations header for ABAQUS .ori files.
Definition: Misc.cpp:157
void WriteValues(std::ostream &Output, T &Values, int iMaxPerLine)
Definition: Misc.h:274
@ SCIRUN_EXPORT
Definition: Misc.h:115
@ INP_EXPORT
Definition: Misc.h:113
@ VTU_EXPORT
Definition: Misc.h:114
WXYZ & Normalise(WXYZ &Quaternion)
Normalise the quaternion and return it.
Definition: mymath.h:609
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
XYZ CrossProduct(const XYZ &left, const XYZ &right)
Get the cross product of two vectors.
Definition: mymath.h:524
@ NO_BOUNDARY_CONDITIONS
Definition: Misc.h:270
Struct for representing points in 3D space.
Definition: mymath.h:56