Line data Source code
1 : /*
2 : * Copyright (c) 2012-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 <cassert>
34 : #include "common/profiler/profiler.h"
35 : #include "global_fractured_media_refiner.h"
36 : #include "lib_grid/algorithms/algorithms.h"
37 : #include "lib_grid/common_attachments.h"
38 : #include "lib_grid/file_io/file_io.h"
39 :
40 : //define PROFILE_GLOBAL_FRACTURED_DOMAIN_REFINER if you want to profile
41 : //the refinement code.
42 : //#define PROFILE_GLOBAL_FRACTURED_DOMAIN_REFINER
43 : #ifdef PROFILE_GLOBAL_FRACTURED_DOMAIN_REFINER
44 : #define GFDR_PROFILE_FUNC() PROFILE_FUNC_GROUP("grid")
45 : #define GFDR_PROFILE(name) PROFILE_BEGIN_GROUP(name, "grid")
46 : #define GFDR_PROFILE_END() PROFILE_END()
47 : #else
48 : #define GFDR_PROFILE_FUNC()
49 : #define GFDR_PROFILE(name)
50 : #define GFDR_PROFILE_END()
51 : #endif
52 :
53 : using namespace std;
54 :
55 : namespace ug
56 : {
57 :
58 : //template <class TAPosition>
59 : //GlobalFracturedMediaRefiner<TAPosition>::
60 0 : GlobalFracturedMediaRefiner::
61 0 : GlobalFracturedMediaRefiner(SPRefinementProjector projector) :
62 : IRefiner(projector),
63 0 : m_pMG(NULL),
64 0 : m_pSH(NULL)
65 : {
66 : // initialize m_aPos to the default position attachment
67 : //m_aPos = GetDefaultPositionAttachment<TAPosition>();
68 0 : }
69 :
70 :
71 : //template <class TAPosition>
72 : //GlobalFracturedMediaRefiner<TAPosition>::
73 0 : GlobalFracturedMediaRefiner::
74 0 : GlobalFracturedMediaRefiner(MultiGrid& mg, SPRefinementProjector projector) :
75 : IRefiner(projector),
76 0 : m_pSH(NULL)
77 : {
78 0 : m_pMG = NULL;
79 0 : assign_grid(mg);
80 : // initialize m_aPos to the default position attachment
81 : //m_aPos = GetDefaultPositionAttachment<TAPosition>();
82 0 : }
83 :
84 :
85 : //template <class TAPosition>
86 : //GlobalFracturedMediaRefiner<TAPosition>::
87 0 : GlobalFracturedMediaRefiner::
88 0 : ~GlobalFracturedMediaRefiner()
89 : {
90 0 : if(m_pMG)
91 0 : m_pMG->unregister_observer(this);
92 0 : }
93 :
94 :
95 : //template <class TAPosition>
96 : //void GlobalFracturedMediaRefiner<TAPosition>::
97 0 : void GlobalFracturedMediaRefiner::
98 : grid_to_be_destroyed(Grid* grid)
99 : {
100 0 : m_pMG = NULL;
101 0 : }
102 :
103 :
104 : //template <class TAPosition>
105 : //void GlobalFracturedMediaRefiner<TAPosition>::
106 0 : void GlobalFracturedMediaRefiner::
107 : assign_grid(MultiGrid& mg)
108 : {
109 0 : assign_grid(&mg);
110 0 : }
111 :
112 :
113 : //template <class TAPosition>
114 : //void GlobalFracturedMediaRefiner<TAPosition>::
115 0 : void GlobalFracturedMediaRefiner::
116 : assign_grid(MultiGrid* mg)
117 : {
118 0 : if(m_pMG){
119 0 : m_pMG->unregister_observer(this);
120 0 : m_marker.assign_grid(NULL);
121 0 : m_pMG = NULL;
122 : }
123 :
124 0 : if(mg){
125 0 : m_pMG = mg;
126 0 : set_message_hub(mg->message_hub());
127 0 : m_pMG->register_observer(this, OT_GRID_OBSERVER);
128 0 : m_marker.assign_grid(mg);
129 : }
130 0 : }
131 :
132 : //template <class TAPosition>
133 : //void GlobalFracturedMediaRefiner<TAPosition>::
134 0 : void GlobalFracturedMediaRefiner::
135 : mark_as_fracture(int subInd, bool isFracture)
136 : {
137 0 : if(subInd < 0)
138 : return;
139 :
140 0 : while(subInd > (int)m_subsetIsFracture.size())
141 0 : m_subsetIsFracture.push_back(false);
142 0 : if(subInd == (int)m_subsetIsFracture.size())
143 0 : m_subsetIsFracture.push_back(isFracture);
144 : else
145 : m_subsetIsFracture[subInd] = isFracture;
146 : }
147 :
148 : //template <class TAPosition>
149 : //bool GlobalFracturedMediaRefiner<TAPosition>::
150 0 : bool GlobalFracturedMediaRefiner::
151 : is_fracture(int subInd)
152 : {
153 0 : if((subInd < 0) || (subInd >= (int)m_subsetIsFracture.size()))
154 : return false;
155 : return m_subsetIsFracture[subInd];
156 : }
157 :
158 0 : void GlobalFracturedMediaRefiner::
159 : num_marked_edges_local(std::vector<int>& numMarkedEdgesOut)
160 : {
161 0 : num_marked_elems<Edge>(numMarkedEdgesOut);
162 0 : }
163 :
164 0 : void GlobalFracturedMediaRefiner::
165 : num_marked_faces_local(std::vector<int>& numMarkedFacesOut)
166 : {
167 0 : num_marked_elems<Face>(numMarkedFacesOut);
168 0 : }
169 :
170 0 : void GlobalFracturedMediaRefiner::
171 : num_marked_volumes_local(std::vector<int>& numMarkedVolsOut)
172 : {
173 0 : num_marked_elems<Volume>(numMarkedVolsOut);
174 0 : }
175 :
176 :
177 : template <class TElem>
178 0 : void GlobalFracturedMediaRefiner::
179 : num_marked_elems(std::vector<int>& numMarkedElemsOut)
180 : {
181 : numMarkedElemsOut.clear();
182 0 : if(!m_pMG)
183 : return;
184 0 : numMarkedElemsOut.resize(m_pMG->num_levels(), 0);
185 0 : if(m_pMG->num_levels() > 0)
186 0 : numMarkedElemsOut.back() = m_pMG->num<TElem>(m_pMG->top_level());
187 : }
188 :
189 : //template <class TAPosition>
190 : //void GlobalFracturedMediaRefiner<TAPosition>::
191 0 : void GlobalFracturedMediaRefiner::
192 : perform_refinement()
193 : {
194 : GFDR_PROFILE_FUNC();
195 : UG_DLOG(LIB_GRID, 1, "GlobalFracturedMediaRefiner\n");
196 :
197 : // get the grid
198 0 : if(!m_pMG)
199 0 : UG_THROW("No grid assigned!");
200 :
201 : MultiGrid& mg = *m_pMG;
202 :
203 : // adjust marks for refinement
204 0 : adjust_marks();
205 :
206 : // if a debug file was specified, we'll now save the marks to that file
207 0 : if(!m_adjustedMarksDebugFilename.empty())
208 0 : save_marks_to_file(m_adjustedMarksDebugFilename.c_str());
209 :
210 : // make sure that the required options are enabled.
211 0 : if(mg.num_volumes() > 0){
212 0 : if(!mg.option_is_enabled(VOLOPT_AUTOGENERATE_FACES))
213 : {
214 : LOG("WARNING in GlobalFracturedMediaRefiner::refine(): auto-enabling VOLOPT_AUTOGENERATE_FACES.\n");
215 0 : mg.enable_options(VOLOPT_AUTOGENERATE_FACES);
216 : }
217 : }
218 :
219 0 : if(mg.num_faces() > 0){
220 0 : if(!mg.option_is_enabled(FACEOPT_AUTOGENERATE_EDGES))
221 : {
222 : LOG("WARNING in GlobalFracturedMediaRefiner::refine(): auto-enabling FACEOPT_AUTOGENERATE_EDGES.\n");
223 0 : mg.enable_options(FACEOPT_AUTOGENERATE_EDGES);
224 : }
225 : }
226 :
227 0 : if(mg.num_levels() == 0)
228 0 : return;
229 :
230 : // the old top level
231 0 : int oldTopLevel = mg.num_levels() - 1;
232 0 : m_messageHub->post_message(GridMessage_Adaption(GMAT_GLOBAL_REFINEMENT_BEGINS,
233 0 : mg.get_grid_objects(oldTopLevel)));
234 :
235 0 : if(projector()->refinement_begins_requires_subgrid()){
236 0 : SubGrid<ConsiderAll> sg(mg.get_grid_objects(), ConsiderAll());
237 0 : projector()->refinement_begins(&sg);
238 : }
239 : else
240 0 : projector()->refinement_begins(NULL);
241 :
242 : UG_DLOG(LIB_GRID, 1, "REFINER: reserving memory...");
243 :
244 : // reserve enough memory to speed up the algo
245 : // todo: give better estimate based on fracture-subsets...
246 : GFDR_PROFILE(GFDR_Reserve);
247 : {
248 : int l = oldTopLevel;
249 :
250 : GFDR_PROFILE(GFDR_ReserveVrtData);
251 0 : mg.reserve<Vertex>(mg.num<Vertex>() +
252 0 : + mg.num<Vertex>(l) + mg.num<Edge>(l)
253 0 : + mg.num<Quadrilateral>(l) + mg.num<Hexahedron>(l));
254 : GFDR_PROFILE_END();
255 :
256 : GFDR_PROFILE(GFDR_ReserveEdgeData);
257 0 : mg.reserve<Edge>(mg.num<Edge>()
258 0 : + 2 * mg.num<Edge>(l) + 3 * mg.num<Triangle>(l)
259 0 : + 4 * mg.num<Quadrilateral>(l) + 3 * mg.num<Prism>(l)
260 0 : + mg.num<Tetrahedron>(l)
261 0 : + 4 * mg.num<Pyramid>(l) + 6 * mg.num<Hexahedron>(l));
262 : GFDR_PROFILE_END();
263 :
264 : GFDR_PROFILE(GFDR_ReserveFaceData);
265 0 : mg.reserve<Face>(mg.num<Face>()
266 0 : + 4 * mg.num<Face>(l) + 10 * mg.num<Prism>(l)
267 0 : + 8 * mg.num<Tetrahedron>(l)
268 0 : + 9 * mg.num<Pyramid>(l) + 12 * mg.num<Hexahedron>(l));
269 : GFDR_PROFILE_END();
270 :
271 : GFDR_PROFILE(GFDR_ReserveVolData);
272 0 : mg.reserve<Volume>(mg.num<Volume>()
273 0 : + 8 * mg.num<Tetrahedron>(l) + 8 * mg.num<Prism>(l)
274 0 : + 6 * mg.num<Pyramid>(l) + 8 * mg.num<Hexahedron>(l));
275 : GFDR_PROFILE_END();
276 : }
277 : GFDR_PROFILE_END();
278 : UG_DLOG(LIB_GRID, 1, " done.\n");
279 :
280 : UG_DLOG(LIB_GRID, 1, " refinement begins.\n");
281 : // notify derivates that refinement begins
282 0 : refinement_step_begins();
283 :
284 : // cout << "num marked edges: " << m_selMarks.num<Edge>() << endl;
285 : // cout << "num marked faces: " << m_selMarks.num<Face>() << endl;
286 :
287 : // we want to add new elements in a new layer.
288 : bool bHierarchicalInsertionWasEnabled = mg.hierarchical_insertion_enabled();
289 0 : if(!bHierarchicalInsertionWasEnabled)
290 0 : mg.enable_hierarchical_insertion(true);
291 :
292 :
293 : // some buffers
294 : vector<Vertex*> vVrts;
295 : vector<Vertex*> vEdgeVrts;
296 : vector<Vertex*> vFaceVrts;
297 : vector<Edge*> vEdges;
298 : vector<Face*> vFaces;
299 : vector<Volume*> vVols;
300 :
301 : // some repeatedly used objects
302 0 : EdgeDescriptor ed;
303 0 : FaceDescriptor fd;
304 0 : VolumeDescriptor vd;
305 :
306 : UG_DLOG(LIB_GRID, 1, " creating new vertices\n");
307 :
308 : // create new vertices from marked vertices
309 0 : for(VertexIterator iter = mg.begin<Vertex>(oldTopLevel);
310 0 : iter != mg.end<Vertex>(oldTopLevel); ++iter)
311 : {
312 0 : if(!refinement_is_allowed(*iter) || !m_marker.is_marked(*iter))
313 0 : continue;
314 :
315 : Vertex* v = *iter;
316 :
317 : // create a new vertex in the next layer.
318 : //GFDR_PROFILE(GFDR_Refine_CreatingVertices);
319 0 : Vertex* nVrt = *mg.create_by_cloning(v, v);
320 :
321 : // allow refCallback to calculate a new position
322 0 : if(m_projector.valid())
323 0 : m_projector->new_vertex(nVrt, v);
324 : //GFDR_PROFILE_END();
325 : }
326 :
327 :
328 : UG_DLOG(LIB_GRID, 1, " creating new edges\n");
329 :
330 : // create new vertices and edges from marked edges
331 0 : for(EdgeIterator iter = mg.begin<Edge>(oldTopLevel);
332 0 : iter != mg.end<Edge>(oldTopLevel); ++iter)
333 : {
334 0 : if(!refinement_is_allowed(*iter))
335 0 : continue;
336 :
337 : Edge* e = *iter;
338 :
339 0 : if(!(refinement_is_allowed(e->vertex(0)) && refinement_is_allowed(e->vertex(1))))
340 : {
341 0 : UG_THROW("Refinement has to be allowed for both corners of an edge, for"
342 : " which refinement is allowed.");
343 : }
344 :
345 : // associated vertices on next level.
346 : Vertex* substituteVrts[2];
347 0 : substituteVrts[0] = mg.get_child_vertex(e->vertex(0));
348 0 : substituteVrts[1] = mg.get_child_vertex(e->vertex(1));
349 :
350 : // if the face is not marked, we'll simply clone it to the next level
351 0 : if(!m_marker.is_marked(e)){
352 0 : ed.set_vertices(substituteVrts[0], substituteVrts[1]);
353 0 : mg.create_by_cloning(e, ed, e);
354 0 : continue;
355 : }
356 :
357 : //GFDR_PROFILE(GFDR_Refine_CreatingEdgeVertices);
358 : // create two new edges by edge-split
359 0 : RegularVertex* nVrt = *mg.create<RegularVertex>(e);
360 :
361 : // allow refCallback to calculate a new position
362 0 : if(m_projector.valid())
363 0 : m_projector->new_vertex(nVrt, e);
364 : //GFDR_PROFILE_END();
365 :
366 : // split the edge
367 : //GFDR_PROFILE(GFDR_Refine_CreatingEdges);
368 0 : e->refine(vEdges, nVrt, substituteVrts);
369 : assert((vEdges.size() == 2) && "RegularEdge refine produced wrong number of edges.");
370 0 : mg.register_element(vEdges[0], e);
371 0 : mg.register_element(vEdges[1], e);
372 : //GFDR_PROFILE_END();
373 : }
374 :
375 :
376 : UG_DLOG(LIB_GRID, 1, " creating new faces\n");
377 :
378 : // create new vertices and faces from marked faces
379 0 : for(FaceIterator iter = mg.begin<Face>(oldTopLevel);
380 0 : iter != mg.end<Face>(oldTopLevel); ++iter)
381 : {
382 0 : if(!refinement_is_allowed(*iter))
383 0 : continue;
384 :
385 : Face* f = *iter;
386 :
387 : // collect child-vertices
388 : vVrts.clear();
389 0 : for(uint j = 0; j < f->num_vertices(); ++j)
390 0 : vVrts.push_back(mg.get_child_vertex(f->vertex(j)));
391 :
392 : // collect the child vertices of marked and associated edges
393 : vEdgeVrts.clear();
394 0 : for(uint j = 0; j < f->num_edges(); ++j){
395 0 : Edge *edge=mg.get_edge(f,j);
396 0 : if(m_marker.is_marked(edge)){
397 0 : vEdgeVrts.push_back(mg.get_child_vertex(edge));
398 : }else{
399 0 : vEdgeVrts.push_back(NULL);
400 : }
401 : }
402 :
403 : // instead of not marked faces, the cloning condition is changed for those, of which no edges will be refined.
404 : // It means, a face will be cloned to the next level if none of its edges should be refined.
405 : // There are no such cases for testing. If such faces exist in your case, please test to ensure if it works for parallel computing.
406 0 : if(vEdgeVrts.size()==0){
407 0 : fd.set_num_vertices(vVrts.size());
408 0 : for(size_t i = 0; i < vVrts.size(); ++i)
409 0 : fd.set_vertex(i, vVrts[i]);
410 0 : mg.create_by_cloning(f, fd, f);
411 0 : continue;
412 0 : }
413 :
414 : //GFDR_PROFILE(GFDR_Refine_CreatingFaces);
415 : Vertex* newVrt;
416 0 : if(f->refine(vFaces, &newVrt, &vEdgeVrts.front(), NULL, &vVrts.front())){
417 : // if a new vertex was generated, we have to register it
418 0 : if(newVrt){
419 : //GFDR_PROFILE(GFDR_Refine_CreatingVertices);
420 0 : mg.register_element(newVrt, f);
421 : // allow refCallback to calculate a new position
422 0 : if(m_projector.valid())
423 0 : m_projector->new_vertex(newVrt, f);
424 : //GFDR_PROFILE_END();
425 : }
426 :
427 : // register the new faces and assign status
428 0 : for(size_t j = 0; j < vFaces.size(); ++j)
429 0 : mg.register_element(vFaces[j], f);
430 : }
431 : else{
432 : LOG(" WARNING in Refine: could not refine face.\n");
433 : }
434 : //GFDR_PROFILE_END();
435 : }
436 :
437 :
438 : UG_DLOG(LIB_GRID, 1, " creating new volumes\n");
439 :
440 : // only used for tetrahedron refinement
441 0 : vector<vector3> corners(4, vector3(0, 0, 0));
442 :
443 : // create new vertices and volumes from marked volumes
444 0 : for(VolumeIterator iter = mg.begin<Volume>(oldTopLevel);
445 0 : iter != mg.end<Volume>(oldTopLevel); ++iter)
446 : {
447 0 : if(!refinement_is_allowed(*iter) || !m_marker.is_marked(*iter))
448 0 : continue;
449 :
450 : Volume* v = *iter;
451 : //GFDR_PROFILE(GFDR_Refining_Volume);
452 :
453 : // collect child-vertices
454 : //GFDR_PROFILE(GFDR_CollectingVolumeVertices);
455 : vVrts.clear();
456 0 : for(uint j = 0; j < v->num_vertices(); ++j)
457 0 : vVrts.push_back(mg.get_child_vertex(v->vertex(j)));
458 : //GFDR_PROFILE_END();
459 :
460 : // collect the associated edges
461 : vEdgeVrts.clear();
462 : //GFDR_PROFILE(GFDR_CollectingVolumeEdgeVertices);
463 : //bool bIrregular = false;
464 0 : for(uint j = 0; j < v->num_edges(); ++j)
465 0 : vEdgeVrts.push_back(mg.get_child_vertex(mg.get_edge(v, j)));
466 : //GFDR_PROFILE_END();
467 :
468 : // collect associated face-vertices
469 : vFaceVrts.clear();
470 : //GFDR_PROFILE(GFDR_CollectingVolumeFaceVertices);
471 0 : for(uint j = 0; j < v->num_faces(); ++j)
472 0 : vFaceVrts.push_back(mg.get_child_vertex(mg.get_face(v, j)));
473 : //GFDR_PROFILE_END();
474 :
475 : // if we're performing tetrahedral refinement, we have to collect
476 : // the corner coordinates, so that the refinement algorithm may choose
477 : // the best interior diagonal.
478 : vector3* pCorners = NULL;
479 0 : if((v->num_vertices() == 4) && m_projector.valid()){
480 0 : for(size_t i = 0; i < 4; ++i){
481 0 : corners[i] = m_projector->geometry()->pos(v->vertex(i));
482 : }
483 : pCorners = &corners.front();
484 : }
485 :
486 : Vertex* newVrt;
487 0 : if(v->refine(vVols, &newVrt, &vEdgeVrts.front(), &vFaceVrts.front(),
488 0 : NULL, RegularVertex(), &vVrts.front(), pCorners)){
489 : // if a new vertex was generated, we have to register it
490 0 : if(newVrt){
491 0 : mg.register_element(newVrt, v);
492 : // allow refCallback to calculate a new position
493 0 : if(m_projector.valid())
494 0 : m_projector->new_vertex(newVrt, v);
495 : }
496 :
497 : // register the new faces and assign status
498 0 : for(size_t j = 0; j < vVols.size(); ++j)
499 0 : mg.register_element(vVols[j], v);
500 : }
501 : else{
502 : LOG(" WARNING in Refine: could not refine volume.\n");
503 : }
504 : //GFDR_PROFILE_END();
505 : }
506 :
507 : // done - clean up
508 0 : if(!bHierarchicalInsertionWasEnabled)
509 0 : mg.enable_hierarchical_insertion(false);
510 :
511 : // notify derivates that refinement ends
512 0 : refinement_step_ends();
513 :
514 0 : projector()->refinement_ends();
515 0 : m_messageHub->post_message(GridMessage_Adaption(GMAT_GLOBAL_REFINEMENT_ENDS,
516 0 : mg.get_grid_objects(oldTopLevel)));
517 :
518 : UG_DLOG(LIB_GRID, 1, " refinement done.");
519 0 : }
520 :
521 : /*
522 : //template <class TAPosition>
523 : template <class TElem>
524 : //void GlobalFracturedMediaRefiner<TAPosition>::
525 : void GlobalFracturedMediaRefiner::
526 : assign_elem_and_side_marks()
527 : {
528 : typedef typename MultiGrid::traits<TElem>::iterator ElemIter;
529 : typedef typename TElem::side Side;
530 :
531 : if(!m_pMG)
532 : UG_THROW("No grid assigned!");
533 :
534 : if(!m_pSH)
535 : UG_THROW("No subset handler assigned");
536 :
537 : if(m_pMG->num_levels() == 0)
538 : UG_THROW("At least one level has to exist in the associated multi-grid.");
539 :
540 : // we might need an attachment accessor for this one.
541 : // Grid::VertexAttachmentAccessor<TAPosition> aaPos;
542 : // if(!aaPos.access(*m_pMG, m_aPos))
543 : // UG_THROW("You have to specify a valid position attachment through 'set_position_attachment'");
544 :
545 : MultiGrid& mg = *m_pMG;
546 :
547 : // iterate over all elements in the top level and adjust their marks
548 : int topLvl = mg.num_levels() - 1;
549 :
550 : // we only clear side marks and element marks here, since only those will be set again.
551 : m_marker.clear();
552 :
553 : // collect associated sides and elements in this vectors
554 : std::vector<Side*> sides, sidesCon;
555 : std::vector<TElem*> elems;
556 : std::queue<TElem*> unclassified;
557 :
558 : for(ElemIter iter = mg.begin<TElem>(topLvl);
559 : iter != mg.end<TElem>(topLvl); ++iter)
560 : {
561 : TElem* e = *iter;
562 :
563 : m_marker.mark(e);
564 :
565 : // if the element is not a fracture element, we'll mark all sides for refinement
566 : // then we're done
567 : if(!is_fracture_element(e)){
568 : CollectAssociated(sides, mg, e);
569 : for(size_t i = 0; i < sides.size(); ++i)
570 : m_marker.mark(sides[i]);
571 : continue;
572 : }
573 :
574 :
575 : // find the side which is connected to a non-fracture element
576 : CollectAssociated(sides, mg, e);
577 : Side* outerSide = NULL;
578 : for(size_t i_side = 0; i_side < sides.size(); ++i_side){
579 : Side* s = sides[i_side];
580 : CollectAssociated(elems, mg, s);
581 : bool gotOne = false;
582 : for(size_t i_elem = 0; i_elem < elems.size(); ++i_elem){
583 : if(!is_fracture_element(elems[i_elem])){
584 : gotOne = true;
585 : break;
586 : }
587 : }
588 :
589 : if(gotOne){
590 : outerSide = s;
591 : break;
592 : }
593 : }
594 :
595 : if(outerSide){
596 : // the side is an interface between the fracture and the surrounding medium.
597 : m_marker.mark(outerSide);
598 :
599 : // we now have to also mark the opposing side in e. Get the descriptor
600 : // of that side
601 : typename geometry_traits<Side>::GeneralDescriptor desc;
602 : if(e->get_opposing_side(outerSide, desc)){
603 : Side* opSide = mg.get_element(desc);
604 : if(opSide)
605 : m_marker.mark(opSide);
606 : }
607 : else{
608 : // this should only happen if we're on a cap-element
609 : // we'll push the element to the list of unclassified elements
610 : unclassified.push(e);
611 : }
612 : }
613 : else{
614 : UG_THROW("Invalid fracture structure. Each fracture element has to"
615 : " be connected to non-fracture element.");
616 : }
617 : }
618 :
619 : // we'll count the number of failed tries in a row, to avoid an infinite loop
620 : size_t numFailedTries = 0;
621 : while(!unclassified.empty()){
622 : TElem* e = unclassified.front();
623 : unclassified.pop();
624 :
625 : // first check if two sides are already marked (e.g. during an earlier iteration)
626 : // if this is the case, we'll ignore the element
627 : CollectAssociated(sides, mg, e);
628 : size_t numMarked = num_marked(sides);
629 :
630 : if(numMarked == 2){
631 : // the additional side has been marked during iteration over an other element
632 : continue;
633 : }
634 : else if(numMarked != 1){
635 : UG_THROW("Unknown fracture topology encountered. "
636 : << numMarked << " sides marked, but only 1 should be marked."
637 : " In element at " << GetGridObjectCenter(mg, e));
638 : }
639 :
640 : // iterate over neighbors which lie in the fracture. If exactly one
641 : // neighbor has only one marked side (!=markedSide), then the side
642 : // connecting both will be marked.
643 : // If no neighbor is found with exactly one marked side, then an error is raised.
644 : // If more then one neighbors are found with exactly one marked side, then
645 : // the element is pushed back to queue.
646 : Side* connectingSide = NULL;
647 : int elemsFound = 0;// set to true, if an element with exactly one marked side is found.
648 : for(size_t i_side = 0; i_side < sides.size(); ++i_side){
649 : Side* s = sides[i_side];
650 : if(m_marker.is_marked(s))
651 : continue;
652 :
653 : // get the connected element of the side
654 : CollectAssociated(elems, mg, s);
655 : if(elems.size() != 2){
656 : UG_THROW("Invalid fracture topology detected. Aborting.");
657 : }
658 :
659 : TElem* eCon;
660 : if(elems[0] == e) eCon = elems[1];
661 : else eCon = elems[0];
662 :
663 : // count the number of marked edges in eCon
664 : CollectAssociated(sidesCon, mg, eCon);
665 :
666 : if(num_marked(sidesCon) == 1){
667 : ++elemsFound;
668 : connectingSide = s;
669 : }
670 : }
671 :
672 : if(elemsFound == 0){
673 : UG_THROW("Couldn't decide which inner edge to refine."
674 : " Please check your fractures topology at "
675 : << GetGridObjectCenter(mg, e));
676 : }
677 : else if(elemsFound == 1){
678 : // we found the side which we have to mark.
679 : m_marker.mark(connectingSide);
680 : numFailedTries = 0;
681 : }
682 : else{
683 : // we have to push the element back to the queue
684 : unclassified.push(e);
685 : ++numFailedTries;
686 : // prevent an infinite loop
687 : if(numFailedTries >= unclassified.size()){
688 : UG_THROW("Couldn't identify fracture structure.");
689 : }
690 : }
691 : }
692 :
693 : // mark sides of marked side elements
694 : mark_sides_of_marked_top_level_elements<Side>();
695 : }
696 : */
697 :
698 :
699 : template <class TElem>
700 0 : void GlobalFracturedMediaRefiner::
701 : assign_elem_and_side_marks()
702 : {
703 : typedef typename MultiGrid::traits<TElem>::iterator ElemIter;
704 : typedef typename TElem::side Side;
705 : typedef typename Side::side SideOfSide;
706 :
707 0 : if(!m_pMG)
708 0 : UG_THROW("No grid assigned!");
709 :
710 0 : if(!m_pSH)
711 0 : UG_THROW("No subset handler assigned");
712 :
713 0 : if(m_pMG->num_levels() == 0)
714 0 : UG_THROW("At least one level has to exist in the associated multi-grid.");
715 :
716 : MultiGrid& mg = *m_pMG;
717 :
718 : // iterate over all elements in the top level and adjust their marks
719 0 : int topLvl = mg.num_levels() - 1;
720 :
721 : // we only clear side marks and element marks here, since only those will be set again.
722 0 : m_marker.clear();
723 :
724 : // we'll also use a temporary marked to mark fixed sides.
725 0 : BoolMarker fixedMarker(mg);
726 :
727 : // collect associated sides and elements in this vectors
728 :
729 : std::vector<Side*> sides, sidesCon;
730 : std::vector<SideOfSide*> sidesOfSides; //only used in 3d for fixed marks
731 : std::vector<TElem*> elems;
732 : std::vector<TElem*> caps;// elements at fracture caps (rim-boundary)
733 :
734 : // first we'll mark all elements which do not lie in a fracture. We'll also mark
735 : // their sides.
736 0 : for(ElemIter iter = mg.begin<TElem>(topLvl);
737 0 : iter != mg.end<TElem>(topLvl); ++iter)
738 : {
739 : TElem* e = *iter;
740 :
741 0 : if(!refinement_is_allowed(e))
742 0 : continue;
743 :
744 0 : if(!is_fracture_element(e)){
745 : m_marker.mark(e);
746 : CollectAssociated(sides, mg, e);
747 0 : for(size_t i = 0; i < sides.size(); ++i){
748 0 : m_marker.mark(sides[i]);
749 : }
750 : }
751 : }
752 :
753 :
754 : // call a callback, which allows to communicate marks in a parallel environment.
755 : // the default implementation does nothing
756 0 : communicate_marks(m_marker);
757 :
758 : // now iterate over all fracture elements and mark interior sides accordingly
759 0 : for(ElemIter iter = mg.begin<TElem>(topLvl);
760 0 : iter != mg.end<TElem>(topLvl); ++iter)
761 : {
762 0 : TElem* e = *iter;
763 :
764 0 : if(!refinement_is_allowed(e))
765 0 : continue;
766 :
767 0 : if(!is_fracture_element(e))
768 0 : continue;
769 :
770 : m_marker.mark(e);
771 :
772 : // find the side which is already marked. If more than one is marked,
773 : // we'll ignore the element. If no side is marked, we'll throw an error,
774 : // since the fracture does not feature the required topology in this case.
775 : CollectAssociated(sides, mg, e);
776 : Side* markedSide = NULL;
777 : int numMarked = 0;
778 0 : for(size_t i_side = 0; i_side < sides.size(); ++i_side){
779 0 : Side* s = sides[i_side];
780 0 : if(m_marker.is_marked(s)){
781 : markedSide = s;
782 0 : ++numMarked;
783 : }
784 : }
785 :
786 0 : if((numMarked == 0) || (numMarked > 2)){
787 0 : UG_THROW("Bad fracture topology encountered in element with center "
788 : << GetGridObjectCenter(mg, e) << "! Aborting.");
789 : }
790 :
791 0 : typename geometry_traits<Side>::GeneralDescriptor desc;
792 0 : if(e->get_opposing_side(markedSide, desc)){
793 : Side* opSide = mg.get_element(desc);
794 0 : if(opSide){
795 0 : if((numMarked == 2) && (!m_marker.is_marked(opSide))){
796 0 : UG_THROW("Bad fracture topology encountered in element with center "
797 : << GetGridObjectCenter(mg, e) << "! Aborting.");
798 : }
799 : m_marker.mark(opSide);
800 : }
801 : else{
802 0 : UG_THROW("Opposing side could not be found in grid. Check grid options.");
803 : }
804 :
805 : // we'll iterate over the sides of the element again and mark currently
806 : // unmarked sides as fixed sides
807 0 : for(size_t i_side = 0; i_side < sides.size(); ++i_side){
808 0 : Side* s = sides[i_side];
809 0 : if(!m_marker.is_marked(s)){
810 : // in 3d we'll also have to mark sides of sides, to guarantee,
811 : // that fixed marks are communicated correctly.
812 : fixedMarker.mark(s);
813 : }
814 :
815 : if(Side::dim == 2){
816 : CollectAssociated(sidesOfSides, mg, s);
817 0 : for(size_t i = 0; i < sidesOfSides.size(); ++i)
818 0 : fixedMarker.mark(sidesOfSides[i]);
819 : }
820 : }
821 : }
822 : else{
823 : // this should only happen if we're on a cap-element. Store the element for later use.
824 0 : caps.push_back(e);
825 : }
826 : }
827 :
828 : // communicate fixed marks. Default implementation does nothing.
829 : // This can be used by a derived class to communicate marks to other processes.
830 0 : communicate_marks(fixedMarker);
831 :
832 : // now adjust marks at cap elements
833 0 : for(typename std::vector<TElem*>::iterator iter = caps.begin();
834 0 : iter != caps.end(); ++iter)
835 : {
836 0 : TElem* e = *iter;
837 : m_marker.mark(e);
838 :
839 : // sides which are marked as fixed or have a fixed side them selfes (3d only)
840 : // may not be refined. All others may be refined.
841 : CollectAssociated(sides, mg, e);
842 0 : for(size_t i_side = 0; i_side < sides.size(); ++i_side){
843 : bool refineSide = true;
844 0 : Side* s = sides[i_side];
845 :
846 0 : if(fixedMarker.is_marked(s)){
847 : refineSide = false;
848 : break;
849 : }
850 : else{
851 : if(Side::dim == 2){
852 : CollectAssociated(sidesOfSides, mg, s);
853 0 : for(size_t i = 0; i < sidesOfSides.size(); ++i){
854 0 : if(fixedMarker.is_marked(sidesOfSides[i])){
855 : refineSide = false;
856 : break;
857 : }
858 : }
859 :
860 0 : if(!refineSide)
861 : break;
862 : }
863 : }
864 :
865 : // if the side shall be refined, we'll mark it
866 : if(refineSide)
867 : m_marker.mark(s);
868 : }
869 :
870 : }
871 :
872 : // mark sides of marked side elements
873 : if(Side::HAS_SIDES)
874 0 : mark_sides_of_marked_top_level_elements<Side>();
875 0 : }
876 :
877 :
878 : template <class TElem>
879 0 : void GlobalFracturedMediaRefiner::
880 : mark_sides_of_marked_top_level_elements()
881 : {
882 : typedef typename MultiGrid::traits<TElem>::iterator ElemIter;
883 : typedef typename TElem::side Side;
884 0 : if(!m_pMG)
885 0 : UG_THROW("No grid assigned!");
886 :
887 0 : if(m_pMG->num_levels() == 0)
888 0 : UG_THROW("At least one level has to exist in the associated multi-grid.");
889 :
890 : MultiGrid& mg = *m_pMG;
891 :
892 : // iterate over all elements in the top level and adjust their marks
893 0 : int topLvl = mg.num_levels() - 1;
894 :
895 : // collect sides in this container
896 : vector<Side*> sides;
897 :
898 0 : for(ElemIter iter = mg.begin<TElem>(topLvl);
899 0 : iter != mg.end<TElem>(topLvl); ++iter)
900 : {
901 : TElem* e = *iter;
902 0 : if(m_marker.is_marked(e)){
903 0 : CollectAssociated(sides, mg, e);
904 0 : for(size_t i = 0; i < sides.size(); ++i){
905 0 : m_marker.mark(sides[i]);
906 : }
907 : }
908 : }
909 :
910 : if(Side::HAS_SIDES)
911 0 : mark_sides_of_marked_top_level_elements<Side>();
912 0 : }
913 :
914 :
915 : //template <class TAPosition>
916 : //void GlobalFracturedMediaRefiner<TAPosition>::
917 0 : void GlobalFracturedMediaRefiner::
918 : adjust_marks()
919 : {
920 0 : if(!m_pMG){
921 0 : UG_THROW("No grid assigned!");
922 : }
923 :
924 : // call the actual assign method
925 0 : if(m_pMG->num<Volume>() > 0)
926 0 : assign_elem_and_side_marks<Volume>();
927 0 : else if(m_pMG->num<Face>() > 0)
928 0 : assign_elem_and_side_marks<Face>();
929 : else{
930 : // simply mark everything
931 0 : m_marker.mark(m_pMG->vertices_begin(), m_pMG->vertices_end());
932 0 : m_marker.mark(m_pMG->edges_begin(), m_pMG->edges_end());
933 : }
934 0 : }
935 :
936 : //template <class TAPosition>
937 : //bool GlobalFracturedMediaRefiner<TAPosition>::
938 0 : bool GlobalFracturedMediaRefiner::
939 : save_marks_to_file(const char* filename)
940 : {
941 0 : if(!m_pMG){
942 0 : UG_THROW("ERROR in GlobalFracturedMediaRefiner::save_marks_to_file: No grid assigned!");
943 : }
944 :
945 : MultiGrid& mg = *m_pMG;
946 0 : SubsetHandler sh(mg);
947 :
948 0 : AssignGridToSubset(mg, sh, 2);
949 0 : int lvl = mg.num_levels() - 1;
950 :
951 0 : for(VolumeIterator iter = mg.begin<Volume>(lvl); iter != mg.end<Volume>(lvl); ++iter){
952 0 : if(m_marker.is_marked(*iter)){
953 0 : if(is_fracture_element(*iter))
954 0 : sh.assign_subset(*iter, 1);
955 : else
956 0 : sh.assign_subset(*iter, 0);
957 : }
958 : }
959 :
960 : vector<Edge*> edges;
961 0 : for(FaceIterator iter = mg.begin<Face>(lvl); iter != mg.end<Face>(lvl); ++iter){
962 : Face* f = *iter;
963 : CollectAssociated(edges, mg, f);
964 0 : if(m_marker.is_marked(f)){
965 : //if(num_marked(edges) != f->num_sides())
966 : //sh.assign_subset(f, 1);
967 : //else
968 0 : sh.assign_subset(f, 0);
969 0 : }else if(num_marked(edges) != f->num_sides()){
970 0 : sh.assign_subset(f, 1);
971 : }
972 : }
973 :
974 0 : for(EdgeIterator iter = mg.begin<Edge>(lvl); iter != mg.end<Edge>(lvl); ++iter){
975 0 : if(m_marker.is_marked(*iter))
976 0 : sh.assign_subset(*iter, 0);
977 : }
978 :
979 0 : for(VertexIterator iter = mg.begin<Vertex>(lvl); iter != mg.end<Vertex>(lvl); ++iter){
980 0 : if(m_marker.is_marked(*iter))
981 0 : sh.assign_subset(*iter, 0);
982 : }
983 :
984 0 : sh.subset_info(0).name = "refine regular";
985 0 : sh.subset_info(1).name = "refine anisotropic";
986 0 : sh.subset_info(2).name = "no marks";
987 :
988 0 : AssignSubsetColors(sh);
989 :
990 0 : return SaveGridToFile(mg, sh, filename);
991 0 : }
992 :
993 : template <class TElem>
994 0 : size_t GlobalFracturedMediaRefiner::
995 : num_marked(const std::vector<TElem*>& elems) const
996 : {
997 : size_t num = 0;
998 0 : for(size_t i = 0; i < elems.size(); ++i){
999 0 : if(m_marker.is_marked(elems[i]))
1000 0 : ++num;
1001 : }
1002 0 : return num;
1003 : }
1004 :
1005 : }// end of namespace
|
