Line data Source code
1 : /*
2 : * Copyright (c) 2009-2015: G-CSC, Goethe University Frankfurt
3 : * Authors: Martin Stepniewski, 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 <vector>
34 : #include <string>
35 : #include <fstream>
36 : #include <sstream>
37 : #include <cstdlib>
38 : #include <cstring>
39 : #include <algorithm>
40 : #include <set>
41 : #include "file_io_ug.h"
42 : #include "lib_grid/lg_base.h"
43 : #include "lib_grid/algorithms/attachment_util.h"
44 : #include "lib_grid/algorithms/geom_obj_util/geom_obj_util.h"
45 : #include "lib_grid/algorithms/orientation_util.h"
46 : #include "lib_grid/algorithms/polychain_util.h"
47 : #include "lib_grid/callbacks/callbacks.h"
48 :
49 :
50 : using namespace std;
51 :
52 : namespace ug
53 : {
54 :
55 : static bool CollectLines(Grid& grid, const SubsetHandler& shFace, EdgeSelector& LineSel);
56 :
57 : static bool CollectInnerVertices(Grid& grid, VertexSelector& InnVrtSel, VertexSelector& SurfVrtSel);
58 :
59 : static bool CollectSurfaceVertices(Grid& grid, const SubsetHandler& shFace,
60 : VertexSelector& SurfVrtSel);
61 :
62 : static bool CollectAllVerticesForNG(Grid& grid, VertexSelector& NgVrtSel,
63 : VertexSelector& SurfVrtSel, VertexSelector& InnVrtSel);
64 :
65 : static bool WriteLGM(Grid& grid,
66 : const char* lgmFilename,
67 : const char* problemName,
68 : const char* lgmName,
69 : int convex,
70 : const SubsetHandler& shFaces,
71 : const SubsetHandler& shVolumes,
72 : EdgeSelector& LineSel,
73 : VertexSelector& SurfVrtSel,
74 : Grid::EdgeAttachmentAccessor<AInt>& aaLineIndex,
75 : Grid::VertexAttachmentAccessor<AInt>& aaSurfVrtIndex,
76 : Grid::VertexAttachmentAccessor<AVector3>& aaPos);
77 :
78 : static bool GetRightLeftUnitIndex(int& rightIndex, int& leftIndex, Grid& grid, Face* face,
79 : const SubsetHandler& shVolume);
80 :
81 : static bool WriteNG(Grid& grid,
82 : const SubsetHandler& shFaces,
83 : const SubsetHandler& shVolumes,
84 : const char* ngFilename,
85 : VertexSelector& SurfVrtSel,
86 : VertexSelector& InnVrtSel,
87 : VertexSelector& NgVrtSel,
88 : EdgeSelector& LineSel,
89 : Grid::EdgeAttachmentAccessor<AInt>& aaLineIndex,
90 : Grid::VertexAttachmentAccessor<AInt>& aaInnVrtIndex,
91 : Grid::VertexAttachmentAccessor<AInt>& aaNgVrtIndex,
92 : Grid::FaceAttachmentAccessor<AInt>& aaFaceIndex,
93 : Grid::VertexAttachmentAccessor<AVector3>& aaPos);
94 :
95 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
96 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
97 : // ConvertTETGENToUG
98 : /// converts tetgen files (*.node, *.face and *.ele) to UG files (*.lgm, *.ng)
99 0 : bool ExportGridToUG(const Grid& g, const SubsetHandler& shFace, const SubsetHandler& shVolume,
100 : const char* fileNamePrefix, const char* lgmName,
101 : const char* problemName, int convex)
102 : {
103 : // the original grid may not be altered
104 0 : Grid grid = g;
105 :
106 : // we need subset-handlers that operate on the local grid
107 0 : SubsetHandler shFaces(grid, SHE_FACE);
108 0 : SubsetHandler shVolumes(grid, SHE_VOLUME);
109 0 : shFaces = shFace;
110 0 : shVolumes = shVolume;
111 :
112 : // fix orientation of faces
113 0 : for(int i = 0; i < shFaces.num_subsets(); ++i)
114 0 : FixFaceOrientation(grid, shFaces.begin<Face>(i), shFaces.end<Face>(i));
115 :
116 : // make sure that all face-subsets are in consecutive order
117 : {
118 : bool foundEmpty = false;
119 0 : for(int i = 0; i < shFaces.num_subsets(); ++i){
120 0 : if(shFaces.num<Face>(i) == 0)
121 : foundEmpty = true;
122 : else{
123 0 : if(foundEmpty){
124 : // there must be no empty face-subsets between filled ones.
125 : UG_LOG("WARNING in ExportGridToUG: Empty face subset found between filled ones. Aborting...\n");
126 : return false;
127 : }
128 : }
129 : }
130 : }
131 :
132 : // make sure that all volume-subsets are in consecutive order
133 : {
134 : bool foundEmpty = false;
135 0 : for(int i = 0; i < shVolumes.num_subsets(); ++i){
136 0 : if(shVolumes.num<Volume>(i) == 0)
137 : foundEmpty = true;
138 : else{
139 0 : if(foundEmpty){
140 : // there must be no empty volume-subsets between filled ones.
141 : UG_LOG("WARNING in ExportGridToUG: Empty volume subset found between filled ones. Aborting...\n");
142 : return false;
143 : }
144 : }
145 : }
146 : }
147 :
148 : // initialization
149 0 : EdgeSelector LineSel(grid);
150 0 : VertexSelector NgVrtSel(grid);
151 0 : VertexSelector SurfVrtSel(grid);
152 0 : VertexSelector InnVrtSel(grid);
153 :
154 0 : grid.enable_options(FACEOPT_AUTOGENERATE_EDGES);
155 0 : grid.enable_options(EDGEOPT_STORE_ASSOCIATED_FACES);
156 : // for write ng
157 0 : grid.enable_options(FACEOPT_STORE_ASSOCIATED_EDGES);
158 0 : grid.enable_options(VRTOPT_STORE_ASSOCIATED_FACES);
159 :
160 : Grid::VertexAttachmentAccessor<AVector3> aaPos(grid, aPosition);
161 :
162 : // selection of lines, surface-vertices, inner vertices
163 0 : CollectLines(grid, shFaces, LineSel);
164 0 : CollectSurfaceVertices(grid, shFaces, SurfVrtSel);
165 0 : CollectInnerVertices(grid, InnVrtSel, SurfVrtSel);
166 0 : CollectAllVerticesForNG(grid, NgVrtSel, SurfVrtSel, InnVrtSel);
167 :
168 :
169 : //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
170 : // INDEX ASSIGNMENT SECTION >>>>>>>>>>>>>>>>>>>>>>>>>
171 : //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
172 : // assign index to every line
173 : AInt aLineIndex;
174 : grid.attach_to_edges(aLineIndex);
175 : Grid::EdgeAttachmentAccessor<AInt> aaLineIndex(grid, aLineIndex);
176 : int counter = 0;
177 0 : for(EdgeIterator EIter = LineSel.begin(); EIter != LineSel.end(); ++EIter, ++counter)
178 0 : aaLineIndex[*EIter] = counter;
179 :
180 : // assign index to every face regarding EACH surface
181 : AInt aFaceIndex;
182 : grid.attach_to_faces(aFaceIndex);
183 : Grid::FaceAttachmentAccessor<AInt> aaFaceIndex(grid, aFaceIndex);
184 0 : for(int i = 0; i < shFaces.num_subsets(); ++i)
185 : {
186 : counter = 0;
187 : // assign triangle indices
188 0 : for(TriangleIterator iter = shFaces.begin<Triangle>(i);
189 0 : iter != shFaces.end<Triangle>(i); ++iter, ++counter)
190 0 : aaFaceIndex[*iter] = counter;
191 :
192 : // assign quadrilaterals. Increase index by 2, since 2 triangles are written for each quad
193 0 : for(QuadrilateralIterator iter = shFaces.begin<Quadrilateral>(i);
194 0 : iter != shFaces.end<Quadrilateral>(i); ++iter, counter += 2)
195 0 : aaFaceIndex[*iter] = counter;
196 : }
197 :
198 : // assign index to every grid-vertex in *.ng file order
199 : AInt aNgVrtIndex;
200 : grid.attach_to_vertices(aNgVrtIndex);
201 : Grid::VertexAttachmentAccessor<AInt> aaNgVrtIndex(grid, aNgVrtIndex);
202 : counter = 0;
203 0 : for(VertexIterator VIter = NgVrtSel.begin(); VIter != NgVrtSel.end(); ++VIter, ++counter)
204 0 : aaNgVrtIndex[*VIter] = counter;
205 :
206 : // assign index to every inner-vertex
207 : AInt aInnVrtIndex;
208 : grid.attach_to_vertices(aInnVrtIndex);
209 : Grid::VertexAttachmentAccessor<AInt> aaInnVrtIndex(grid, aInnVrtIndex);
210 : counter = 0;
211 0 : for(VertexIterator VIter = InnVrtSel.begin(); VIter != InnVrtSel.end(); ++VIter, ++counter)
212 0 : aaInnVrtIndex[*VIter] = counter;
213 :
214 : // assign index to every surface-vertex
215 : AInt aSurfVrtIndex;
216 : grid.attach_to_vertices(aSurfVrtIndex);
217 : Grid::VertexAttachmentAccessor<AInt> aaSurfVrtIndex(grid, aSurfVrtIndex);
218 : counter = 0;
219 0 : for(VertexIterator SVIter = SurfVrtSel.begin(); SVIter != SurfVrtSel.end(); ++SVIter, ++counter)
220 0 : aaSurfVrtIndex[*SVIter] = counter;
221 : //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
222 : //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
223 0 : string lgmFilename(fileNamePrefix);
224 0 : lgmFilename.append(".lgm");
225 : // write *.lgm file
226 0 : WriteLGM(grid, lgmFilename.c_str(), problemName, lgmName,
227 : convex, shFaces, shVolumes, LineSel, SurfVrtSel, aaLineIndex, aaSurfVrtIndex, aaPos);
228 :
229 0 : string ngFilename(fileNamePrefix);
230 0 : ngFilename.append(".ng");
231 :
232 : // write *.ng file
233 0 : WriteNG(grid, shFaces, shVolumes, ngFilename.c_str(), SurfVrtSel, InnVrtSel, NgVrtSel, LineSel,
234 : aaLineIndex, aaInnVrtIndex, aaNgVrtIndex, aaFaceIndex, aaPos);
235 :
236 : return true;
237 0 : }
238 :
239 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
240 : // CollectLines
241 : /// collects lines in a selector
242 0 : static bool CollectLines(Grid& grid, const SubsetHandler& shFace, EdgeSelector& LineSel)
243 : {
244 : // store associated faces in this vector
245 : vector<Face*> vFaces;
246 :
247 : // iterate through all edges in the grid and identify lines by comparing the subset-membership
248 : // of the associated faces
249 0 : for(EdgeIterator EIter = grid.edges_begin(); EIter != grid.edges_end(); ++EIter)
250 : {
251 0 : CollectFaces(vFaces, grid, *EIter);
252 :
253 : //if(vFaces.size() > 1)
254 : {
255 : size_t i = 0;
256 : Face* f1 = NULL;
257 : // find the first face that is assigned to a subset
258 0 : for(; i < vFaces.size(); ++i)
259 : {
260 0 : if(shFace.get_subset_index(vFaces[i]) != -1){
261 0 : f1 = vFaces[i];
262 : // if the iteration leaves with a break, we have to increase i manually
263 0 : ++i;
264 0 : break;
265 : }
266 : }
267 :
268 : // compare with others. only check the ones that are assigned to a subset.
269 : // if only one face is in a subset, we'll need a line anyway
270 : bool gotTwoSubsetFaces = false;
271 0 : for(; i < vFaces.size(); ++i)
272 : {
273 0 : if(shFace.get_subset_index(vFaces[i]) == -1)
274 0 : continue;
275 :
276 : gotTwoSubsetFaces = true;
277 0 : if(shFace.get_subset_index(f1) != shFace.get_subset_index(vFaces[i]))
278 : {
279 0 : LineSel.select(*EIter);
280 : break;
281 : }
282 : }
283 :
284 0 : if((!gotTwoSubsetFaces) && f1){
285 : // the edge is adjacent to only one subset-face.
286 : // it thus has to be a line
287 0 : LineSel.select(*EIter);
288 : }
289 : }
290 : }
291 :
292 0 : return true;
293 0 : }
294 :
295 :
296 :
297 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
298 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
299 : // CollectInnerVertices
300 : /// collects inner-vertices in a selector
301 0 : static bool CollectInnerVertices(Grid& grid, VertexSelector& InnVrtSel, VertexSelector& SurfVrtSel)
302 : {
303 : // iterate through all grid-vertices and select them with VertexSelector InnVrtSel
304 0 : for(VertexIterator VIter = grid.vertices_begin(); VIter != grid.vertices_end(); ++VIter)
305 : {
306 0 : InnVrtSel.select(*VIter);
307 : }
308 :
309 : // iterate through all surface-vertices and deselect them in VertexSelector InnVrtSel, so that there only remain
310 : // the inner vertices
311 0 : for(VertexIterator VIter = SurfVrtSel.begin(); VIter != SurfVrtSel.end(); ++VIter)
312 : {
313 0 : InnVrtSel.deselect(*VIter);
314 : }
315 :
316 0 : return true;
317 : }
318 :
319 :
320 :
321 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
322 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
323 : // CollectSurfaceVertices
324 : /// collects surface-vertices in a selector
325 0 : static bool CollectSurfaceVertices(Grid& grid, const SubsetHandler& shFace,
326 : VertexSelector& SurfVrtSel)
327 : {
328 : // iterate through all faces that are assigned to a subset and select them
329 0 : for(FaceIterator iter = grid.faces_begin(); iter != grid.faces_end(); ++iter)
330 : {
331 : Face* f = *iter;
332 0 : if(shFace.get_subset_index(f) != -1)
333 : {
334 0 : for(uint i = 0; i < f->num_vertices(); ++i)
335 0 : SurfVrtSel.select(f->vertex(i));
336 : }
337 : }
338 :
339 0 : return true;
340 : }
341 :
342 :
343 :
344 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
345 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
346 : // CollectAllVerticesForNG
347 : /// collects all vertices in *.ng file order in a selector
348 0 : static bool CollectAllVerticesForNG(Grid& grid, VertexSelector& NgVrtSel,
349 : VertexSelector& SurfVrtSel, VertexSelector& InnVrtSel)
350 : {
351 : // collecting all vertices for an *.ng file requires first to select the surface-vertices and then the inner ones
352 :
353 : // iterate through all surface-vertices and select them
354 0 : for(VertexIterator VIter = SurfVrtSel.begin(); VIter != SurfVrtSel.end(); ++VIter)
355 : {
356 0 : NgVrtSel.select(*VIter);
357 : }
358 :
359 : // iterate through all inner vertices and select them
360 0 : for(VertexIterator VIter = InnVrtSel.begin(); VIter != InnVrtSel.end(); ++VIter)
361 : {
362 0 : NgVrtSel.select(*VIter);
363 : }
364 :
365 0 : return true;
366 : }
367 :
368 : /** THIS METHOD USES Grid::mark
369 : * This method is required since double-triangles have to be avoided
370 : * during triangulation of quadrilaterals which share two edges whith
371 : * another quadrilateral (this e.g. occurs in fracture geometries with neumann
372 : * boundaries).
373 : * Note that only faces in the same subset as q are regarded.
374 : */
375 0 : static int GetFirstRegularVertex(Grid& grid, const SubsetHandler& sh, Face* q)
376 : {
377 0 : grid.begin_marking();
378 : vector<Edge*> edges;
379 : vector<Face*> faces;
380 0 : int si = sh.get_subset_index(q);
381 :
382 : // check whether a connected face exists, which shares two sides with q.
383 : // iterate over all edges, collect associated faces and mark them.
384 : // if a face was already marked, then it is shared by at least
385 : // two sides. We will then immediately mark all associated vertices of
386 : // that face.
387 : CollectAssociated(edges, grid, q);
388 0 : for(size_t i_edge = 0; i_edge != edges.size(); ++i_edge){
389 0 : CollectAssociated(faces, grid, edges[i_edge]);
390 : // add all faces to nbrs. Check whether it is already contained.
391 0 : for(size_t i_face = 0; i_face < faces.size(); ++i_face){
392 0 : Face* f = faces[i_face];
393 0 : if(sh.get_subset_index(f) != si)
394 0 : continue;
395 0 : if(f == q)
396 0 : continue;
397 :
398 : // check whether the face was already encountered
399 0 : if(grid.is_marked(f)){
400 : // it was already there!
401 : // mark all of its vertices
402 0 : for(size_t i_vrt = 0; i_vrt < f->num_vertices(); ++i_vrt)
403 0 : grid.mark(f->vertex(i_vrt));
404 : }
405 : else{
406 : // its there for the first time. mark it.
407 : grid.mark(f);
408 : }
409 : }
410 : }
411 :
412 : // find the first unmarked vertex in q
413 : int firstUnmarked = -1;
414 0 : for(int i_vrt = 0; i_vrt < (int)q->num_vertices(); ++i_vrt){
415 0 : if(!grid.is_marked(q->vertex(i_vrt))){
416 : firstUnmarked = i_vrt;
417 : break;
418 : }
419 : }
420 0 : grid.end_marking();
421 0 : return firstUnmarked;
422 0 : }
423 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
424 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
425 : // WriteLGM
426 : /// writes an *lgm file
427 0 : static bool WriteLGM(Grid& grid,
428 : const char* lgmFilename,
429 : const char* problemName,
430 : const char* lgmName,
431 : int convex,
432 : const SubsetHandler& shFaces,
433 : const SubsetHandler& shVolumes,
434 : EdgeSelector& LineSel,
435 : VertexSelector& SurfVrtSel,
436 : Grid::EdgeAttachmentAccessor<AInt>& aaLineIndex,
437 : Grid::VertexAttachmentAccessor<AInt>& aaSurfVrtIndex,
438 : Grid::VertexAttachmentAccessor<AVector3>& aaPos)
439 : {
440 : // initialization
441 0 : ofstream out(lgmFilename);
442 0 : if(!out)
443 : return false;
444 : out.setf(ios::scientific);
445 : out.precision(24);
446 :
447 : // write the header
448 : out << "#Domain-Info" << endl;
449 0 : out << "name = " << lgmName << endl;
450 0 : out << "problemname = " << problemName << endl;
451 0 : out << "convex = " << convex << endl << endl;
452 :
453 : // write the units
454 : out << "#Unit-Info" << endl;
455 0 : for(int i = 0; i < shVolumes.num_subsets(); ++i)
456 : {
457 : // if the name is empty then write a standard name
458 : // be cautious with the unit indices: id = 0 is reserved for outer space ONLY by UG
459 0 : if(shVolumes.subset_info(i).name.size() > 0)
460 0 : out << "unit " << i + 1 << " " << shVolumes.subset_info(i).name << endl;
461 : else {
462 0 : out << "unit " << i + 1 << " " << "unit_" << i + 1 << endl;
463 : }
464 :
465 : }
466 : out << endl;
467 :
468 : // write the lines
469 : out << "#Line-Info" << endl;
470 0 : for(EdgeIterator iter = LineSel.begin(); iter != LineSel.end(); ++iter)
471 : {
472 : Edge* e = *iter;
473 0 : out << "line " << aaLineIndex[e] << ": points: ";
474 0 : out << aaSurfVrtIndex[e->vertex(0)] << " " << aaSurfVrtIndex[e->vertex(1)] << ";" << endl;
475 : }
476 : out << endl;
477 :
478 : // write the surfaces
479 : out << "#Surface-Info" << endl;
480 : {
481 : // used to find vertices that belong to the i-th subset.
482 0 : VertexSelector tmpSurfVrtSel(grid);
483 0 : for(int i = 0; i < shFaces.num_subsets(); ++i)
484 : {
485 0 : tmpSurfVrtSel.clear();
486 :
487 : // identify left and right unit index of each surface
488 : int tmpLeft, tmpRight;
489 0 : if(!GetRightLeftUnitIndex(tmpRight, tmpLeft, grid, *shFaces.begin<Face>(i), shVolumes))
490 : {
491 : LOG("- GetRightLeftUnitIndex failed during lgm-write.\n");
492 0 : LOG("- In surface " << i << " face 0\n");
493 : LOG("- This can happen due to volume-elements with bad orinentation.\n");
494 : LOG("- IMPLEMENT a geometrical method for fallback!\n");
495 : }
496 :
497 0 : out << "surface " << i << ": left=" << tmpLeft << "; right=" << tmpRight << "; points:";
498 : for(ConstFaceIterator FIter = shFaces.begin<Face>(i);
499 0 : FIter != shFaces.end<Face>(i); ++FIter)
500 : {
501 : Face* f = *FIter;
502 0 : for(uint j = 0; j < f->num_vertices(); ++j)
503 : {
504 0 : Vertex* v = f->vertex(j);
505 0 : if(!tmpSurfVrtSel.is_selected(v))
506 : {
507 : tmpSurfVrtSel.select(v);
508 0 : out << " " << aaSurfVrtIndex[v];
509 : }
510 : }
511 : }
512 :
513 : // write lines
514 : // we have to avoid to write lines twice which lie inside a subset.
515 0 : grid.begin_marking();
516 0 : out << "; lines:";
517 : {
518 : vector<Edge*> vEdges;
519 : for(ConstFaceIterator FIter = shFaces.begin<Face>(i);
520 0 : FIter != shFaces.end<Face>(i); ++FIter)
521 : {
522 0 : CollectEdges(vEdges, grid, *FIter);
523 0 : for(vector<Edge*>::iterator EIter = vEdges.begin();
524 0 : EIter != vEdges.end(); ++EIter)
525 : {
526 0 : Edge* e = *EIter;
527 0 : if(LineSel.is_selected(e) && (!grid.is_marked(e)))
528 : {
529 : grid.mark(e);
530 0 : out << " " << aaLineIndex[e];
531 : }
532 : }
533 : }
534 0 : }
535 0 : grid.end_marking();
536 :
537 : // write triangles
538 0 : out << "; triangles:";
539 : {
540 : vector<Vertex*> vVertices;
541 0 : for(ConstTriangleIterator TIter = shFaces.begin<Triangle>(i);
542 0 : TIter != shFaces.end<Triangle>(i); ++TIter)
543 : {
544 : Triangle* t = *TIter;
545 0 : for(uint j = 0; j < t->num_vertices(); ++j)
546 : {
547 0 : out << " " << aaSurfVrtIndex[t->vertex(j)];
548 : }
549 0 : out << ";";
550 : }
551 0 : }
552 :
553 : // write quadrilaterals as triangles
554 : // special care has to be taken since quadrliaterals can possibly
555 : // be connected at two sides. If one isn't careful one could thus produce
556 : // two identical triangles during triangulation.
557 : {
558 : // we need Grid::mark to check which point is not contained in a neighbour
559 0 : for(ConstQuadrilateralIterator iter = shFaces.begin<Quadrilateral>(i);
560 0 : iter != shFaces.end<Quadrilateral>(i); ++iter)
561 : {
562 : Quadrilateral* q = *iter;
563 :
564 0 : int firstRegular = GetFirstRegularVertex(grid, shFaces, q);
565 :
566 : // if we didn't find one, we can't split the face. Schedule a warning.
567 0 : if(firstRegular == -1){
568 : UG_LOG("Can't split quadrilateral due to too many associated degenerated faces.\n");
569 : }
570 : else{
571 0 : for(int j = firstRegular; j < firstRegular + 3; ++j)
572 0 : out << " " << aaSurfVrtIndex[q->vertex(j%4)];
573 :
574 0 : out << ";";
575 :
576 0 : for(int j = firstRegular + 2; j < firstRegular + 5; ++j)
577 0 : out << " " << aaSurfVrtIndex[q->vertex(j%4)];
578 :
579 0 : out << ";";
580 :
581 : /*
582 : for(uint j = 0; j < 3; ++j)
583 : out << " " << aaSurfVrtIndex[q->vertex(j)];
584 :
585 : out << ";";
586 :
587 : for(uint j = 2; j < 5; ++j)
588 : out << " " << aaSurfVrtIndex[q->vertex(j%4)];
589 :
590 : out << ";";
591 : */
592 : }
593 : }
594 : }
595 : // surface written...
596 : out << endl;
597 : }
598 : }
599 :
600 : // write the vertices position data
601 : out << endl << "#Point-Info" << endl;
602 0 : for(VertexIterator iter = SurfVrtSel.begin(); iter != SurfVrtSel.end(); ++iter)
603 : {
604 0 : out << aaPos[*iter].x() << " " << aaPos[*iter].y() << " " << aaPos[*iter].z() << ";" << endl;
605 : }
606 :
607 : // close out file
608 0 : out.close();
609 :
610 : return true;
611 0 : }
612 :
613 :
614 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
615 : // GetRightLeftUnitIndices
616 : /// identifies the right and left unit index to a given surface
617 0 : static bool GetRightLeftUnitIndex(int& rightIndex, int& leftIndex, Grid& grid, Face* face,
618 : const SubsetHandler& shVolume)
619 : {
620 : // initialization
621 : vector<Volume*> vVolumes;
622 0 : FaceDescriptor fd;
623 :
624 0 : rightIndex = 0;
625 0 : leftIndex = 0;
626 :
627 : // collect all volumes which are adjacent to face into vector<Volume*> vVolumes
628 0 : CollectVolumes(vVolumes, grid, face);
629 : // iterate through all volumes adjacent to the face and identify left and right unit index
630 0 : for(uint j = 0; j < vVolumes.size(); ++j)
631 : {
632 0 : Volume* v = vVolumes[j];
633 :
634 : // find the face in the volume that matches the face
635 0 : for(uint i = 0; i < v->num_faces(); ++i)
636 : {
637 0 : v->face_desc(i, fd);
638 0 : if(CompareVertices(face, &fd))
639 : {
640 : // we found a matching face.
641 : // check whether the orientation is the same as in the face or not
642 0 : int i0 = GetVertexIndex(face, fd.vertex(0));
643 0 : int i1 = GetVertexIndex(face, fd.vertex(1));
644 0 : if(i1 == (i0 + 1) % (int)fd.num_vertices())
645 : {
646 0 : if(rightIndex)
647 : return false;
648 : // the orientation is the same. the volume is on the right.
649 0 : rightIndex = shVolume.get_subset_index(v) + 1;
650 : }
651 : else
652 : {
653 0 : if(leftIndex)
654 : return false;
655 : // the orientation is not the same. the volume is on the left.
656 0 : leftIndex = shVolume.get_subset_index(v) + 1;
657 : }
658 :
659 : break;
660 : }
661 : }
662 :
663 : // stop iterating through all adjacent volumes, if indices have already been assigned
664 0 : if(leftIndex && rightIndex)
665 : break;
666 : }
667 :
668 : return true;
669 0 : }
670 :
671 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
672 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
673 : // WriteNG
674 : /// writes an *ng file
675 0 : static bool WriteNG(Grid& grid,
676 : const SubsetHandler& shFaces,
677 : const SubsetHandler& shVolumes,
678 : const char* ngFilename,
679 : VertexSelector& SurfVrtSel,
680 : VertexSelector& InnVrtSel,
681 : VertexSelector& NgVrtSel,
682 : EdgeSelector& LineSel,
683 : Grid::EdgeAttachmentAccessor<AInt>& aaLineIndex,
684 : Grid::VertexAttachmentAccessor<AInt>& aaInnVrtIndex,
685 : Grid::VertexAttachmentAccessor<AInt>& aaNgVrtIndex,
686 : Grid::FaceAttachmentAccessor<AInt>& aaFaceIndex,
687 : Grid::VertexAttachmentAccessor<AVector3>& aaPos)
688 : {
689 : // initialization
690 : vector<Edge*> vEdges;
691 : vector<Face*> vFaces;
692 0 : vector<bool> faceFlags(shFaces.num_subsets());
693 :
694 0 : ofstream out(ngFilename);
695 0 : if(!out)
696 : return false;
697 : out.setf(ios::scientific);
698 : out.precision(24);
699 : //UG_LOG("INFO: WRITING DEBUG INFO TO .ng FILE!\n");//ONLY FOR DEBUG
700 : //size_t nodeCount = 0;//ONLY FOR DEBUG
701 :
702 : // write the boundary nodes section
703 : out << "# boundary nodes" << endl;
704 0 : for(VertexIterator VIter = SurfVrtSel.begin(); VIter != SurfVrtSel.end(); ++VIter)
705 : {
706 : Vertex* v = *VIter;
707 :
708 0 : for(int i = 0; i < shFaces.num_subsets(); ++i)
709 : {
710 : faceFlags[i] = false;
711 : }
712 :
713 : //out << "# node-id: " << nodeCount << endl;//ONLY FOR DEBUG
714 : //++nodeCount;//ONLY FOR DEBUG
715 :
716 : // write the boundary-vertex position data
717 0 : out << "B " << aaPos[v].x() << " " << aaPos[v].y() << " " << aaPos[v].z() << endl;
718 0 : out << " ";
719 :
720 : // iterate through all lines
721 0 : for(EdgeIterator EIter = LineSel.begin(); EIter != LineSel.end(); ++EIter)
722 : {
723 : Edge* e = *EIter;
724 :
725 : // check if the current boundary-vertex is also a line-vertex
726 0 : if(GetVertexIndex(e, v) != -1)
727 : {
728 : // write the line-index and line-vertex-index (corresponds to local position)
729 0 : out << " L " << aaLineIndex[e] << " " << GetVertexIndex(e, v);
730 : }
731 : }
732 :
733 : // surface data
734 : // iterate through all associated boundary-faces of vertex v
735 0 : for(Grid::AssociatedFaceIterator FIter = grid.associated_faces_begin(v);
736 0 : FIter != grid.associated_faces_end(v); ++FIter)
737 : {
738 0 : Face* f = *FIter;
739 :
740 : // write the surface section including its index and the triangle index of
741 : // ONE triangle-representative (per face-subset) associated to vertex v
742 0 : if(shFaces.get_subset_index(f) != -1)
743 : {
744 0 : if(faceFlags[shFaces.get_subset_index(f)] == false)
745 : {
746 : // write local coordinates of vertex v on the triangle-representative
747 0 : int faceInd = aaFaceIndex[f];
748 0 : int vrtInd = GetVertexIndex(f, v);
749 : vector2 vCoord(0, 0);
750 0 : if(f->num_vertices() == 4){
751 0 : int firstRegular = GetFirstRegularVertex(grid, shFaces, f);
752 0 : if(firstRegular == -1){
753 : UG_LOG("Can't split quadrilateral due to too many associated degenerated faces.\n");
754 : continue;
755 : }
756 0 : vrtInd = (vrtInd + 4 - firstRegular)%4;
757 : }
758 :
759 0 : switch(vrtInd)
760 : {
761 0 : case 0: vCoord = vector2(0, 0); break;
762 0 : case 1: vCoord = vector2(1.0, 0); break;
763 0 : case 2: vCoord = vector2(0, 1.0); break;
764 0 : case 3: vCoord = vector2(1.0, 0); faceInd++; break; //index into the second sub-triangle
765 : }
766 :
767 0 : out << " S " << shFaces.get_subset_index(f) << " " << faceInd;
768 : out << " " << vCoord.x() << " " << vCoord.y();
769 :
770 : faceFlags[shFaces.get_subset_index(f)] = true;
771 : }
772 : }
773 : }
774 :
775 : out << ";" << endl;
776 : }
777 :
778 : // write inner nodes
779 : out << endl;
780 : out << "# inner nodes" << endl;
781 :
782 0 : if(InnVrtSel.num() > 0)
783 : {
784 : // if there are inner vertices, iterate through all of them
785 0 : for(VertexIterator VIter = InnVrtSel.begin(); VIter != InnVrtSel.end(); ++VIter)
786 : {
787 : Vertex* v = *VIter;
788 : //out << "# node-id: " << nodeCount << endl;//ONLY FOR DEBUG
789 : //++nodeCount;//ONLY FOR DEBUG
790 : // write the position data of the inner vertices
791 0 : out << "I " << aaPos[v].x() << " " << aaPos[v].y() << " " << aaPos[v].z() << ";" << endl;
792 : }
793 : }
794 :
795 : // write elements
796 : out << endl;
797 : out << "# elements" << endl;
798 :
799 0 : for(int i = 0; i < shVolumes.num_subsets(); ++i)
800 : {
801 : // iterate through all volumes
802 : for(ConstVolumeIterator VIter = shVolumes.begin<Volume>(i);
803 0 : VIter != shVolumes.end<Volume>(i); ++VIter)
804 : {
805 : Volume* v = *VIter;
806 :
807 : // be cautious with the unit indices: id = 0 is reserved for outer space ONLY by UG
808 0 : out << " E " << i + 1;
809 :
810 : // iterate through all element vertices
811 0 : for(uint j = 0; j < v->num_vertices(); ++j)
812 : {
813 : // write volume-vertex-index (in *.ng file order)
814 0 : out << " " << aaNgVrtIndex[v->vertex(j)];
815 : }
816 :
817 : // collect all associated boundary-faces of volume v and write their *.ng file indices
818 : vFaces.clear();
819 0 : CollectFaces(vFaces, grid, v);
820 0 : for(vector<Face*>::iterator iter = vFaces.begin(); iter != vFaces.end(); ++iter)
821 : {
822 0 : Face* f = *iter;
823 0 : if(shFaces.get_subset_index(f) != -1)
824 : {
825 0 : out << " F " << aaNgVrtIndex[f->vertex(0)];
826 0 : for(uint j = 1; j < f->num_vertices(); ++j)
827 0 : out << " " << aaNgVrtIndex[f->vertex(j)];
828 : }
829 : }
830 :
831 : out << ";" << endl;
832 : }
833 : }
834 :
835 0 : out.close();
836 : return true;
837 0 : }
838 :
839 :
840 : ////////////////////////////////////////////////////////////////////////////////////////////////
841 : ////////////////////////////////////////////////////////////////////////////////////////////////
842 : // 2D EXPORT
843 : ////////////////////////////////////////////////////////////////////////////////////////////////
844 : // ExportMeshToUG_2D and helper functions
845 0 : static bool FaceIsOnRightSide(Face* f, Edge* e)
846 : {
847 : // If the vertices in e are in the same order as the vertices in f,
848 : // then f is on the left side of e (since vertices are specified counter-clockwise).
849 0 : size_t ind1 = GetVertexIndex(f, e->vertex(0));
850 0 : size_t ind2 = GetVertexIndex(f, e->vertex(1));
851 :
852 0 : if(ind2 == (ind1 + 1) % f->num_vertices())
853 0 : return false;
854 :
855 : return true;
856 : }
857 :
858 0 : bool ExportGridToUG_2D(Grid& grid, const char* fileName, const char* lgmName,
859 : const char* problemName, int convex,
860 : SubsetHandler* psh)
861 : {
862 0 : string lgmFileName = fileName;
863 0 : lgmFileName.append(".lgm");
864 :
865 0 : string ngFileName = fileName;
866 0 : ngFileName.append(".ng");
867 :
868 : // open the file
869 0 : ofstream out(lgmFileName.c_str());
870 0 : if(!out)
871 : {
872 : LOG("Failure in ExportGridToUG_2D: couldn't open " << lgmFileName << " for write" << endl);
873 : return false;
874 : }
875 :
876 : out.precision(12);
877 :
878 : // vectors are used to collect associated elements
879 : vector<Face*> vFaces;
880 :
881 : // Position accessor
882 : Grid::VertexAttachmentAccessor<APosition> aaPos(grid, aPosition);
883 :
884 : // initialise all indices with -1
885 : AInt aInt;
886 : grid.attach_to_vertices(aInt);
887 : Grid::VertexAttachmentAccessor<AInt> aaInt(grid, aInt);
888 : SetAttachmentValues(aaInt, grid.vertices_begin(), grid.vertices_end(), -1);
889 :
890 : // write the header
891 : out << "#Domain-Info" << endl;
892 0 : out << "name = " << lgmName << endl;
893 0 : out << "problemname = " << problemName << endl;
894 0 : out << "convex = " << convex << endl << endl;
895 :
896 : bool bUnitsSupplied;
897 : out << "#Unit-Info" << endl;
898 :
899 0 : if(psh)
900 : {
901 : // there are units
902 : bUnitsSupplied = true;
903 :
904 : // write the units
905 : bool bGotEmpty = false;// helps us to print a warning if required.
906 0 : for(int i = 0; i < psh->num_subsets(); i++)
907 : {
908 : // units that do not contain faces are ignored.
909 0 : if(psh->num<Face>(i) == 0){
910 : bGotEmpty = true;
911 0 : continue;
912 : }
913 :
914 0 : if(bGotEmpty){
915 : // schedule a warning
916 : UG_LOG("WARNING in ExportGridToUG_2D: Empty units between filled ones.\n");
917 : }
918 :
919 0 : SubsetInfo& sh = psh->subset_info(i);
920 : string unitName = sh.name;
921 0 : size_t k = sh.name.find(".");
922 0 : if(k != string::npos)
923 0 : unitName.erase(k, unitName.size() - k);
924 0 : if(unitName.size() > 0)
925 0 : out << "unit " << i+1 << " " << unitName.c_str() << endl;
926 : else
927 0 : out << "unit " << i+1 << " unit_" << i+1 << endl;
928 : }
929 : }
930 : else
931 : {
932 : // no subsets are attached. This means we'll have to create one unit for all the triangles
933 : bUnitsSupplied = false;
934 : out << "unit 1 unit_1" << endl;
935 : }
936 : out << endl;
937 :
938 : // determine the vertices that go into the lgm file and assign indices to them.
939 : // This are all vertices wich lie on a border edge or wich are connected to triangles or quads of different subsets.
940 : {
941 : int pointIndexCounter = 0;
942 :
943 : for(VertexIterator iter = grid.vertices_begin();
944 0 : iter != grid.vertices_end(); iter++)
945 : {
946 0 : if(IsBoundaryVertex2D(grid, *iter))
947 : {
948 0 : aaInt[*iter] = pointIndexCounter++;
949 : }
950 0 : else if(bUnitsSupplied)
951 : {
952 : // check if the point is connected to triangles and quads of different subsets
953 : int subsetIndex = -1;
954 :
955 0 : CollectFaces(vFaces, grid, *iter);
956 :
957 0 : if(!vFaces.empty())
958 0 : subsetIndex = psh->get_subset_index(vFaces[0]);
959 :
960 0 : for(size_t i = 1; i < vFaces.size(); ++i){
961 0 : if(psh->get_subset_index(vFaces[i]) != subsetIndex){
962 0 : aaInt[*iter] = pointIndexCounter++;
963 0 : break;
964 : }
965 : }
966 : }
967 : }
968 : }
969 :
970 : // write the lines
971 : out << "#Line-Info" << endl;
972 : // we attach an int value to the edges wich describes their line-index. -1 if the edge is not a line
973 : // initialise line indices with -1
974 : grid.attach_to_edges(aInt);
975 : Grid::EdgeAttachmentAccessor<AInt> aaLineInt(grid, aInt);
976 : SetAttachmentValues(aaLineInt, grid.edges_begin(), grid.edges_end(), -1);
977 :
978 : int numLines = 0;
979 : vector<vector<Vertex*> > lines;
980 : {
981 : // each edge which is connected to two different subsets or which is
982 : // a boundary edge has to be written as a line
983 0 : if(bUnitsSupplied)
984 : {
985 :
986 : // write the edge-subsets as lines
987 0 : for(int i = 0; i < psh->num_subsets(); ++i)
988 : {
989 : // at this point we assume that lines are oriented
990 0 : if(!psh->empty<Edge>(i)){
991 : // we have to follow the polygonal chain in subset i.
992 : // We assume that the chain is regular
993 : // create a callback for this subset, since we use it several times
994 0 : IsInSubset cbIsInSubset(*psh, i);
995 : std::pair<Vertex*, Edge*> curSec =
996 0 : GetFirstSectionOfPolyChain(grid, psh->begin<Edge>(i),
997 0 : psh->end<Edge>(i),
998 : cbIsInSubset);
999 :
1000 0 : if(curSec.first == NULL)
1001 0 : continue;
1002 :
1003 : // find left and right unit
1004 : // curSec.second contains the first edge of the polychain.
1005 0 : CollectFaces(vFaces, grid, curSec.second);
1006 :
1007 0 : if(vFaces.size() < 1)
1008 0 : continue;
1009 :
1010 : // we now have to check which subset lies on which side of the line
1011 0 : int subLeft = psh->get_subset_index(vFaces[0]) + 1;
1012 : int subRight = 0;
1013 0 : if(vFaces.size() > 1)
1014 0 : subRight = psh->get_subset_index(vFaces[1]) + 1;
1015 :
1016 0 : if(FaceIsOnRightSide(vFaces[0], curSec.second))
1017 : swap(subLeft, subRight);
1018 :
1019 : // since the edge and the first line segment can be flipped,
1020 : // we eventually have to swap again
1021 0 : if(curSec.first != curSec.second->vertex(0))
1022 : swap(subLeft, subRight);
1023 :
1024 0 : out << "line " << numLines << ": left="<< subLeft
1025 0 : << "; right=" << subRight << "; points:";
1026 :
1027 : size_t numEdgesInLine = 0;
1028 0 : lines.push_back(vector<Vertex*>());
1029 : vector<Vertex*>& curLine = lines.back();
1030 :
1031 0 : while(curSec.first)
1032 : {
1033 0 : Vertex* curVrt = curSec.first;
1034 : Edge* curEdge = curSec.second;
1035 :
1036 : // write the vertex index
1037 0 : if(aaInt[curVrt] == -1){
1038 0 : throw(UGError("RegularEdge subsets have not been assigned correctly."));
1039 : }
1040 :
1041 0 : if(aaLineInt[curEdge] != -1){
1042 0 : throw(UGError("Implementation error: One edge is in different lines."));
1043 : }
1044 :
1045 0 : aaLineInt[curEdge] = numLines;
1046 0 : ++numEdgesInLine;
1047 :
1048 0 : out << " " << aaInt[curVrt];
1049 0 : curLine.push_back(curVrt);
1050 :
1051 : std::pair<Vertex*, Edge*> nextSec =
1052 0 : GetNextSectionOfPolyChain(grid, curSec, cbIsInSubset);
1053 :
1054 : // check whether the next section is still valid.
1055 : // if not, write the last vertex and exit.
1056 0 : if(!nextSec.first){
1057 0 : out << " " << aaInt[GetConnectedVertex(curEdge, curVrt)];
1058 0 : curLine.push_back(GetConnectedVertex(curEdge, curVrt));
1059 0 : break;
1060 : }
1061 0 : else if(aaLineInt[nextSec.second] != -1){
1062 : // the edge has already been considered
1063 0 : out << " " << aaInt[nextSec.first];
1064 0 : curLine.push_back(nextSec.first);
1065 : break;
1066 : }
1067 :
1068 : curSec = nextSec;
1069 : }
1070 : out << ";" << endl;
1071 0 : ++numLines;
1072 :
1073 : // make sure that all edges in the subset have been written to the line
1074 0 : if(numEdgesInLine != psh->num<Edge>(i)){
1075 0 : stringstream ss;
1076 : ss << "Couldn't write all edges of subset " << i << " to a line. "
1077 : << "This can happen if the subset is not a regular polychain. "
1078 0 : << "It is either separated into different parts or irregular (contains junctions).";
1079 :
1080 0 : throw(UGError(ss.str().c_str()));
1081 0 : }
1082 : }
1083 : }
1084 :
1085 : /*
1086 : for(EdgeIterator iter = grid.edges_begin(); iter != grid.edges_end(); iter++)
1087 : {
1088 : Edge* e = *iter;
1089 : // check for adjacent faces of different subset indices
1090 : CollectFaces(vFaces, grid, e);
1091 :
1092 : if(vFaces.size() == 2)
1093 : {
1094 : if(psh->get_subset_index(vFaces[0]) != psh->get_subset_index(vFaces[1]))
1095 : {
1096 : // e is a line.
1097 : // assign the index
1098 : aaLineInt[e] = numLines;
1099 :
1100 : // write the line
1101 : int subLeft = psh->get_subset_index(vFaces[0]) + 1;
1102 : int subRight = psh->get_subset_index(vFaces[1]) + 1;
1103 : if(!FaceIsOnRightSide(vFaces[0], e))
1104 : swap(subLeft, subRight);
1105 :
1106 : out << "line " << numLines << ": left="<< subLeft
1107 : << "; right=" << subRight << "; points: "
1108 : << aaInt[e->vertex(0)] << " " << aaInt[e->vertex(1)] << ";" << endl;
1109 :
1110 : numLines++;
1111 :
1112 : // endvertices of lines have to be marked as boundary-vertices
1113 : assert((aaInt[e->vertex(0)] != -1) && "This point should be marked as boundary vertex!");
1114 : assert((aaInt[e->vertex(1)] != -1) && "This point should be marked as boundary vertex!");
1115 : }
1116 : }
1117 : }
1118 : */
1119 : }
1120 : else{
1121 : // check for boundary lines
1122 0 : for(EdgeIterator iter = grid.edges_begin(); iter != grid.edges_end(); iter++)
1123 : {
1124 : Edge* e = *iter;
1125 0 : CollectFaces(vFaces, grid, e);
1126 : // check if it is a boundary edge
1127 0 : if(vFaces.size() == 1)
1128 : {
1129 : // assign the index
1130 0 : aaLineInt[e] = numLines;
1131 :
1132 : int unitIndex = 1;
1133 : if(bUnitsSupplied)
1134 : unitIndex = psh->get_subset_index(vFaces[0]) + 1;
1135 :
1136 : int subLeft = 0;
1137 : int subRight = 0;
1138 0 : if(FaceIsOnRightSide(vFaces[0], e))
1139 : subRight = unitIndex;
1140 : else
1141 : subLeft = unitIndex;
1142 :
1143 0 : out << "line " << numLines << ": left="<< subLeft
1144 0 : << "; right=" << subRight << "; points: "
1145 0 : << aaInt[e->vertex(0)] << " " << aaInt[e->vertex(1)] << ";" << endl;
1146 :
1147 0 : numLines++;
1148 :
1149 : // endvertices of lines have to be marked as boundary-vertices
1150 : assert((aaInt[e->vertex(0)] != -1) && "This point should be marked as boundary vertex!");
1151 : assert((aaInt[e->vertex(1)] != -1) && "This point should be marked as boundary vertex!");
1152 : }
1153 : }
1154 : }
1155 : }
1156 :
1157 : // write the vertices position data
1158 : {
1159 : out << endl << "#Point-Info" << endl;
1160 0 : for(VertexIterator iter = grid.vertices_begin(); iter != grid.vertices_end(); iter++)
1161 : {
1162 : // only write the point if it lies on a line
1163 0 : if(aaInt[*iter] != -1)
1164 0 : out << aaPos[*iter].x() << " " << aaPos[*iter].y() << ";" << endl;
1165 : }
1166 : }
1167 :
1168 : // lgm-write done
1169 0 : out.close();
1170 :
1171 : // now we have to write the ng file
1172 : // open the file
1173 0 : out.open(ngFileName.c_str());
1174 :
1175 : // enable scientific number format�
1176 : //out.setf(ios::scientific);
1177 :
1178 0 : if(!out)
1179 : {
1180 : LOG("Failure in ExportGridToUG_2D: couldn't open " << ngFileName << " for write" << endl);
1181 : grid.detach_from_vertices(aInt);
1182 : grid.detach_from_edges(aInt);
1183 : return false;
1184 : }
1185 :
1186 : // write the nodes
1187 : {
1188 : int numNGVertexs = 0;
1189 : // first we'll write the boundary vertices
1190 : {
1191 0 : for(VertexIterator iter = grid.vertices_begin(); iter != grid.vertices_end(); iter++)
1192 : {
1193 : Vertex* p = *iter;
1194 : // check if p is a boundary node
1195 0 : if(aaInt[p] != -1)
1196 : {
1197 : // it is. write it to the file
1198 0 : out << "B ";
1199 : // position:
1200 0 : out << aaPos[p].x() << " " << aaPos[p].y();
1201 : // connected lines
1202 : // make sure that each is only used once
1203 : set<int> encounteredLines;
1204 0 : for(Grid::AssociatedEdgeIterator eIter = grid.associated_edges_begin(p);
1205 0 : eIter != grid.associated_edges_end(p); ++eIter)
1206 : {
1207 0 : if(aaLineInt[(*eIter)] != -1){
1208 : // get the local coordinate of the vertex on its associated line
1209 0 : int lineIndex = aaLineInt[(*eIter)];
1210 0 : if(encounteredLines.count(lineIndex) > 0)
1211 0 : continue;// the line was already encountered.
1212 :
1213 0 : vector<Vertex*>& curLine = lines.at(lineIndex);
1214 : assert(curLine.size() > 1 && "A line has to contain at least 2 vertices.");
1215 : int localCoord = -1;
1216 0 : for(size_t i = 0; i < curLine.size(); ++i){
1217 0 : if(curLine[i] == p)
1218 0 : localCoord = i;
1219 : }
1220 0 : if(localCoord >= 0){
1221 0 : out << " L " << lineIndex << " " << localCoord;
1222 : encounteredLines.insert(lineIndex);
1223 : }
1224 : else{
1225 0 : UG_LOG("WARNING in ExportGridToUG_2D: "
1226 : << "Couldn't find local coordinate of boundary vertex at"
1227 : << aaPos[p] << endl);
1228 : }
1229 : }
1230 : }
1231 : // done
1232 : out << ";" << endl;
1233 0 : numNGVertexs++;
1234 : }
1235 : }
1236 : }
1237 :
1238 : // write the inner vertices
1239 : {
1240 0 : for(VertexIterator iter = grid.vertices_begin(); iter != grid.vertices_end(); iter++)
1241 : {
1242 : Vertex* p = *iter;
1243 0 : if(aaInt[p] == -1)
1244 : {
1245 : // write the point
1246 0 : out << "I " << aaPos[p].x() << " " << aaPos[p].y() << ";" << endl;
1247 0 : aaInt[p] = numNGVertexs++;
1248 : }
1249 : }
1250 : }
1251 : }
1252 :
1253 : // write the elements
1254 : {
1255 : // loop through all the triangles
1256 0 : for(FaceIterator iter = grid.faces_begin(); iter != grid.faces_end(); iter++)
1257 : {
1258 : Face* f = *iter;
1259 :
1260 : int unitIndex = 1;
1261 0 : if(bUnitsSupplied)
1262 0 : unitIndex = psh->get_subset_index(f) + 1;
1263 :
1264 0 : out << "E " << unitIndex;
1265 0 : for(size_t i = 0; i < f->num_vertices(); ++i)
1266 0 : out << " " << aaInt[f->vertex(i)];
1267 :
1268 : // check if one of its edges is a line.
1269 0 : for(size_t i = 0; i < f->num_edges(); ++i)
1270 : {
1271 0 : Edge* e = grid.get_edge(f, i);
1272 0 : if(aaLineInt[e] != -1)
1273 : {
1274 0 : out << " S " << aaInt[e->vertex(0)] << " " << aaInt[e->vertex(1)];
1275 : }
1276 : }
1277 : // done
1278 : out << ";" << endl;
1279 : }
1280 : }
1281 :
1282 : // write complete
1283 0 : out.close();
1284 :
1285 : // clean up
1286 : grid.detach_from_vertices(aInt);
1287 : grid.detach_from_edges(aInt);
1288 :
1289 : return true;
1290 0 : }
1291 :
1292 : }// end of namespace
|
