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 <vector>
35 : #include <iostream>
36 : #include <cstring>
37 : #include "file_io_art.h"
38 : #include "lib_grid/lg_base.h"
39 : #include "lib_grid/algorithms/geom_obj_util/geom_obj_util.h"
40 : #include "lib_grid/algorithms/attachment_util.h"
41 : using namespace std;
42 :
43 : namespace ug
44 : {
45 :
46 : /**
47 : * If you havn't already used strtok on the given buffer,
48 : * pass true to newTokBuff (default == true).
49 : *
50 : * Make sure that the buffer only contains index values!
51 : */
52 0 : static void ReadIndices(vector<int>& vIndsOut, char* buffer,
53 : const char* delim, bool newTokBuff = true)
54 : {
55 : vIndsOut.clear();
56 :
57 : char* tok;
58 0 : if(newTokBuff)
59 0 : tok = strtok(buffer, delim);
60 : else
61 0 : tok = strtok(NULL, delim);
62 :
63 0 : while(tok)
64 : {
65 0 : vIndsOut.push_back(atoi(tok));
66 0 : tok = strtok(NULL, delim);
67 : }
68 0 : }
69 :
70 : template <class TType>
71 :
72 0 : static void CollectUniqueObjects(vector<TType>& vecOut,
73 : const vector<TType>& vecIn)
74 : {
75 : // copy indices to avoid problems when vIndsOut and vIndsInd operate on
76 : // the same storage.
77 : bool identicContainers = (&vecOut == &vecIn);
78 :
79 : vector<TType> vTmp;
80 : vector<TType>* pDest = &vecOut;
81 0 : if(identicContainers)
82 : pDest = &vTmp;
83 :
84 : // we're operating on destInds
85 : vector<TType>& vecDest = *pDest;
86 : vecDest.clear();
87 :
88 : // iterate over all elements of vInds
89 0 : for(size_t i = 0; i < vecIn.size(); ++i){
90 : const TType& val = vecIn[i];
91 : // check whether val is already contained in vecDest
92 : bool gotOne = false;
93 0 : for(size_t j = 0; j < vecDest.size(); ++j){
94 0 : if(vecDest[j] == val){
95 : gotOne = true;
96 : break;
97 : }
98 : }
99 :
100 : // if we havn't found one, we'll insert val into vecDest
101 0 : if(!gotOne)
102 0 : vecDest.push_back(val);
103 : }
104 :
105 : // swap storage if required
106 0 : if(identicContainers)
107 : vecOut.swap(vTmp);
108 0 : }
109 :
110 : // uses grid::mark
111 : static Tetrahedron*
112 0 : CreateTetrahedron(Grid& grid, vector<Triangle*>& vTris,
113 : Grid::VertexAttachmentAccessor<APosition>& aaPos)
114 : {
115 : const char* errorMsg = "ERROR in file_io_art.cpp CreateTetrahedron. ";
116 :
117 0 : if(vTris.size() != 4){
118 0 : UG_LOG(errorMsg << "Bad number of triangles: " << vTris.size() << endl);
119 0 : return NULL;
120 : }
121 :
122 : Vertex* vrts[4];
123 : int vrtCount = 0;
124 :
125 : // get the 4 points of the tetrahedron
126 :
127 0 : grid.begin_marking();
128 :
129 0 : for(size_t i = 0; i < 4; ++i){
130 0 : for(size_t j = 0; j < 3; ++j){
131 0 : Vertex* vrt = vTris[i]->vertex(j);
132 0 : if(!grid.is_marked(vrt)){
133 : // make sure that we won't collect too many vertices.
134 0 : if(vrtCount == 4){
135 0 : UG_LOG(errorMsg << "Triangles do not build a tetrahedron.\n");
136 0 : return NULL;
137 : }
138 :
139 : grid.mark(vrt);
140 0 : vrts[vrtCount++] = vrt;
141 : }
142 : }
143 : }
144 :
145 0 : grid.end_marking();
146 :
147 0 : if(vrtCount < 4){
148 0 : UG_LOG(errorMsg << "Triangles have less than 4 distinct vertices.\n");
149 0 : return NULL;
150 : }
151 :
152 : // we have to check the orientation of the tetrahedron
153 : // compare the normal of the first triangle to the top-vertex
154 0 : if(PointFaceTest(aaPos[vrts[3]], vTris[0], aaPos) < 0){
155 : // we have to change the order of the first three vrts
156 : swap(vrts[0], vrts[1]);
157 : }
158 :
159 : // create the tetrahedron
160 0 : return *grid.create<Tetrahedron>(
161 0 : TetrahedronDescriptor(vrts[0], vrts[1],
162 : vrts[2], vrts[3]));
163 :
164 : }
165 :
166 : // uses grid::mark
167 : static Pyramid*
168 0 : CreatePyramid(Grid& grid, vector<Triangle*>& vTris,
169 : vector<Quadrilateral*>& vQuads,
170 : Grid::VertexAttachmentAccessor<APosition>& aaPos)
171 : {
172 : const char* errorMsg = "ERROR in file_io_art.cpp CreatePyramid. ";
173 :
174 0 : if(vTris.size() != 4){
175 0 : UG_LOG(errorMsg << "Bad number of triangles: " << vTris.size() << endl);
176 0 : return NULL;
177 : }
178 :
179 0 : if(vQuads.size() != 1){
180 0 : UG_LOG(errorMsg << "Bad number of quadrilaterals: " << vQuads.size() << endl);
181 0 : return NULL;
182 : }
183 :
184 : Vertex* vrts[5];
185 : int vrtCount = 0;
186 :
187 0 : grid.begin_marking();
188 :
189 : // get the 5 points of the pyramid
190 : // take the vertices of the base-quad first
191 0 : for(size_t j = 0; j < 4; ++j){
192 0 : Vertex* vrt = vQuads[0]->vertex(j);
193 0 : if(!grid.is_marked(vrt)){
194 : grid.mark(vrt);
195 0 : vrts[vrtCount++] = vrt;
196 : }
197 : }
198 :
199 : // now find the top vertex by checking the first triangle
200 0 : for(size_t j = 0; j < 3; ++j){
201 0 : Vertex* vrt = vTris[0]->vertex(j);
202 0 : if(!grid.is_marked(vrt)){
203 : grid.mark(vrt);
204 0 : vrts[vrtCount++] = vrt;
205 0 : break;
206 : }
207 : }
208 :
209 0 : grid.end_marking();
210 :
211 : // make sure that all vertices have been found
212 0 : if(vrtCount != 5){
213 0 : UG_LOG(errorMsg << "Faces do not build a pyramid.\n");
214 0 : return NULL;
215 : }
216 :
217 : // we have to check the orientation of the prism
218 : // compare the normal of the base-quad to the top
219 0 : if(PointFaceTest(aaPos[vrts[4]], vQuads[0], aaPos) < 0){
220 : // we have to change the order of the first four vrts
221 : swap(vrts[1], vrts[2]);
222 : swap(vrts[0], vrts[3]);
223 : }
224 :
225 : // create the pyramid
226 0 : return *grid.create<Pyramid>(
227 0 : PyramidDescriptor(vrts[0], vrts[1],
228 : vrts[2], vrts[3], vrts[4]));
229 : }
230 :
231 :
232 : // this method assumes that the given faces already exist in the grid
233 : // together with their associated edges
234 : static Prism*
235 0 : CreatePrism(Grid& grid, vector<Triangle*>& vTris,
236 : vector<Quadrilateral*>& vQuads,
237 : Grid::VertexAttachmentAccessor<APosition>& aaPos)
238 : {
239 : const char* errorMsg = "ERROR in file_io_art.cpp CreatePrism. ";
240 :
241 0 : if(vTris.size() != 2){
242 0 : UG_LOG(errorMsg << "Bad number of triangles: " << vTris.size() << endl);
243 0 : return NULL;
244 : }
245 :
246 0 : if(vQuads.size() != 3){
247 0 : UG_LOG(errorMsg << "Bad number of quadrilaterals: " << vQuads.size() << endl);
248 0 : return NULL;
249 : }
250 :
251 : Vertex* vrts[6];
252 : int vrtCount = 0;
253 :
254 : // get the 6 points of the prism
255 : // calculate the center on the fly
256 : vector3 vCenter(0, 0, 0);
257 :
258 0 : for(size_t i = 0; i < 2; ++i){
259 0 : for(size_t j = 0; j < 3; ++j){
260 0 : vrts[vrtCount++] = vTris[i]->vertex(j);
261 0 : VecAdd(vCenter, vCenter, aaPos[vTris[i]->vertex(j)]);
262 : }
263 : }
264 :
265 : VecScale(vCenter, vCenter, 1. / 6.);
266 :
267 : // we have to check the orientation of the prism
268 : // compare the normal of the first triangle to the center
269 0 : if(PointFaceTest(vCenter, vTris[0], aaPos) < 0){
270 : // we have to change the order of the first three vrts
271 : swap(vrts[0], vrts[1]);
272 : }
273 :
274 : // compare the normal of the second triangle to the center
275 0 : if(PointFaceTest(vCenter, vTris[1], aaPos) > 0){
276 : // we have to change the order of the last three vrts
277 : swap(vrts[3], vrts[4]);
278 : }
279 :
280 : // we have to find the vertex of the second triangle that
281 : // is connected to the first vertex of the first triangle
282 : int index = 0;
283 0 : for(; index < 3; ++index){
284 0 : if(grid.get_edge(vrts[0], vrts[3 + index]))
285 : break;
286 : }
287 :
288 0 : if(index == 3)
289 : return NULL;
290 :
291 : // create the tetrahedron
292 0 : return *grid.create<Prism>(
293 0 : PrismDescriptor(vrts[0], vrts[1], vrts[2],
294 0 : vrts[3 + index],
295 0 : vrts[3 + (index + 1) % 3],
296 0 : vrts[3 + (index + 2) % 3]));
297 : }
298 :
299 :
300 0 : bool LoadGridFromART(Grid& grid, const char* filename,
301 : ISubsetHandler* pSH,
302 : AVector3& aPos)
303 : {
304 : // open the stream
305 0 : ifstream in(filename);
306 :
307 0 : if(!in)
308 : return false;
309 :
310 : // The subset handler:
311 0 : SubsetHandler shTmp;
312 0 : if(!pSH)
313 : {
314 0 : shTmp.assign_grid(grid);
315 : pSH = &shTmp;
316 : }
317 : ISubsetHandler& sh = *pSH;
318 :
319 : // make sure that the position attachment is attached properly
320 0 : if(!grid.has_vertex_attachment(aPos))
321 0 : grid.attach_to_vertices(aPos);
322 : Grid::VertexAttachmentAccessor<AVector3> aaPos(grid, aPos);
323 :
324 : // this buffer will be used to store all the lines.
325 0 : char* buf = new char[512];
326 : const char* delim = " :";
327 :
328 : // those vectors will be reused at different sections
329 : vector<int> vInds;
330 : vector<Vertex*> vVrtDump;
331 : vector<Vertex*> vLocalVrts;
332 :
333 : // read the vertices
334 : vector<RegularVertex*> vVrts;
335 :
336 0 : while((!in.eof()))
337 : {
338 0 : in.getline(buf, 512);
339 0 : char* tok = strtok(buf, delim);
340 0 : if(!tok)
341 0 : continue;
342 :
343 : // ignore all lines that start with a %
344 0 : if(*tok == '%')
345 0 : continue;
346 : // if we found a $ this element-type is done
347 0 : if(*tok == '$')
348 : break;
349 :
350 : // create a new vertex
351 0 : RegularVertex* v = *grid.create<RegularVertex>();
352 : // read the coordinates
353 : //TODO: make sure that everything is ok.
354 0 : aaPos[v].x() = atof(tok);
355 0 : tok = strtok(NULL, delim);
356 0 : aaPos[v].y() = atof(tok);
357 0 : tok = strtok(NULL, delim);
358 0 : aaPos[v].z() = atof(tok);
359 :
360 : // store the vertex in an array
361 0 : vVrts.push_back(v);
362 : }
363 :
364 : // read the edges
365 : vector<RegularEdge*> vEdges;
366 :
367 0 : while((!in.eof()))
368 : {
369 0 : in.getline(buf, 512);
370 0 : char* tok = strtok(buf, delim);
371 0 : if(!tok)
372 0 : continue;
373 :
374 : // ignore all lines that start with a %
375 0 : if(*tok == '%')
376 0 : continue;
377 : // if we found a $ this element-type is done
378 0 : if(*tok == '$')
379 : break;
380 :
381 : // read the indices
382 : //TODO: make sure that everything is ok.
383 : int i1, i2;
384 : //int si = atoi(tok);
385 0 : tok = strtok(NULL, delim);
386 : i1 = atoi(tok);
387 0 : tok = strtok(NULL, delim);
388 : i2 = atoi(tok);
389 :
390 : // create a new edge
391 0 : RegularEdge* e = *grid.create<RegularEdge>(EdgeDescriptor(vVrts[i1], vVrts[i2]));
392 : // edges won't be assigned to subsets in the moment, since things are a little chaotic in the files...
393 : /*
394 : if(si < 1)
395 : si = -1;
396 : else if(si == 10000)
397 : si = 0;
398 :
399 : if(si != -1)
400 : sh.assign_subset(e, si);
401 : */
402 : // store the edge in an array
403 0 : vEdges.push_back(e);
404 : }
405 :
406 : // read the faces
407 : vector<Face*> vFaces;
408 :
409 0 : while((!in.eof()))
410 : {
411 0 : in.getline(buf, 512);
412 0 : char* tok = strtok(buf, delim);
413 0 : if(!tok)
414 0 : continue;
415 :
416 : // ignore all lines that start with a %
417 0 : if(*tok == '%')
418 0 : continue;
419 : // if we found a $ this element-type is done
420 0 : if(*tok == '$')
421 : break;
422 :
423 : int si = atoi(tok);
424 :
425 : // read the indices
426 0 : ReadIndices(vInds, buf, delim, false);
427 :
428 : // there have to be at least three edges
429 0 : if(vInds.size() < 3)
430 0 : continue;
431 :
432 : // copy all vertices of the edges to vVrtDump
433 : //TODO: make sure that all indices are in the correct ranges
434 : vVrtDump.clear();
435 0 : for(size_t i = 0; i < vInds.size(); ++i){
436 0 : Edge* e = vEdges[vInds[i]];
437 0 : vVrtDump.push_back(e->vertex(0));
438 0 : vVrtDump.push_back(e->vertex(1));
439 : }
440 :
441 : // now collect the unique vertices in vVrtDump
442 0 : CollectUniqueObjects(vLocalVrts, vVrtDump);
443 :
444 : //TODO: check orientation
445 : // create the faces
446 0 : Face* f = NULL;
447 : size_t numVrts = vLocalVrts.size();
448 0 : switch(numVrts)
449 : {
450 : case 3:
451 0 : f = *grid.create<Triangle>(TriangleDescriptor(vLocalVrts[0], vLocalVrts[1],
452 : vLocalVrts[2]));
453 0 : break;
454 : case 4:
455 0 : f = *grid.create<Quadrilateral>(QuadrilateralDescriptor(vLocalVrts[0], vLocalVrts[1],
456 : vLocalVrts[2], vLocalVrts[3]));
457 0 : break;
458 : default:
459 : LOG(" LoadGridFromART: bad number of vertices in face: " << numVrts << " (3 or 4 supported).\n");
460 0 : continue;
461 : };
462 :
463 0 : if(si < 1)
464 : si = -1;
465 0 : else if(si == 10000)
466 : si = 0;
467 :
468 : if(si != -1)
469 0 : sh.assign_subset(f, si);
470 :
471 : // store the face in an array
472 0 : vFaces.push_back(f);
473 : }
474 :
475 : // read the volumes
476 0 : vector<Triangle*> vTris;
477 0 : vector<Quadrilateral*> vQuads;
478 :
479 0 : while((!in.eof()))
480 : {
481 0 : in.getline(buf, 512);
482 0 : char* tok = strtok(buf, delim);
483 0 : if(!tok)
484 0 : continue;
485 :
486 : // ignore all lines that start with a %
487 0 : if(*tok == '%')
488 0 : continue;
489 : // if we found a $ this element-type is done
490 0 : if(*tok == '$')
491 : break;
492 :
493 : int si = atoi(tok);
494 :
495 : // read the indices
496 0 : ReadIndices(vInds, buf, delim, false);
497 :
498 : // there have to be at least three faces
499 0 : if(vInds.size() < 3)
500 0 : continue;
501 :
502 : // collect faces in vTris and vQuads
503 : vTris.clear();
504 : vQuads.clear();
505 0 : for(size_t i = 0; i < vInds.size(); ++i){
506 0 : Face* f = vFaces[vInds[i]];
507 0 : if(f->num_vertices() == 3){
508 0 : Triangle* t = dynamic_cast<Triangle*>(f);
509 0 : if(t)
510 0 : vTris.push_back(t);
511 : else{
512 : UG_LOG(" bad triangle in LoadGridFromART!");
513 : }
514 : }
515 0 : else if(f->num_vertices() == 4){
516 0 : Quadrilateral* q = dynamic_cast<Quadrilateral*>(f);
517 0 : if(q)
518 0 : vQuads.push_back(q);
519 : else{
520 : UG_LOG(" bad quadrilateral in LoadGridFromART!");
521 : }
522 : }
523 : else {
524 : UG_LOG(" Bad face in LoadGridFromART! Aborting...");
525 : return false;
526 : }
527 : }
528 :
529 : // get the type of volume
530 : int volType = ROID_UNKNOWN;
531 0 : if(vTris.size() == 4 && vQuads.size() == 0)
532 : volType = ROID_TETRAHEDRON;
533 0 : else if(vTris.size() == 4 && vQuads.size() == 1)
534 : volType = ROID_PYRAMID;
535 0 : else if(vTris.size() == 2 && vQuads.size() == 3)
536 : volType = ROID_PRISM;
537 0 : else if(vTris.size() == 0 && vQuads.size() == 6)
538 : volType = ROID_HEXAHEDRON;
539 :
540 :
541 : // create the volume
542 : Volume* v = NULL;
543 :
544 0 : switch(volType)
545 : {
546 : case ROID_TETRAHEDRON:
547 0 : v = CreateTetrahedron(grid, vTris, aaPos);
548 : break;
549 : case ROID_PYRAMID:
550 0 : v = CreatePyramid(grid, vTris, vQuads, aaPos);
551 : break;
552 : case ROID_PRISM:
553 0 : v = CreatePrism(grid, vTris, vQuads, aaPos);
554 : break;
555 : case ROID_HEXAHEDRON:
556 : //TODO: create hexahedron
557 : break;
558 : default:
559 : LOG(" LoadGridFromART: bad volume type. volume has "
560 : << vTris.size() << " triangles and "
561 : << vQuads.size() << " quadrilaterals.\n");
562 0 : continue;
563 : };
564 :
565 0 : si -= 1;
566 0 : if(v)
567 0 : sh.assign_subset(v, si);
568 : else {
569 : LOG(" LoadGridFromART: could not create volume element.\n");
570 : }
571 : }
572 :
573 : // delete the buffer
574 0 : delete[] buf;
575 :
576 0 : return true;
577 0 : }
578 :
579 : ////////////////////////////////////////////////////////////////////////
580 0 : bool SaveGridToART(Grid& srcGrid, const char* filename,
581 : ISubsetHandler* pSH, AVector3& aPos)
582 : {
583 0 : if(!srcGrid.has_vertex_attachment(aPos)){
584 : LOG(" Aborting SaveGridToART: position attachment is missing.\n");
585 0 : return false;
586 : }
587 :
588 : // open the file
589 0 : ofstream out(filename);
590 : out.precision(20);
591 : out.setf(ios_base::scientific);
592 :
593 0 : if(!out)
594 : return false;
595 :
596 : // Options FACEOPT_AUTOGENERATE_EDGES and VOLOPT_AUTOGENERATE_FACES
597 : // have to be enabled. If they are not, we'll create a copy of the
598 : // grid and enable those options there.
599 : // note that we need a second subset handler too
600 0 : Grid tGrid;
601 0 : SubsetHandler tSH(tGrid);
602 : Grid* pGrid; // Used to initialise the reference Grid& grid later on.
603 0 : if(!srcGrid.option_is_enabled(FACEOPT_AUTOGENERATE_EDGES
604 : | VOLOPT_AUTOGENERATE_FACES))
605 : {
606 : // copy grid and subset-handler
607 0 : tGrid = srcGrid;
608 0 : if(pSH){
609 0 : tSH = *pSH;
610 : pSH = &tSH;
611 : }
612 :
613 0 : tGrid.enable_options(VOLOPT_AUTOGENERATE_FACES);
614 0 : tGrid.enable_options(FACEOPT_AUTOGENERATE_EDGES);
615 : pGrid = &tGrid;
616 :
617 : // make sure that the position attachment has been copied
618 0 : if(!tGrid.has_vertex_attachment(aPos)){
619 : // copy it manually
620 0 : if(!CopyAttachments<Vertex>(srcGrid, aPos, tGrid, aPos))
621 : return false;
622 : }
623 : }
624 : else {
625 : pGrid = &srcGrid;
626 : }
627 :
628 : // We'll work on this reference.
629 : Grid& grid = *pGrid;
630 :
631 : // write the header
632 : out << "%% Version 3.0" << endl;
633 : out << "%% VertexNumber: " << grid.num<Vertex>() << endl;
634 : out << "%% EdgeNumber: " << grid.num<Edge>() << endl;
635 : out << "%% FaceNumber: " << grid.num<Face>() << endl;
636 : out << "%% ElementNumber: " << grid.num<Volume>() << endl;
637 : out << "%% Translation: (0, 0, 0)" << endl;
638 : out << "%% Cornermark: 10001" << endl;
639 : out << "%% DO NOT CHANGE LINES ABOVE !!!!" << endl;
640 : out << "% NET: Vertices <-> Edges <-> Faces <-> Elements" << endl;
641 :
642 :
643 : // write vertices
644 : {
645 : out << "% Vertices: x y z" << endl;
646 : Grid::VertexAttachmentAccessor<AVector3> aaPos(grid, aPos);
647 : for(VertexIterator iter = grid.begin<Vertex>();
648 0 : iter != grid.end<Vertex>(); ++iter)
649 : {
650 0 : out << aaPos[*iter].x() << " ";
651 0 : out << aaPos[*iter].y() << " ";
652 0 : out << aaPos[*iter].z() << endl;
653 : }
654 : // done - write end sign
655 : out << "$" << endl;
656 : }
657 :
658 : // write edges
659 : {
660 : // vertices have to be associated with indices
661 : AInt aInt;
662 : grid.attach_to_vertices(aInt);
663 : Grid::VertexAttachmentAccessor<AInt> aaInt(grid, aInt);
664 : AssignIndices<Vertex>(grid.begin<Vertex>(), grid.end<Vertex>(), aaInt);
665 :
666 : // the subset index
667 : int si = 0;
668 :
669 : // begin writing
670 : out << "% Edges (Indices to List of Points):" << endl;
671 : // iterate over edges
672 : for(EdgeIterator iter = grid.begin<Edge>();
673 0 : iter != grid.end<Edge>(); ++iter)
674 : {
675 : Edge* e = *iter;
676 :
677 0 : if(pSH){
678 : si = pSH->get_subset_index(e);
679 0 : if(si == 0)
680 : si = 10000;
681 0 : else if(si == -1)
682 : si = 0;
683 : }
684 :
685 0 : out << si << ": " << aaInt[e->vertex(0)] << " " << aaInt[e->vertex(1)] << endl;
686 : }
687 :
688 : // done - write end sign
689 : out << "$" << endl;
690 :
691 : // detach indices
692 : grid.detach_from_vertices(aInt);
693 : }
694 :
695 : // write faces
696 : {
697 :
698 : // edges have to be associated with indices
699 : AInt aInt;
700 : grid.attach_to_edges(aInt);
701 : Grid::EdgeAttachmentAccessor<AInt> aaInt(grid, aInt);
702 : AssignIndices<Edge>(grid.begin<Edge>(), grid.end<Edge>(), aaInt);
703 :
704 : // the subset index
705 : int si = 0;
706 :
707 : // we'll collect edges in here
708 : vector<Edge*> vEdges;
709 :
710 : // begin writing
711 : out << "% Faces (Indices to List of Edges):" << endl;
712 : // iterate over faces
713 : for(FaceIterator iter = grid.begin<Face>();
714 0 : iter != grid.end<Face>(); ++iter)
715 : {
716 : Face* f = *iter;
717 0 : if(pSH){
718 : si = pSH->get_subset_index(f);
719 0 : if(si == 0)
720 : si = 10000;
721 0 : else if(si == -1)
722 : si = 0;
723 : }
724 :
725 : // collect edges
726 0 : CollectEdges(vEdges, grid, f);
727 :
728 : // write edges
729 0 : out << si << ":";
730 0 : for(size_t i = 0; i < vEdges.size(); ++i)
731 0 : out << " " << aaInt[vEdges[i]];
732 : out << endl;
733 : }
734 :
735 : // done - write end sign
736 : out << "$" << endl;
737 :
738 : // detach indices
739 : grid.detach_from_edges(aInt);
740 0 : }
741 :
742 : // write volumes
743 : {
744 :
745 : // faces have to be associated with indices
746 : AInt aInt;
747 : grid.attach_to_faces(aInt);
748 : Grid::FaceAttachmentAccessor<AInt> aaInt(grid, aInt);
749 : AssignIndices<Face>(grid.begin<Face>(), grid.end<Face>(), aaInt);
750 :
751 : // the subset index
752 : int si = 1;
753 :
754 : // we'll collect faces in here
755 : vector<Face*> vFaces;
756 :
757 : // begin writing
758 : out << "% Elements (Indices to List of Faces):" << endl;
759 : // iterate over volumes
760 : for(VolumeIterator iter = grid.begin<Volume>();
761 0 : iter != grid.end<Volume>(); ++iter)
762 : {
763 : Volume* v = *iter;
764 0 : if(pSH){
765 0 : si = pSH->get_subset_index(v) + 1;
766 : }
767 :
768 : // collect faces
769 0 : CollectFaces(vFaces, grid, v);
770 :
771 : // write faces
772 0 : out << si << ":";
773 0 : for(size_t i = 0; i < vFaces.size(); ++i)
774 0 : out << " " << aaInt[vFaces[i]];
775 : out << endl;
776 : }
777 :
778 : // done - write end sign
779 : out << "$" << endl;
780 :
781 : // detach indices
782 : grid.detach_from_faces(aInt);
783 0 : }
784 :
785 0 : out.close();
786 : return true;
787 0 : }
788 :
789 : }// end of namespace
|
