Line data Source code
1 : /*
2 : * Copyright (c) 2011-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 <boost/function.hpp>
34 : #include <stack>
35 : #include <vector>
36 : #include "expand_layers.h"
37 : #include "lib_grid/algorithms/geom_obj_util/geom_obj_util.h"
38 : #include "lib_grid/callbacks/callbacks.h"
39 : #include "lib_grid/grid/grid_util.h"
40 : //#include "lib_grid/util/simple_algebra/least_squares_solver.h"
41 :
42 : using namespace std;
43 :
44 : namespace ug{
45 :
46 :
47 : /// This class can be used in Element callbacks.
48 : /** Returns true, if the attachmed value in the given element matches a predefined value.*/
49 : template <class TElem, class TAttachmentAccessor>
50 0 : class AttachmentUnequal{
51 : public:
52 : AttachmentUnequal(const TAttachmentAccessor& aa,
53 : const typename TAttachmentAccessor::ValueType& val) :
54 : m_aa(aa), m_val(val) {}
55 :
56 0 : bool operator() (TElem* e) {return m_aa[e] != m_val;}
57 :
58 : private:
59 : TAttachmentAccessor m_aa;
60 : typename TAttachmentAccessor::ValueType m_val;
61 :
62 : };
63 :
64 : /// returns true if the vertex lies on a surface.
65 : /** This method uses Grid::mark.
66 : *
67 : * This method tries to find a closed surface around the vertex.
68 : * Note that it returns true, even if the vertex lies on a
69 : * boundary edge at the same time (this can happen if surfaces
70 : * intersect each other).
71 : *
72 : * Requires the option FACEOPT_AUTOGENERATE_EDGES.
73 : */
74 : /*
75 : static bool VertexLiesOnSurface(Grid& grid, Vertex* vrt,
76 : CB_ConsiderFace funcIsSurfFace)
77 : {
78 : if(!grid.option_is_enabled(FACEOPT_AUTOGENERATE_EDGES)){
79 : UG_LOG("WARNING in VertexLiesOnSurface: autoenabling FACEOPT_AUTOGENERATE_EDGES.\n");
80 : grid.enable_options(FACEOPT_AUTOGENERATE_EDGES);
81 : }
82 :
83 : if(grid.associated_edges_begin(vrt) == grid.associated_edges_end(vrt))
84 : {
85 : return false;
86 : }
87 :
88 : grid.begin_marking();
89 :
90 : stack<Edge*> stk;
91 : vector<Edge*> edges;
92 : vector<Face*> faces;
93 :
94 : // collect associated faces of vrt, which lie on the surface
95 : for(Grid::AssociatedFaceIterator iter = grid.associated_faces_begin(vrt);
96 : iter != grid.associated_faces_end(vrt); ++iter)
97 : {
98 : if(funcIsSurfFace(*iter))
99 : faces.push_back(*iter);
100 : }
101 :
102 : // start with an associated edge of vrt
103 : stk.push(*grid.associated_edges_begin(vrt));
104 : grid.mark(stk.top());
105 :
106 : while(!stk.empty()){
107 : Edge* e = stk.top();
108 :
109 : // find an unmarked associated face
110 : Face* f = NULL;
111 : for(size_t i = 0; i < faces.size(); ++i){
112 : if(!grid.is_marked(faces[i])){
113 : if(FaceContains(faces[i], e)){
114 : f = faces[i];
115 : break;
116 : }
117 : }
118 : }
119 :
120 : if(!f){
121 : stk.pop();
122 : continue;
123 : }
124 :
125 : grid.mark(f);
126 :
127 : // collect associated edges of f
128 : CollectEdges(edges, grid, f);
129 :
130 : // find the edge that is not e an which is connected to vrt
131 : for(size_t i = 0; i < edges.size(); ++i){
132 : Edge* ne = edges[i];
133 : if(ne != e){
134 : if(EdgeContains(ne, vrt)){
135 : // if the edge is marked, then we found a surface
136 : if(grid.is_marked(ne)){
137 : grid.end_marking();
138 : return true;
139 : }
140 : else{
141 : grid.mark(ne);
142 : stk.push(ne);
143 : break;
144 : }
145 : }
146 : }
147 : }
148 : }
149 :
150 : grid.end_marking();
151 : return false;
152 : }
153 : */
154 : /// calculates the normal of the crease vertex vrt on the side of f
155 : /**
156 : * This algorithm uses grid::mark.
157 : * f has to contain vrt. vrt has to have at least two associated
158 : * crease edges.
159 : *
160 : * Note that the resulting normal is not normalized. This is important,
161 : * since rounding errors could else lead to problems with normals which
162 : * have a length of nearly 0.
163 : *
164 : * This algorithm requires the option FACEOPT_AUTOGENERATE_EDGES.
165 : * The option is automatically enabled if required.
166 : */
167 : template <class TAAPosVRT>
168 : typename TAAPosVRT::ValueType
169 0 : CalculateCreaseNormal(Grid& grid, Face* f, Vertex* vrt,
170 : Grid::edge_traits::callback funcIsCreaseEdge,
171 : TAAPosVRT& aaPos)
172 : {
173 0 : if(!grid.option_is_enabled(FACEOPT_AUTOGENERATE_EDGES)){
174 : UG_LOG("WARNING in CalculateCreaseNormal: grid option FACEOPT_AUTOGENERATE_EDGES autoenabled.\n");
175 0 : grid.enable_options(FACEOPT_AUTOGENERATE_EDGES);
176 : }
177 :
178 : typedef typename TAAPosVRT::ValueType vector_t;
179 : vector_t n;
180 : VecSet(n, 0);
181 :
182 0 : grid.begin_marking();
183 :
184 : // we'll use a stack to find all faces on this side of the crease.
185 : // each face in the stack has to be marked.
186 : stack<Face*> stk;
187 : stk.push(f);
188 0 : grid.mark(f);
189 :
190 : // objects for temporary results
191 : vector<Edge*> edges;
192 : vector<Face*> faces;
193 :
194 : // we'll loop while there are faces in the stack
195 0 : while(!stk.empty()){
196 0 : Face* face = stk.top();
197 : stk.pop();
198 :
199 : // the center might be required later on
200 0 : vector_t center = CalculateCenter(face, aaPos);
201 :
202 : // iterate over the edges of face
203 0 : CollectEdges(edges, grid, face);
204 0 : for(size_t i_e = 0; i_e < edges.size(); ++i_e){
205 0 : Edge* e = edges[i_e];
206 0 : if(EdgeContains(e, vrt)){
207 : // check whether e is a crease
208 0 : if(funcIsCreaseEdge(e)){
209 : // we have to add the edges normal to n.
210 : // to make sure that the algorithm works for manifolds too,
211 : // we'll perform a more complicated calculation.
212 :
213 : // project the center onto the edge
214 : vector_t p;
215 0 : DropAPerpendicular(p, center, aaPos[e->vertex(0)], aaPos[e->vertex(1)]);
216 :
217 : // vector from projection to center is the unnormalized normal
218 : vector_t tmpN;
219 : VecSubtract(tmpN, center, p);
220 0 : VecNormalize(tmpN, tmpN);
221 : VecAdd(n, n, tmpN);
222 : }
223 : else{
224 : // since the edge is not a crease edge, we have to add associated unmarked
225 : // faces to the stack
226 0 : CollectFaces(faces, grid, e);
227 0 : for(size_t i = 0; i < faces.size(); ++i){
228 0 : if(!grid.is_marked(faces[i])){
229 : grid.mark(faces[i]);
230 : stk.push(faces[i]);
231 : }
232 : }
233 : }
234 : }
235 : }
236 : }
237 :
238 0 : grid.end_marking();
239 :
240 : //VecNormalize(n, n);
241 0 : return n;
242 0 : }
243 :
244 : /// calculates the normal of the crease vertex vrt on the side of vol
245 : /**
246 : * This algorithm uses grid::mark.
247 : * vol has to contain vrt. vrt has to have at least two associated
248 : * crease faces.
249 : *
250 : * Note that the resulting normal is not normalized. This is important,
251 : * since rounding errors could else lead to problems with normals which
252 : * have a length of nearly 0.
253 : *
254 : * This algorithm requires the option VOLOPT_AUTOGENERATE_FACES.
255 : * The option is automatically enabled if required.
256 : */
257 : template <class TAAPosVRT>
258 : typename TAAPosVRT::ValueType
259 0 : CalculateCreaseNormal(Grid& grid, Volume* vol, Vertex* vrt,
260 : Grid::face_traits::callback funcIsCreaseFace,
261 : TAAPosVRT& aaPos)
262 : {
263 0 : if(!grid.option_is_enabled(VOLOPT_AUTOGENERATE_FACES)){
264 : UG_LOG("WARNING in CalculateCreaseNormal: grid option VOLOPT_AUTOGENERATE_FACES autoenabled.\n");
265 0 : grid.enable_options(VOLOPT_AUTOGENERATE_FACES);
266 : }
267 :
268 : typedef typename TAAPosVRT::ValueType vector_t;
269 : vector_t n;
270 : VecSet(n, 0);
271 :
272 0 : grid.begin_marking();
273 :
274 : // we'll use a stack to find all volumes on this side of the crease.
275 : // each volume in the stack has to be marked.
276 : stack<Volume*> stk;
277 : stk.push(vol);
278 0 : grid.mark(vol);
279 :
280 : // objects for temporary results
281 0 : FaceDescriptor fd;
282 : vector<Volume*> vols;
283 :
284 : // we'll loop while there are faces in the stack
285 0 : while(!stk.empty()){
286 0 : Volume* curVol = stk.top();
287 : stk.pop();
288 :
289 : // iterate over all sides of curVol
290 0 : for(size_t i_side = 0; i_side < curVol->num_sides(); ++i_side){
291 0 : Face* f = grid.get_side(curVol, i_side);
292 0 : if(f){
293 : // only proceed if the face contains vrt
294 0 : if(FaceContains(f, vrt)){
295 : // check whether f is a crease
296 0 : if(funcIsCreaseFace(f)){
297 : // calculate the normal of the side
298 : vector_t tmpN;
299 0 : curVol->face_desc(i_side, fd);
300 0 : CalculateNormal(tmpN, &fd, aaPos);
301 :
302 : // the normal points away from the volume, so
303 : // we have to subtract it from n
304 : VecSubtract(n, n, tmpN);
305 : }
306 : else{
307 : // since the face is not a crease face, we have to add associated unmarked
308 : // volumes to the stack
309 0 : CollectVolumes(vols, grid, f);
310 0 : for(size_t i = 0; i < vols.size(); ++i){
311 0 : if(!grid.is_marked(vols[i])){
312 : grid.mark(vols[i]);
313 : stk.push(vols[i]);
314 : }
315 : }
316 : }
317 : }
318 : }
319 : }
320 : }
321 :
322 0 : grid.end_marking();
323 :
324 : //VecNormalize(n, n);
325 0 : return n;
326 0 : }
327 :
328 : /**
329 : * This algorithm indirectly uses Grid::mark.
330 : *
331 : * 1 dimensional fractures specified in fracInfos are expanded to 2 dimensional subsets.
332 : * the resulting fractures will then consist of 2 layers of quadrilaterals. On the
333 : * boundaries triangles are inserted.
334 : *
335 : * Through expandFracBoundaries you can tell the algorithm whether inner fracture
336 : * boundaries shall be expanded. Note that this means that an additional node is
337 : * introduced at each inner fracture boundary vertex and that the associated
338 : * fracture elements are connected at two sides.
339 : * Note that fractures are always expanded at boundaries which lie on the geometries
340 : * boundary.
341 : *
342 : * This algorithm requires the option FACEOPT_AUTOGENERATE_EDGES.
343 : * The option is automatically enabled if required.
344 : */
345 0 : bool ExpandFractures2d(Grid& grid, SubsetHandler& sh, const vector<FractureInfo>& fracInfos,
346 : bool expandInnerFracBnds, bool expandOuterFracBnds)
347 : {
348 : // access position attachment
349 0 : if(!grid.has_vertex_attachment(aPosition)){
350 : UG_LOG("Error in ExpandFractures: Missing position attachment");
351 0 : return false;
352 : }
353 : Grid::VertexAttachmentAccessor<APosition> aaPos(grid, aPosition);
354 :
355 0 : if(!grid.option_is_enabled(FACEOPT_AUTOGENERATE_EDGES)){
356 : UG_LOG("WARNING in CalculateCreaseNormal: grid option FACEOPT_AUTOGENERATE_EDGES autoenabled.\n");
357 0 : grid.enable_options(FACEOPT_AUTOGENERATE_EDGES);
358 : }
359 :
360 : // objects for temporary results
361 0 : FaceDescriptor fd;
362 : vector<Edge*> edges; // used for temporary results.
363 : vector<Face*> faces; // used for temporary results.
364 :
365 : // vectors that allow to access fracture properties by subset index
366 0 : vector<FractureInfo> fracInfosBySubset(sh.num_subsets(), FractureInfo(-1, -1, 0));
367 0 : for(size_t i = 0; i < fracInfos.size(); ++i){
368 0 : if(fracInfos[i].subsetIndex >= sh.num_subsets()){
369 0 : throw(UGError("Bad subsetIndex in given fracInfos."));
370 : }
371 :
372 0 : fracInfosBySubset[fracInfos[i].subsetIndex] = fracInfos[i];
373 : }
374 :
375 : ////////////////////////////////
376 : // Collect surrounding faces of all fractures in a selector
377 : // and select edges and vertices too.
378 0 : Selector sel(grid);
379 0 : sel.enable_autoselection(false);
380 0 : sel.enable_selection_inheritance(false);
381 :
382 : AInt aAdjMarker; // used to mark how many adjacent fractures a vertex has.
383 : // 0: no frac, 1: frac-boundary, >1: inner frac vertex
384 0 : grid.attach_to_vertices_dv(aAdjMarker, 0);
385 : Grid::VertexAttachmentAccessor<AInt> aaMarkVRT(grid, aAdjMarker);
386 0 : grid.attach_to_edges_dv(aAdjMarker, 0);
387 : Grid::EdgeAttachmentAccessor<AInt> aaMarkEDGE(grid, aAdjMarker);
388 :
389 : // iterate over the given fracture infos and select all fracture edges
390 : // and fracture vertices.
391 0 : for(size_t i_fi = 0; i_fi < fracInfos.size(); ++i_fi){
392 0 : int fracInd = fracInfos[i_fi].subsetIndex;
393 : for(EdgeIterator iter = sh.begin<Edge>(fracInd);
394 0 : iter != sh.end<Edge>(fracInd); ++iter)
395 : {
396 : // mark edge and vertices
397 : sel.select(*iter);
398 0 : aaMarkEDGE[*iter] = 1;
399 :
400 : // select associated vertices
401 0 : for(size_t i = 0; i < 2; ++i){
402 0 : Vertex* v = (*iter)->vertex(i);
403 : sel.select(v);
404 :
405 : // if fracture boundaries are expanded, we'll regard all fracture vertices
406 : // as inner vertices
407 0 : if(expandInnerFracBnds){
408 0 : if(!expandOuterFracBnds){
409 0 : if(IsBoundaryVertex2D(grid, v))
410 0 : aaMarkVRT[v]++;
411 : else
412 0 : aaMarkVRT[v] = 2;
413 : }
414 : else
415 0 : aaMarkVRT[v] = 2;
416 : }
417 : else
418 0 : aaMarkVRT[v]++;
419 : }
420 : }
421 : }
422 :
423 : // Make sure that selected vertices that lie on the boundary of the geometry
424 : // are treated as inner fracture vertices.
425 : // This is only required if frac-boundaries are not expanded anyways.
426 0 : if(expandOuterFracBnds && !expandInnerFracBnds){
427 : for(VertexIterator iter = sel.vertices_begin();
428 0 : iter != sel.vertices_end(); ++iter)
429 : {
430 : Vertex* v = *iter;
431 0 : if(aaMarkVRT[v] == 1){
432 0 : if(IsBoundaryVertex2D(grid, v))
433 0 : aaMarkVRT[v] = 2;
434 : }
435 : }
436 : }
437 :
438 : // Select all edges and faces which are connected to inner fracture vertices
439 : for(VertexIterator iter = sel.begin<Vertex>();
440 0 : iter != sel.end<Vertex>(); ++iter)
441 : {
442 0 : if(aaMarkVRT[*iter] > 1){
443 0 : sel.select(grid.associated_edges_begin(*iter),
444 : grid.associated_edges_end(*iter));
445 0 : sel.select(grid.associated_faces_begin(*iter),
446 : grid.associated_faces_end(*iter));
447 : }
448 : }
449 :
450 : ////////////////////////////////
451 : // create new vertices
452 :
453 : // we have to associate a vector of vertices with each node in the fracture.
454 : // since an empty vector is quite small, we can associate one with each vertex in
455 : // the whole grid. This could be optimized if required, by using subset attachments.
456 : typedef Attachment<vector<Vertex*> > AVrtVec;
457 : AVrtVec aVrtVec;
458 : grid.attach_to_vertices(aVrtVec);
459 : Grid::VertexAttachmentAccessor<AVrtVec> aaVrtVecVRT(grid, aVrtVec);
460 :
461 : // we also have to associate a vector of vertices for each face adjacent to the frac.
462 : // it will store the a second set of vertices. An entry contains the new vertex, if the
463 : // corresponding vertex is an inner fracture vertex, and NULL if not.
464 : grid.attach_to_faces(aVrtVec);
465 : Grid::FaceAttachmentAccessor<AVrtVec> aaVrtVecFACE(grid, aVrtVec);
466 :
467 : // a callback that returns true if the edge is a fracture edge
468 : AttachmentUnequal<Edge, Grid::EdgeAttachmentAccessor<AInt> > isFracEdge(aaMarkEDGE, 0);
469 :
470 : // iterate over all surrounding faces and create new vertices.
471 0 : for(FaceIterator iter_sf = sel.faces_begin(); iter_sf != sel.faces_end(); ++iter_sf)
472 : {
473 : Face* sf = *iter_sf;
474 :
475 : vector<Vertex*>& newVrts = aaVrtVecFACE[sf];
476 0 : newVrts.resize(sf->num_vertices());
477 :
478 : // check for each vertex whether it lies in the fracture
479 : // (aaMarkVRT > 1 in this case)
480 : // if so, we have to copy or create a vertex from/in aaVrtVec[vrt] which is
481 : // associated with the crease normal on the side of sf.
482 0 : for(size_t i_vrt = 0; i_vrt < sf->num_vertices(); ++i_vrt)
483 : {
484 0 : newVrts[i_vrt] = NULL;
485 0 : Vertex* vrt = sf->vertex(i_vrt);
486 0 : if(aaMarkVRT[vrt] > 1){
487 : // calculate the normal on this side of the frac
488 0 : vector3 n = CalculateCreaseNormal(grid, sf, vrt, isFracEdge, aaPos);
489 : //UG_LOG("calculated crease normal: " << n << endl);
490 :
491 : // check whether aaVrtVecs already contains a vertex associated with n.
492 : // the normal of new vrts is stored in their position attachment
493 : vector<Vertex*>& vrtVec = aaVrtVecVRT[vrt];
494 0 : for(size_t i = 0; i < vrtVec.size(); ++i){
495 : //UG_LOG("comparing to: " << aaPos[vrtVec[i]] << endl);
496 0 : if(VecDistanceSq(aaPos[vrtVec[i]], n) < SMALL){
497 : // got one
498 0 : newVrts[i_vrt] = vrtVec[i];
499 0 : break;
500 : }
501 : }
502 :
503 : // if we didn't find one then create and associate one.
504 : // store the normal in the position attachment of the new vertex
505 0 : if(!newVrts[i_vrt]){
506 0 : newVrts[i_vrt] = *grid.create<RegularVertex>();
507 : aaPos[newVrts[i_vrt]] = n;
508 0 : aaVrtVecVRT[vrt].push_back(newVrts[i_vrt]);
509 : }
510 : }
511 : }
512 : }
513 :
514 : // assign the new positions
515 : for(VertexIterator iter = sel.vertices_begin();
516 0 : iter != sel.vertices_end(); ++iter)
517 : {
518 : Vertex* vrt = *iter;
519 :
520 : // calculate the width as the maximum of associated fracture widths
521 0 : CollectEdges(edges, grid, vrt);
522 :
523 0 : number width = 0;
524 0 : for(size_t i = 0; i < edges.size(); ++i){
525 0 : if(aaMarkEDGE[edges[i]])
526 0 : width = max<number>(width, fracInfosBySubset.at(sh.get_subset_index(edges[i])).width);
527 : }
528 :
529 : // iterate over associated vertices
530 : vector<Vertex*>& vrtVec = aaVrtVecVRT[vrt];
531 :
532 : // note that the position attachment of new vertices holds their normal.
533 0 : for(size_t i = 0; i < vrtVec.size(); ++i){
534 0 : Vertex* nVrt = vrtVec[i];
535 0 : if(width > 0){
536 : vector3 n = aaPos[nVrt];
537 0 : if(VecLengthSq(n) > SMALL)
538 0 : VecNormalize(n, n);
539 :
540 0 : VecScale(n, n, width / 2.);
541 :
542 0 : UG_LOG("n: " << n << endl);
543 :
544 : // n now holds the offset for nVrt relative to vrt.
545 : // if width is higher than 0, we'll have to adjust the offset at
546 : // boundary vertices.
547 0 : if(IsBoundaryVertex2D(grid, vrt)){
548 : // First determine the normal pointing outwards
549 : vector3 nOut;
550 0 : CalculateBoundaryVertexNormal2D(nOut, grid, vrt, aaPos);
551 :
552 : // flip it by 90 degrees
553 0 : number tmp = nOut.x();
554 0 : nOut.x() = -nOut.y();
555 0 : nOut.y() = tmp;
556 :
557 0 : UG_LOG("nOut: " << nOut << endl);
558 :
559 : // now project the offset onto this vector
560 : VecScale(nOut, nOut, VecDot(nOut, n));
561 :
562 : // and now scale the new offset so that we receive the final offset.
563 : number dot = VecDot(n, nOut);
564 0 : if(dot > SMALL)
565 0 : VecScale(n, nOut, VecLengthSq(n) / dot);
566 : }
567 :
568 0 : UG_LOG("nFinal: " << n << endl);
569 : VecAdd(aaPos[nVrt], n, aaPos[vrt]);
570 : UG_LOG("\n");
571 : }
572 : else
573 : aaPos[nVrt] = aaPos[vrt];
574 : }
575 :
576 : // the current position is only a guess. Especially vertices where
577 : // fractures cross, this is not yet optimal.
578 : //todo: create an iterative spring system to find the new position.
579 : }
580 :
581 : ////////////////////////////////
582 : // create new elements
583 :
584 : // first we create new edges from selected ones which are connected to
585 : // inner vertices. This allows to preserve old subsets.
586 : // Since we have to make sure that we use the right vertices,
587 : // we have to iterate over the selected faces and perform all actions on the edges
588 : // of those faces.
589 0 : for(FaceIterator iter_sf = sel.faces_begin(); iter_sf != sel.faces_end(); ++iter_sf)
590 : {
591 : Face* sf = *iter_sf;
592 : // check for each edge whether it has to be copied.
593 0 : for(size_t i_edge = 0; i_edge < sf->num_edges(); ++i_edge){
594 0 : Edge* e = grid.get_edge(sf, i_edge);
595 0 : if(sel.is_selected(e)){
596 : // check the associated vertices through the volumes aaVrtVecVol attachment.
597 : // If at least one has an associated new vertex and if no edge between the
598 : // new vertices already exists, we'll create the new edge.
599 : size_t ind0 = i_edge;
600 0 : size_t ind1 = (i_edge + 1) % sf->num_edges();
601 :
602 0 : Vertex* nv0 = (aaVrtVecFACE[sf])[ind0];
603 0 : Vertex* nv1 = (aaVrtVecFACE[sf])[ind1];
604 :
605 0 : if(nv0 || nv1){
606 : // if one vertex has no associated new one, then we use the vertex itself
607 0 : if(!nv0)
608 0 : nv0 = sf->vertex(ind0);
609 0 : if(!nv1)
610 0 : nv1 = sf->vertex(ind1);
611 :
612 : // create the new edge if it not already exists.
613 0 : if(!grid.get_edge(nv0, nv1))
614 0 : grid.create_by_cloning(e, EdgeDescriptor(nv0, nv1), e);
615 : }
616 : }
617 : }
618 : }
619 :
620 : // iterate over all surrounding faces and create new vertices.
621 : // Since faces are replaced on the fly, we have to take care with the iterator.
622 0 : for(FaceIterator iter_sf = sel.faces_begin(); iter_sf != sel.faces_end();)
623 : {
624 : Face* sf = *iter_sf;
625 : ++iter_sf;
626 :
627 0 : vector<Vertex*> newVrts = aaVrtVecFACE[sf];
628 :
629 : // all new vertices have been assigned to newVrts.
630 : // Note that if newVrts[i] == NULL, then we have to take the
631 : // old vertex sf->vertex(i).
632 : // now expand the fracture edges of sf to faces.
633 0 : for(size_t i_vrt = 0; i_vrt < sf->num_vertices(); ++i_vrt){
634 : size_t iv1 = i_vrt;
635 0 : size_t iv2 = (i_vrt + 1) % sf->num_vertices();
636 0 : Edge* tEdge = grid.get_edge(sf->vertex(iv1), sf->vertex(iv2));
637 0 : if(tEdge){
638 0 : if(aaMarkEDGE[tEdge]){
639 : Face* expFace = NULL;
640 0 : if(newVrts[iv1] && newVrts[iv2]){
641 : // create a new quadrilateral
642 0 : expFace = *grid.create<Quadrilateral>(
643 0 : QuadrilateralDescriptor(sf->vertex(iv1), sf->vertex(iv2),
644 : newVrts[iv2], newVrts[iv1]));
645 : }
646 0 : else if(newVrts[iv1]){
647 : // create a new triangle
648 0 : expFace = *grid.create<Triangle>(
649 0 : TriangleDescriptor(sf->vertex(iv1), sf->vertex(iv2), newVrts[iv1]));
650 : }
651 0 : else if(newVrts[iv2]){
652 : // create a new triangle
653 0 : expFace = *grid.create<Triangle>(
654 0 : TriangleDescriptor(sf->vertex(iv1), sf->vertex(iv2), newVrts[iv2]));
655 : }
656 : else{
657 : // this code-block should never be entered. If it is entered then
658 : // we selected the wrong faces. This is would be a BUG!!!
659 : // remove the temporary attachments and throw an error
660 : grid.detach_from_vertices(aVrtVec);
661 : grid.detach_from_faces(aVrtVec);
662 : grid.detach_from_vertices(aAdjMarker);
663 : grid.detach_from_edges(aAdjMarker);
664 0 : throw(UGError("Implementation error in ExpandFractures2d."));
665 : }
666 :
667 0 : sh.assign_subset(expFace, fracInfosBySubset.at(sh.get_subset_index(tEdge)).newSubsetIndex);
668 : }
669 : }
670 : }
671 :
672 :
673 : // now set up a new face descriptor and replace the face.
674 0 : if(fd.num_vertices() != sf->num_vertices())
675 0 : fd.set_num_vertices(sf->num_vertices());
676 :
677 0 : for(size_t i_vrt = 0; i_vrt < sf->num_vertices(); ++i_vrt){
678 0 : if(newVrts[i_vrt])
679 0 : fd.set_vertex(i_vrt, newVrts[i_vrt]);
680 : else
681 0 : fd.set_vertex(i_vrt, sf->vertex(i_vrt));
682 : }
683 0 : grid.create_by_cloning(sf, fd, sf);
684 0 : grid.erase(sf);
685 0 : }
686 :
687 : // we have to clean up unused edges.
688 : // All selected edges with mark 0 have to be deleted.
689 0 : for(EdgeIterator iter = sel.edges_begin(); iter != sel.edges_end();)
690 : {
691 : // be careful with the iterator
692 : Edge* e = *iter;
693 : ++iter;
694 :
695 0 : if(!aaMarkEDGE[e])
696 0 : grid.erase(e);
697 : }
698 :
699 : // remove the temporary attachments
700 : grid.detach_from_vertices(aVrtVec);
701 : grid.detach_from_faces(aVrtVec);
702 : grid.detach_from_vertices(aAdjMarker);
703 : grid.detach_from_edges(aAdjMarker);
704 :
705 : return true;
706 0 : }
707 :
708 : /** Selects all involved geometic objects and assigns marks to them.
709 : * If required, som edges may be split, so that we always operate
710 : * on a fully expandable fracture.
711 : *
712 : * Make sure that selection_inheritance is enabled and that
713 : * strict_inheritance is disabled in sel.*/
714 0 : static void DistributeExpansionMarks3D(Grid& grid, SubsetHandler& sh, Selector& sel,
715 : const vector<FractureInfo>& fracInfos,
716 : bool expandInnerFracBnds,
717 : bool expandOuterFracBnds,
718 : Grid::VertexAttachmentAccessor<AInt>& aaMarkVRT,
719 : Grid::EdgeAttachmentAccessor<AInt>& aaMarkEDGE,
720 : Grid::FaceAttachmentAccessor<AInt>& aaMarkFACE)
721 : {
722 : Grid::VertexAttachmentAccessor<APosition> aaPos(grid, aPosition);
723 :
724 : // objects for temporary results
725 : //VolumeDescriptor vd;
726 : vector<Edge*> edges; // used for temporary results.
727 : vector<Face*> faces; // used for temporary results.
728 : vector<Volume*> vols; // used for temporary results.
729 :
730 : // iterate over the given fracture infos and select all fracture faces,
731 : // fracture edges and fracture vertices
732 0 : for(size_t i_fi = 0; i_fi < fracInfos.size(); ++i_fi){
733 0 : int fracInd = fracInfos[i_fi].subsetIndex;
734 : for(FaceIterator iter = sh.begin<Face>(fracInd);
735 0 : iter != sh.end<Face>(fracInd); ++iter)
736 : {
737 : Face* f = *iter;
738 :
739 : // mark face and vertices
740 0 : sel.select(f);
741 :
742 : // collect associated volumes of the vertices and add them to surroundingFaces
743 0 : for(size_t i = 0; i < f->num_vertices(); ++i){
744 0 : Vertex* v = f->vertex(i);
745 : sel.select(v);
746 : // CollectVolumes(vols, grid, v);
747 : // sel.select(vols.begin(), vols.end());
748 : }
749 :
750 : // collect associated edges of the faces and
751 : // increase their adjacency counter too, since this helps
752 : // to identify whether they lie on the selection boundary.
753 0 : CollectEdges(edges, grid, f);
754 :
755 0 : for(size_t i = 0; i < edges.size(); ++i){
756 0 : aaMarkEDGE[edges[i]]++;
757 : sel.select(edges[i]);
758 : }
759 : }
760 : }
761 :
762 : // all edges that lie on the geometries boundary have to be regarded as inner edges
763 0 : if(expandOuterFracBnds){
764 0 : for(EdgeIterator iter = sel.edges_begin(); iter != sel.edges_end(); ++iter)
765 : {
766 : Edge* e = *iter;
767 0 : if(aaMarkEDGE[e] == 1){
768 0 : if(IsBoundaryEdge3D(grid, e))
769 0 : aaMarkEDGE[e] = 2;
770 : }
771 : }
772 : }
773 :
774 : // iterate over selected vertices and check the adjacency status of associated edges.
775 : // the vertex can either be a boundary vertex (1) or a surface vertex(2) or both (3).
776 : // Note that vertices lying on the geometries boundary are regarded as surface vertices.
777 : // At this point it is important to mark all vertices that lie on any boundary as
778 : // a boundary vertex, since we have to split inner edges connecting boundary-vertices
779 : // in the next step.
780 : // Note that this is important even for the degenerate case, since otherwise identical
781 : // elements may be created.
782 : //todo: Currently only surface or boundary vertices are regarded. The mix is
783 : // treated as a boundary vertex for the non-degenerated case.
784 : // This is because regarding the mix of both would result in problematic element shapes.
785 : // For the degenerated case we regard the mix as an inner vertex.
786 : for(VertexIterator iter = sel.vertices_begin();
787 0 : iter != sel.vertices_end(); ++iter)
788 : {
789 0 : CollectEdges(edges, grid, *iter);
790 0 : aaMarkVRT[*iter] = 2;
791 0 : for(size_t i = 0; i < edges.size(); ++i){
792 0 : if(aaMarkEDGE[edges[i]] == 1){
793 : // if(VertexLiesOnSurface(grid, *iter, IsSelected(sel))){
794 : // UG_LOG("found mixed\n");
795 : // aaMarkVRT[*iter] = 3;
796 : // }
797 0 : aaMarkVRT[*iter] = 1;
798 0 : break;
799 : }
800 : }
801 : }
802 :
803 : // todo: Quadrilaterals which have more than 2 boundary vertices or have two boundary
804 : // vertices, which are not adjacent to each other, have to be transformed to
805 : // triangles.
806 :
807 :
808 : // now make sure that no inner edge is associated with two
809 : // boundary vertices (referring to the selection)
810 : edges.clear();
811 : for(EdgeIterator iter = sel.begin<Edge>();
812 0 : iter != sel.end<Edge>(); ++iter)
813 : {
814 0 : Edge* e = *iter;
815 0 : if(aaMarkVRT[e->vertex(0)] != 2 &&
816 0 : aaMarkVRT[e->vertex(1)] != 2 &&
817 0 : aaMarkEDGE[e] > 1)
818 : {
819 0 : edges.push_back(e);
820 : }
821 : }
822 :
823 0 : for(size_t i = 0; i < edges.size(); ++i){
824 0 : vector3 center = CalculateCenter(edges[i], aaPos);
825 0 : RegularVertex* v = SplitEdge<RegularVertex>(grid, edges[i], false);
826 : aaPos[v] = center;
827 0 : aaMarkVRT[v] = 2;
828 0 : sel.select(v);
829 : // assign adjacency values for associated selected edges (2 to each)
830 0 : for(Grid::AssociatedEdgeIterator iter = grid.associated_edges_begin(v);
831 0 : iter != grid.associated_edges_end(v); ++iter)
832 : {
833 0 : if(sel.is_selected(*iter))
834 0 : aaMarkEDGE[*iter] = 2;
835 : }
836 : }
837 :
838 : // if fracture boundaries shall be extended, then we have to regard all
839 : // boundary vertices as inner vertices
840 0 : if(expandInnerFracBnds && expandOuterFracBnds){
841 0 : for(VertexIterator iter = sel.vertices_begin(); iter != sel.vertices_end(); ++iter)
842 : {
843 : Vertex* v = *iter;
844 0 : if(aaMarkVRT[v] == 1)
845 0 : aaMarkVRT[v] = 2;
846 : }
847 : }
848 0 : else if(expandInnerFracBnds){
849 0 : for(VertexIterator iter = sel.vertices_begin(); iter != sel.vertices_end(); ++iter)
850 : {
851 : Vertex* v = *iter;
852 0 : if(!IsBoundaryVertex3D(grid, v)){
853 0 : if(aaMarkVRT[v] == 1)
854 0 : aaMarkVRT[v] = 2;
855 : }
856 : }
857 : }
858 0 : else if(expandOuterFracBnds){
859 0 : for(VertexIterator iter = sel.vertices_begin(); iter != sel.vertices_end(); ++iter)
860 : {
861 : Vertex* v = *iter;
862 0 : if(IsBoundaryVertex3D(grid, v)){
863 : // get state from marked associated boundary edges
864 0 : aaMarkVRT[v] = 0;
865 : CollectAssociated(edges, grid, v);
866 :
867 0 : for(size_t i = 0; i < edges.size(); ++i){
868 0 : Edge* e = edges[i];
869 0 : if(aaMarkEDGE[e] > 1){
870 0 : if(IsBoundaryEdge3D(grid, e))
871 0 : aaMarkVRT[v]++;
872 : }
873 : }
874 : }
875 : }
876 : }
877 :
878 :
879 : // All fracture faces shall be marked with 1. We do this here, since new
880 : // faces may have been created during the edge-splits.
881 0 : for(FaceIterator iter = sel.begin<Face>(); iter != sel.end<Face>(); ++iter)
882 0 : aaMarkFACE[*iter] = 1;
883 :
884 : // select all non-fracture edges, faces and volumes, which are connected to an
885 : // inner fracture vertex.
886 : for(VertexIterator iter = sel.begin<Vertex>();
887 0 : iter != sel.end<Vertex>(); ++iter)
888 : {
889 : Vertex* vrt = *iter;
890 0 : if(aaMarkVRT[vrt] > 1){
891 : // select all associated edges, faces and volumes
892 0 : sel.select(grid.associated_edges_begin(vrt),
893 : grid.associated_edges_end(vrt));
894 0 : sel.select(grid.associated_faces_begin(vrt),
895 : grid.associated_faces_end(vrt));
896 0 : sel.select(grid.associated_volumes_begin(vrt),
897 : grid.associated_volumes_end(vrt));
898 : }
899 : }
900 0 : }
901 :
902 :
903 : /**
904 : * This algorithm indirectly uses Grid::mark.
905 : *
906 : * 2 dimensional fractures specified in fracInfos are expanded to 3 dimensional subsets.
907 : * the resulting fractures will then consist of 2 layers of hexahedrons. On the
908 : * boundaries tetrahedrons, prisms and pyramids are inserted.
909 : *
910 : * Through expandFracBoundaries you can tell the algorithm whether inner fracture
911 : * boundaries shall be expanded. Note that this means that an additional node is
912 : * introduced at each inner fracture boundary vertex and that the associated
913 : * fracture elements are connected at two sides.
914 : * Note that fractures are always expanded at boundaries which lie on the geometries
915 : * boundary.
916 : *
917 : * This algorithm requires the option FACEOPT_AUTOGENERATE_EDGES.
918 : * The option is automatically enabled if required.
919 : *
920 : * This algorithm requires the option VOLOPT_AUTOGENERATE_FACES.
921 : * The option is automatically enabled if required.
922 : */
923 0 : bool ExpandFractures3d(Grid& grid, SubsetHandler& sh, const vector<FractureInfo>& fracInfos,
924 : bool expandInnerFracBnds, bool expandOuterFracBnds)
925 : {
926 : // access position attachment
927 0 : if(!grid.has_vertex_attachment(aPosition)){
928 : UG_LOG("Error in ExpandFractures: Missing position attachment");
929 0 : return false;
930 : }
931 : Grid::VertexAttachmentAccessor<APosition> aaPos(grid, aPosition);
932 :
933 : // make sure that the required options are enabled.
934 0 : if(!grid.option_is_enabled(VOLOPT_AUTOGENERATE_FACES)){
935 : UG_LOG("WARNING in CalculateCreaseNormal: grid option VOLOPT_AUTOGENERATE_FACES autoenabled.\n");
936 0 : grid.enable_options(VOLOPT_AUTOGENERATE_FACES);
937 : }
938 :
939 0 : if(!grid.option_is_enabled(FACEOPT_AUTOGENERATE_EDGES)){
940 : UG_LOG("WARNING in CalculateCreaseNormal: grid option FACEOPT_AUTOGENERATE_EDGES autoenabled.\n");
941 0 : grid.enable_options(FACEOPT_AUTOGENERATE_EDGES);
942 : }
943 :
944 : // objects for temporary results
945 0 : FaceDescriptor fd;
946 0 : VolumeDescriptor vd;
947 : vector<Edge*> edges; // used for temporary results.
948 : vector<Face*> faces; // used for temporary results.
949 : vector<Volume*> vols; // used for temporary results.
950 :
951 : ////////////////////////////////
952 : // Collect surrounding volumes, faces and edges of all fractures in a selector
953 : // and select fracture faces, edges and vertices too.
954 0 : Selector sel(grid);
955 : //Selector& sel = tmpSel;
956 0 : sel.enable_autoselection(false);
957 0 : sel.enable_selection_inheritance(true); //required for DistributeExpansionMarks3D. disabled later on.
958 0 : sel.enable_strict_inheritance(false);
959 :
960 : AInt aAdjMarker; // used to mark where in a fracture an edge or a vertex lies.
961 : // 0: no frac, 1: frac-boundary, >1: inner frac
962 0 : grid.attach_to_vertices_dv(aAdjMarker, 0);
963 : Grid::VertexAttachmentAccessor<AInt> aaMarkVRT(grid, aAdjMarker);
964 0 : grid.attach_to_edges_dv(aAdjMarker, 0);
965 : Grid::EdgeAttachmentAccessor<AInt> aaMarkEDGE(grid, aAdjMarker);
966 0 : grid.attach_to_faces_dv(aAdjMarker, 0);
967 : Grid::FaceAttachmentAccessor<AInt> aaMarkFACE(grid, aAdjMarker);
968 :
969 : // Distribute marks and select elements.
970 :
971 0 : DistributeExpansionMarks3D(grid, sh, sel, fracInfos, expandInnerFracBnds,
972 : expandOuterFracBnds, aaMarkVRT, aaMarkEDGE, aaMarkFACE);
973 :
974 : // We'll now store all fracture faces in a vector.
975 : // They will help to adjust positions of new vertices later on.
976 : std::vector<Face*> originalFractureFaces;
977 0 : for(FaceIterator iter = sel.begin<Face>(); iter != sel.end<Face>(); ++iter){
978 0 : if(aaMarkFACE[*iter] == 1)
979 0 : originalFractureFaces.push_back(*iter);
980 : }
981 :
982 : // vectors that allow to access fracture properties by subset index
983 0 : vector<FractureInfo> fracInfosBySubset(sh.num_subsets(), FractureInfo(-1, -1, 0));
984 0 : for(size_t i = 0; i < fracInfos.size(); ++i){
985 0 : if(fracInfos[i].subsetIndex >= sh.num_subsets()){
986 0 : throw(UGError("Bad subsetIndex in given fracInfos."));
987 : }
988 :
989 0 : fracInfosBySubset[fracInfos[i].subsetIndex] = fracInfos[i];
990 : }
991 :
992 : // disable selection inheritance to avoid infinite recursion.
993 0 : sel.enable_selection_inheritance(false);
994 : // clear buffers for later use
995 : edges.clear();
996 :
997 : ////////////////////////////////
998 : // create new vertices
999 : // we have to associate a vector of vertices with each node in the fracture.
1000 : // since an empty vector is quite small, we can associate one with each vertex in
1001 : // the whole grid. This could be optimized if required, by using subset attachments.
1002 : typedef Attachment<vector<Vertex*> > AVrtVec;
1003 : AVrtVec aVrtVec;
1004 : grid.attach_to_vertices(aVrtVec);
1005 : Grid::VertexAttachmentAccessor<AVrtVec> aaVrtVecVRT(grid, aVrtVec);
1006 :
1007 : // we also have to associate a vector of vertices for each volume adjacent to the frac.
1008 : // it will store a second set of vertices. An entry contains the new vertex, if the
1009 : // corresponding vertex is an inner fracture vertex, and NULL if not.
1010 : grid.attach_to_volumes(aVrtVec);
1011 : Grid::VolumeAttachmentAccessor<AVrtVec> aaVrtVecVOL(grid, aVrtVec);
1012 :
1013 : // a callback which tells whether a face is inside the fracture or not
1014 : AttachmentUnequal<Face, Grid::FaceAttachmentAccessor<AInt> > faceIsInFrac(aaMarkFACE, 0);
1015 :
1016 : // iterate over all surrounding volumes and create new vertices.
1017 0 : for(VolumeIterator iter_sv = sel.volumes_begin(); iter_sv != sel.volumes_end(); ++iter_sv)
1018 : {
1019 : Volume* sv = *iter_sv;
1020 :
1021 : vector<Vertex*>& newVrts = aaVrtVecVOL[sv];
1022 0 : newVrts.resize(sv->num_vertices());
1023 :
1024 : // check for each vertex whether it lies in the fracture
1025 : // (aaMarkVRT > 1 in this case)
1026 : // if so, we have to copy or create a vertex from/in aaVrtVec[vrt] which is
1027 : // associated with the crease normal on the side of sf.
1028 0 : for(size_t i_vrt = 0; i_vrt < sv->num_vertices(); ++i_vrt)
1029 : {
1030 0 : newVrts[i_vrt] = NULL;
1031 0 : Vertex* vrt = sv->vertex(i_vrt);
1032 0 : if(aaMarkVRT[vrt] > 1){
1033 : // calculate the normal on this side of the frac
1034 0 : vector3 n = CalculateCreaseNormal(grid, sv, vrt, faceIsInFrac, aaPos);
1035 : //UG_LOG("calculated crease normal: " << n << endl);
1036 :
1037 : // check whether aaVrtVecs already contains a vertex associated with n.
1038 : // the normal of new vrts is stored in their position attachment
1039 : vector<Vertex*>& vrtVec = aaVrtVecVRT[vrt];
1040 0 : for(size_t i = 0; i < vrtVec.size(); ++i){
1041 : //UG_LOG("comparing to: " << aaPos[vrtVec[i]] << endl);
1042 0 : if(VecDistanceSq(aaPos[vrtVec[i]], n) < SMALL){
1043 : // got one
1044 0 : newVrts[i_vrt] = vrtVec[i];
1045 0 : break;
1046 : }
1047 : }
1048 :
1049 : // if we didn't find one then create and associate one.
1050 : // store the normal in the position attachment of the new vertex
1051 0 : if(!newVrts[i_vrt]){
1052 0 : newVrts[i_vrt] = *grid.create<RegularVertex>();
1053 : aaPos[newVrts[i_vrt]] = n;
1054 0 : aaVrtVecVRT[vrt].push_back(newVrts[i_vrt]);
1055 : }
1056 : }
1057 : }
1058 : }
1059 :
1060 : ////////////////////////////////
1061 : // assign positions
1062 : for(VertexIterator iter = sel.vertices_begin();
1063 0 : iter != sel.vertices_end(); ++iter)
1064 : {
1065 : Vertex* vrt = *iter;
1066 :
1067 : // calculate the width as the maximum of associated fracture widths
1068 0 : CollectFaces(faces, grid, vrt);
1069 :
1070 0 : number width = 0;
1071 0 : for(size_t i = 0; i < faces.size(); ++i){
1072 0 : if(aaMarkFACE[faces[i]])
1073 0 : width = max<number>(width, fracInfosBySubset.at(sh.get_subset_index(faces[i])).width);
1074 : }
1075 :
1076 : // iterate over associated vertices
1077 : vector<Vertex*>& vrtVec = aaVrtVecVRT[vrt];
1078 :
1079 : // note that the position attachment of new vertices holds their normal.
1080 0 : for(size_t i = 0; i < vrtVec.size(); ++i){
1081 0 : Vertex* nVrt = vrtVec[i];
1082 0 : if(width > 0){
1083 : vector3 n = aaPos[nVrt];
1084 0 : if(VecLengthSq(n) > SMALL)
1085 0 : VecNormalize(n, n);
1086 :
1087 0 : VecScale(n, n, width / 2.);
1088 :
1089 0 : if(IsBoundaryVertex3D(grid, vrt)){
1090 : // First determine the normal pointing outwards
1091 : vector3 nOut;
1092 0 : CalculateBoundaryVertexNormal3D(nOut, grid, vrt, aaPos);
1093 :
1094 : // project the normal into the plane with the normal nOut
1095 : vector3 nNew;
1096 0 : ProjectPointToPlane(nNew, n, vector3(0, 0, 0), nOut);
1097 :
1098 : // and now scale the new offset so that we receive the final offset.
1099 : number dot = VecDot(n, nNew);
1100 0 : if(dot > SMALL)
1101 0 : VecScale(n, nNew, VecLengthSq(n) / dot);
1102 :
1103 : }
1104 :
1105 : VecAdd(aaPos[nVrt], n, aaPos[vrt]);
1106 : }
1107 : else
1108 : aaPos[nVrt] = aaPos[vrt];
1109 : }
1110 :
1111 : // the current position is only a guess. Especially at vertices where
1112 : // fractures cross, this is not yet optimal.
1113 : //todo: create an iterative spring system to find the new position.
1114 : }
1115 :
1116 : ////////////////////////////////
1117 : // create new elements
1118 :
1119 : // holds local side vertex indices
1120 : vector<size_t> locVrtInds;
1121 : // all new vertices have been assigned to newVrts.
1122 : // Note that if newVrts[i] == NULL, then we have to take the
1123 : // old vertex sf->vertex(i).
1124 :
1125 : // first we create new edges from selected ones which are connected to
1126 : // inner vertices. This allows to preserve old subsets.
1127 : // Since we have to make sure that we use the right vertices,
1128 : // we have to iterate over the selected volumes and perform all actions on the edges
1129 : // of those volumes.
1130 0 : for(VolumeIterator iter_sv = sel.volumes_begin(); iter_sv != sel.volumes_end(); ++iter_sv)
1131 : {
1132 : Volume* sv = *iter_sv;
1133 : // check for each edge whether it has to be copied.
1134 0 : for(size_t i_edge = 0; i_edge < sv->num_edges(); ++i_edge){
1135 0 : Edge* e = grid.get_edge(sv, i_edge);
1136 0 : if(sel.is_selected(e)){
1137 : // check the associated vertices through the volumes aaVrtVecVol attachment.
1138 : // If at least one has an associated new vertex and if no edge between the
1139 : // new vertices already exists, we'll create the new edge.
1140 : size_t ind0, ind1;
1141 0 : sv->get_vertex_indices_of_edge(ind0, ind1, i_edge);
1142 0 : Vertex* nv0 = (aaVrtVecVOL[sv])[ind0];
1143 0 : Vertex* nv1 = (aaVrtVecVOL[sv])[ind1];
1144 :
1145 0 : if(nv0 || nv1){
1146 : // if one vertex has no associated new one, then we use the vertex itself
1147 0 : if(!nv0)
1148 0 : nv0 = sv->vertex(ind0);
1149 0 : if(!nv1)
1150 0 : nv1 = sv->vertex(ind1);
1151 :
1152 : // create the new edge if it not already exists.
1153 0 : if(!grid.get_edge(nv0, nv1))
1154 0 : grid.create_by_cloning(e, EdgeDescriptor(nv0, nv1), e);
1155 : }
1156 : }
1157 : }
1158 : }
1159 :
1160 : // now we create new faces from selected ones which are connected to
1161 : // inner vertices. This allows to preserve old subsets.
1162 : // Since we have to make sure that we use the right vertices,
1163 : // we have to iterate over the selected volumes and perform all actions on the side-faces
1164 : // of those volumes.
1165 0 : for(VolumeIterator iter_sv = sel.volumes_begin(); iter_sv != sel.volumes_end(); ++iter_sv)
1166 : {
1167 : Volume* sv = *iter_sv;
1168 : // check for each face whether it has to be copied.
1169 0 : for(size_t i_face = 0; i_face < sv->num_faces(); ++i_face){
1170 0 : Face* sf = grid.get_face(sv, i_face);
1171 0 : if(sel.is_selected(sf)){
1172 : // check the associated vertices through the volumes aaVrtVecVol attachment.
1173 : // If no face between the new vertices already exists, we'll create the new face.
1174 0 : sv->get_vertex_indices_of_face(locVrtInds, i_face);
1175 0 : fd.set_num_vertices(sf->num_vertices());
1176 :
1177 0 : for(size_t i = 0; i < sf->num_vertices(); ++i){
1178 0 : Vertex* nVrt = (aaVrtVecVOL[sv])[locVrtInds[i]];
1179 0 : if(nVrt)
1180 0 : fd.set_vertex(i, nVrt);
1181 : else
1182 0 : fd.set_vertex(i, sv->vertex(locVrtInds[i]));
1183 : }
1184 :
1185 : // if the new face does not already exist, we'll create it
1186 0 : if(!grid.get_face(fd))
1187 0 : grid.create_by_cloning(sf, fd, sf);
1188 : }
1189 : }
1190 : }
1191 :
1192 : // Expand all faces.
1193 : // Since volumes are replaced on the fly, we have to take care with the iterator.
1194 : // record all new volumes in a vector. This will help to adjust positions later on.
1195 : vector<Volume*> newFractureVolumes;
1196 0 : for(VolumeIterator iter_sv = sel.volumes_begin(); iter_sv != sel.volumes_end();)
1197 : {
1198 : Volume* sv = *iter_sv;
1199 : ++iter_sv;
1200 :
1201 : // now expand the fracture faces of sv to volumes.
1202 0 : for(size_t i_side = 0; i_side < sv->num_sides(); ++i_side){
1203 : // get the local vertex indices of the side of the volume
1204 0 : sv->get_vertex_indices_of_face(locVrtInds, i_side);
1205 :
1206 0 : Face* tFace = grid.get_side(sv, i_side);
1207 :
1208 0 : if(tFace){
1209 0 : if(aaMarkFACE[tFace]){
1210 0 : Volume* expVol = NULL;
1211 0 : if(locVrtInds.size() == 3){
1212 0 : size_t iv0 = locVrtInds[0];
1213 0 : size_t iv1 = locVrtInds[1];
1214 0 : size_t iv2 = locVrtInds[2];
1215 :
1216 0 : if((aaVrtVecVOL[sv])[iv0] && (aaVrtVecVOL[sv])[iv1] && (aaVrtVecVOL[sv])[iv2]){
1217 : // create a new prism
1218 0 : expVol = *grid.create<Prism>(
1219 0 : PrismDescriptor(sv->vertex(iv2), sv->vertex(iv1), sv->vertex(iv0),
1220 : (aaVrtVecVOL[sv])[iv2], (aaVrtVecVOL[sv])[iv1], (aaVrtVecVOL[sv])[iv0]));
1221 : }
1222 0 : else if((aaVrtVecVOL[sv])[iv0] && (aaVrtVecVOL[sv])[iv1]){
1223 : // create a new Pyramid
1224 0 : expVol = *grid.create<Pyramid>(
1225 0 : PyramidDescriptor(sv->vertex(iv0), sv->vertex(iv1),
1226 0 : (aaVrtVecVOL[sv])[iv1], (aaVrtVecVOL[sv])[iv0], sv->vertex(iv2)));
1227 : }
1228 0 : else if((aaVrtVecVOL[sv])[iv1] && (aaVrtVecVOL[sv])[iv2]){
1229 : // create a new Pyramid
1230 0 : expVol = *grid.create<Pyramid>(
1231 0 : PyramidDescriptor(sv->vertex(iv1), sv->vertex(iv2),
1232 0 : (aaVrtVecVOL[sv])[iv2], (aaVrtVecVOL[sv])[iv1], sv->vertex(iv0)));
1233 : }
1234 0 : else if((aaVrtVecVOL[sv])[iv0] && (aaVrtVecVOL[sv])[iv2]){
1235 : // create a new Pyramid
1236 0 : expVol = *grid.create<Pyramid>(
1237 0 : PyramidDescriptor(sv->vertex(iv2), sv->vertex(iv0),
1238 0 : (aaVrtVecVOL[sv])[iv0], (aaVrtVecVOL[sv])[iv2], sv->vertex(iv1)));
1239 : }
1240 0 : else if((aaVrtVecVOL[sv])[iv0]){
1241 : // create a new Tetrahedron
1242 0 : expVol = *grid.create<Tetrahedron>(
1243 0 : TetrahedronDescriptor(sv->vertex(iv2), sv->vertex(iv1), sv->vertex(iv0),
1244 : (aaVrtVecVOL[sv])[iv0]));
1245 : }
1246 0 : else if((aaVrtVecVOL[sv])[iv1]){
1247 : // create a new Tetrahedron
1248 0 : expVol = *grid.create<Tetrahedron>(
1249 0 : TetrahedronDescriptor(sv->vertex(iv2), sv->vertex(iv1), sv->vertex(iv0),
1250 : (aaVrtVecVOL[sv])[iv1]));
1251 : }
1252 0 : else if((aaVrtVecVOL[sv])[iv2]){
1253 : // create a new Tetrahedron
1254 0 : expVol = *grid.create<Tetrahedron>(
1255 0 : TetrahedronDescriptor(sv->vertex(iv2), sv->vertex(iv1), sv->vertex(iv0),
1256 : (aaVrtVecVOL[sv])[iv2]));
1257 : }
1258 : else{
1259 : // this code-block should never be entered. If it is entered then
1260 : // we either selected the wrong faces (this shouldn't happen), or there
1261 : // are selected faces, which have fracture-boundary-vertices only.
1262 : // This is the same is if inner fracture edges exists, which are
1263 : // connected to two boundary vertices.
1264 : // Since we tried to remove those edges above, something went wrong.
1265 : // remove the temporary attachments and throw an error
1266 : grid.detach_from_vertices(aVrtVec);
1267 : grid.detach_from_volumes(aVrtVec);
1268 : grid.detach_from_vertices(aAdjMarker);
1269 : grid.detach_from_edges(aAdjMarker);
1270 0 : throw(UGError("Error in ExpandFractures3d. Implementation Error."));
1271 : }
1272 : }
1273 : else{
1274 : // currently only tetrahedrons are supported. This section thus raises an error
1275 : grid.detach_from_vertices(aVrtVec);
1276 : grid.detach_from_volumes(aVrtVec);
1277 : grid.detach_from_vertices(aAdjMarker);
1278 : grid.detach_from_edges(aAdjMarker);
1279 0 : throw(UGError("Incomplete implementation error in ExpandFractures3d: Only tetrahedrons are supported in the current implementation."));
1280 : }
1281 0 : if(expVol){
1282 0 : sh.assign_subset(expVol, fracInfosBySubset.at(sh.get_subset_index(tFace)).newSubsetIndex);
1283 0 : newFractureVolumes.push_back(expVol);
1284 : }
1285 : }
1286 : }
1287 : }
1288 :
1289 : // now set up a new volume descriptor and replace the volume.
1290 0 : if(vd.num_vertices() != sv->num_vertices())
1291 0 : vd.set_num_vertices(sv->num_vertices());
1292 :
1293 0 : for(size_t i_vrt = 0; i_vrt < sv->num_vertices(); ++i_vrt){
1294 0 : if((aaVrtVecVOL[sv])[i_vrt])
1295 0 : vd.set_vertex(i_vrt, (aaVrtVecVOL[sv])[i_vrt]);
1296 : else
1297 0 : vd.set_vertex(i_vrt, sv->vertex(i_vrt));
1298 : }
1299 :
1300 0 : grid.create_by_cloning(sv, vd, sv);
1301 0 : grid.erase(sv);
1302 : }
1303 :
1304 : // we have to clean up unused faces and edges.
1305 : // note that all selected edges with mark 0 may safley be deleted.
1306 : for(EdgeIterator iter = sel.begin<Edge>();
1307 0 : iter!= sel.end<Edge>();)
1308 : {
1309 : // take care of the iterator
1310 : Edge* e = *iter;
1311 : ++iter;
1312 :
1313 0 : if(aaMarkEDGE[e] == 0)
1314 0 : grid.erase(e);
1315 : }
1316 :
1317 : // make sure that no unused faces linger around (This should never happen!)
1318 : bool foundUnusedFaces = false;
1319 : for(FaceIterator iter = sel.begin<Face>();
1320 0 : iter != sel.end<Face>();)
1321 : {
1322 : Face* f = *iter;
1323 : ++iter;
1324 :
1325 0 : if(aaMarkFACE[f] == 0){
1326 : foundUnusedFaces = true;
1327 0 : grid.erase(f);
1328 : }
1329 : }
1330 :
1331 0 : if(foundUnusedFaces){
1332 : UG_LOG("WARNING in ExpandFractures3D: Unused faces encountered during cleanup. Removing them...\n");
1333 : }
1334 : /*
1335 : // finally assign the new positions
1336 : // find the maximal fracture width. If it is 0, we only have to copy positions.
1337 : number maxFracWidth = 0;
1338 : for(size_t i = 0; i < fracInfos.size(); ++i)
1339 : maxFracWidth = max(fracInfos[i].width, maxFracWidth);
1340 :
1341 : // in this case equality with 0 is desired.
1342 : if(maxFracWidth == 0){
1343 : // set all positions of new vertices to the positions of their parents
1344 : for(VertexIterator iter = sel.vertices_begin();
1345 : iter != sel.vertices_end(); ++iter)
1346 : {
1347 : Vertex* vrt = *iter;
1348 : // iterate over associated vertices
1349 : vector<Vertex*>& vrtVec = aaVrtVecVRT[vrt];
1350 : for(size_t i = 0; i < vrtVec.size(); ++i){
1351 : Vertex* nVrt = vrtVec[i];
1352 : aaPos[nVrt] = aaPos[vrt];
1353 : }
1354 : }
1355 : }
1356 : else{
1357 : // we will find the new positions by solving a minimal least squares problem.
1358 : //todo: Special treatment has to be to boundary vertices.
1359 :
1360 : // currently all original fracture vertices are selected.
1361 : }*/
1362 :
1363 : // remove the temporary attachments
1364 : grid.detach_from_vertices(aVrtVec);
1365 : grid.detach_from_faces(aVrtVec);
1366 : grid.detach_from_vertices(aAdjMarker);
1367 : grid.detach_from_edges(aAdjMarker);
1368 : return true;
1369 0 : }
1370 :
1371 : }// end of namespace
1372 :
1373 :
1374 : // This method is unused. And I can't really remeber what its use was.
1375 : // It looks quite complicated for a rather simple task. Probably some
1376 : // deeper thought was involved...
1377 : /// this method returns true if the face has to be treated as
1378 : /// an inner fracture boundary face.
1379 : /** This method uses Grid::mark.
1380 : *
1381 : * This method starts at all open boundary edges and tries
1382 : * to find f by checking faces which are connected to one
1383 : * of the boundary-edges endpoints and are reachable from
1384 : * the initial selection by traversing faces over regular surface
1385 : * edges.
1386 : * If it can be found, the face is a fracture boundary face,
1387 : * if not, it is a inner fracture face.
1388 : *
1389 : * This method only makes sense, if funcIsSurfFace(f) evaluates to true.
1390 : */
1391 : /*
1392 : bool FractureBoundaryFace(Grid& grid, Face* f,
1393 : CB_ConsiderFace funcIsSurfFace)
1394 : {
1395 : if(!funcIsSurfFace(f))
1396 : return false;
1397 :
1398 : if(!grid.option_is_enabled(FACEOPT_AUTOGENERATE_EDGES)){
1399 : UG_LOG("WARNING in FractureBoundaryFace: autoenabling FACEOPT_AUTOGENERATE_EDGES.\n");
1400 : grid.enable_options(FACEOPT_AUTOGENERATE_EDGES);
1401 : }
1402 :
1403 : grid.begin_marking();
1404 :
1405 : stack<Edge*> stk;
1406 : vector<Edge*> edges;
1407 : vector<Face*> allFaces;
1408 : vector<Face*> faces;
1409 :
1410 : // collect all associated fracture-boundary-edges of vertices of e
1411 : for(size_t i = 0; i < f->num_vertices(); ++i){
1412 : for(Grid::AssociatedEdgeIterator iter = grid.associated_edges_begin(f->vertex(i));
1413 : iter != grid.associated_edges_end(f->vertex(i)); ++iter)
1414 : {
1415 : Edge* e = *iter;
1416 : if(!grid.is_marked(e)){
1417 : if(IsBoundaryEdge(grid, e, funcIsSurfFace)){
1418 : grid.mark(e);
1419 : grid.mark(e->vertex(0));
1420 : grid.mark(e->vertex(1));
1421 : stk.push(e);
1422 : }
1423 : }
1424 : }
1425 : }
1426 :
1427 : while(!stk.empty()){
1428 : Edge* e = stk.top();
1429 : // we can safely pop e from the stack
1430 : stk.pop();
1431 :
1432 : // collect unmarked surface-faces associated with e
1433 : CollectFaces(allFaces, grid, e);
1434 : faces.clear();
1435 : for(size_t i = 0; i < allFaces.size(); ++i){
1436 : if(!grid.is_marked(allFaces[i])){
1437 : if(funcIsSurfFace(allFaces[i]))
1438 : faces.push_back(allFaces[i]);
1439 : }
1440 : }
1441 :
1442 : // make sure that we found exactly one face. Otherwise
1443 : // we would not cross a regular surface edge
1444 : if(faces.size() != 1)
1445 : continue;
1446 :
1447 : Face* nf = faces[0];
1448 :
1449 : // if we found f, it is clear that f is a boundary face.
1450 : if(nf == f){
1451 : grid.end_marking();
1452 : return true;
1453 : }
1454 :
1455 : grid.mark(nf);
1456 :
1457 : // add unmarked edges to the stack and mark them, if they
1458 : // are connected to a marked vertex
1459 : CollectEdges(edges, grid, nf);
1460 :
1461 : for(size_t i = 0; i < edges.size(); ++i){
1462 : Edge* e = edges[i];
1463 : if(!grid.is_marked(e)){
1464 : grid.mark(e);
1465 : if(grid.is_marked(e->vertex(0)) || grid.is_marked(e->vertex(1))){
1466 : stk.push(e);
1467 : }
1468 : }
1469 : }
1470 : }
1471 :
1472 : grid.end_marking();
1473 : return false;
1474 : }
1475 : */
|
