Line data Source code
1 : /*
2 : * Copyright (c) 2009-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #include <vector>
34 : #include "edge_util.h"
35 : #include "lib_grid/grid/grid_util.h"
36 : #include "vertex_util.h"
37 : #include "face_util.h"
38 : #include "lib_grid/algorithms/debug_util.h"
39 : #include "lib_grid/refinement/regular_refinement.h"
40 : #include "lib_grid/refinement/projectors/plane_cut_projector.h"
41 : #include "common/util/vec_for_each.h"
42 : #include "common/util/vec_for_each.h"
43 :
44 : using namespace std;
45 :
46 : namespace ug
47 : {
48 :
49 : ////////////////////////////////////////////////////////////////////////
50 0 : int GetEdgeIndex(Face* f, Edge* e)
51 : {
52 : uint numEdges = f->num_edges();
53 0 : EdgeDescriptor ed;
54 0 : for(uint i = 0; i < numEdges; ++i)
55 : {
56 0 : f->edge_desc(i, ed);
57 0 : if(CompareVertices(e, &ed))
58 0 : return (int)i;
59 : }
60 : return -1;
61 : }
62 :
63 : ////////////////////////////////////////////////////////////////////////
64 0 : int GetEdgeIndex(Volume* vol, Edge* e)
65 : {
66 0 : uint numEdges = vol->num_edges();
67 0 : EdgeDescriptor ed;
68 0 : for(uint i = 0; i < numEdges; ++i)
69 : {
70 0 : vol->edge_desc(i, ed);
71 0 : if(CompareVertices(e, &ed))
72 0 : return (int)i;
73 : }
74 : return -1;
75 : }
76 :
77 : ////////////////////////////////////////////////////////////////////////
78 0 : bool IsBoundaryEdge(Grid& grid, Edge* e,
79 : Grid::face_traits::callback funcIsSurfFace)
80 : {
81 : int counter = 0;
82 0 : if(grid.option_is_enabled(EDGEOPT_STORE_ASSOCIATED_FACES))
83 : {
84 0 : for(Grid::AssociatedFaceIterator iter = grid.associated_faces_begin(e);
85 0 : iter != grid.associated_faces_end(e); ++iter)
86 : {
87 0 : if(funcIsSurfFace(*iter))
88 0 : ++counter;
89 0 : if(counter > 1)
90 : return false;
91 : }
92 : }
93 : else
94 : {
95 : // fill a vector using a helper function
96 : vector<Face*> faces;
97 0 : CollectFaces(faces, grid, e, false);
98 0 : for(size_t i = 0; i < faces.size(); ++i){
99 0 : if(funcIsSurfFace(faces[i]))
100 0 : ++counter;
101 0 : if(counter > 1)
102 : return false;
103 : }
104 0 : }
105 :
106 0 : if(counter == 1)
107 : return true;
108 : return false;
109 : }
110 :
111 : ////////////////////////////////////////////////////////////////////////
112 0 : bool IsBoundaryEdge2D(Grid& grid, Edge* e)
113 : {
114 : // get the number of connected faces. if only one face is connected then
115 : // the edge is considered to be a boundary edge.
116 : /* int counter = 0;
117 : if(grid.option_is_enabled(EDGEOPT_STORE_ASSOCIATED_FACES))
118 : {
119 : for(Grid::AssociatedFaceIterator iter = grid.associated_faces_begin(e);
120 : iter != grid.associated_faces_end(e); ++iter)
121 : ++counter;
122 :
123 : if(counter == 1)
124 : return true;
125 : }
126 : else
127 : {
128 : // fill a vector using a helper function
129 : vector<Face*> vFaces;
130 : CollectFaces(vFaces, grid, e, false);
131 : if(vFaces.size() == 1)
132 : return true;
133 : }
134 : return false;*/
135 0 : return NumAssociatedFaces(grid, e) == 1;
136 : }
137 :
138 : ////////////////////////////////////////////////////////////////////////
139 0 : bool IsBoundaryEdge3D(Grid& grid, Edge* e)
140 : {
141 0 : if(!grid.option_is_enabled(VOLOPT_AUTOGENERATE_FACES)){
142 : UG_LOG("WARNING in IsBoundaryEdge3D: Autoenabling VOLOPT_AUTOGENERATE_FACES.\n");
143 0 : grid.enable_options(VOLOPT_AUTOGENERATE_FACES);
144 : }
145 :
146 : // check whether an associated face is a boundary face
147 0 : if(grid.option_is_enabled(EDGEOPT_STORE_ASSOCIATED_FACES))
148 : {
149 0 : for(Grid::AssociatedFaceIterator iter = grid.associated_faces_begin(e);
150 0 : iter != grid.associated_faces_end(e); ++iter)
151 : {
152 0 : if(IsBoundaryFace3D(grid, *iter))
153 : return true;
154 : }
155 : }
156 : else
157 : {
158 : // fill a vector using a helper function
159 : vector<Face*> vFaces;
160 0 : CollectFaces(vFaces, grid, e, false);
161 0 : for(size_t i = 0; i < vFaces.size(); ++i){
162 0 : if(IsBoundaryFace3D(grid, vFaces[i]))
163 : return true;
164 : }
165 0 : }
166 : return false;
167 : }
168 :
169 0 : bool LiesOnBoundary(Grid& grid, Edge* e)
170 : {
171 : // first check whether the edge is a 2d boundary element
172 0 : if((grid.num<Face>() > 0) && IsBoundaryEdge2D(grid, e)){
173 : return true;
174 : }
175 :
176 : // since it isn't a 2d boundary element, it might be a 3d boundary element
177 0 : if((grid.num<Volume>() > 0) && IsBoundaryEdge3D(grid, e))
178 : return true;
179 :
180 : // ok - it isn't a boundary element
181 : return false;
182 : }
183 :
184 : ////////////////////////////////////////////////////////////////////////
185 : // GetAssociatedFaces
186 0 : int GetAssociatedFaces(Face** facesOut, Grid& grid,
187 : Edge* e, int maxNumFaces)
188 : {
189 0 : if(grid.option_is_enabled(EDGEOPT_STORE_ASSOCIATED_FACES))
190 : {
191 : int counter = 0;
192 0 : Grid::AssociatedFaceIterator iterEnd = grid.associated_faces_end(e);
193 0 : for(Grid::AssociatedFaceIterator iter = grid.associated_faces_begin(e);
194 0 : iter != iterEnd; ++iter)
195 : {
196 0 : if(counter < maxNumFaces)
197 0 : facesOut[counter] = *iter;
198 :
199 0 : counter++;
200 : }
201 : return counter;
202 : }
203 : else
204 : {
205 : // we're using grid::mark for maximal speed.
206 : //grid.begin_marking();
207 : // mark the end-points of the edge
208 : //grid.mark(e->vertex(0));
209 : //grid.mark(e->vertex(1));
210 :
211 0 : Vertex* v0 = e->vertex(0);
212 0 : Vertex* v1 = e->vertex(1);
213 :
214 : // we have to find the triangles 'by hand'
215 : // iterate over all associated faces of vertex 0
216 : int counter = 0;
217 0 : Grid::AssociatedFaceIterator iterEnd = grid.associated_faces_end(v0);
218 0 : for(Grid::AssociatedFaceIterator iter = grid.associated_faces_begin(v0);
219 0 : iter != iterEnd; ++iter)
220 : {
221 0 : Face* tf = *iter;
222 0 : uint numVrts = tf->num_vertices();
223 : int numMarked = 0;
224 0 : for(uint i = 0; i < numVrts; ++i){
225 0 : Vertex* v = tf->vertex(i);
226 0 : if((v == v0) || (v == v1))
227 0 : numMarked++;
228 : }
229 :
230 0 : if(numMarked > 1)
231 : {
232 : // the face is connected with the edge
233 0 : if(counter < maxNumFaces)
234 0 : facesOut[counter] = tf;
235 0 : counter++;
236 : }
237 : }
238 : // done with marking
239 : //grid.end_marking();
240 :
241 : return counter;
242 : }
243 : }
244 :
245 : ////////////////////////////////////////////////////////////////////////
246 : // NumAssociatedFaces
247 0 : int NumAssociatedFaces(Grid& grid, Edge* e)
248 : {
249 0 : if(grid.option_is_enabled(EDGEOPT_STORE_ASSOCIATED_FACES))
250 : {
251 : int counter = 0;
252 0 : Grid::AssociatedFaceIterator iterEnd = grid.associated_faces_end(e);
253 0 : for(Grid::AssociatedFaceIterator iter = grid.associated_faces_begin(e);
254 0 : iter != iterEnd; ++iter)
255 : {
256 0 : counter++;
257 : }
258 : return counter;
259 : }
260 : else
261 : {
262 : // we're using grid::mark for maximal speed.
263 0 : grid.begin_marking();
264 : // mark the end-points of the edge
265 0 : grid.mark(e->vertex(0));
266 0 : grid.mark(e->vertex(1));
267 :
268 : // we have to find the triangles 'by hand'
269 : // iterate over all associated faces of vertex 0
270 : int counter = 0;
271 0 : Vertex* v = e->vertex(0);
272 0 : Grid::AssociatedFaceIterator iterEnd = grid.associated_faces_end(v);
273 0 : for(Grid::AssociatedFaceIterator iter = grid.associated_faces_begin(v);
274 0 : iter != iterEnd; ++iter)
275 : {
276 0 : Face* tf = *iter;
277 0 : uint numVrts = tf->num_vertices();
278 : int numMarked = 0;
279 0 : for(uint i = 0; i < numVrts; ++i){
280 0 : if(grid.is_marked(tf->vertex(i)))
281 0 : numMarked++;
282 : }
283 :
284 0 : if(numMarked > 1)
285 : {
286 : // the face is connected with the edge
287 0 : counter++;
288 : }
289 : }
290 : // done with marking
291 0 : grid.end_marking();
292 :
293 : return counter;
294 : }
295 : }
296 :
297 : ////////////////////////////////////////////////////////////////////////
298 0 : Edge* GetConnectingEdge(Grid& grid, Face* f1, Face* f2)
299 : {
300 : Grid::edge_traits::secure_container edges1, edges2;
301 : grid.associated_elements(edges1, f1);
302 : grid.associated_elements(edges2, f2);
303 :
304 0 : for(size_t i = 0; i < edges1.size(); ++i){
305 0 : for(size_t j = 0; j < edges2.size(); ++j){
306 0 : if(edges1[i] == edges2[j])
307 : return edges1[i];
308 : }
309 : }
310 : return NULL;
311 : }
312 :
313 : ////////////////////////////////////////////////////////////////////////
314 0 : int CalculateNormal(vector3& vNormOut, Grid& grid, Edge* e,
315 : Grid::AttachmentAccessor<Vertex, APosition>& aaPos,
316 : Grid::AttachmentAccessor<Face, ANormal>* paaNormFACE)
317 : {
318 : Face* f[2];
319 :
320 0 : int numFaces = GetAssociatedFaces(f, grid, e, 2);
321 :
322 0 : switch(numFaces){
323 :
324 : case 0:{ // if there are no associated faces.
325 : // we'll assume that the edge lies in the xy plane and return its normal
326 : vector3 dir;
327 0 : VecSubtract(dir, aaPos[e->vertex(1)], aaPos[e->vertex(0)]);
328 0 : vNormOut.x() = dir.y();
329 0 : vNormOut.y() = -dir.x();
330 0 : vNormOut.z() = 0;
331 0 : VecNormalize(vNormOut, vNormOut);
332 0 : }break;
333 :
334 0 : case 1: // if there is one face, the normal will be set to the faces normal
335 0 : if(paaNormFACE)
336 0 : vNormOut = (*paaNormFACE)[f[0]];
337 : else{
338 0 : CalculateNormal(vNormOut, f[0], aaPos);
339 : }
340 : break;
341 :
342 0 : default: // there are at least 2 associated faces
343 0 : if(paaNormFACE)
344 0 : VecAdd(vNormOut, (*paaNormFACE)[f[0]], (*paaNormFACE)[f[1]]);
345 : else{
346 : vector3 fn0, fn1;
347 0 : CalculateNormalNoNormalize(fn0, f[0], aaPos);
348 0 : CalculateNormalNoNormalize(fn1, f[1], aaPos);
349 : VecAdd(vNormOut, fn0, fn1);
350 : }
351 0 : VecNormalize(vNormOut, vNormOut);
352 : break;
353 : }
354 :
355 0 : return numFaces;
356 : }
357 :
358 0 : int CalculateNormalNoNormalize(vector3& vNormOut, Grid& grid, Edge* e,
359 : Grid::VertexAttachmentAccessor<APosition>& aaPos,
360 : Grid::FaceAttachmentAccessor<ANormal>* paaNormFACE)
361 : {
362 : Face* f[2];
363 :
364 0 : int numFaces = GetAssociatedFaces(f, grid, e, 2);
365 :
366 0 : switch(numFaces){
367 :
368 : case 0:{ // if there are no associated faces.
369 : // we'll assume that the edge lies in the xy plane and return its normal
370 : vector3 dir;
371 0 : VecSubtract(dir, aaPos[e->vertex(1)], aaPos[e->vertex(0)]);
372 0 : vNormOut.x() = dir.y();
373 0 : vNormOut.y() = -dir.x();
374 0 : vNormOut.z() = 0;
375 0 : }break;
376 :
377 0 : case 1: // if there is one face, the normal will be set to the faces normal
378 0 : if(paaNormFACE)
379 0 : vNormOut = (*paaNormFACE)[f[0]];
380 : else{
381 0 : CalculateNormalNoNormalize(vNormOut, f[0], aaPos);
382 : }
383 : break;
384 :
385 0 : default: // there are at least 2 associated faces
386 0 : if(paaNormFACE)
387 0 : VecAdd(vNormOut, (*paaNormFACE)[f[0]], (*paaNormFACE)[f[1]]);
388 : else{
389 : vector3 fn0, fn1;
390 0 : CalculateNormalNoNormalize(fn0, f[0], aaPos);
391 0 : CalculateNormalNoNormalize(fn1, f[1], aaPos);
392 : VecAdd(vNormOut, fn0, fn1);
393 : VecScale(vNormOut, vNormOut, 0.5);
394 : }
395 0 : VecNormalize(vNormOut, vNormOut);
396 : break;
397 : }
398 :
399 0 : return numFaces;
400 : }
401 :
402 : ////////////////////////////////////////////////////////////////////////
403 0 : bool CollapseEdge(Grid& grid, Edge* e, Vertex* newVrt)
404 : {
405 : // prepare the grid, so that we may perform Grid::replace_vertex.
406 : // create collapse geometries first and delete old geometries.
407 : Grid::face_traits::secure_container assFaces;
408 0 : if(grid.num_faces() > 0){
409 : // collect adjacent faces
410 : grid.associated_elements(assFaces, e);
411 :
412 : // a vector thats is used to store the collapse-geometries.
413 : vector<Face*> vNewFaces;
414 :
415 0 : for(uint i = 0; i < assFaces.size(); ++i)
416 : {
417 :
418 : Face* f = assFaces[i];
419 0 : int eInd = GetEdgeIndex(f, e);
420 :
421 : // create the collapse-geometry
422 0 : f->collapse_edge(vNewFaces, eInd, newVrt);
423 :
424 : // register the new faces
425 0 : for(uint j = 0; j < vNewFaces.size(); ++j)
426 0 : grid.register_element(vNewFaces[j], f);
427 :
428 : // before we erase the face, we first have to notify whether
429 : // edges were merged
430 0 : if(f->num_edges() == 3){
431 : // two edges will be merged. we have to inform the grid.
432 0 : Vertex* conVrt = GetConnectedVertex(e, f);
433 : // now get the edge between conVrt and newVrt
434 0 : Edge* target = grid.get_edge(conVrt, newVrt);
435 : // now get the two old edges
436 0 : Edge* e1 = grid.get_edge(f, (eInd + 1) % 3);
437 0 : Edge* e2 = grid.get_edge(f, (eInd + 2) % 3);
438 : grid.objects_will_be_merged(target, e1, e2);
439 : }
440 : }
441 0 : }
442 :
443 : Grid::volume_traits::secure_container assVols;
444 0 : if(grid.num_volumes() > 0){
445 : // used to identify merge-faces
446 : Grid::edge_traits::secure_container assEdges;
447 :
448 : // collect adjacent volumes
449 : grid.associated_elements(assVols, e);
450 :
451 : // a vector thats used to store the collapse-geometries.
452 : vector<Volume*> vNewVolumes;
453 :
454 0 : for(uint i = 0; i < assVols.size(); ++i)
455 : {
456 : Volume* v = assVols[i];
457 : // create the collapse-geometry
458 0 : int eInd = GetEdgeIndex(v, e);
459 0 : v->collapse_edge(vNewVolumes, eInd, newVrt);
460 :
461 : // register the new volumes
462 0 : for(uint j = 0; j < vNewVolumes.size(); ++j)
463 0 : grid.register_element(vNewVolumes[j], v);
464 :
465 : // if v is a tetrahedron, two faces will be merged
466 0 : if(v->num_vertices() == 4){
467 : // find the opposing edge of e
468 : Edge* oppEdge = NULL;
469 : grid.associated_elements(assEdges, v);
470 0 : for_each_in_vec(Edge* te, assEdges){
471 0 : if(!GetSharedVertex(e, te)){
472 : oppEdge = te;
473 : break;
474 : }
475 : }end_for;
476 :
477 0 : if(oppEdge){
478 0 : Face* f0 = grid.get_face(FaceDescriptor(
479 0 : e->vertex(0),
480 0 : oppEdge->vertex(0),
481 0 : oppEdge->vertex(1)));
482 0 : Face* f1 = grid.get_face(FaceDescriptor(
483 0 : e->vertex(1),
484 0 : oppEdge->vertex(0),
485 0 : oppEdge->vertex(1)));
486 0 : Face* fNew = grid.get_face(FaceDescriptor(
487 : newVrt,
488 0 : oppEdge->vertex(0),
489 0 : oppEdge->vertex(1)));
490 0 : if(f0 && f1 && fNew)
491 : grid.objects_will_be_merged(fNew, f0, f1);
492 : }
493 : }
494 : }
495 0 : }
496 :
497 : // store the end-points of e
498 : Vertex* v[2];
499 0 : v[0] = e->vertex(0);
500 0 : v[1] = e->vertex(1);
501 :
502 : // notify the grid that the endpoints will be merged
503 : grid.objects_will_be_merged(newVrt, v[0], v[1]);
504 :
505 : // erase e, associated faces and associated volumes
506 0 : for_each_in_vec(Volume* v, assVols){
507 0 : grid.erase(v);
508 : }end_for;
509 :
510 0 : for_each_in_vec(Face* f, assFaces){
511 0 : grid.erase(f);
512 : }end_for;
513 :
514 0 : grid.erase(e);
515 :
516 : // perform replace_vertex for both endpoints
517 0 : for(uint i = 0; i < 2; ++i)
518 : {
519 0 : if(v[i] != newVrt)
520 0 : grid.replace_vertex(v[i], newVrt);
521 : }
522 :
523 0 : return true;
524 : }
525 :
526 : ////////////////////////////////////////////////////////////////////////
527 0 : bool EdgeCollapseIsValid(Grid& grid, Edge* e)
528 : {
529 : //TODO: Make sure that this is sufficient for geometries including Quads.
530 : //TODO: Check validity for volumes.
531 : // in order to not destroy the topology of the grid, a collapse is
532 : // only valid if no three edges build a triangle that does not exist
533 : // in the grid.
534 :
535 : // first we need all vertices that are connected with end-points of e.
536 : vector<Vertex*> vVertices1;
537 : vector<Vertex*> vVertices2;
538 :
539 0 : CollectNeighbors(vVertices1, grid, e->vertex(0)); // e->vertex(0) is not contained in vVertices1!
540 0 : CollectNeighbors(vVertices2, grid, e->vertex(1)); // e->vertex(1) is not contained in vVertices2!
541 :
542 : // we need access to the triangles connected with e.
543 : vector<Face*> vFaces;
544 0 : CollectFaces(vFaces, grid, e);
545 :
546 : // this face descriptor will be needed during the algorithm.
547 0 : FaceDescriptor fd(3);
548 0 : fd.set_vertex(0, e->vertex(0));
549 0 : fd.set_vertex(1, e->vertex(1));
550 :
551 : // now check for each vertex in vVertices1 if it also exists in vVertices2.
552 0 : for(uint i = 0; i < vVertices1.size(); ++i)
553 : {
554 0 : Vertex* v1 = vVertices1[i];
555 :
556 0 : if(v1 != e->vertex(1))
557 : {
558 0 : for(uint j = 0; j < vVertices2.size(); ++j)
559 : {
560 0 : Vertex* v2 = vVertices2[j];
561 :
562 0 : if(v1 == v2)
563 : {
564 : // v1 and v2 exist in both arrays.
565 : // check if a triangle exists that connects e with v1.
566 : // it is sufficient to search in vFaces.
567 : bool bGotOne = false;
568 : fd.set_vertex(2, v1);
569 0 : for(uint k = 0; k < vFaces.size(); ++k)
570 : {
571 0 : if(CompareVertices(vFaces[k], &fd))
572 : {
573 : bGotOne = true;
574 : break;
575 : }
576 : }
577 :
578 : // if there was none, the collapse would be illegal.
579 0 : if(!bGotOne)
580 : return false;
581 : }
582 : }
583 : }
584 : }
585 :
586 : // everything is fine.
587 : return true;
588 0 : }
589 :
590 : ////////////////////////////////////////////////////////////////////////
591 0 : bool CreateEdgeSplitGeometry(Grid& destGrid, Grid& srcGrid, Edge* e,
592 : Vertex* newVertex, AVertex* paAssociatedVertices)
593 : {
594 :
595 0 : if((paAssociatedVertices == NULL) && (&destGrid != &srcGrid))
596 : return false;
597 :
598 : vector<Vertex*> vVrts;
599 :
600 : // If paAssociatedVertices is specified, we'll have to find the vertices
601 : // in destGrid for each element that is adjacent to e. We then have
602 : // to store these vertices in vVrts and pass them to the
603 : // split-routine of the element.
604 :
605 : // the attachment-accessor
606 : Grid::VertexAttachmentAccessor<AVertex> aaAssociatedVertices;
607 :
608 0 : if(paAssociatedVertices != NULL)
609 : {
610 : AVertex& aAssociatedVertices = *paAssociatedVertices;
611 :
612 : // check if aVertex is properly attached.
613 0 : if(!srcGrid.has_vertex_attachment(aAssociatedVertices))
614 : // attach it and initialize its values.
615 0 : srcGrid.attach_to_vertices_dv(aAssociatedVertices, NULL, false);
616 :
617 : // initialize the attachment-accessor
618 : aaAssociatedVertices.access(srcGrid, aAssociatedVertices);
619 : }
620 :
621 : // we will store the substitute-vertices in this vector - if they are needed.
622 : vector<Vertex*> vSubstituteVertices;
623 :
624 : // split the edge
625 : {
626 : // simply create two new edges
627 : Edge* parent = e;
628 : // the grids do not match then we can't pass e as a parent
629 0 : if(&srcGrid != &destGrid)
630 : parent = NULL;
631 :
632 0 : if(paAssociatedVertices){
633 0 : destGrid.create<RegularEdge>(EdgeDescriptor(aaAssociatedVertices[e->vertex(0)], newVertex), parent);
634 0 : destGrid.create<RegularEdge>(EdgeDescriptor(newVertex, aaAssociatedVertices[e->vertex(1)]), parent);
635 : }
636 : else{
637 0 : destGrid.create<RegularEdge>(EdgeDescriptor(e->vertex(0), newVertex), parent);
638 0 : destGrid.create<RegularEdge>(EdgeDescriptor(newVertex, e->vertex(1)), parent);
639 : }
640 : }
641 :
642 : // split faces
643 : {
644 : // collect all faces which are associated with e
645 : vector<Face*> vOldFaces;
646 0 : CollectFaces(vOldFaces, srcGrid, e, false);
647 :
648 : // we will collect new faces in this vector
649 : vector<Face*> vNewFaces;
650 :
651 : // iterate through those faces and split each.
652 : // If vertices are missing in destGrid, create them first.
653 0 : for(vector<Face*>::iterator oldFaceIter = vOldFaces.begin();
654 0 : oldFaceIter != vOldFaces.end(); ++oldFaceIter)
655 : {
656 0 : Face* oldFace = *oldFaceIter;
657 0 : uint numVrts = oldFace->num_vertices();
658 :
659 : // get the index of e in oldFace
660 0 : int edgeIndex = GetEdgeIndex(oldFace, e);
661 :
662 : // clear the new faces-vec
663 : vNewFaces.clear();
664 :
665 : // get the substitute-vertices if they are required
666 0 : if(paAssociatedVertices != NULL)
667 : {
668 0 : if(numVrts > vSubstituteVertices.size())
669 0 : vSubstituteVertices.resize(numVrts);
670 :
671 : // check for each vertex in oldFace, if an associated vertex exists in destGrid.
672 : // if not create a new one by cloning the associated one in srcGrid.
673 : // store the vertices in vSubstituteVertices
674 : {
675 0 : for(uint i = 0; i < numVrts; ++i)
676 : {
677 0 : if(aaAssociatedVertices[oldFace->vertex(i)] == NULL)
678 0 : aaAssociatedVertices[oldFace->vertex(i)] = *destGrid.create_by_cloning(oldFace->vertex(i));
679 0 : vSubstituteVertices[i] = aaAssociatedVertices[oldFace->vertex(i)];
680 : }
681 : }
682 :
683 : // create the new faces by splitting the old face. use substitutes.
684 0 : oldFace->create_faces_by_edge_split(edgeIndex, newVertex, vNewFaces, &vSubstituteVertices.front());
685 : }
686 : else
687 : {
688 : // create the new faces by splitting the old face.
689 : // no substitutes required
690 0 : oldFace->create_faces_by_edge_split(edgeIndex, newVertex, vNewFaces);
691 : }
692 :
693 : // register all new faces at destGrid
694 : {
695 : Face* pParent = NULL;
696 0 : if(&srcGrid == &destGrid)
697 : pParent = oldFace;
698 :
699 0 : for(vector<Face*>::iterator iter = vNewFaces.begin();
700 0 : iter != vNewFaces.end(); ++iter)
701 : {
702 0 : destGrid.register_element(*iter, pParent);
703 0 : if(pParent)
704 0 : destGrid.pass_on_values(pParent, *iter);
705 : }
706 : }
707 : }
708 0 : }
709 :
710 : // split volumes
711 0 : if(srcGrid.num<Volume>() > 0)
712 : {
713 : // this vector will be used to specify on which edge a vertex has to be inserted
714 : vector<Vertex*> edgeVrts;
715 :
716 : // collect all volumes associated with the edge
717 : vector<Volume*> vols, newVols;
718 0 : CollectVolumes(vols, srcGrid, e);
719 :
720 0 : for(size_t i_vols = 0; i_vols < vols.size(); ++i_vols)
721 : {
722 0 : Volume* oldVol = vols[i_vols];
723 0 : uint numVrts = oldVol->num_vertices();
724 :
725 : newVols.clear();
726 :
727 : // get the index of e in oldFace and fill the edgeVrts array
728 : edgeVrts.clear();
729 0 : for(size_t i_edge = 0; i_edge < oldVol->num_edges(); ++i_edge)
730 : {
731 0 : if(srcGrid.get_edge(oldVol, i_edge) == e)
732 0 : edgeVrts.push_back(newVertex);
733 : else
734 0 : edgeVrts.push_back(NULL);
735 : }
736 :
737 : // if refine creates a new vertex in the center of the volume,
738 : // it will be stored in this var.
739 0 : Vertex* newVolVrt = NULL;
740 :
741 : // get the substitute-vertices if they are required
742 0 : if(paAssociatedVertices != NULL)
743 : {
744 0 : if(numVrts > vSubstituteVertices.size())
745 0 : vSubstituteVertices.resize(numVrts);
746 :
747 : // check for each vertex in oldFace, if an associated vertex exists in destGrid.
748 : // if not create a new one by cloning the associated one in srcGrid.
749 : // store the vertices in vSubstituteVertices
750 : {
751 0 : for(uint i = 0; i < numVrts; ++i)
752 : {
753 0 : if(aaAssociatedVertices[oldVol->vertex(i)] == NULL)
754 0 : aaAssociatedVertices[oldVol->vertex(i)] = *destGrid.create_by_cloning(oldVol->vertex(i));
755 0 : vSubstituteVertices[i] = aaAssociatedVertices[oldVol->vertex(i)];
756 : }
757 : }
758 :
759 : // create the new faces by splitting the old face. use substitutes.
760 0 : oldVol->refine(newVols, &newVolVrt, &edgeVrts.front(), NULL, NULL,
761 0 : RegularVertex(), &vSubstituteVertices.front());
762 : }
763 : else
764 : {
765 : // create the new faces by splitting the old face.
766 : // no substitutes required
767 0 : oldVol->refine(newVols, &newVolVrt, &edgeVrts.front(), NULL, NULL,
768 0 : RegularVertex(), NULL);
769 : }
770 :
771 : // register all new vertices and volumes at destGrid
772 : {
773 : Volume* pParent = NULL;
774 0 : if(&srcGrid == &destGrid)
775 : pParent = oldVol;
776 :
777 0 : if(newVolVrt)
778 : destGrid.register_element(newVolVrt, pParent);
779 :
780 0 : for(vector<Volume*>::iterator iter = newVols.begin();
781 0 : iter != newVols.end(); ++iter)
782 : {
783 0 : destGrid.register_element(*iter, pParent);
784 0 : if(pParent)
785 0 : destGrid.pass_on_values(pParent, *iter);
786 : }
787 : }
788 : }
789 0 : }
790 : return true;
791 0 : }
792 :
793 0 : Edge* SwapEdge(Grid& grid, Edge* e)
794 : {
795 : // get the associated faces.
796 : // Only 2 associated faces are allowed.
797 : Face* f[2];
798 0 : if(GetAssociatedFaces(f, grid, e, 2) != 2){
799 : UG_LOG("Swap Failed: #neighbor-faces != 2.\n");
800 0 : return NULL;
801 : }
802 :
803 : // make sure that both faces are triangles
804 0 : if((f[0]->num_vertices() != 3) || (f[1]->num_vertices() != 3)){
805 : UG_LOG("Swap Failed: At least one neighbor-face is not a triangle.\n");
806 0 : return NULL;
807 : }
808 :
809 : // begin marking
810 0 : grid.begin_marking();
811 :
812 : // mark the end-points of the edge
813 0 : grid.mark(e->vertex(0));
814 0 : grid.mark(e->vertex(1));
815 :
816 : // get the two vertices that will be connected by the new edge
817 : Vertex* v[2];
818 : int vrtInd[2];
819 0 : for(int j = 0; j < 2; ++j){
820 0 : v[j] = NULL;
821 0 : for(int i = 0; i < 3; ++i){
822 0 : Vertex* vrt = f[j]->vertex(i);
823 0 : if(!grid.is_marked(vrt)){
824 0 : vrtInd[j] = i;
825 0 : v[j] = vrt;
826 0 : break;
827 : }
828 : }
829 : }
830 :
831 : // we're done with marking
832 0 : grid.end_marking();
833 :
834 : // make sure that both vertices have been found.
835 0 : if(!(v[0] && v[1])){
836 : UG_LOG("Swap Failed: connected vertices were not found correctly.\n");
837 0 : return NULL;
838 : }
839 :
840 : // make sure that no edge exists between v[0] and v[1]
841 0 : if(grid.get_edge(v[0], v[1])){
842 : UG_LOG("Swap Failed: New edge already exists in the grid.\n");
843 0 : return NULL;
844 : }
845 :
846 : // the indices of the marked points in the first triangle
847 0 : int ind1 = (vrtInd[0] + 1) % 3;
848 0 : int ind2 = (vrtInd[0] + 2) % 3;
849 :
850 : // create the new edge
851 0 : Edge* nEdge = *grid.create_by_cloning(e, EdgeDescriptor(v[0], v[1]), e);
852 :
853 : // create the new faces
854 0 : grid.create<Triangle>(TriangleDescriptor(v[0], f[0]->vertex(ind1), v[1]), f[0]);
855 0 : grid.create<Triangle>(TriangleDescriptor(f[0]->vertex(ind2), v[0], v[1]), f[1]);
856 :
857 : // erase the old faces
858 0 : grid.erase(f[0]);
859 0 : grid.erase(f[1]);
860 :
861 : // erase the old edge
862 0 : grid.erase(e);
863 :
864 : return nEdge;
865 : }
866 :
867 : ////////////////////////////////////////////////////////////////////////
868 0 : bool CutEdgesWithPlane(Selector& sel, const vector3& p, const vector3& n,
869 : APosition& aPos)
870 : {
871 0 : if(!sel.grid()){
872 : UG_LOG("ERROR in CutEdgesWithPlane: sel has to be assigned to a grid.\n");
873 0 : return false;
874 : }
875 :
876 : Grid& grid = *sel.grid();
877 :
878 0 : if(!grid.has_vertex_attachment(aPos)){
879 : UG_LOG("ERROR in CutEdgesWithPlane: aPos has to be attached to the vertices of the grid.\n");
880 0 : return false;
881 : }
882 :
883 : Grid::VertexAttachmentAccessor<APosition> aaPos(grid, aPos);
884 :
885 : // used for plane-intersection
886 : number t;
887 : vector3 v;
888 :
889 : // iterate through all edges and deselect all that do not intersect the plane
890 : // deselect all vertices, faces and volumes, too.
891 0 : sel.clear<Vertex>();
892 0 : sel.clear<Face>();
893 0 : sel.clear<Volume>();
894 :
895 : EdgeIterator iter = sel.begin<Edge>();
896 0 : while(iter != sel.end<Edge>()){
897 : Edge* e = *iter;
898 : iter++;
899 :
900 : // check whether the edge intersects the plane
901 : vector3 rayDir;
902 0 : VecSubtract(rayDir, aaPos[e->vertex(1)], aaPos[e->vertex(0)]);
903 :
904 0 : bool bIntersect = RayPlaneIntersection(v, t, aaPos[e->vertex(0)],
905 : rayDir, p, n);
906 0 : if(!bIntersect || t < SMALL || t > 1. - SMALL)
907 : {
908 : // it doesn't. Remove it from the selector
909 0 : sel.deselect(e);
910 : }
911 : }
912 :
913 : // refine all selected edges. RefinementCallbackEdgePlaneCut will insert
914 : // new vertices on the plane.
915 0 : PlaneCutProjector planeCutProjector(MakeGeometry3d(grid, aPos), p, n);
916 0 : const bool success = Refine(grid, sel, &planeCutProjector);
917 :
918 : // deselect all elements and edges which do not connect two selected vertices
919 0 : if(success){
920 : // deselect all vertices which are not very close to the plane
921 : VertexIterator vrtIter = sel.begin<Vertex>();
922 0 : while(vrtIter != sel.end<Vertex>()){
923 : Vertex* vrt = *vrtIter;
924 : ++vrtIter;
925 0 : if(DistancePointToPlane(aaPos[vrt], p, n) > SMALL)
926 0 : sel.deselect(vrt);
927 : }
928 :
929 : EdgeIterator iter = sel.begin<Edge>();
930 0 : while(iter != sel.end<Edge>()){
931 : Edge* e = *iter;
932 : iter++;
933 0 : if(!(sel.is_selected(e->vertex(0)) && sel.is_selected(e->vertex(1)))){
934 0 : sel.deselect(e);
935 : }
936 : }
937 :
938 0 : sel.clear<Face>();
939 0 : sel.clear<Volume>();
940 : }
941 :
942 : return success;
943 : }
944 :
945 : }// end of namespace
|
