Line data Source code
1 : /*
2 : * Copyright (c) 2009-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #include <fstream>
34 : #include "file_io_tetgen.h"
35 : #include "common/util/string_util.h"
36 : #include "../lg_base.h"
37 :
38 : using namespace std;
39 :
40 : namespace ug
41 : {
42 :
43 : ////////////////////////////////////////////////////////////////////////
44 0 : bool LoadGridFromELE(Grid& grid, const char* filename, ISubsetHandler* pSH,
45 : APosition& aPos)
46 : {
47 : // build the correct filenames
48 0 : string sElements = filename;
49 : string sNodes, sFaces;
50 0 : if(sElements.find(".ele", sElements.size() - 4) != string::npos)
51 : {
52 : // the filename ends as we would expect
53 0 : sNodes = sFaces = sElements.substr(0, sElements.size() - 4);
54 0 : sNodes.append(".node");
55 0 : sFaces.append(".face");
56 : }
57 : else
58 : {
59 0 : LOG("Problem in LoadGridFromELE with file " << filename << ". filename has to end with .ele\n");
60 : return false;
61 : }
62 :
63 0 : return ImportGridFromTETGEN(grid, sNodes.c_str(), sFaces.c_str(), sElements.c_str(),
64 : aPos, pSH);
65 : }
66 :
67 : ////////////////////////////////////////////////////////////////////////
68 0 : bool SaveGridToELE(Grid& grid, const char* filename, ISubsetHandler* pSH,
69 : APosition& aPos, ANumber* paVolumeConstraint)
70 : {
71 : // give a warning if the grid doesn't consist of tetrahedrons only.
72 0 : if(grid.num<Tetrahedron>() < grid.num<Volume>()){
73 : UG_LOG(" INFO in SaveGridToELE: Non-tetrahedral elements will be skipped.\n");
74 : }
75 :
76 : // if a subset-handler was supplied we have to copy its subset-indices
77 : // to a simple int-attachment.
78 0 : if(pSH)
79 : {
80 : AInt aInt;
81 : grid.attach_to_volumes(aInt);
82 : Grid::VolumeAttachmentAccessor<AInt> aaInt(grid, aInt);
83 :
84 0 : for(VolumeIterator iter = grid.begin<Volume>(); iter != grid.end<Volume>(); ++iter)
85 0 : aaInt[*iter] = pSH->get_subset_index(*iter);
86 :
87 : // face subset-indices
88 : grid.attach_to_faces(aInt);
89 : Grid::FaceAttachmentAccessor<AInt> aaIntFace(grid, aInt);
90 0 : for(FaceIterator iter = grid.begin<Face>(); iter != grid.end<Face>(); ++iter)
91 0 : aaIntFace[*iter] = pSH->get_subset_index(*iter);
92 :
93 0 : bool retVal = ExportGridToTETGEN(grid, filename, aPos, NULL, NULL, &aInt, &aInt,
94 : paVolumeConstraint);
95 :
96 : grid.detach_from_faces(aInt);
97 : grid.detach_from_volumes(aInt);
98 :
99 : return retVal;
100 : }
101 :
102 0 : return ExportGridToTETGEN(grid, filename, aPos,
103 0 : NULL, NULL, NULL, NULL, paVolumeConstraint);
104 :
105 : }
106 :
107 : ////////////////////////////////////////////////////////////////////////
108 : // ImportGridFromTETGEN
109 0 : bool ImportGridFromTETGEN(Grid& grid,
110 : const char* nodesFilename, const char* facesFilename,
111 : const char* elemsFilename, AVector3& aPos,
112 : std::vector<AFloat>* pvNodeAttributes,
113 : AInt* paNodeBoundaryMarker,
114 : AInt* paFaceBoundaryMarker,
115 : AInt* paElementAttribute)
116 : {
117 : // read nodes and store them in an array for index access
118 : vector<RegularVertex*> vVertices;
119 :
120 : {
121 0 : ifstream in(nodesFilename);
122 0 : if(!in)
123 : {
124 0 : LOG("WARNING in ImportGridFromTETGEN: nodes file not found: " << nodesFilename << endl);
125 : return false;
126 : }
127 :
128 : uint numNodes, dim, numAttribs, numBoundaryMarkers;
129 : in >> numNodes;
130 : in >> dim;
131 : in >> numAttribs;
132 : in >> numBoundaryMarkers;
133 :
134 0 : vVertices.reserve(numNodes + 1);
135 :
136 : // set up attachment accessors
137 0 : if(!grid.has_vertex_attachment(aPos))
138 0 : grid.attach_to_vertices(aPos);
139 : Grid::VertexAttachmentAccessor<AVector3> aaPosVRT(grid, aPos);
140 :
141 : Grid::VertexAttachmentAccessor<AInt> aaBMVRT;
142 0 : if(paNodeBoundaryMarker != NULL)
143 : aaBMVRT.access(grid, *paNodeBoundaryMarker);
144 :
145 : vector<Grid::VertexAttachmentAccessor<AFloat> > vaaAttributesVRT;
146 0 : if(pvNodeAttributes != NULL)
147 : {
148 0 : vaaAttributesVRT.resize(pvNodeAttributes->size());
149 0 : for(uint i = 0; i < pvNodeAttributes->size(); ++i)
150 : vaaAttributesVRT[i].access(grid, (*pvNodeAttributes)[i]);
151 : }
152 :
153 : // read the vertices
154 : int index;
155 0 : for(uint i = 0; i < numNodes; ++i)
156 : {
157 0 : RegularVertex* v = *grid.create<RegularVertex>();
158 0 : vVertices.push_back(v);
159 :
160 : // read index and coords
161 0 : in >> index;
162 0 : if(index > (int) vVertices.size())
163 0 : vVertices.resize(index, NULL);
164 : in >> aaPosVRT[v].x();
165 : in >> aaPosVRT[v].y();
166 : in >> aaPosVRT[v].z();
167 :
168 : // read attributes
169 0 : if(numAttribs > 0)
170 : {
171 0 : for(uint j = 0; j < numAttribs; ++j)
172 : {
173 : float tmp;
174 : in >> tmp;
175 0 : if(j < vaaAttributesVRT.size())
176 0 : (vaaAttributesVRT[j])[v] = tmp;
177 : }
178 : }
179 :
180 : // read boundary marker
181 0 : if(numBoundaryMarkers > 0)
182 : {
183 : int bm;
184 0 : in >> bm;
185 0 : if(paNodeBoundaryMarker != NULL)
186 0 : aaBMVRT[v] = bm;
187 : }
188 : }
189 :
190 0 : in.close();
191 0 : }
192 :
193 : // read faces
194 0 : if(facesFilename != NULL)
195 : {
196 0 : ifstream in(facesFilename);
197 0 : if(in)
198 : {
199 : int numFaces, numBoundaryMarkers;
200 0 : in >> numFaces;
201 0 : in >> numBoundaryMarkers;
202 :
203 : Grid::FaceAttachmentAccessor<AInt> aaBMFACE;
204 0 : if(paFaceBoundaryMarker != NULL)
205 : aaBMFACE.access(grid, *paFaceBoundaryMarker);
206 :
207 0 : for(int i = 0; i < numFaces; ++i)
208 : {
209 : int index, i1, i2, i3;
210 0 : in >> index;
211 0 : in >> i1;
212 0 : in >> i2;
213 0 : in >> i3;
214 :
215 0 : Triangle* t = *grid.create<Triangle>(TriangleDescriptor(vVertices[i1], vVertices[i2], vVertices[i3]));
216 :
217 :
218 :
219 0 : if(numBoundaryMarkers > 0)
220 : {
221 : int bm;
222 0 : in >> bm;
223 0 : if(aaBMFACE.valid())
224 0 : aaBMFACE[t] = bm;
225 : }
226 : else
227 : {
228 0 : if(aaBMFACE.valid())
229 0 : aaBMFACE[t] = -1;
230 : }
231 : }
232 : }
233 : else
234 0 : LOG("WARNING in ImportGridFromTETGEN: faces file not found: " << facesFilename << endl);
235 :
236 0 : in.close();
237 0 : }
238 :
239 : // read volumes
240 0 : if(elemsFilename != NULL)
241 : {
242 0 : ifstream in(elemsFilename);
243 0 : if(in)
244 : {
245 : int numTets, numNodesPerTet, numAttribs;
246 0 : in >> numTets;
247 0 : in >> numNodesPerTet;
248 0 : in >> numAttribs;
249 :
250 : // attachment accessors:
251 : Grid::VolumeAttachmentAccessor<AInt> aaAttributeVOL;
252 0 : if(paElementAttribute)
253 : aaAttributeVOL.access(grid, *paElementAttribute);
254 :
255 0 : vector<int> vTetNodes(numNodesPerTet);
256 0 : for(int i = 0; i < numTets; ++i)
257 : {
258 : int index;
259 0 : in >> index;
260 0 : for(int j = 0; j < numNodesPerTet; ++j)
261 0 : in >> vTetNodes[j];
262 :
263 0 : Tetrahedron* t = *grid.create<Tetrahedron>(TetrahedronDescriptor(
264 0 : vVertices[vTetNodes[0]], vVertices[vTetNodes[1]],
265 0 : vVertices[vTetNodes[2]], vVertices[vTetNodes[3]]));
266 :
267 0 : if(numAttribs > 0)
268 : {
269 : int a;
270 0 : in >> a;
271 0 : if(paElementAttribute != NULL)
272 0 : aaAttributeVOL[t] = a;
273 : }
274 : }
275 0 : }
276 : else
277 0 : LOG("WARNING in ImportGridFromTETGEN: elems file not found: " << elemsFilename << endl);
278 :
279 0 : in.close();
280 0 : }
281 : return true;
282 0 : }
283 :
284 : ////////////////////////////////////////////////////////////////////////
285 : // ImportGridFromTETGEN
286 0 : bool ImportGridFromTETGEN(Grid& grid,
287 : const char* nodesFilename, const char* facesFilename,
288 : const char* elemsFilename, AVector3& aPos,
289 : ISubsetHandler* psh,
290 : std::vector<AFloat>* pvNodeAttributes)
291 : {
292 : PROFILE_FUNC();
293 : // read nodes and store them in an array for index access
294 : vector<RegularVertex*> vVertices;
295 :
296 : {
297 : PROFILE_BEGIN(read_vertices);
298 0 : ifstream in(nodesFilename);
299 0 : if(!in)
300 : {
301 0 : LOG("WARNING in ImportGridFromTETGEN: nodes file not found: " << nodesFilename << endl);
302 : return false;
303 : }
304 :
305 : uint numNodes, dim, numAttribs, numBoundaryMarkers;
306 : in >> numNodes;
307 : in >> dim;
308 : in >> numAttribs;
309 : in >> numBoundaryMarkers;
310 :
311 0 : vVertices.reserve(numNodes + 1);
312 0 : grid.reserve<Vertex>(numNodes);
313 :
314 : // set up attachment accessors
315 0 : if(!grid.has_vertex_attachment(aPos))
316 0 : grid.attach_to_vertices(aPos);
317 : Grid::VertexAttachmentAccessor<AVector3> aaPosVRT(grid, aPos);
318 :
319 : vector<Grid::VertexAttachmentAccessor<AFloat> > vaaAttributesVRT;
320 0 : if(pvNodeAttributes != NULL)
321 : {
322 0 : vaaAttributesVRT.resize(pvNodeAttributes->size());
323 0 : for(uint i = 0; i < pvNodeAttributes->size(); ++i)
324 : vaaAttributesVRT[i].access(grid, (*pvNodeAttributes)[i]);
325 : }
326 :
327 : // read the vertices
328 : int index;
329 0 : for(uint i = 0; i < numNodes; ++i)
330 : {
331 0 : RegularVertex* v = *grid.create<RegularVertex>();
332 :
333 : // read index and coords
334 0 : in >> index;
335 0 : if(index > (int) vVertices.size())
336 0 : vVertices.resize(index, NULL);
337 0 : vVertices.push_back(v);
338 : in >> aaPosVRT[v].x();
339 : in >> aaPosVRT[v].y();
340 : in >> aaPosVRT[v].z();
341 :
342 : // read attributes
343 0 : if(numAttribs > 0)
344 : {
345 0 : for(uint j = 0; j < numAttribs; ++j)
346 : {
347 : float tmp;
348 : in >> tmp;
349 0 : if(j < vaaAttributesVRT.size())
350 0 : (vaaAttributesVRT[j])[v] = tmp;
351 : }
352 : }
353 :
354 : // read boundary marker
355 0 : if(numBoundaryMarkers > 0)
356 : {
357 : int bm;
358 0 : in >> bm;
359 0 : if(psh != NULL)
360 0 : psh->assign_subset(v, abs(bm));
361 : }
362 : }
363 :
364 0 : in.close();
365 0 : }
366 :
367 : // read faces
368 0 : if(facesFilename != NULL)
369 : {
370 : PROFILE_BEGIN(read_faces);
371 0 : ifstream in(facesFilename);
372 0 : if(in)
373 : {
374 : int numFaces, numBoundaryMarkers;
375 0 : in >> numFaces;
376 0 : in >> numBoundaryMarkers;
377 :
378 0 : grid.reserve<Face>(numFaces);
379 :
380 0 : for(int i = 0; i < numFaces; ++i)
381 : {
382 : int index, i1, i2, i3;
383 0 : in >> index;
384 0 : in >> i1;
385 0 : in >> i2;
386 0 : in >> i3;
387 :
388 0 : Triangle* t = *grid.create<Triangle>(TriangleDescriptor(vVertices[i1], vVertices[i2], vVertices[i3]));
389 :
390 0 : if(numBoundaryMarkers > 0){
391 : int bm;
392 0 : in >> bm;
393 0 : if(psh != NULL){
394 0 : psh->assign_subset(t, abs(bm));
395 : }
396 : }
397 0 : else if(psh != NULL){
398 0 : psh->assign_subset(t, 0);
399 : }
400 : }
401 : }
402 : else
403 0 : LOG("WARNING in ImportGridFromTETGEN: faces file not found: " << facesFilename << endl);
404 :
405 0 : in.close();
406 0 : }
407 :
408 : // read volumes
409 0 : if(elemsFilename != NULL)
410 : {
411 : PROFILE_BEGIN(read_volumes);
412 0 : ifstream in(elemsFilename);
413 0 : if(in)
414 : {
415 : int numTets, numNodesPerTet, numAttribs;
416 0 : in >> numTets;
417 0 : in >> numNodesPerTet;
418 0 : in >> numAttribs;
419 :
420 0 : grid.reserve<Volume>(numTets);
421 :
422 0 : vector<int> vTetNodes(numNodesPerTet);
423 0 : for(int i = 0; i < numTets; ++i)
424 : {
425 : int index;
426 0 : in >> index;
427 : // UG_LOG(" index = " << index << std::endl);
428 0 : for(int j = 0; j < numNodesPerTet; ++j)
429 0 : in >> vTetNodes[j];
430 :
431 0 : Tetrahedron* t = *grid.create<Tetrahedron>(TetrahedronDescriptor(
432 0 : vVertices[vTetNodes[0]], vVertices[vTetNodes[1]],
433 0 : vVertices[vTetNodes[2]], vVertices[vTetNodes[3]]));
434 :
435 0 : int a = 0;
436 0 : if(numAttribs > 0)
437 0 : in >> a;
438 0 : if(psh != NULL)
439 0 : psh->assign_subset(t, a);
440 : }
441 0 : }
442 : else
443 0 : LOG("WARNING in ImportGridFromTETGEN: elems file not found: " << elemsFilename << endl);
444 :
445 0 : in.close();
446 0 : }
447 :
448 : return true;
449 0 : }
450 :
451 : ////////////////////////////////////////////////////////////////////////
452 : // ExportGridToSMESH
453 0 : bool ExportGridToSMESH(Grid& grid, const char* filename, AVector3& aPos,
454 : std::vector<AFloat>* pvNodeAttributes,
455 : AInt* paNodeBoundaryMarker,
456 : AInt* paFaceBoundaryMarker,
457 : std::vector<vector3>* pvHoles,
458 : std::vector<vector3>* pvRegionPositions,
459 : std::vector<int>* pvRegionAttributes,
460 : std::vector<float>* pvRegionVolumeConstraints)
461 :
462 : {
463 0 : if(!grid.has_vertex_attachment(aPos))
464 : return false;
465 :
466 0 : ofstream out(filename);
467 :
468 0 : if(!out)
469 : return false;
470 :
471 : // check if regions are specified in the correct way.
472 0 : if(pvRegionPositions != NULL)
473 : {
474 0 : if(pvRegionAttributes != NULL)
475 0 : if(pvRegionAttributes->size() != pvRegionPositions->size())
476 : return false;
477 :
478 0 : if(pvRegionVolumeConstraints != NULL)
479 0 : if(pvRegionVolumeConstraints->size() != pvRegionPositions->size())
480 : return false;
481 : }
482 :
483 : vector<int> vTmpRegionAttributes;
484 0 : if((pvRegionAttributes == NULL) && (pvRegionPositions != NULL) && (pvRegionVolumeConstraints != NULL))
485 : {
486 : // attributes have to be supplied if constraints are given.
487 : // since no attributes were passed, generate your own.
488 : pvRegionAttributes = &vTmpRegionAttributes;
489 0 : for(uint i = 0; i < pvRegionPositions->size(); ++i)
490 0 : pvRegionAttributes->push_back(0);
491 : }
492 :
493 : // write points. store indices with each point for later array-like access.
494 : AInt aInt;
495 : grid.attach_to_vertices(aInt);
496 : Grid::VertexAttachmentAccessor<AInt> aaIntVRT(grid, aInt);
497 : Grid::VertexAttachmentAccessor<AVector3> aaPos(grid, aPos);
498 :
499 : {
500 : int numAttribs = 0;
501 : int numBoundaryMarkers = 0;
502 :
503 : // attachment-accessors for the nodes attributes
504 : vector<Grid::VertexAttachmentAccessor<AFloat> > vaaFloatVRT;
505 0 : if(pvNodeAttributes != NULL)
506 : {
507 0 : numAttribs = pvNodeAttributes->size();
508 0 : vaaFloatVRT.resize(pvNodeAttributes->size());
509 0 : for(uint i = 0; i < pvNodeAttributes->size(); ++i)
510 : vaaFloatVRT[i].access(grid, (*pvNodeAttributes)[i]);
511 : }
512 :
513 : // attachment-accessor for the nodes boundary-marker
514 : Grid::VertexAttachmentAccessor<AInt> aaBMVRT;
515 0 : if(paNodeBoundaryMarker != NULL)
516 : {
517 : numBoundaryMarkers = 1;
518 : aaBMVRT.access(grid, *paNodeBoundaryMarker);
519 : }
520 :
521 : // write number of nodes, dimension, number of attributes, boundary markers (0 or 1)
522 0 : out << grid.num_vertices() << " 3 " << numAttribs << " " << numBoundaryMarkers << endl;
523 :
524 : int counter = 0;
525 0 : for(VertexIterator iter = grid.vertices_begin(); iter != grid.vertices_end(); iter++, counter++)
526 : {
527 0 : aaIntVRT[*iter] = counter;
528 0 : out << counter << " " << aaPos[*iter].x() << " " <<
529 0 : aaPos[*iter].y() << " " <<
530 0 : aaPos[*iter].z();
531 :
532 : // write attributes:
533 0 : for(uint i = 0; i < vaaFloatVRT.size(); ++i)
534 : {
535 0 : out << " " << (vaaFloatVRT[i])[*iter];
536 : }
537 :
538 : // write boundary markers:
539 0 : if(paNodeBoundaryMarker != NULL)
540 0 : out << " " << aaBMVRT[*iter];
541 :
542 : out << endl;
543 :
544 : }
545 0 : }
546 :
547 : out << endl;
548 :
549 : // write facets
550 : {
551 : Grid::FaceAttachmentAccessor<AInt> aaBMFACE;
552 0 : if(paFaceBoundaryMarker != NULL)
553 : {
554 : out << grid.num_faces() << " 1" << endl;
555 : aaBMFACE.access(grid, *paFaceBoundaryMarker);
556 : }
557 : else
558 : out << grid.num_faces() << " 0" << endl;
559 :
560 0 : for(FaceIterator iter = grid.faces_begin(); iter != grid.faces_end(); ++iter)
561 : {
562 : Face* f = *iter;
563 0 : out << f->num_vertices();
564 :
565 0 : for(uint i = 0; i < f->num_vertices(); ++i)
566 0 : out << " " << aaIntVRT[f->vertex(i)];
567 :
568 0 : if(paFaceBoundaryMarker != NULL)
569 0 : out << " " << aaBMFACE[f];
570 :
571 : out << endl;
572 : }
573 : }
574 :
575 : out << endl;
576 :
577 : // write holes
578 0 : if(pvHoles != NULL)
579 : {
580 : vector<vector3>& vHoles = *pvHoles;
581 : out << vHoles.size() << endl;
582 0 : for(uint i = 0; i < vHoles.size(); ++i)
583 : {
584 0 : out << i << " " << vHoles[i].x() << " " <<
585 0 : vHoles[i].y() << " " <<
586 0 : vHoles[i].z() << endl;
587 : }
588 : }
589 : else
590 : {
591 : out << "0" << endl;
592 : }
593 :
594 : out << endl;
595 :
596 : // write regions (optional)
597 0 : if(pvRegionPositions != NULL)
598 : {
599 : vector<vector3>& vRegionPositions = *pvRegionPositions;
600 : out << vRegionPositions.size() << endl;
601 0 : for(uint i = 0; i < vRegionPositions.size(); ++i)
602 : {
603 0 : out << i << " " << vRegionPositions[i].x() << " " <<
604 0 : vRegionPositions[i].y() << " " <<
605 0 : vRegionPositions[i].z();
606 :
607 0 : if(pvRegionAttributes != NULL)
608 0 : out << " " << (*pvRegionAttributes)[i];
609 :
610 0 : if(pvRegionVolumeConstraints != NULL)
611 0 : out << " " << (*pvRegionVolumeConstraints)[i];
612 :
613 : out << endl;
614 : }
615 : }
616 :
617 0 : out.close();
618 :
619 : grid.detach_from_vertices(aInt);
620 :
621 : return true;
622 0 : }
623 :
624 :
625 : ////////////////////////////////////////////////////////////////////////
626 : // ExportGridToSMESH
627 0 : bool ExportGridToSMESH(Grid& grid, const char* filename, AVector3& aPos,
628 : ISubsetHandler* psh,
629 : std::vector<AFloat>* pvNodeAttributes,
630 : std::vector<vector3>* pvHoles,
631 : std::vector<vector3>* pvRegionPositions,
632 : std::vector<int>* pvRegionAttributes,
633 : std::vector<float>* pvRegionVolumeConstraints)
634 :
635 : {
636 0 : if(!grid.has_vertex_attachment(aPos))
637 : return false;
638 :
639 0 : ofstream out(filename);
640 :
641 0 : if(!out)
642 : return false;
643 :
644 : // check if regions are specified in the correct way.
645 0 : if(pvRegionPositions != NULL)
646 : {
647 0 : if(pvRegionAttributes != NULL)
648 0 : if(pvRegionAttributes->size() != pvRegionPositions->size())
649 : return false;
650 :
651 0 : if(pvRegionVolumeConstraints != NULL)
652 0 : if(pvRegionVolumeConstraints->size() != pvRegionPositions->size())
653 : return false;
654 : }
655 :
656 : vector<int> vTmpRegionAttributes;
657 0 : if((pvRegionAttributes == NULL) && (pvRegionPositions != NULL) && (pvRegionVolumeConstraints != NULL))
658 : {
659 : // attributes have to be supplied if constraints are given.
660 : // since no attributes were passed, generate your own.
661 : pvRegionAttributes = &vTmpRegionAttributes;
662 0 : for(uint i = 0; i < pvRegionPositions->size(); ++i)
663 0 : pvRegionAttributes->push_back(0);
664 : }
665 :
666 : // write points. store indices with each point for later array-like access.
667 : AInt aInt;
668 : grid.attach_to_vertices(aInt);
669 : Grid::VertexAttachmentAccessor<AInt> aaIntVRT(grid, aInt);
670 : Grid::VertexAttachmentAccessor<AVector3> aaPos(grid, aPos);
671 :
672 : {
673 : int numAttribs = 0;
674 : int numBoundaryMarkers = 0;
675 :
676 : // attachment-accessors for the nodes attributes
677 : vector<Grid::VertexAttachmentAccessor<AFloat> > vaaFloatVRT;
678 0 : if(pvNodeAttributes != NULL)
679 : {
680 0 : numAttribs = pvNodeAttributes->size();
681 0 : vaaFloatVRT.resize(pvNodeAttributes->size());
682 0 : for(uint i = 0; i < pvNodeAttributes->size(); ++i)
683 : vaaFloatVRT[i].access(grid, (*pvNodeAttributes)[i]);
684 : }
685 :
686 : // attachment-accessor for the nodes boundary-marker
687 0 : if(psh != NULL)
688 : numBoundaryMarkers = 1;
689 :
690 : // write number of nodes, dimension, number of attributes, boundary markers (0 or 1)
691 0 : out << grid.num_vertices() << " 3 " << numAttribs << " " << numBoundaryMarkers << endl;
692 :
693 : int counter = 0;
694 0 : for(VertexIterator iter = grid.vertices_begin(); iter != grid.vertices_end(); iter++, counter++)
695 : {
696 0 : aaIntVRT[*iter] = counter;
697 0 : out << counter << " " << aaPos[*iter].x() << " " <<
698 0 : aaPos[*iter].y() << " " <<
699 0 : aaPos[*iter].z();
700 :
701 : // write attributes:
702 0 : for(uint i = 0; i < vaaFloatVRT.size(); ++i)
703 0 : out << " " << (vaaFloatVRT[i])[*iter];
704 :
705 : // write boundary markers:
706 0 : if(psh != NULL)
707 0 : out << " " << psh->get_subset_index(*iter);
708 :
709 : out << endl;
710 :
711 : }
712 0 : }
713 :
714 : out << endl;
715 :
716 : // write facets
717 : {
718 0 : if(psh != NULL)
719 : out << grid.num_faces() << " 1" << endl;
720 : else
721 : out << grid.num_faces() << " 0" << endl;
722 :
723 0 : for(FaceIterator iter = grid.faces_begin(); iter != grid.faces_end(); ++iter)
724 : {
725 : Face* f = *iter;
726 0 : out << f->num_vertices();
727 :
728 0 : for(uint i = 0; i < f->num_vertices(); ++i)
729 0 : out << " " << aaIntVRT[f->vertex(i)];
730 :
731 0 : if(psh != NULL)
732 0 : out << " " << psh->get_subset_index(f);
733 :
734 : out << endl;
735 : }
736 : }
737 :
738 : out << endl;
739 :
740 : // write holes
741 0 : if(pvHoles != NULL)
742 : {
743 : vector<vector3>& vHoles = *pvHoles;
744 : out << vHoles.size() << endl;
745 0 : for(uint i = 0; i < vHoles.size(); ++i)
746 : {
747 0 : out << i << " " << vHoles[i].x() << " " <<
748 0 : vHoles[i].y() << " " <<
749 0 : vHoles[i].z() << endl;
750 : }
751 : }
752 : else
753 : {
754 : out << "0" << endl;
755 : }
756 :
757 : out << endl;
758 :
759 : // write regions (optional)
760 0 : if(pvRegionPositions != NULL)
761 : {
762 : vector<vector3>& vRegionPositions = *pvRegionPositions;
763 : out << vRegionPositions.size() << endl;
764 0 : for(uint i = 0; i < vRegionPositions.size(); ++i)
765 : {
766 0 : out << i << " " << vRegionPositions[i].x() << " " <<
767 0 : vRegionPositions[i].y() << " " <<
768 0 : vRegionPositions[i].z();
769 :
770 0 : if(pvRegionAttributes != NULL)
771 0 : out << " " << (*pvRegionAttributes)[i];
772 :
773 0 : if(pvRegionVolumeConstraints != NULL)
774 0 : out << " " << (*pvRegionVolumeConstraints)[i];
775 :
776 : out << endl;
777 : }
778 : }
779 :
780 0 : out.close();
781 :
782 : grid.detach_from_vertices(aInt);
783 :
784 : return true;
785 0 : }
786 :
787 : ////////////////////////////////////////////////////////////////////////
788 0 : bool LoadGridFromSMESH(Grid& grid, const char* filename, AVector3& aPos,
789 : ISubsetHandler* psh)
790 : {
791 0 : ifstream file(filename);
792 0 : UG_COND_THROW(!file, "Couldn't open " << filename << " for reading.");
793 :
794 0 : if(!grid.has_vertex_attachment(aPos))
795 0 : grid.attach_to_vertices(aPos);
796 : Grid::VertexAttachmentAccessor<AVector3> aaPos(grid, aPos);
797 :
798 : enum ReadState{
799 : READ_VRT_HEADER,
800 : READ_VERTICES,
801 : READ_FACE_HEADER,
802 : READ_FACES
803 : };
804 :
805 : ReadState readState = READ_VRT_HEADER;
806 :
807 0 : size_t numVrts = 0;
808 0 : size_t dim = 0;
809 0 : size_t numAttribs = 0;
810 0 : size_t bndMarker = 0;
811 0 : size_t numFaces = 0;
812 : size_t numFacesRead = 0;
813 : vector<Vertex*> vrts;
814 :
815 : size_t curLine = 0;
816 :
817 0 : while(!file.eof()){
818 0 : ++curLine;
819 :
820 : string line;
821 0 : getline(file, line);
822 0 : string trimmedLine = TrimString(line);
823 :
824 0 : if(trimmedLine.empty())
825 0 : continue;
826 :
827 0 : if(trimmedLine[0] == '#')
828 0 : continue;
829 :
830 0 : stringstream in(trimmedLine);
831 :
832 : try{
833 0 : switch(readState){
834 : case READ_VRT_HEADER:{
835 : in >> numVrts >> dim >> numAttribs >> bndMarker;
836 0 : UG_COND_THROW(!in, "Couldn't read vertex header");
837 0 : if(numVrts == 0)
838 : return true;
839 : readState = READ_VERTICES;
840 : }break;
841 :
842 :
843 0 : case READ_VERTICES:{
844 : size_t vrtInd;
845 : int tmp;
846 0 : Vertex* vrt = *grid.create<RegularVertex>();
847 : in >> vrtInd >> aaPos[vrt].x() >> aaPos[vrt].y() >> aaPos[vrt].z();
848 0 : for(size_t i = 0; i < numAttribs; ++i)
849 0 : in >> tmp;
850 :
851 0 : if(psh && bndMarker){
852 0 : in >> tmp;
853 0 : psh->assign_subset(vrt, tmp);
854 : }
855 :
856 :
857 0 : UG_COND_THROW(!in, "Couldn't read vertex");
858 :
859 0 : vrts.push_back(vrt);
860 0 : if(vrts.size() >= numVrts)
861 : readState = READ_FACE_HEADER;
862 0 : }break;
863 :
864 :
865 : case READ_FACE_HEADER:{
866 : in >> numFaces >> bndMarker;
867 0 : UG_COND_THROW(!in, "Couldn't read face header");
868 0 : if(numFaces == 0)
869 : return true;
870 : readState = READ_FACES;
871 : }break;
872 :
873 :
874 : case READ_FACES:{
875 : size_t numCorners;
876 : size_t inds[4];
877 :
878 : in >> numCorners;
879 0 : for(size_t i = 0; i < numCorners; ++i)
880 0 : in >> inds[i];
881 :
882 : Face* f = NULL;
883 0 : switch(numCorners){
884 0 : case 3:{
885 0 : f = *grid.create<Triangle>(
886 0 : TriangleDescriptor(
887 0 : vrts[inds[0]], vrts[inds[1]],
888 0 : vrts[inds[2]]));
889 0 : }break;
890 :
891 0 : case 4:{
892 0 : f = *grid.create<Quadrilateral>(
893 0 : QuadrilateralDescriptor(
894 0 : vrts[inds[0]], vrts[inds[1]],
895 0 : vrts[inds[2]], vrts[inds[3]]));
896 0 : }break;
897 :
898 0 : default:
899 0 : UG_THROW("Faces with " << numCorners << " corners are currently not supported");
900 : }
901 :
902 0 : if(psh && bndMarker){
903 : int si;
904 0 : in >> si;
905 0 : psh->assign_subset(f, si);
906 : }
907 :
908 0 : ++numFacesRead;
909 0 : if(numFacesRead >= numFaces)
910 0 : return true;
911 0 : }break;
912 : }
913 : }
914 0 : UG_CATCH_THROW("LoadGridFromSMESH: Failed to interprete data in line " << curLine
915 : << " of file '" << filename << "'");
916 0 : }
917 :
918 : return false;
919 0 : }
920 :
921 :
922 : ////////////////////////////////////////////////////////////////////////
923 : // ExportGridToTETGEN
924 0 : bool ExportGridToTETGEN(Grid& grid, const char* filename,
925 : AVector3& aPos, std::vector<AFloat>* pvNodeAttributes,
926 : AInt* paNodeBoundaryMarker,
927 : AInt* paFaceBoundaryMarker,
928 : AInt* paElementAttribute,
929 : ANumber* paVolumeConstraint)
930 : {
931 0 : if(!grid.has_vertex_attachment(aPos))
932 : return false;
933 :
934 : // set up filenames
935 0 : string eleName = filename;
936 0 : if(eleName.find(".ele", eleName.size() - 4) == string::npos)
937 0 : eleName.append(".ele");
938 0 : string baseName = eleName.substr(0, eleName.size() - 4);
939 0 : string nodeName = baseName + string(".node");
940 0 : string edgeName = baseName + string(".edge");
941 0 : string faceName = baseName + string(".face");
942 0 : string volName = baseName + string(".vol");
943 :
944 :
945 : // write nodes
946 : Grid::VertexAttachmentAccessor<AInt> aaIntVRT;
947 : {
948 0 : ofstream out(nodeName.c_str());
949 :
950 0 : if(!out)
951 : return false;
952 :
953 : // write points. store indices with each point for later array-like access.
954 : AInt aInt;
955 : grid.attach_to_vertices(aInt);
956 : Grid::VertexAttachmentAccessor<AVector3> aaPos(grid, aPos);
957 : aaIntVRT.access(grid, aInt);
958 :
959 :
960 : int numAttribs = 0;
961 : int numBoundaryMarkers = 0;
962 :
963 : // attachment-accessors for the nodes attributes
964 : vector<Grid::VertexAttachmentAccessor<AFloat> > vaaFloatVRT;
965 0 : if(pvNodeAttributes != NULL)
966 : {
967 0 : numAttribs = pvNodeAttributes->size();
968 0 : vaaFloatVRT.resize(pvNodeAttributes->size());
969 0 : for(uint i = 0; i < pvNodeAttributes->size(); ++i)
970 : vaaFloatVRT[i].access(grid, (*pvNodeAttributes)[i]);
971 : }
972 :
973 : // attachment-accessor for the nodes boundary-marker
974 : Grid::VertexAttachmentAccessor<AInt> aaBMVRT;
975 0 : if(paNodeBoundaryMarker != NULL)
976 : {
977 : numBoundaryMarkers = 1;
978 : aaBMVRT.access(grid, *paNodeBoundaryMarker);
979 : }
980 :
981 : // write number of nodes, dimension, number of attributes, boundary markers (0 or 1)
982 0 : out << grid.num_vertices() << " 3 " << numAttribs << " " << numBoundaryMarkers << endl;
983 :
984 : int counter = 0;
985 0 : for(VertexIterator iter = grid.vertices_begin(); iter != grid.vertices_end(); iter++, counter++)
986 : {
987 0 : aaIntVRT[*iter] = counter;
988 0 : out << counter << " " << aaPos[*iter].x() << " " <<
989 0 : aaPos[*iter].y() << " " <<
990 0 : aaPos[*iter].z();
991 :
992 : // write attributes:
993 0 : for(uint i = 0; i < vaaFloatVRT.size(); ++i)
994 : {
995 0 : out << " " << (vaaFloatVRT[i])[*iter];
996 : }
997 :
998 : // write boundary markers:
999 0 : if(paNodeBoundaryMarker != NULL)
1000 0 : out << " " << aaBMVRT[*iter];
1001 :
1002 : out << endl;
1003 : }
1004 :
1005 0 : out.close();
1006 0 : }
1007 :
1008 : // write facets
1009 : {
1010 0 : ofstream out(faceName.c_str());
1011 0 : if(out)
1012 : {
1013 : Grid::FaceAttachmentAccessor<AInt> aaBMFACE;
1014 0 : if(paFaceBoundaryMarker != NULL)
1015 : {
1016 : out << grid.num_faces() << " 1" << endl;
1017 : aaBMFACE.access(grid, *paFaceBoundaryMarker);
1018 : }
1019 : else
1020 : out << grid.num_faces() << " 0" << endl;
1021 :
1022 0 : for(FaceIterator iter = grid.faces_begin(); iter != grid.faces_end(); ++iter)
1023 : {
1024 : Face* f = *iter;
1025 0 : out << f->num_vertices();
1026 :
1027 0 : for(uint i = 0; i < f->num_vertices(); ++i)
1028 0 : out << " " << aaIntVRT[f->vertex(i)];
1029 :
1030 0 : if(paFaceBoundaryMarker != NULL)
1031 0 : out << " " << aaBMFACE[f];
1032 :
1033 : out << endl;
1034 : }
1035 : }
1036 0 : out.close();
1037 0 : }
1038 :
1039 : // write tetrahedrons
1040 : {
1041 0 : ofstream out(eleName.c_str());
1042 0 : if(out)
1043 : {
1044 : Grid::VolumeAttachmentAccessor<AInt> aaElementAttributeVOL;
1045 0 : if(paElementAttribute != NULL)
1046 : {
1047 : aaElementAttributeVOL.access(grid, *paElementAttribute);
1048 : out << grid.num<Tetrahedron>() << " 4 1" << endl;
1049 : }
1050 : else
1051 : out << grid.num<Tetrahedron>() << " 4 0" << endl;
1052 :
1053 : int counter = 0;
1054 : for(TetrahedronIterator iter = grid.begin<Tetrahedron>();
1055 0 : iter != grid.end<Tetrahedron>(); iter++, counter++)
1056 : {
1057 : Tetrahedron* tet = *iter;
1058 0 : out << counter << " " << aaIntVRT[tet->vertex(0)] << " " <<
1059 0 : aaIntVRT[tet->vertex(1)] << " " <<
1060 0 : aaIntVRT[tet->vertex(2)] << " " <<
1061 0 : aaIntVRT[tet->vertex(3)];
1062 :
1063 0 : if(paElementAttribute != NULL)
1064 0 : out << " " << aaElementAttributeVOL[tet];
1065 :
1066 : out << endl;
1067 : }
1068 : }
1069 :
1070 0 : out.close();
1071 0 : }
1072 :
1073 : // write volume constraints
1074 0 : if(paVolumeConstraint)
1075 : {
1076 0 : ofstream out(volName.c_str());
1077 0 : if(out)
1078 : {
1079 : Grid::VolumeAttachmentAccessor<ANumber> aaVolCon(grid, *paVolumeConstraint);
1080 : out << grid.num<Tetrahedron>() << endl;
1081 : int counter = 0;
1082 : for(TetrahedronIterator iter = grid.begin<Tetrahedron>();
1083 0 : iter != grid.end<Tetrahedron>(); iter++, counter++)
1084 : {
1085 0 : out << counter << " " << aaVolCon[*iter] << endl;
1086 : }
1087 : }
1088 :
1089 0 : out.close();
1090 0 : }
1091 :
1092 : return true;
1093 : }
1094 :
1095 : }// end of namespace
|
