Line data Source code
1 : /*
2 : * Copyright (c) 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 <algorithm>
34 : #include <vector>
35 : #include "horizontal_layers_mesher.h"
36 : #include "lib_grid/callbacks/callbacks.h"
37 : #include "lib_grid/algorithms/extrusion/extrude.h"
38 : #include "lib_grid/algorithms/geom_obj_util/face_util.h"
39 : #include "lib_grid/algorithms/geom_obj_util/vertex_util.h"
40 : #include "lib_grid/file_io/file_io_asc.h"
41 : #include "lib_grid/iterators/associated_elements_iterator.h"
42 : #include "lib_grid/iterators/lg_for_each.h"
43 :
44 : using namespace std;
45 :
46 : namespace ug{
47 :
48 : ////////////////////////////////////////////////////////////////////////////////
49 0 : void MeshLayerBoundaries(Grid& grid, const RasterLayers& layers,
50 : Grid::VertexAttachmentAccessor<AVector3> aaPos,
51 : ISubsetHandler* pSH)
52 : {
53 : int defSubInd = -1;
54 0 : if(pSH)
55 : defSubInd = pSH->get_default_subset_index();
56 :
57 0 : for(size_t lvl = 0; lvl < layers.size(); ++lvl){
58 0 : if(pSH)
59 0 : pSH->set_default_subset_index((int)lvl);
60 0 : CreateGridFromFieldBoundary(grid, layers.heightfield(lvl), aaPos);
61 : }
62 :
63 0 : if(pSH)
64 0 : pSH->set_default_subset_index(defSubInd);
65 0 : }
66 :
67 : ////////////////////////////////////////////////////////////////////////////////
68 0 : void MeshLayers(Grid& grid, const RasterLayers& layers,
69 : Grid::VertexAttachmentAccessor<AVector3> aaPos,
70 : ISubsetHandler* pSH)
71 : {
72 : int defSubInd = -1;
73 0 : if(pSH)
74 : defSubInd = pSH->get_default_subset_index();
75 :
76 0 : for(size_t lvl = 0; lvl < layers.size(); ++lvl){
77 0 : if(pSH)
78 0 : pSH->set_default_subset_index((int)lvl);
79 0 : CreateGridFromField(grid, layers.heightfield(lvl), aaPos);
80 : }
81 :
82 0 : if(pSH)
83 0 : pSH->set_default_subset_index(defSubInd);
84 0 : }
85 :
86 :
87 : ////////////////////////////////////////////////////////////////////////////////
88 : ////////////////////////////////////////////////////////////////////////////////
89 : // EXTRUDE LAYERS
90 : struct ConnectedToOneMarkedVrt{
91 : ConnectedToOneMarkedVrt(Grid& grid) : m_grid(grid) {}
92 0 : bool operator() (Edge* e) const{
93 0 : return (m_grid.is_marked(e->vertex(0)) || m_grid.is_marked(e->vertex(1)))
94 0 : && !(m_grid.is_marked(e->vertex(0)) && m_grid.is_marked(e->vertex(1)));
95 : }
96 : Grid& m_grid;
97 : };
98 :
99 0 : void ExtrudeLayers (
100 : Grid& grid,
101 : const RasterLayers& layers,
102 : Grid::VertexAttachmentAccessor<AVector3> aaPos,
103 : ISubsetHandler& sh,
104 : bool allowForTetsAndPyras,
105 : const ANumber* aRelZOut)
106 : {
107 0 : if(allowForTetsAndPyras){
108 0 : ExtrudeLayersMixed(grid, layers, aaPos, sh, aRelZOut);
109 0 : return;
110 : }
111 :
112 0 : UG_COND_THROW(layers.size() < 2, "At least 2 layers are required to perform extrusion!");
113 :
114 0 : grid.begin_marking();
115 :
116 : vector<Vertex*> curVrts; // list of vertices that are considered for extrusion
117 : vector<Vertex*> tmpVrts; // used to determine the set of vertices that can be extruded
118 : vector<Vertex*> smoothVrts;
119 : vector<Face*> curFaces; // list of faces that are considered for extrusion
120 : vector<Face*> tmpFaces; // used to determine the set of faces that can be extruded
121 : vector<number> vrtHeightVals; // here we'll record height-values at which new vertices will be placed
122 : vector<number> volHeightVals; // here we'll record height-values at the center of new volumes
123 : vector<number> curUpperLayerHeight;
124 : vector<int> volSubsetInds;
125 : vector<Volume*> newVols;
126 : queue<Volume*> volCandidates;
127 :
128 : Grid::edge_traits::callback cbIsMarked = IsMarked(grid);
129 : Grid::edge_traits::callback cbConnectedToOneMarkedVrt = ConnectedToOneMarkedVrt(grid);
130 0 : AssocElemIter<Vertex, Edge> assocVrtEdgeIterMarkedEdge(cbIsMarked);
131 0 : AssocElemIter<Vertex, Edge> assocVrtEdgeIterOneMarked(cbConnectedToOneMarkedVrt);
132 0 : AssocElemIter<Face, Edge> assocFaceEdgeIter;
133 0 : AssocElemIter<Volume, Edge> assocVolEdgeIter(cbConnectedToOneMarkedVrt);
134 :
135 0 : grid.reserve<Vertex>(grid.num_vertices() * layers.size());
136 0 : grid.reserve<Volume>(grid.num_faces() * layers.size() - 1);
137 :
138 : // todo: this accessor is primarily used during smoothing. Maybe it can be removed?
139 : ANumber aHeight;
140 : grid.attach_to_vertices(aHeight);
141 : Grid::VertexAttachmentAccessor<ANumber> aaHeight(grid, aHeight);
142 :
143 : Grid::VertexAttachmentAccessor<ANumber> aaRelZOut;
144 0 : if(aRelZOut){
145 : ANumber aRelZ = *aRelZOut;
146 0 : if(!grid.has_vertex_attachment(aRelZ))
147 : grid.attach_to_vertices(aRelZ);
148 : aaRelZOut.access(grid, aRelZ);
149 : }
150 :
151 : // we have to determine the vertices that can be projected onto the top of the
152 : // given layers-stack. Only those will be used during extrusion
153 : // all considered vertices will be marked.
154 : const RasterLayers::layer_t& top = layers.top();
155 0 : const int topLayerInd = (int)layers.num_layers() - 1;
156 0 : for(VertexIterator i = grid.begin<Vertex>(); i != grid.end<Vertex>(); ++i){
157 0 : Vertex* v = *i;
158 0 : vector2 c(aaPos[v].x(), aaPos[v].y());
159 : // number val = top.heightfield.interpolate(c);
160 0 : number val = layers.relative_to_global_height(c, topLayerInd);
161 0 : if(val != top.heightfield.no_data_value()){
162 0 : aaPos[v].z() = val;
163 0 : curVrts.push_back(v);
164 0 : curUpperLayerHeight.push_back(val);
165 : grid.mark(v);
166 0 : sh.assign_subset(v, topLayerInd);
167 : }
168 : }
169 :
170 0 : if(curVrts.size() < 3){
171 : UG_LOG("Not enough vertices lie in the region of the surface layer\n");
172 : return;
173 : }
174 :
175 0 : if(aaRelZOut.valid()){
176 0 : for(size_t ivrt = 0; ivrt < curVrts.size(); ++ivrt)
177 0 : aaRelZOut[curVrts[ivrt]] = topLayerInd;
178 : }
179 :
180 : // all faces of the initial triangulation that are only connected to marked
181 : // vertices are considered for extrusion
182 0 : for(FaceIterator fi = grid.begin<Face>(); fi != grid.end<Face>(); ++fi){
183 0 : Face* f = *fi;
184 : bool allMarked = true;
185 0 : Face::ConstVertexArray vrts = f->vertices();
186 0 : const size_t numVrts = f->num_vertices();
187 0 : for(size_t i = 0; i < numVrts; ++i){
188 0 : if(!grid.is_marked(vrts[i])){
189 : allMarked = false;
190 : break;
191 : }
192 : }
193 0 : if(allMarked)
194 0 : curFaces.push_back(f);
195 : }
196 :
197 0 : if(curFaces.empty()){
198 : UG_LOG("Not enough faces lie in the region of the surface layer\n");
199 : return;
200 : }
201 :
202 :
203 0 : tmpVrts.reserve(curVrts.size());
204 0 : smoothVrts.reserve(curVrts.size());
205 0 : curUpperLayerHeight.reserve(curVrts.size());
206 0 : tmpFaces.reserve(curFaces.size());
207 0 : newVols.reserve(curFaces.size());
208 0 : volHeightVals.reserve(curFaces.size());
209 0 : volSubsetInds.reserve(curFaces.size());
210 0 : const int invalidSub = max<int>(sh.num_subsets(), layers.size() + 1);
211 : vector<Volume*> invalidVols;
212 :
213 0 : for(int ilayer = (int)layers.size() - 2; ilayer >= 0; --ilayer){
214 :
215 : tmpVrts.clear();
216 : vrtHeightVals.clear();
217 : tmpFaces.clear();
218 : newVols.clear();
219 : volHeightVals.clear();
220 : volSubsetInds.clear();
221 0 : grid.clear_marks();
222 : // trace rays from the current vertices down through the layers until the
223 : // next valid entry is found. If none is found, the vertex will either be ignored
224 : // from then on (allowForTetsAndPyras == false) or a dummy vertex will be inserted.
225 0 : for(size_t icur = 0; icur < curVrts.size(); ++icur){
226 0 : Vertex* v = curVrts[icur];
227 0 : vector2 c(aaPos[v].x(), aaPos[v].y());
228 0 : number upperVal = curUpperLayerHeight[icur];
229 0 : number height = layers.relative_to_global_height(c, (number) ilayer);
230 :
231 : // bool searching = true;
232 : int curLayer = ilayer;
233 :
234 : // while(searching && curLayer >= 0){
235 : // searching = false;
236 0 : pair<int, number> val = layers.trace_line_down(c, curLayer);
237 : number lowerVal = 0;
238 0 : if(val.first >= 0)
239 0 : lowerVal = layers.relative_to_global_height(c, (number) val.first);
240 :
241 : // if(curLayer == 0 && val.first >= 0
242 : // && upperVal - lowerVal < layers.min_height(curLayer))
243 : // {
244 : // val.first = -1;
245 : // }
246 :
247 0 : if(val.first >= 0){
248 : // if(upperVal - lowerVal < layers.min_height(curLayer)){
249 : // searching = true;
250 : // --curLayer;
251 : // continue;
252 : // }
253 : // if val.first == curLayer height will equal lowerVal. If not,
254 : // a linear interpolation is performed, considering the height-val
255 : // of the current vertex, the layer distance and the target value.
256 0 : tmpVrts.push_back(v);
257 0 : vrtHeightVals.push_back(height);
258 0 : aaHeight[v] = upperVal - lowerVal;//total height of curLayer
259 0 : sh.assign_subset(v, val.first);
260 : grid.mark(v);
261 0 : if(val.first == ilayer){
262 : // UG_LOG("DBG: setting upper value: " << lowerVal << "(old: " << upperVal << ")\n");
263 0 : curUpperLayerHeight[icur] = lowerVal;
264 : }
265 : }
266 : else if(allowForTetsAndPyras){
267 : // we insert a dummy-vertex which will later on allow for easier
268 : // edge-collapses of inner vertical rim edges
269 : tmpVrts.push_back(v);
270 : vrtHeightVals.push_back(height);
271 : aaHeight[v] = upperVal - val.second;//total height of ilayer
272 : sh.assign_subset(v, invalidSub);
273 : grid.mark(v);
274 : }
275 : // }
276 : }
277 : // now find the faces which connect those vertices
278 0 : for(size_t iface = 0; iface < curFaces.size(); ++iface){
279 0 : Face* f = curFaces[iface];
280 :
281 0 : vector3 center = CalculateCenter(f, aaPos);
282 0 : vector2 c(center.x(), center.y());
283 0 : pair<int, number> val = layers.trace_line_down(c, ilayer);
284 0 : if((val.first < 0) && !allowForTetsAndPyras)
285 0 : continue;
286 :
287 : // now check whether all vertices are marked
288 : bool allMarked = true;
289 0 : Face::ConstVertexArray vrts = f->vertices();
290 0 : size_t numVrts = f->num_vertices();
291 0 : for(size_t i = 0; i < numVrts; ++i){
292 0 : if(!grid.is_marked(vrts[i])){
293 : allMarked = false;
294 : break;
295 : }
296 : }
297 :
298 0 : if(allMarked){
299 0 : int subsetIndex = invalidSub;
300 :
301 : // if(val.first >= 0) {
302 : // // subset from center trace down
303 : // subsetIndex = val.first;
304 : // }
305 : // else {
306 : // the highest valid side of the volume determines its subset.
307 : number maxHeight = 0;
308 0 : for(size_t i = 0; i < numVrts; ++i){
309 0 : int si = sh.get_subset_index(vrts[i]);
310 0 : if(si == invalidSub)
311 0 : continue;
312 :
313 0 : number vh = aaHeight[vrts[i]];
314 0 : if(vh > maxHeight){
315 : maxHeight = vh;
316 0 : subsetIndex = si;
317 : }
318 : }
319 : // }
320 :
321 :
322 :
323 0 : pair<int, number> upVal = layers.trace_line_up(c, ilayer+1);
324 0 : if(val.first < 0 || upVal.first < 0){
325 : if(allowForTetsAndPyras){
326 : tmpFaces.push_back(f);
327 : // volSubsetInds.push_back(invalidSub);
328 : volSubsetInds.push_back(subsetIndex);
329 : volHeightVals.push_back(layers.min_height(ilayer));
330 : }
331 : }
332 : else{
333 0 : tmpFaces.push_back(f);
334 : // volSubsetInds.push_back(val.first);
335 0 : volSubsetInds.push_back(subsetIndex);
336 0 : volHeightVals.push_back(upVal.second - val.second);
337 : }
338 : }
339 : }
340 :
341 0 : if(tmpFaces.empty())
342 : break;
343 :
344 : // swap containers and perform extrusion
345 : curVrts.swap(tmpVrts);
346 : curFaces.swap(tmpFaces);
347 0 : Extrude(grid, &curVrts, NULL, &curFaces, vector3(0, 0, 0), aaPos,
348 : EO_DEFAULT, &newVols);
349 :
350 : // assign pre-determined subsets
351 0 : for(size_t ivol = 0; ivol < newVols.size(); ++ivol){
352 0 : sh.assign_subset(newVols[ivol], volSubsetInds[ivol]);
353 0 : if(volSubsetInds[ivol] == invalidSub)
354 0 : invalidVols.push_back(newVols[ivol]);
355 : }
356 :
357 : // set the precalculated height of new vertices
358 0 : for(size_t ivrt = 0; ivrt < curVrts.size(); ++ivrt){
359 0 : aaPos[curVrts[ivrt]].z() = vrtHeightVals[ivrt];
360 : }
361 :
362 0 : if(aaRelZOut.valid()){
363 0 : for(size_t ivrt = 0; ivrt < curVrts.size(); ++ivrt)
364 0 : aaRelZOut[curVrts[ivrt]] = ilayer;
365 : }
366 : }
367 :
368 : // remove unnecessary prisms through edge-collapses and thus introduce pyramids and tetrahedra
369 : if(allowForTetsAndPyras){
370 : ABool aInterface;
371 : grid.attach_to_vertices_dv(aInterface, false, true);
372 : Grid::VertexAttachmentAccessor<ABool> aaInterface(grid, aInterface);
373 :
374 : Grid::volume_traits::secure_container assVols;
375 :
376 : // all triangle-interface-elements shall store 'true' in aaIsInterface
377 : lg_for_each(Face, f, grid){
378 : if(f->num_vertices() != 3)
379 : continue;
380 :
381 : grid.associated_elements(assVols, f);
382 : if(assVols.size() == 1){
383 : lg_for_each_vertex_in_elem(vrt, f){
384 : aaInterface[vrt] = true;
385 : }lg_end_for;
386 : }
387 : else{
388 : int si = -1;
389 : for_each_in_vec(Volume* v, assVols){
390 : if(si == -1)
391 : si = sh.get_subset_index(v);
392 : else if(sh.get_subset_index(v) != si){
393 : lg_for_each_vertex_in_elem(vrt, f){
394 : aaInterface[vrt] = true;
395 : }lg_end_for;
396 : break;
397 : }
398 : }end_for;
399 : }
400 : }lg_end_for;
401 :
402 : // all lower vertices of elements in invalidSub are collapse candidates and
403 : // are thus not considered to be interface elements
404 : Grid::edge_traits::secure_container assEdges;
405 : for_each_in_vec(Volume* vol, invalidVols){
406 : grid.associated_elements(assEdges, vol);
407 : for_each_in_vec(Edge* e, assEdges){
408 : Vertex* v0 = e->vertex(0);
409 : Vertex* v1 = e->vertex(1);
410 :
411 : vector3 dir;
412 : VecSubtract(dir, aaPos[v1], aaPos[v0]);
413 : if((dir.z() > SMALL) && (fabs(dir.x()) < SMALL) && (fabs(dir.y()) < SMALL)){
414 : aaInterface[v0] = false;
415 : }
416 : else if((-dir.z() > SMALL) && (fabs(dir.x()) < SMALL) && (fabs(dir.y()) < SMALL)){
417 : aaInterface[v1] = false;
418 : }
419 : }end_for;
420 : }end_for;
421 :
422 : // all unmarked vertices are collapse candidates
423 : vector<Vertex*> candidates;
424 : lg_for_each(Vertex, vrt, grid){
425 : if(!aaInterface[vrt])
426 : candidates.push_back(vrt);
427 : }lg_end_for;
428 :
429 : // merge each candidate with the next upper vertex.
430 : for_each_in_vec(Vertex* vrt, candidates){
431 : grid.associated_elements(assEdges, vrt);
432 : for_each_in_vec(Edge* e, assEdges){
433 : Vertex* cv = GetConnectedVertex(e, vrt);
434 : if(cv == vrt)
435 : continue;
436 :
437 : vector3 dir;
438 : VecSubtract(dir, aaPos[cv], aaPos[vrt]);
439 : if((dir.z() > SMALL) && (fabs(dir.x()) < SMALL) && (fabs(dir.y()) < SMALL)){
440 : CollapseEdge(grid, e, cv);
441 : break;
442 : }
443 : }end_for;
444 : }end_for;
445 :
446 : // clean up
447 : grid.detach_from_vertices(aInterface);
448 :
449 : // delete invalidSub from the subset handler
450 : if(invalidSub < sh.num_subsets())
451 : sh.erase_subset(invalidSub);
452 : }
453 0 : grid.end_marking();
454 : grid.detach_from_vertices(aHeight);
455 0 : }
456 :
457 :
458 0 : void ExtrudeLayersMixed (
459 : Grid& grid,
460 : const RasterLayers& layers,
461 : Grid::VertexAttachmentAccessor<AVector3> aaPos,
462 : ISubsetHandler& sh,
463 : const ANumber* aRelZOut)
464 : {
465 0 : UG_COND_THROW(layers.size() < 2, "At least 2 layers are required to perform extrusion!");
466 :
467 : vector<Vertex*> curVrts; // list of vertices that are considered for extrusion
468 : vector<Vertex*> nextVrts; // used to determine the set of vertices that can be extruded
469 : vector<Face*> curFaces; // list of faces that are considered for extrusion
470 : vector<Face*> nextFaces; // used to determine the set of faces that can be extruded
471 :
472 0 : grid.reserve<Vertex>(grid.num_vertices() * layers.size());
473 0 : grid.reserve<Volume>(grid.num_faces() * layers.size() - 1);
474 :
475 : AInt aVrtInd;
476 : grid.attach_to_vertices(aVrtInd);
477 : Grid::VertexAttachmentAccessor<AInt> aaVrtInd(grid, aVrtInd);
478 :
479 : Grid::VertexAttachmentAccessor<ANumber> aaRelZOut;
480 0 : if(aRelZOut){
481 : ANumber aRelZ = *aRelZOut;
482 0 : if(!grid.has_vertex_attachment(aRelZ))
483 : grid.attach_to_vertices(aRelZ);
484 : aaRelZOut.access(grid, aRelZ);
485 : }
486 :
487 : // we have to determine the vertices that can be projected onto the top of the
488 : // given layers-stack. Only those will be used during extrusion
489 : // all considered vertices will be marked.
490 : const RasterLayers::layer_t& top = layers.top();
491 0 : const int topLayerInd = (int)layers.num_layers() - 1;
492 0 : for(VertexIterator i = grid.begin<Vertex>(); i != grid.end<Vertex>(); ++i){
493 0 : Vertex* v = *i;
494 0 : vector2 c(aaPos[v].x(), aaPos[v].y());
495 : // number val = layers.relative_to_global_height(c, topLayerInd);
496 0 : number val = layers.heightfield(topLayerInd).interpolate(c, 1);
497 0 : if(val != top.heightfield.no_data_value()){
498 0 : aaPos[v].z() = val;
499 0 : aaVrtInd[v] = (int)curVrts.size();
500 0 : curVrts.push_back(v);
501 0 : sh.assign_subset(v, topLayerInd);
502 : }
503 : else
504 0 : aaVrtInd[v] = -1;
505 : }
506 :
507 0 : if(curVrts.size() < 3){
508 : UG_LOG("Not enough vertices lie in the region of the surface layer\n");
509 : return;
510 : }
511 :
512 0 : if(aaRelZOut.valid()){
513 0 : for(size_t ivrt = 0; ivrt < curVrts.size(); ++ivrt)
514 0 : aaRelZOut[curVrts[ivrt]] = topLayerInd;
515 : }
516 :
517 : // make sure that only triangles are used
518 0 : if(grid.num<Quadrilateral>() > 0){
519 : UG_LOG("WARNING: Converting all quadrilaterals to triangles\n");
520 0 : Triangulate(grid, grid.begin<Quadrilateral>(), grid.end<Quadrilateral>(), &aaPos);
521 : }
522 :
523 : // all faces of the initial triangulation that are only connected to marked
524 : // vertices are considered for extrusion
525 0 : for(FaceIterator fi = grid.begin<Face>(); fi != grid.end<Face>(); ++fi){
526 0 : Face* f = *fi;
527 : bool allConsidered = true;
528 0 : Face::ConstVertexArray vrts = f->vertices();
529 0 : const size_t numVrts = f->num_vertices();
530 0 : for(size_t i = 0; i < numVrts; ++i){
531 0 : if(aaVrtInd[vrts[i]] == -1){
532 : allConsidered = false;
533 : break;
534 : }
535 : }
536 0 : if(allConsidered)
537 0 : curFaces.push_back(f);
538 : }
539 :
540 0 : if(curFaces.empty()){
541 : UG_LOG("Not enough faces lie in the region of the surface layer\n");
542 : return;
543 : }
544 :
545 :
546 0 : nextVrts.reserve(curVrts.size());
547 0 : nextFaces.reserve(curFaces.size());
548 :
549 :
550 0 : for(int ilayer = (int)layers.size() - 1; ilayer > 0; --ilayer){
551 :
552 : nextVrts.clear();
553 : nextFaces.clear();
554 :
555 0 : const int nextLayer = ilayer - 1;
556 :
557 : // trace rays from the current vertices down through the layers until the
558 : // next valid entry is found. If none is found, the vertex will either be ignored
559 : // from then on (allowForTetsAndPyras == false) or a dummy vertex will be inserted.
560 0 : for(size_t icur = 0; icur < curVrts.size(); ++icur){
561 0 : Vertex* v = curVrts[icur];
562 0 : vector2 c(aaPos[v].x(), aaPos[v].y());
563 0 : number curH = aaPos[v].z();
564 : // number nextH = layers.relative_to_global_height(c, (number) nextLayer);
565 0 : number nextH = layers.heightfield(nextLayer).interpolate(c, 1);
566 0 : if(curH - nextH < layers.min_height(nextLayer)){
567 0 : nextVrts.push_back(v);
568 : }
569 : else{
570 0 : Vertex* nextVrt = *grid.create<RegularVertex>(v);
571 0 : aaVrtInd[nextVrt] = aaVrtInd[v];
572 : aaPos[nextVrt] = aaPos[v];
573 0 : aaPos[nextVrt].z() = nextH;
574 0 : sh.assign_subset(nextVrt, nextLayer);
575 0 : nextVrts.push_back(nextVrt);
576 :
577 0 : if(aaRelZOut.valid()){
578 0 : aaRelZOut[nextVrt] = nextLayer;
579 : }
580 : }
581 : }
582 :
583 : // now find the faces which connect those vertices
584 0 : for(size_t iface = 0; iface < curFaces.size(); ++iface){
585 0 : Face* f = curFaces[iface];
586 0 : Face::ConstVertexArray vrts = f->vertices();
587 : int numNew = 0;
588 : int firstSame = -1;
589 : int firstNew = -1;
590 : TriangleDescriptor nextTriDesc;
591 0 : for(int i = 0; i < 3; ++i){
592 0 : int vrtInd = aaVrtInd[vrts[i]];
593 0 : nextTriDesc.set_vertex(i, nextVrts[vrtInd]);
594 0 : if(nextVrts[vrtInd] == curVrts[vrtInd]){
595 0 : if(firstSame == -1)
596 : firstSame = i;
597 : }
598 : else{
599 0 : ++numNew;
600 0 : if(firstNew == -1)
601 : firstNew = i;
602 : }
603 : }
604 :
605 0 : Face* nextFace = f;
606 0 : if(numNew > 0){
607 0 : nextFace = *grid.create<Triangle>(nextTriDesc, f);
608 0 : sh.assign_subset (nextFace, nextLayer);
609 : }
610 0 : nextFaces.push_back(nextFace);
611 :
612 : Volume* newVol = NULL;
613 :
614 0 : switch(numNew){
615 : case 0:// no new volume to be built.
616 : break;
617 0 : case 1: {// build a tetrahedron
618 : TetrahedronDescriptor tetDesc (nextTriDesc.vertex(0),
619 : nextTriDesc.vertex(1),
620 : nextTriDesc.vertex(2),
621 0 : vrts[firstNew]);
622 0 : newVol = *grid.create<Tetrahedron>(tetDesc);
623 0 : } break;
624 :
625 0 : case 2: {// build a pyramid
626 0 : int i0 = (firstSame + 1) % 3;
627 0 : int i1 = (firstSame + 2) % 3;
628 : PyramidDescriptor pyrDesc (nextTriDesc.vertex(i0),
629 : nextTriDesc.vertex(i1),
630 0 : vrts[i1],
631 0 : vrts[i0],
632 0 : vrts[firstSame]);
633 0 : newVol = *grid.create<Pyramid>(pyrDesc);
634 0 : } break;
635 :
636 0 : case 3: {// build a prism
637 : PrismDescriptor prismDesc (nextTriDesc.vertex(0),
638 : nextTriDesc.vertex(1),
639 : nextTriDesc.vertex(2),
640 : vrts[0],
641 : vrts[1],
642 0 : vrts[2]);
643 0 : newVol = *grid.create<Prism>(prismDesc);
644 0 : } break;
645 : }
646 :
647 0 : if(newVol)
648 0 : sh.assign_subset(newVol, nextLayer);
649 : }
650 :
651 : // swap containers and perform extrusion
652 : curVrts.swap(nextVrts);
653 : curFaces.swap(nextFaces);
654 : }
655 :
656 : // clean up
657 : grid.detach_from_vertices(aVrtInd);
658 0 : }
659 :
660 : /// projects the given (surface-) grid to the specified raster
661 0 : void ProjectToLayer(
662 : Grid& grid,
663 : const RasterLayers& layers,
664 : int layerIndex,
665 : Grid::VertexAttachmentAccessor<AVector3> aaPos)
666 : {
667 0 : UG_COND_THROW(layerIndex < 0 || layerIndex >= (int)layers.size(),
668 : "Bad layerIndex in ProjectToLayer");
669 :
670 0 : const RasterLayers::layer_t& layer = layers.layer(layerIndex);
671 0 : for(VertexIterator i = grid.begin<Vertex>(); i != grid.end<Vertex>(); ++i){
672 : Vertex* v = *i;
673 0 : vector2 c(aaPos[v].x(), aaPos[v].y());
674 0 : number val = layers.relative_to_global_height(c, layerIndex);
675 :
676 0 : if(val != layer.heightfield.no_data_value()){
677 0 : aaPos[v].z() = val;
678 : }
679 : }
680 0 : }
681 :
682 :
683 0 : void SnapToHorizontalRaster(
684 : Grid& grid,
685 : const RasterLayers& layers,
686 : Grid::VertexAttachmentAccessor<AVector3> aaPos)
687 : {
688 0 : UG_COND_THROW(layers.empty(), "Can't snap to empty raster!");
689 :
690 0 : const Heightfield& field = layers.top().heightfield;
691 :
692 : for(VertexIterator ivrt = grid.vertices_begin();
693 0 : ivrt != grid.vertices_end(); ++ivrt)
694 : {
695 : Vertex* vrt = *ivrt;
696 : vector3& v = aaPos[vrt];
697 0 : pair<int, int> ci = field.coordinate_to_index(v.x(), v.y());
698 0 : vector2 nv = field.index_to_coordinate(ci.first, ci.second);
699 0 : v.x() = nv.x();
700 0 : v.y() = nv.y();
701 : }
702 0 : }
703 :
704 : }// end of namespace
|
