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 <sstream>
34 : #include "common/assert.h"
35 : #include "lib_grid/algorithms/multi_grid_util.h"
36 : #include "lib_grid/algorithms/selection_util.h"
37 : #include "lib_grid/algorithms/debug_util.h"
38 : #include "hanging_node_refiner_multi_grid.h"
39 : #include "ref_mark_adjusters/mg_hnode_adjuster.h"
40 :
41 : // Only required for debug saves
42 : #include "lib_grid/file_io/file_io.h"
43 :
44 : #ifdef UG_PARALLEL
45 : #include "pcl/pcl.h"
46 : #include "lib_grid/parallelization/distributed_grid.h"
47 : #include "lib_grid/parallelization/parallelization_util.h"
48 : #endif
49 :
50 : //define PROFILE_HANGING_NODE_REFINER_MG if you want to profile
51 : //the refinement code.
52 : //#define PROFILE_HANGING_NODE_REFINER_MG
53 : #ifdef PROFILE_HANGING_NODE_REFINER_MG
54 : #define MGHNODE_PROFILE_FUNC() PROFILE_FUNC()
55 : #define MGHNODE_PROFILE_BEGIN(name) PROFILE_BEGIN(name)
56 : #define MGHNODE_PROFILE_END() PROFILE_END()
57 : #else
58 : #define MGHNODE_PROFILE_FUNC()
59 : #define MGHNODE_PROFILE_BEGIN(name)
60 : #define MGHNODE_PROFILE_END()
61 : #endif
62 :
63 : using namespace std;
64 :
65 : namespace ug{
66 :
67 0 : HangingNodeRefiner_MultiGrid::
68 0 : HangingNodeRefiner_MultiGrid(SPRefinementProjector projector) :
69 : BaseClass(projector),
70 0 : m_pMG(NULL)
71 : {
72 0 : add_ref_mark_adjuster(MGHNodeAdjuster::create());
73 0 : }
74 :
75 0 : HangingNodeRefiner_MultiGrid::
76 : HangingNodeRefiner_MultiGrid(MultiGrid& mg,
77 0 : SPRefinementProjector projector) :
78 : BaseClass(projector),
79 0 : m_pMG(NULL)
80 : {
81 0 : add_ref_mark_adjuster(MGHNodeAdjuster::create());
82 0 : set_grid(&mg);
83 0 : }
84 :
85 0 : HangingNodeRefiner_MultiGrid::
86 0 : ~HangingNodeRefiner_MultiGrid()
87 : {
88 0 : set_grid(NULL);
89 0 : }
90 :
91 0 : void HangingNodeRefiner_MultiGrid::
92 : grid_to_be_destroyed(Grid* grid)
93 : {
94 0 : set_grid(NULL);
95 0 : }
96 :
97 0 : void HangingNodeRefiner_MultiGrid::
98 : assign_grid(MultiGrid& mg)
99 : {
100 0 : set_grid(&mg);
101 0 : }
102 :
103 0 : void HangingNodeRefiner_MultiGrid::
104 : num_marked_edges_local(std::vector<int>& numMarkedEdgesOut)
105 : {
106 0 : num_marked_elems<Edge>(numMarkedEdgesOut);
107 0 : }
108 :
109 0 : void HangingNodeRefiner_MultiGrid::
110 : num_marked_faces_local(std::vector<int>& numMarkedFacesOut)
111 : {
112 0 : num_marked_elems<Face>(numMarkedFacesOut);
113 0 : }
114 :
115 0 : void HangingNodeRefiner_MultiGrid::
116 : num_marked_volumes_local(std::vector<int>& numMarkedVolsOut)
117 : {
118 0 : num_marked_elems<Volume>(numMarkedVolsOut);
119 0 : }
120 :
121 :
122 : template <class TElem>
123 0 : void HangingNodeRefiner_MultiGrid::
124 : num_marked_elems(std::vector<int>& numMarkedElemsOut)
125 : {
126 : MGSelector& sel = get_refmark_selector();
127 : numMarkedElemsOut.clear();
128 0 : numMarkedElemsOut.resize(sel.num_levels(), 0);
129 0 : for(size_t i = 0; i < sel.num_levels(); ++i){
130 0 : numMarkedElemsOut[i] = sel.num<TElem>(i);
131 : }
132 0 : }
133 :
134 0 : void HangingNodeRefiner_MultiGrid::
135 : set_grid(MultiGrid* mg)
136 : {
137 : // call base class implementation
138 0 : BaseClass::set_grid(mg);
139 :
140 0 : m_pMG = mg;
141 0 : }
142 :
143 : bool
144 0 : HangingNodeRefiner_MultiGrid::
145 : refinement_is_allowed(Vertex* elem)
146 : {
147 0 : return !m_pMG->has_children(elem);
148 : }
149 :
150 : bool
151 0 : HangingNodeRefiner_MultiGrid::
152 : refinement_is_allowed(Edge* elem)
153 : {
154 0 : return (!m_pMG->has_children(elem))
155 0 : || elem->is_constraining();
156 : }
157 :
158 : bool
159 0 : HangingNodeRefiner_MultiGrid::
160 : refinement_is_allowed(Face* elem)
161 : {
162 0 : return (!m_pMG->has_children(elem))
163 0 : || elem->is_constraining();
164 : }
165 :
166 : bool
167 0 : HangingNodeRefiner_MultiGrid::
168 : refinement_is_allowed(Volume* elem)
169 : {
170 0 : return !m_pMG->has_children(elem);
171 : }
172 :
173 0 : void HangingNodeRefiner_MultiGrid::assign_hnode_marks()
174 : {
175 : // first call the base implementation
176 0 : BaseClass::assign_hnode_marks();
177 :
178 : // now use multigrid info to do another important set of hnode assignments
179 0 : if (!m_pMG)
180 0 : throw(UGError("HangingNodeRefiner_MultiGrid::refine: No grid assigned."));
181 :
182 : MultiGrid& mg = *m_pMG;
183 :
184 : // TODO: This had better be done in the base implementation. But how?
185 : // Also take care of (unmarkable) SHADOW_COPY edges
186 : // whose children are marked for refinement:
187 : // These children need to be marked HNRM_TO_CONSTRAINING!
188 : Grid::face_traits::secure_container fl;
189 : EdgeIterator it = mg.begin<Edge>();
190 : EdgeIterator itEnd = mg.end<Edge>();
191 0 : for (; it != itEnd; ++it)
192 : {
193 : Edge* e = *it;
194 :
195 : // find out whether this is a SHADOW-COPY edge;
196 : // which is the case iff it has exactly one child
197 : // and there exists an associated face that does not have any
198 0 : if (mg.num_child_edges(e) != 1)
199 0 : continue;
200 :
201 : bool faceWithoutChildExists = false;
202 0 : mg.associated_elements(fl, e);
203 : const size_t flSz = fl.size();
204 0 : for (size_t f = 0; f < flSz; ++f)
205 : {
206 0 : if (refinement_is_allowed(fl[f]) && !mg.num_child_faces(fl[f]))
207 : {
208 : faceWithoutChildExists = true;
209 : break;
210 : }
211 : }
212 :
213 0 : if (!faceWithoutChildExists)
214 0 : continue;
215 :
216 : Edge* child = mg.get_child_edge(e, 0);
217 :
218 0 : if (marked_refine(child))
219 0 : add_hmark(child, HNRM_TO_CONSTRAINING);
220 : }
221 0 : }
222 :
223 0 : void HangingNodeRefiner_MultiGrid::
224 : pre_refine()
225 : {
226 0 : if(!m_pMG)
227 0 : throw(UGError("HangingNodeRefiner_MultiGrid::refine: No grid assigned."));
228 :
229 : MultiGrid& mg = *m_pMG;
230 :
231 : // Resize the attachment containers
232 :
233 : {
234 : selector_t& sel = get_refmark_selector();
235 :
236 : MGHNODE_PROFILE_BEGIN(MGHNode_ReserveAttachmentMemory);
237 :
238 : MGHNODE_PROFILE_BEGIN(MGHNODE_ReserveVrtData);
239 0 : mg.reserve<Vertex>(mg.num<Vertex>() +
240 : + sel.num<Vertex>()
241 0 : + sel.num<Edge>() - sel.num<ConstrainingEdge>()
242 0 : + sel.num<Quadrilateral>()
243 0 : + sel.num<ConstrainedQuadrilateral>()
244 0 : + sel.num<Hexahedron>());
245 : MGHNODE_PROFILE_END();
246 :
247 : MGHNODE_PROFILE_BEGIN(MGHNODE_ReserveEdgeData);
248 0 : mg.reserve<Edge>(mg.num<Edge>()
249 0 : + 2 * (sel.num<Edge>() - sel.num<ConstrainingEdge>())
250 0 : + 3 * (sel.num<Triangle>() + sel.num<ConstrainedTriangle>())
251 0 : + 4 * (sel.num<Quadrilateral>() + sel.num<ConstrainedQuadrilateral>())
252 0 : + 3 * sel.num<Prism>() + sel.num<Tetrahedron>()
253 0 : + 4 * sel.num<Pyramid>() + 6 * sel.num<Hexahedron>());
254 : MGHNODE_PROFILE_END();
255 :
256 : MGHNODE_PROFILE_BEGIN(MGHNODE_ReserveFaceData);
257 0 : mg.reserve<Face>(mg.num<Face>()
258 0 : + 4 * (sel.num<Face>() - sel.num<ConstrainingTriangle>()
259 0 : - sel.num<ConstrainingQuadrilateral>())
260 0 : + 10 * sel.num<Prism>()
261 0 : + 8 * sel.num<Tetrahedron>()
262 0 : + 9 * sel.num<Pyramid>() + 12 * sel.num<Hexahedron>());
263 : MGHNODE_PROFILE_END();
264 :
265 : MGHNODE_PROFILE_BEGIN(MGHNODE_ReserveVolData);
266 0 : mg.reserve<Volume>(mg.num<Volume>()
267 0 : + 8 * sel.num<Tetrahedron>() + 8 * sel.num<Prism>()
268 0 : + 6 * sel.num<Pyramid>() + 8 * sel.num<Hexahedron>());
269 : MGHNODE_PROFILE_END();
270 :
271 : MGHNODE_PROFILE_END();
272 : }
273 :
274 :
275 : // create a child vertex for each marked vertex
276 0 : for(sel_vrt_iter iter = m_selMarkedElements.begin<Vertex>();
277 0 : iter != m_selMarkedElements.end<Vertex>(); ++iter)
278 : {
279 0 : if(marked_refine(*iter) && refinement_is_allowed(*iter)){
280 0 : Vertex* vrt = *mg.create<RegularVertex>(*iter);
281 0 : if(m_projector.valid())
282 0 : m_projector->new_vertex(vrt, *iter);
283 : }
284 : }
285 0 : }
286 :
287 0 : void HangingNodeRefiner_MultiGrid::
288 : process_constraining_edge(ConstrainingEdge* cge)
289 : {
290 : // call base implementation to perform refinement
291 0 : BaseClass::process_constraining_edge(cge);
292 :
293 : // the constrained edge is now a normal edge
294 : UG_ASSERT(marked_to_normal(cge), "A constraining has to be converted to a"
295 : " normal edge when it is refined.");
296 0 : Vertex* centerVrt = get_center_vertex(cge);
297 0 : RegularEdge* e = *m_pMG->create_and_replace<RegularEdge>(cge);
298 0 : if(centerVrt)
299 0 : set_center_vertex(e, centerVrt);
300 0 : }
301 :
302 0 : void HangingNodeRefiner_MultiGrid::
303 : refine_edge_with_normal_vertex(Edge* e, Vertex** newCornerVrts)
304 : {
305 : // collect child corners
306 : std::vector<Vertex*> childCorners;
307 0 : collect_child_corners(childCorners, e);
308 :
309 : // call original implementation
310 0 : BaseClass::refine_edge_with_normal_vertex(e, &childCorners.front());
311 0 : }
312 :
313 0 : void HangingNodeRefiner_MultiGrid::
314 : refine_edge_with_hanging_vertex(Edge* e, Vertex** newCornerVrts)
315 : {
316 : // collect child corners
317 : std::vector<Vertex*> childCorners;
318 0 : collect_child_corners(childCorners, e);
319 :
320 : // call original implementation
321 0 : BaseClass::refine_edge_with_hanging_vertex(e, &childCorners.front());
322 0 : }
323 :
324 0 : void HangingNodeRefiner_MultiGrid::
325 : refine_face_with_normal_vertex(Face* f, Vertex** newCornerVrts)
326 : {
327 : // collect child corners
328 : std::vector<Vertex*> childCorners;
329 0 : collect_child_corners(childCorners, f);
330 :
331 : // call original implementation
332 0 : BaseClass::refine_face_with_normal_vertex(f, &childCorners.front());
333 0 : }
334 :
335 0 : void HangingNodeRefiner_MultiGrid::
336 : refine_face_with_hanging_vertex(Face* f, Vertex** newCornerVrts)
337 : {
338 : // collect child corners
339 : std::vector<Vertex*> childCorners;
340 0 : collect_child_corners(childCorners, f);
341 :
342 : // call original implementation
343 0 : BaseClass::refine_face_with_hanging_vertex(f, &childCorners.front());
344 0 : }
345 :
346 0 : void HangingNodeRefiner_MultiGrid::
347 : refine_volume_with_normal_vertex(Volume* v, Vertex** newCornerVrts)
348 : {
349 : // collect child corners
350 : std::vector<Vertex*> childCorners;
351 0 : collect_child_corners(childCorners, v);
352 :
353 : // call original implementation
354 0 : BaseClass::refine_volume_with_normal_vertex(v, &childCorners.front());
355 0 : }
356 :
357 0 : Vertex* HangingNodeRefiner_MultiGrid::
358 : get_center_vertex(Edge* e)
359 : {
360 0 : return m_pMG->get_child_vertex(e);
361 : }
362 :
363 :
364 0 : void HangingNodeRefiner_MultiGrid::
365 : set_center_vertex(Edge* e, Vertex* v)
366 : {
367 : // the child vertex should already have been associated since the
368 : // multigrid is an observer itself.
369 : UG_ASSERT(m_pMG->get_child_vertex(e) == v, "child/vertex mismatch.");
370 0 : }
371 :
372 :
373 0 : Vertex* HangingNodeRefiner_MultiGrid::
374 : get_center_vertex(Face* f)
375 : {
376 0 : return m_pMG->get_child_vertex(f);
377 : }
378 :
379 :
380 0 : void HangingNodeRefiner_MultiGrid::
381 : set_center_vertex(Face* f, Vertex* v)
382 : {
383 : // the child vertex should already have been associated since the
384 : // multigrid is an observer itself.
385 : UG_ASSERT(m_pMG->get_child_vertex(f) == v, "child/vertex mismatch.");
386 0 : }
387 :
388 :
389 :
390 0 : void HangingNodeRefiner_MultiGrid::
391 : restrict_selection_to_surface_coarsen_elements()
392 : {
393 0 : restrict_selection_to_surface_coarsen_elements<Vertex>();
394 0 : restrict_selection_to_surface_coarsen_elements<Edge>();
395 0 : restrict_selection_to_surface_coarsen_elements<Face>();
396 0 : restrict_selection_to_surface_coarsen_elements<Volume>();
397 0 : }
398 :
399 : template <class TElem>
400 0 : void HangingNodeRefiner_MultiGrid::
401 : restrict_selection_to_surface_coarsen_elements()
402 : {
403 0 : MultiGrid& mg = *m_pMG;
404 : selector_t& sel = get_refmark_selector();
405 :
406 0 : for(typename selector_t::template traits<TElem>::iterator iter = sel.template begin<TElem>();
407 0 : iter != sel.template end<TElem>();)
408 : {
409 : TElem* e = *iter;
410 : ++iter;
411 :
412 : // make sure that only coarsen-marks are applied
413 0 : if(get_mark(e) != RM_COARSEN){
414 0 : sel.deselect(e);
415 0 : continue;
416 : }
417 :
418 : // make sure that the element is a surface element
419 0 : if(mg.has_children(e))
420 0 : sel.deselect(e);
421 : }
422 0 : }
423 :
424 :
425 0 : void HangingNodeRefiner_MultiGrid::
426 : restrict_selection_to_coarsen_families()
427 : {
428 0 : restrict_selection_to_coarsen_families<Vertex>();
429 0 : restrict_selection_to_coarsen_families<Edge>();
430 0 : restrict_selection_to_coarsen_families<Face>();
431 0 : restrict_selection_to_coarsen_families<Volume>();
432 0 : }
433 :
434 :
435 : template <class TElem>
436 0 : void HangingNodeRefiner_MultiGrid::
437 : restrict_selection_to_coarsen_families()
438 : {
439 : typedef typename TElem::grid_base_object TBaseElem;
440 0 : MultiGrid& mg = *m_pMG;
441 : selector_t& sel = get_refmark_selector();
442 :
443 : // coarsening is only performed for complete families. If some siblings aren't
444 : // selected for coarsening, the whole family may not be coarsened. All siblings
445 : // have to be deselected in this case.
446 :
447 : // mark parents, which have been processed and mark siblings, which are valid
448 : // candidates. If thus a surface element is encountered, which has a marked
449 : // parent, but is not marked itself, then it is clear that it is not a valid
450 : // candidate.
451 0 : mg.begin_marking();
452 0 : for(typename selector_t::traits<TElem>::iterator iter = sel.begin<TElem>();
453 0 : iter != sel.end<TElem>();)
454 : {
455 : TElem* e = *iter;
456 : ++iter;
457 :
458 : // make sure that only surface elements are selected
459 : UG_ASSERT(mg.num_children<TBaseElem>(e) == 0,
460 : "Only surface elements may be passed to this method.");
461 :
462 : // make sure that only RM_COARSEN marks are used
463 : UG_ASSERT(get_mark(e) == RM_COARSEN,
464 : "Only RM_COARSEN marks may be used in this method.");
465 :
466 : // if the element is marked, we'll continue (the family is complete)
467 0 : if(mg.is_marked(e))
468 0 : continue;
469 :
470 : // get the parent
471 0 : TBaseElem* parent = dynamic_cast<TBaseElem*>(mg.get_parent(e));
472 0 : if(parent){
473 0 : if(mg.is_marked(parent)){
474 : // the parent is marked and e is not. We thus have to deselect e
475 0 : sel.deselect(e);
476 0 : continue;
477 : }
478 :
479 : mg.mark(parent);
480 :
481 : // check whether all children of e of type TBaseElem are marked
482 : bool allMarked = true;
483 : size_t numChildren = mg.num_children<TBaseElem>(parent);
484 0 : for(size_t i = 0; i < numChildren; ++i){
485 0 : if(get_mark(mg.get_child<TBaseElem>(parent, i)) != RM_COARSEN){
486 : allMarked = false;
487 : break;
488 : }
489 : }
490 :
491 0 : if(allMarked){
492 : // mark all children of parent, so that all it is clear that they
493 : // belong to a complete family
494 0 : for(size_t i = 0; i < numChildren; ++i){
495 : mg.mark(mg.get_child<TBaseElem>(parent, i));
496 : }
497 : }
498 : else
499 0 : sel.deselect(e);
500 : }
501 : else{
502 : // the parent of this element has a different type or it doesn't exist
503 : // at all. The element is thus not regarded as a valid candidate.
504 0 : sel.deselect(e);
505 : }
506 : }
507 0 : mg.end_marking();
508 0 : }
509 :
510 0 : static void DeselectFamily(ISelector& sel, MultiGrid& mg, Vertex* elem)
511 : {
512 0 : sel.deselect(elem);
513 : size_t numVertexChildren = mg.num_children<Vertex>(elem);
514 :
515 0 : for(size_t i = 0; i < numVertexChildren; ++i)
516 0 : DeselectFamily(sel, mg, mg.get_child<Vertex>(elem, i));
517 0 : }
518 :
519 0 : static void DeselectFamily(ISelector& sel, MultiGrid& mg, Edge* elem)
520 : {
521 0 : sel.deselect(elem);
522 : size_t numVertexChildren = mg.num_children<Vertex>(elem);
523 : size_t numEdgeChildren = mg.num_children<Edge>(elem);
524 :
525 0 : for(size_t i = 0; i < numVertexChildren; ++i)
526 0 : DeselectFamily(sel, mg, mg.get_child<Vertex>(elem, i));
527 :
528 0 : for(size_t i = 0; i < numEdgeChildren; ++i)
529 0 : DeselectFamily(sel, mg, mg.get_child<Edge>(elem, i));
530 0 : }
531 :
532 0 : static void DeselectFamily(ISelector& sel, MultiGrid& mg, Face* elem)
533 : {
534 0 : sel.deselect(elem);
535 : size_t numVertexChildren = mg.num_children<Vertex>(elem);
536 : size_t numEdgeChildren = mg.num_children<Edge>(elem);
537 : size_t numFaceChildren = mg.num_children<Face>(elem);
538 :
539 0 : for(size_t i = 0; i < numVertexChildren; ++i)
540 0 : DeselectFamily(sel, mg, mg.get_child<Vertex>(elem, i));
541 :
542 0 : for(size_t i = 0; i < numEdgeChildren; ++i)
543 0 : DeselectFamily(sel, mg, mg.get_child<Edge>(elem, i));
544 :
545 0 : for(size_t i = 0; i < numFaceChildren; ++i)
546 0 : DeselectFamily(sel, mg, mg.get_child<Face>(elem, i));
547 0 : }
548 :
549 0 : static void DeselectFamily(ISelector& sel, MultiGrid& mg, Volume* elem)
550 : {
551 0 : sel.deselect(elem);
552 : size_t numVertexChildren = mg.num_children<Vertex>(elem);
553 : size_t numEdgeChildren = mg.num_children<Edge>(elem);
554 : size_t numFaceChildren = mg.num_children<Face>(elem);
555 : size_t numVolumeChildren = mg.num_children<Volume>(elem);
556 :
557 0 : for(size_t i = 0; i < numVertexChildren; ++i)
558 0 : DeselectFamily(sel, mg, mg.get_child<Vertex>(elem, i));
559 :
560 0 : for(size_t i = 0; i < numEdgeChildren; ++i)
561 0 : DeselectFamily(sel, mg, mg.get_child<Edge>(elem, i));
562 :
563 0 : for(size_t i = 0; i < numFaceChildren; ++i)
564 0 : DeselectFamily(sel, mg, mg.get_child<Face>(elem, i));
565 :
566 0 : for(size_t i = 0; i < numVolumeChildren; ++i)
567 0 : DeselectFamily(sel, mg, mg.get_child<Volume>(elem, i));
568 0 : }
569 :
570 0 : static void DeselectFamily(ISelector& sel, MultiGrid& mg, GridObject* elem)
571 : {
572 0 : switch(elem->base_object_id()){
573 0 : case VERTEX: return DeselectFamily(sel, mg, static_cast<Vertex*>(elem));
574 0 : case EDGE: return DeselectFamily(sel, mg, static_cast<Edge*>(elem));
575 0 : case FACE: return DeselectFamily(sel, mg, static_cast<Face*>(elem));
576 0 : case VOLUME: return DeselectFamily(sel, mg, static_cast<Volume*>(elem));
577 : }
578 : }
579 :
580 :
581 0 : static void SaveCoarsenMarksToFile(ISelector& sel, const char* filename)
582 : {
583 : assert(sel.grid());
584 : Grid& g = *sel.grid();
585 :
586 0 : SubsetHandler sh(g);
587 :
588 0 : AssignGridToSubset(g, sh, 8);
589 0 : for(VertexIterator iter = g.vertices_begin(); iter != g.vertices_end(); ++iter){
590 : ISelector::status_t status = sel.get_selection_status(*iter);
591 0 : switch(status){
592 0 : case RM_COARSEN: sh.assign_subset(*iter, 0); break;
593 0 : case HangingNodeRefiner_MultiGrid::HNCM_ALL: sh.assign_subset(*iter, 1); break;
594 0 : case HangingNodeRefiner_MultiGrid::HNCM_PARTIAL: sh.assign_subset(*iter, 2); break;
595 0 : case HangingNodeRefiner_MultiGrid::HNCM_NONE: sh.assign_subset(*iter, 3); break;
596 0 : case HangingNodeRefiner_MultiGrid::HNCM_UNKNOWN: sh.assign_subset(*iter, 4); break;
597 0 : case HangingNodeRefiner_MultiGrid::HNCM_INVALID: sh.assign_subset(*iter, 5); break;
598 0 : case HangingNodeRefiner_MultiGrid::HNCM_REPLACE: sh.assign_subset(*iter, 6); break;
599 0 : case HangingNodeRefiner_MultiGrid::HNCM_NO_NBRS: sh.assign_subset(*iter, 7); break;
600 : }
601 : }
602 :
603 0 : for(EdgeIterator iter = g.edges_begin(); iter != g.edges_end(); ++iter){
604 : ISelector::status_t status = sel.get_selection_status(*iter);
605 0 : switch(status){
606 0 : case RM_COARSEN: sh.assign_subset(*iter, 0); break;
607 0 : case HangingNodeRefiner_MultiGrid::HNCM_ALL: sh.assign_subset(*iter, 1); break;
608 0 : case HangingNodeRefiner_MultiGrid::HNCM_PARTIAL: sh.assign_subset(*iter, 2); break;
609 0 : case HangingNodeRefiner_MultiGrid::HNCM_NONE: sh.assign_subset(*iter, 3); break;
610 0 : case HangingNodeRefiner_MultiGrid::HNCM_UNKNOWN: sh.assign_subset(*iter, 4); break;
611 0 : case HangingNodeRefiner_MultiGrid::HNCM_INVALID: sh.assign_subset(*iter, 5); break;
612 0 : case HangingNodeRefiner_MultiGrid::HNCM_REPLACE: sh.assign_subset(*iter, 6); break;
613 0 : case HangingNodeRefiner_MultiGrid::HNCM_NO_NBRS: sh.assign_subset(*iter, 7); break;
614 : }
615 : }
616 :
617 0 : for(FaceIterator iter = g.faces_begin(); iter != g.faces_end(); ++iter){
618 : ISelector::status_t status = sel.get_selection_status(*iter);
619 0 : switch(status){
620 0 : case RM_COARSEN: sh.assign_subset(*iter, 0); break;
621 0 : case HangingNodeRefiner_MultiGrid::HNCM_ALL: sh.assign_subset(*iter, 1); break;
622 0 : case HangingNodeRefiner_MultiGrid::HNCM_PARTIAL: sh.assign_subset(*iter, 2); break;
623 0 : case HangingNodeRefiner_MultiGrid::HNCM_NONE: sh.assign_subset(*iter, 3); break;
624 0 : case HangingNodeRefiner_MultiGrid::HNCM_UNKNOWN: sh.assign_subset(*iter, 4); break;
625 0 : case HangingNodeRefiner_MultiGrid::HNCM_INVALID: sh.assign_subset(*iter, 5); break;
626 0 : case HangingNodeRefiner_MultiGrid::HNCM_REPLACE: sh.assign_subset(*iter, 6); break;
627 0 : case HangingNodeRefiner_MultiGrid::HNCM_NO_NBRS: sh.assign_subset(*iter, 7); break;
628 : }
629 : }
630 :
631 0 : for(VolumeIterator iter = g.volumes_begin(); iter != g.volumes_end(); ++iter){
632 : ISelector::status_t status = sel.get_selection_status(*iter);
633 0 : switch(status){
634 0 : case RM_COARSEN: sh.assign_subset(*iter, 0); break;
635 0 : case HangingNodeRefiner_MultiGrid::HNCM_ALL: sh.assign_subset(*iter, 1); break;
636 0 : case HangingNodeRefiner_MultiGrid::HNCM_PARTIAL: sh.assign_subset(*iter, 2); break;
637 0 : case HangingNodeRefiner_MultiGrid::HNCM_NONE: sh.assign_subset(*iter, 3); break;
638 0 : case HangingNodeRefiner_MultiGrid::HNCM_UNKNOWN: sh.assign_subset(*iter, 4); break;
639 0 : case HangingNodeRefiner_MultiGrid::HNCM_INVALID: sh.assign_subset(*iter, 5); break;
640 0 : case HangingNodeRefiner_MultiGrid::HNCM_REPLACE: sh.assign_subset(*iter, 6); break;
641 0 : case HangingNodeRefiner_MultiGrid::HNCM_NO_NBRS: sh.assign_subset(*iter, 7); break;
642 : }
643 : }
644 :
645 0 : sh.subset_info(0).name = "coarsen";
646 0 : sh.subset_info(1).name = "all";
647 0 : sh.subset_info(2).name = "partial";
648 0 : sh.subset_info(3).name = "none";
649 0 : sh.subset_info(4).name = "unknown";
650 0 : sh.subset_info(5).name = "invalid";
651 0 : sh.subset_info(6).name = "replace";
652 0 : sh.subset_info(7).name = "no_nbrs";
653 0 : sh.subset_info(8).name = "unassigned";
654 0 : AssignSubsetColors(sh);
655 :
656 :
657 0 : if(MultiGrid* mg = dynamic_cast<MultiGrid*>(&g)){
658 0 : if(mg->num<Volume>() > 0)
659 0 : SaveGridHierarchyTransformed(*mg, sh, filename, 3);
660 : else
661 0 : SaveGridHierarchyTransformed(*mg, sh, filename, 0.1);
662 : }
663 : else
664 0 : SaveGridToFile(g, sh, filename);
665 0 : }
666 :
667 0 : void HangingNodeRefiner_MultiGrid::
668 : save_coarsen_marks_to_file(ISelector& sel, const char* filename)
669 : {
670 0 : SaveCoarsenMarksToFile(sel, filename);
671 0 : }
672 :
673 :
674 : /// temporary method, which will be removed when debugging is done.
675 0 : void HangingNodeRefiner_MultiGrid::
676 : debug_save(ISelector& sel, const char* filePrefix)
677 : {
678 0 : if(debugging_enabled()){
679 0 : stringstream ss;
680 0 : ss << filePrefix << "_p";
681 : #ifdef UG_PARALLEL
682 : ss << pcl::ProcRank();
683 : #else
684 0 : ss << "0";
685 : #endif
686 0 : ss << ".ugx";
687 : //UG_LOG("Saving coarsen debug marks: " << ss.str() << endl);
688 0 : save_coarsen_marks_to_file(sel, ss.str().c_str());
689 0 : }
690 0 : }
691 :
692 : /// temporary method, which will be removed when debugging is done.
693 0 : static void ContinuousDebugSave(ISelector& sel)
694 : {
695 : static int counter = 0;
696 0 : stringstream ss;
697 0 : ss << "coarsen_debug_" << counter << "_p";
698 : #ifdef UG_PARALLEL
699 : ss << pcl::ProcRank();
700 : #else
701 0 : ss << "0";
702 : #endif
703 0 : ss << ".ugx";
704 0 : UG_LOG("Saving coarsen debug marks: " << ss.str() << endl);
705 0 : SaveCoarsenMarksToFile(sel, ss.str().c_str());
706 0 : ++counter;
707 0 : }
708 :
709 : /// temporary method, which will be removed when debugging is done.
710 0 : static void ParallelLayoutDebugSave(MultiGrid& mg)
711 : {
712 : static int counter = 0;
713 0 : stringstream ss;
714 0 : ss << "coarsen_parallel_layout_" << counter << "_p";
715 : #ifdef UG_PARALLEL
716 : ss << pcl::ProcRank();
717 : #else
718 0 : ss << "0";
719 : #endif
720 0 : ss << ".ugx";
721 0 : UG_LOG("Saving coarsen parallel layout: " << ss.str() << endl);
722 0 : if(mg.num<Volume>() > 0)
723 0 : SaveParallelGridLayout(mg, ss.str().c_str(), 3);
724 : else
725 0 : SaveParallelGridLayout(mg, ss.str().c_str(), 0.1);
726 0 : ++counter;
727 0 : }
728 :
729 :
730 : //void HangingNodeRefiner_MultiGrid::
731 : //collect_objects_for_coarsen(bool scheduleCoarseningBeginsMessage)
732 : //{
733 : // MultiGrid& mg = *m_pMG;
734 : // selector_t& sel = get_refmark_selector();
735 : //
736 : // bool gotVols = contains_volumes();
737 : // bool gotFaces = contains_faces();
738 : // bool gotEdges = contains_edges();
739 : //
740 : //debug_save(sel, "coarsen_marks_01_start");
741 : //// deselect all elements of lower dimension, since we currently can't handle them correctly...
742 : ////todo: A more sophisticated approach could be used. E.g. only deselect elements if they
743 : //// are connected to an element of higher dimension.
744 : // if(gotVols)
745 : // sel.clear<Face>();
746 : //
747 : // if(gotFaces)
748 : // sel.clear<Edge>();
749 : //
750 : // if(gotEdges)
751 : // sel.clear<Vertex>();
752 : //
753 : //debug_save(sel, "coarsen_marks_01.5_only_full_dim_selected");
754 : // restrict_selection_to_surface_coarsen_elements();
755 : // copy_marks_to_vmasters(true, true, true, true);
756 : // restrict_selection_to_coarsen_families();
757 : // copy_marks_to_vslaves(true, true, true, true);
758 : //
759 : // SelectAssociatedEdges(sel, sel.begin<Face>(), sel.end<Face>(), HNCM_UNKNOWN);
760 : // SelectAssociatedEdges(sel, sel.begin<Volume>(), sel.end<Volume>(), HNCM_UNKNOWN);
761 : // broadcast_marks_horizontally(false, true, false);
762 : // copy_marks_to_vmasters(false, true, false, false);
763 : //
764 : //debug_save(sel, "coarsen_marks_02_restricted_to_surface_families");
765 : //
766 : // typedef vector<Edge*> EdgeVec;
767 : // EdgeVec vedges;
768 : // vedges.assign(sel.begin<Edge>(), sel.end<Edge>());
769 : //
770 : // Grid::edge_traits::secure_container edges;
771 : // Grid::face_traits::secure_container faces;
772 : // Grid::volume_traits::secure_container vols;
773 : //
774 : //// make sure that coarsening won't generate constraining edges with more than
775 : //// one constrained vertex.
776 : // bool running = gotVols || gotFaces;
777 : // while(running){
778 : // // classify unknown edges
779 : // for(EdgeVec::iterator iter = vedges.begin(); iter != vedges.end(); ++iter)
780 : // {
781 : // Edge* e = *iter;
782 : //
783 : // if(sel.get_selection_status(e) == HNCM_UNKNOWN){
784 : // size_t numSel = 0;
785 : // size_t numNbrs = 0;
786 : // if(gotVols){
787 : // mg.associated_elements(vols, e);
788 : // numNbrs = vols.size();
789 : // for(size_t i = 0; i < numNbrs; ++i){
790 : // if(sel.is_selected(vols[i]))
791 : // ++numSel;
792 : // }
793 : // }
794 : // else if(gotFaces){
795 : // mg.associated_elements(faces, e);
796 : // numNbrs = faces.size();
797 : // for(size_t i = 0; i < numNbrs; ++i){
798 : // if(sel.is_selected(faces[i]))
799 : // ++numSel;
800 : // }
801 : // }
802 : //
803 : // if(numNbrs > 0){
804 : // if(numSel == 0)
805 : // sel.select(e, HNCM_NONE);
806 : // else if(numSel < numNbrs)
807 : // sel.select(e, HNCM_PARTIAL);
808 : // else
809 : // sel.select(e, HNCM_ALL);
810 : // }
811 : // else
812 : // sel.select(e, HNCM_NO_NBRS);
813 : // }
814 : // }
815 : //
816 : // broadcast_marks_horizontally(false, true, false);
817 : // copy_marks_to_vmasters(false, true, false, false);
818 : //
819 : // // clear all edges which are marked with HNCM_NONE from the selection
820 : // for(selector_t::traits<Edge>::iterator iter = sel.begin<Edge>();
821 : // iter != sel.end<Edge>();)
822 : // {
823 : // Edge* e = *iter;
824 : // ++iter;
825 : // if(sel.get_selection_status(e) == HNCM_NONE)
826 : // sel.deselect(e);
827 : // }
828 : //
829 : // // if a constrained edge of a constraining edge won't be fully coarsened,
830 : // // then the constraining edge may not be coarsened at all.
831 : // bool foundInvalid = false;
832 : // for(selector_t::traits<ConstrainingEdge>::iterator iter = sel.begin<ConstrainingEdge>();
833 : // iter != sel.end<ConstrainingEdge>(); ++iter)
834 : // {
835 : // ConstrainingEdge* cge = *iter;
836 : // for(size_t i = 0; i < cge->num_constrained_edges(); ++i){
837 : // if(sel.get_selection_status(cge->constrained_edge(i)) != HNCM_ALL){
838 : // foundInvalid = true;
839 : // sel.select(cge, HNCM_INVALID);
840 : // break;
841 : // }
842 : // }
843 : // }
844 : //
845 : // // has anybody found an invalid parent? If not, exit adjustment.
846 : // if(!one_proc_true(foundInvalid)){
847 : // running = false;
848 : // break;
849 : // }
850 : //
851 : // // clear the current candidate vector
852 : // vedges.clear();
853 : //
854 : //
855 : // // exchange invalid marks
856 : // // we have to copy them to vmasters, since constraining ghost edges don't
857 : // // have an associated constrained edge from which they could obtain their
858 : // // invalid mark.
859 : // copy_marks_to_vmasters(false, true, false, false);
860 : //
861 : // for(selector_t::traits<Edge>::iterator iter = sel.begin<Edge>();
862 : // iter != sel.end<Edge>();)
863 : // {
864 : // Edge* e = *iter;
865 : // ++iter;
866 : // if(sel.get_selection_status(e) == HNCM_INVALID){
867 : // if(gotVols){
868 : // mg.associated_elements(vols, e);
869 : // for(size_t i_vol = 0; i_vol < vols.size(); ++i_vol){
870 : // if(!sel.is_selected(vols[i_vol]))
871 : // continue;
872 : // GridObject* parent = mg.get_parent(vols[i_vol]);
873 : // if(parent){
874 : // DeselectFamily(sel, mg, parent);
875 : // mg.associated_elements(edges, parent);
876 : // for(size_t i_edge = 0; i_edge < edges.size(); ++i_edge){
877 : // Edge* e = edges[i_edge];
878 : // if(sel.is_selected(e))
879 : // sel.select(e, HNCM_UNKNOWN);
880 : // for(size_t i = 0; i < mg.num_children<Edge>(e); ++i){
881 : // Edge* child = mg.get_child<Edge>(e, i);
882 : // if(sel.is_selected(child)){
883 : // sel.select(child, HNCM_UNKNOWN);
884 : // }
885 : // }
886 : // }
887 : //
888 : // // we also have to mark edge-children of side-faces as unknown
889 : // mg.associated_elements(faces, parent);
890 : // for(size_t i_face = 0; i_face < faces.size(); ++i_face){
891 : // Face* f = faces[i_face];
892 : // if(sel.is_selected(f))
893 : // sel.select(f, HNCM_UNKNOWN);
894 : // for(size_t i = 0; i < mg.num_children<Edge>(f); ++i){
895 : // Edge* child = mg.get_child<Edge>(f, i);
896 : // if(sel.is_selected(child)){
897 : // sel.select(child, HNCM_UNKNOWN);
898 : // }
899 : // }
900 : // }
901 : //
902 : // }
903 : // }
904 : // }
905 : // else if(gotFaces){
906 : // mg.associated_elements(faces, e);
907 : // for(size_t i_face = 0; i_face < faces.size(); ++i_face){
908 : // if(!sel.is_selected(faces[i_face]))
909 : // continue;
910 : // GridObject* parent = mg.get_parent(faces[i_face]);
911 : // if(parent){
912 : // DeselectFamily(sel, mg, parent);
913 : // mg.associated_elements(edges, parent);
914 : // for(size_t i_edge = 0; i_edge < edges.size(); ++i_edge){
915 : // if(sel.is_selected(edges[i_edge]))
916 : // sel.select(edges[i_edge], HNCM_UNKNOWN);
917 : // for(size_t i = 0; i < mg.num_children<Edge>(edges[i_edge]); ++i){
918 : // Edge* child = mg.get_child<Edge>(edges[i_edge], i);
919 : // if(sel.is_selected(child)){
920 : // sel.select(child, HNCM_UNKNOWN);
921 : // }
922 : // }
923 : // }
924 : // }
925 : // }
926 : // }
927 : //
928 : // // mark the formerly invalid edge as unknown
929 : // sel.select(e, HNCM_UNKNOWN);
930 : // }
931 : // }
932 : //
933 : // if(debugging_enabled()){
934 : // ContinuousDebugSave(sel);
935 : // }
936 : //
937 : // // since h-interfaces do not always exist between v-master copies,
938 : // // we have to perform an additional broadcast, to make sure, that all copies
939 : // // have the correct marks.
940 : // broadcast_marks_vertically(false, true, true, true, true);
941 : // if(debugging_enabled()){
942 : // ContinuousDebugSave(sel);
943 : // }
944 : //
945 : //
946 : //
947 : // // find edges which were marked as unknown and prepare qedges for the next run
948 : // for(selector_t::traits<Edge>::iterator iter = sel.begin<Edge>();
949 : // iter != sel.end<Edge>(); ++iter)
950 : // {
951 : // if(sel.get_selection_status(*iter) == HNCM_UNKNOWN)
952 : // vedges.push_back(*iter);
953 : // }
954 : // }
955 : //
956 : //debug_save(sel, "coarsen_marks_03_irregularities_resolved");
957 : //
958 : //// finally classify faces and vertices
959 : // SelectAssociatedFaces(sel, sel.begin<Volume>(), sel.end<Volume>(), HNCM_ALL);
960 : // SelectAssociatedVertices(sel, sel.begin<Edge>(), sel.end<Edge>(), HNCM_ALL);
961 : //
962 : // broadcast_marks_horizontally(true, false, true);
963 : // copy_marks_to_vmasters(true, false, true, false);
964 : //
965 : // if(gotVols){
966 : // for(selector_t::traits<Face>::iterator iter = sel.begin<Face>();
967 : // iter != sel.end<Face>();)
968 : // {
969 : // Face* f = *iter;
970 : // ++iter;
971 : //
972 : // mg.associated_elements(vols, f);
973 : // for(size_t i = 0; i < vols.size(); ++i){
974 : // if(!sel.is_selected(vols[i])){
975 : // sel.select(f, HNCM_PARTIAL);
976 : // break;
977 : // }
978 : // }
979 : // }
980 : // }
981 : //
982 : // for(selector_t::traits<Vertex>::iterator iter = sel.begin<Vertex>();
983 : // iter != sel.end<Vertex>();)
984 : // {
985 : // Vertex* v = *iter;
986 : // ++iter;
987 : //
988 : // mg.associated_elements(edges, v);
989 : // for(size_t i = 0; i < edges.size(); ++i){
990 : // if(sel.get_selection_status(edges[i]) != HNCM_ALL){
991 : // sel.select(v, HNCM_PARTIAL);
992 : // break;
993 : // }
994 : // }
995 : // }
996 : //
997 : // broadcast_marks_horizontally(true, false, true);
998 : // copy_marks_to_vmasters(true, false, true, false);
999 : //
1000 : //debug_save(sel, "coarsen_marks_04_faces_and_vertices_classified");
1001 : //
1002 : //
1003 : //// ATTENTION: We post the message that coarsening begins at this point, since
1004 : //// we don't want REPLACE elements to be selected and thus contained
1005 : //// in the specified geometric-object-collection.
1006 : // if(scheduleCoarseningBeginsMessage){
1007 : // m_messageHub->post_message(GridMessage_Adaption(GMAT_HNODE_COARSENING_BEGINS,
1008 : // m_selMarkedElements.get_grid_objects()));
1009 : // }
1010 : //
1011 : //
1012 : //// select unselected constraining edges and faces of selected constrained
1013 : //// edges and faces with HNCM_REPLACE to indicate, that they have to be
1014 : //// transformed to a normal edge
1015 : //// mark parents of normal and constraining edges and faces which were marked for
1016 : //// partial refinement with a replace mark
1017 : // for(selector_t::traits<Edge>::iterator
1018 : // iter = sel.begin<Edge>(); iter != sel.end<Edge>(); ++iter)
1019 : // {
1020 : // Edge* e = *iter;
1021 : // if(GridObject* parent = mg.get_parent(e)){
1022 : // bool isConstrained = e->is_constrained();
1023 : // if((isConstrained && (sel.get_selection_status(e) == HNCM_ALL)) ||
1024 : // ((!isConstrained) && (sel.get_selection_status(*iter) == HNCM_PARTIAL)))
1025 : // {
1026 : // if(!sel.is_selected(parent))
1027 : // sel.select(parent, HNCM_REPLACE);
1028 : // }
1029 : // }
1030 : // }
1031 : //
1032 : // for(selector_t::traits<Face>::iterator
1033 : // iter = sel.begin<Face>(); iter != sel.end<Face>(); ++iter)
1034 : // {
1035 : // Face* e = *iter;
1036 : // if(GridObject* parent = mg.get_parent(e)){
1037 : // bool isConstrained = e->is_constrained();
1038 : // if((isConstrained && (sel.get_selection_status(e) == HNCM_ALL)) ||
1039 : // ((!isConstrained) && (sel.get_selection_status(*iter) == HNCM_PARTIAL)))
1040 : // {
1041 : // if(!sel.is_selected(parent))
1042 : // sel.select(parent, HNCM_REPLACE);
1043 : // }
1044 : // }
1045 : // }
1046 : //
1047 : // broadcast_marks_horizontally(false, true, true);
1048 : // copy_marks_to_vmasters(false, true, true, false);
1049 : //
1050 : //debug_save(sel, "coarsen_marks_05_adjustment_done");
1051 : //}
1052 :
1053 0 : void HangingNodeRefiner_MultiGrid::
1054 : collect_objects_for_coarsen(bool scheduleCoarseningBeginsMessage)
1055 : {
1056 : PROFILE_FUNC();
1057 :
1058 0 : MultiGrid& mg = *m_pMG;
1059 : selector_t& sel = get_refmark_selector();
1060 :
1061 0 : bool gotVols = contains_volumes();
1062 0 : bool gotFaces = contains_faces();
1063 0 : bool gotEdges = contains_edges();
1064 :
1065 0 : debug_save(sel, "coarsen_marks_01_start");
1066 : // deselect all elements of lower dimension, since we currently can't handle them correctly...
1067 : //todo: A more sophisticated approach could be used. E.g. only deselect elements if they
1068 : // are connected to an element of higher dimension.
1069 0 : if(gotVols)
1070 0 : sel.clear<Face>();
1071 :
1072 0 : if(gotFaces)
1073 0 : sel.clear<Edge>();
1074 :
1075 0 : if(gotEdges)
1076 0 : sel.clear<Vertex>();
1077 :
1078 0 : restrict_selection_to_surface_coarsen_elements();
1079 0 : copy_marks_to_vmasters(true, true, true, true);
1080 0 : restrict_selection_to_coarsen_families();
1081 0 : copy_marks_to_vslaves(true, true, true, true);
1082 :
1083 0 : SelectAssociatedVertices(sel, sel.begin<Edge>(), sel.end<Edge>(), HNCM_UNKNOWN);
1084 0 : SelectAssociatedVertices(sel, sel.begin<Face>(), sel.end<Face>(), HNCM_UNKNOWN);
1085 0 : SelectAssociatedVertices(sel, sel.begin<Volume>(), sel.end<Volume>(), HNCM_UNKNOWN);
1086 0 : broadcast_marks_horizontally(true, false, false);
1087 0 : copy_marks_to_vmasters(true, false, false, false);
1088 :
1089 0 : debug_save(sel, "coarsen_marks_02_restricted_to_surface_families");
1090 :
1091 : typedef vector<Vertex*> VrtVec;
1092 : VrtVec vvrts;
1093 0 : vvrts.assign(sel.begin<Vertex>(), sel.end<Vertex>());
1094 :
1095 : Grid::vertex_traits::secure_container vrts;
1096 : Grid::edge_traits::secure_container edges;
1097 : Grid::face_traits::secure_container faces;
1098 : Grid::volume_traits::secure_container vols;
1099 :
1100 : // make sure that coarsening won't generate situations in which surface elements
1101 : // which meet at a common vertex have a level difference of more than 1.
1102 0 : bool running = gotVols || gotFaces;
1103 0 : while(running){
1104 : // classify unknown vertices
1105 0 : for(VrtVec::iterator iter = vvrts.begin(); iter != vvrts.end(); ++iter)
1106 : {
1107 0 : Vertex* e = *iter;
1108 :
1109 0 : if(sel.get_selection_status(e) == HNCM_UNKNOWN){
1110 : size_t numSel = 0;
1111 : size_t numNbrs = 0;
1112 0 : if(gotVols){
1113 0 : mg.associated_elements(vols, e);
1114 : numNbrs = vols.size();
1115 0 : for(size_t i = 0; i < numNbrs; ++i){
1116 0 : if(sel.is_selected(vols[i]))
1117 0 : ++numSel;
1118 : }
1119 : }
1120 0 : else if(gotFaces){
1121 0 : mg.associated_elements(faces, e);
1122 : numNbrs = faces.size();
1123 0 : for(size_t i = 0; i < numNbrs; ++i){
1124 0 : if(sel.is_selected(faces[i]))
1125 0 : ++numSel;
1126 : }
1127 : }
1128 :
1129 0 : if(numNbrs > 0){
1130 0 : if(numSel == 0)
1131 0 : sel.select(e, HNCM_NONE);
1132 0 : else if(numSel < numNbrs)
1133 0 : sel.select(e, HNCM_PARTIAL);
1134 : else
1135 0 : sel.select(e, HNCM_ALL);
1136 : }
1137 : else
1138 0 : sel.select(e, HNCM_NO_NBRS);
1139 : }
1140 : }
1141 :
1142 0 : broadcast_marks_horizontally(true, false, false);
1143 0 : copy_marks_to_vmasters(true, false, false, false);
1144 :
1145 : // clear all vertices which are marked with HNCM_NONE from the selection
1146 0 : for(selector_t::traits<Vertex>::iterator iter = sel.begin<Vertex>();
1147 0 : iter != sel.end<Vertex>();)
1148 : {
1149 : Vertex* e = *iter;
1150 : ++iter;
1151 0 : if(sel.get_selection_status(e) == HNCM_NONE)
1152 0 : sel.deselect(e);
1153 : }
1154 :
1155 : // if a vertex has a child which will not be fully coarsed then the
1156 : // vertex itself may not be coarsened at all (nor neighboring elements)
1157 : bool foundInvalid = false;
1158 0 : for(selector_t::traits<Vertex>::iterator iter = sel.begin<Vertex>();
1159 0 : iter != sel.end<Vertex>(); ++iter)
1160 : {
1161 : Vertex* e = *iter;
1162 : Vertex* c = mg.get_child_vertex(e);
1163 0 : if(c && sel.get_selection_status(c) != HNCM_ALL){
1164 : foundInvalid = true;
1165 0 : sel.select(e, HNCM_INVALID);
1166 : }
1167 : }
1168 :
1169 : // has anybody marked an element as invalid? If not, exit adjustment.
1170 0 : if(!one_proc_true(foundInvalid)){
1171 : break;
1172 : }
1173 :
1174 : // clear the current candidate vector
1175 : vvrts.clear();
1176 :
1177 : // exchange invalid marks
1178 : // we have to copy them to vmasters, since in some situations there may be
1179 : // ghost-vertices which are shadows of vertices which lie on the rim of a level...
1180 0 : copy_marks_to_vmasters(true, false, false, false);
1181 :
1182 : // we have to deselect all families which are connected to invalid vertices
1183 0 : for(selector_t::traits<Vertex>::iterator iter = sel.begin<Vertex>();
1184 0 : iter != sel.end<Vertex>();)
1185 : {
1186 : Vertex* e = *iter;
1187 : ++iter;
1188 0 : if(sel.get_selection_status(e) == HNCM_INVALID){
1189 0 : if(gotVols){
1190 0 : mg.associated_elements(vols, e);
1191 0 : for(size_t i_vol = 0; i_vol < vols.size(); ++i_vol){
1192 0 : if(!sel.is_selected(vols[i_vol]))
1193 0 : continue;
1194 : GridObject* parent = mg.get_parent(vols[i_vol]);
1195 0 : if(parent){
1196 0 : DeselectFamily(sel, mg, parent);
1197 : mg.associated_elements(vrts, parent);
1198 0 : for(size_t i_vrt = 0; i_vrt < vrts.size(); ++i_vrt){
1199 : Vertex* vrt = vrts[i_vrt];
1200 0 : if(sel.is_selected(vrt))
1201 0 : sel.select(vrt, HNCM_UNKNOWN);
1202 :
1203 : Vertex* c = mg.get_child_vertex(vrt);
1204 0 : if(c && sel.is_selected(c)){
1205 0 : sel.select(c, HNCM_UNKNOWN);
1206 : }
1207 : }
1208 :
1209 : // we also have to mark vertex-children of side-faces
1210 : // and side-edges as unknown
1211 : mg.associated_elements(edges, parent);
1212 0 : for(size_t i_edge = 0; i_edge < edges.size(); ++i_edge){
1213 : Edge* edge = edges[i_edge];
1214 0 : if(sel.is_selected(edge))
1215 0 : sel.select(edge, HNCM_UNKNOWN);
1216 :
1217 : Vertex* c = mg.get_child_vertex(edge);
1218 0 : if(c && sel.is_selected(c)){
1219 0 : sel.select(c, HNCM_UNKNOWN);
1220 : }
1221 : }
1222 :
1223 : mg.associated_elements(faces, parent);
1224 0 : for(size_t i_face = 0; i_face < faces.size(); ++i_face){
1225 : Face* f = faces[i_face];
1226 0 : if(sel.is_selected(f))
1227 0 : sel.select(f, HNCM_UNKNOWN);
1228 :
1229 : Vertex* c = mg.get_child_vertex(f);
1230 0 : if(c && sel.is_selected(c)){
1231 0 : sel.select(c, HNCM_UNKNOWN);
1232 : }
1233 : }
1234 :
1235 : }
1236 : }
1237 : }
1238 0 : else if(gotFaces){
1239 0 : mg.associated_elements(faces, e);
1240 0 : for(size_t i_face = 0; i_face < faces.size(); ++i_face){
1241 0 : if(!sel.is_selected(faces[i_face]))
1242 0 : continue;
1243 : GridObject* parent = mg.get_parent(faces[i_face]);
1244 0 : if(parent){
1245 0 : DeselectFamily(sel, mg, parent);
1246 : mg.associated_elements(vrts, parent);
1247 0 : for(size_t i_vrt = 0; i_vrt < vrts.size(); ++i_vrt){
1248 : Vertex* vrt = vrts[i_vrt];
1249 0 : if(sel.is_selected(vrt))
1250 0 : sel.select(vrt, HNCM_UNKNOWN);
1251 :
1252 : Vertex* c = mg.get_child_vertex(vrt);
1253 0 : if(c && sel.is_selected(c)){
1254 0 : sel.select(c, HNCM_UNKNOWN);
1255 : }
1256 : }
1257 :
1258 : // we also have to mark vertex-children of side-edges
1259 : mg.associated_elements(edges, parent);
1260 0 : for(size_t i_edge = 0; i_edge < edges.size(); ++i_edge){
1261 : Edge* edge = edges[i_edge];
1262 0 : if(sel.is_selected(edge))
1263 0 : sel.select(edge, HNCM_UNKNOWN);
1264 :
1265 : Vertex* c = mg.get_child_vertex(edge);
1266 0 : if(c && sel.is_selected(c)){
1267 0 : sel.select(c, HNCM_UNKNOWN);
1268 : }
1269 : }
1270 : }
1271 : }
1272 : }
1273 :
1274 : // mark the formerly invalid vertex as unknown
1275 0 : sel.select(e, HNCM_UNKNOWN);
1276 : }
1277 : }
1278 :
1279 0 : if(debugging_enabled()){
1280 0 : ContinuousDebugSave(sel);
1281 : }
1282 :
1283 : // since h-interfaces do not always exist between v-master copies,
1284 : // we have to perform an additional broadcast, to make sure, that all copies
1285 : // have the correct marks.
1286 0 : broadcast_marks_vertically(true, true, true, true, true);
1287 0 : if(debugging_enabled()){
1288 0 : ContinuousDebugSave(sel);
1289 : }
1290 :
1291 :
1292 :
1293 : // find vertices which were marked as unknown and prepare the candidates array for the next run
1294 0 : for(selector_t::traits<Vertex>::iterator iter = sel.begin<Vertex>();
1295 0 : iter != sel.end<Vertex>(); ++iter)
1296 : {
1297 0 : if(sel.get_selection_status(*iter) == HNCM_UNKNOWN)
1298 0 : vvrts.push_back(*iter);
1299 : }
1300 : }
1301 :
1302 0 : debug_save(sel, "coarsen_marks_03_irregularities_resolved");
1303 :
1304 : // finally classify faces and edges
1305 0 : SelectAssociatedFaces(sel, sel.begin<Volume>(), sel.end<Volume>(), HNCM_ALL);
1306 0 : SelectAssociatedEdges(sel, sel.begin<Face>(), sel.end<Face>(), HNCM_ALL);
1307 :
1308 : // we only have to communicate face marks if volumes exist
1309 0 : if(gotVols){
1310 0 : broadcast_marks_horizontally(false, true, true);
1311 0 : copy_marks_to_vmasters(false, true, true, false);
1312 : }
1313 : else{
1314 0 : broadcast_marks_horizontally(false, true, false);
1315 0 : copy_marks_to_vmasters(false, true, false, false);
1316 : }
1317 :
1318 0 : if(gotVols){
1319 0 : for(selector_t::traits<Face>::iterator iter = sel.begin<Face>();
1320 0 : iter != sel.end<Face>();)
1321 : {
1322 : Face* f = *iter;
1323 : ++iter;
1324 :
1325 0 : mg.associated_elements(vols, f);
1326 0 : for(size_t i = 0; i < vols.size(); ++i){
1327 0 : if(!sel.is_selected(vols[i])){
1328 0 : sel.select(f, HNCM_PARTIAL);
1329 : break;
1330 : }
1331 : }
1332 : }
1333 :
1334 0 : for(selector_t::traits<Edge>::iterator iter = sel.begin<Edge>();
1335 0 : iter != sel.end<Edge>();)
1336 : {
1337 : Edge* e = *iter;
1338 : ++iter;
1339 :
1340 0 : mg.associated_elements(vols, e);
1341 0 : for(size_t i = 0; i < vols.size(); ++i){
1342 0 : if(!sel.is_selected(vols[i])){
1343 0 : sel.select(e, HNCM_PARTIAL);
1344 : break;
1345 : }
1346 : }
1347 : }
1348 : }
1349 0 : else if(gotFaces){
1350 0 : for(selector_t::traits<Edge>::iterator iter = sel.begin<Edge>();
1351 0 : iter != sel.end<Edge>();)
1352 : {
1353 : Edge* e = *iter;
1354 : ++iter;
1355 :
1356 0 : mg.associated_elements(faces, e);
1357 0 : for(size_t i = 0; i < faces.size(); ++i){
1358 0 : if(!sel.is_selected(faces[i])){
1359 0 : sel.select(e, HNCM_PARTIAL);
1360 : break;
1361 : }
1362 : }
1363 : }
1364 : }
1365 :
1366 0 : if(gotVols){
1367 0 : broadcast_marks_horizontally(false, true, true);
1368 0 : copy_marks_to_vmasters(false, true, true, false);
1369 : }
1370 : else{
1371 0 : broadcast_marks_horizontally(false, true, false);
1372 0 : copy_marks_to_vmasters(false, true, false, false);
1373 : }
1374 :
1375 0 : debug_save(sel, "coarsen_marks_04_faces_and_vertices_classified");
1376 :
1377 :
1378 : // ATTENTION: We post the message that coarsening begins at this point, since
1379 : // we don't want REPLACE elements to be selected and thus contained
1380 : // in the specified geometric-object-collection.
1381 0 : if(scheduleCoarseningBeginsMessage){
1382 0 : m_messageHub->post_message(GridMessage_Adaption(GMAT_HNODE_COARSENING_BEGINS,
1383 0 : m_selMarkedElements.get_grid_objects()));
1384 : }
1385 :
1386 :
1387 : // select unselected constraining edges and faces of selected constrained
1388 : // edges and faces with HNCM_REPLACE to indicate, that they have to be
1389 : // transformed to a normal edge
1390 : // mark parents of normal and constraining edges and faces which were marked for
1391 : // partial refinement with a replace mark
1392 : for(selector_t::traits<Edge>::iterator
1393 0 : iter = sel.begin<Edge>(); iter != sel.end<Edge>(); ++iter)
1394 : {
1395 : Edge* e = *iter;
1396 0 : if(GridObject* parent = mg.get_parent(e)){
1397 0 : bool isConstrained = e->is_constrained();
1398 0 : if((isConstrained && (sel.get_selection_status(e) == HNCM_ALL)) ||
1399 0 : ((!isConstrained) && (sel.get_selection_status(*iter) == HNCM_PARTIAL)))
1400 : {
1401 0 : if(!sel.is_selected(parent))
1402 0 : sel.select(parent, HNCM_REPLACE);
1403 : }
1404 : }
1405 : }
1406 :
1407 0 : if(gotVols){
1408 : for(selector_t::traits<Face>::iterator
1409 0 : iter = sel.begin<Face>(); iter != sel.end<Face>(); ++iter)
1410 : {
1411 : Face* e = *iter;
1412 0 : if(GridObject* parent = mg.get_parent(e)){
1413 0 : bool isConstrained = e->is_constrained();
1414 0 : if((isConstrained && (sel.get_selection_status(e) == HNCM_ALL)) ||
1415 0 : ((!isConstrained) && (sel.get_selection_status(*iter) == HNCM_PARTIAL)))
1416 : {
1417 0 : if(!sel.is_selected(parent))
1418 0 : sel.select(parent, HNCM_REPLACE);
1419 : }
1420 : }
1421 : }
1422 : }
1423 :
1424 : if(gotVols){
1425 0 : broadcast_marks_horizontally(false, true, true);
1426 0 : copy_marks_to_vmasters(false, true, true, false);
1427 : }
1428 : else{
1429 0 : broadcast_marks_horizontally(false, true, false);
1430 0 : copy_marks_to_vmasters(false, true, false, false);
1431 : }
1432 :
1433 0 : debug_save(sel, "coarsen_marks_05_adjustment_done");
1434 0 : }
1435 :
1436 0 : bool HangingNodeRefiner_MultiGrid::
1437 : perform_coarsening()
1438 : {
1439 : PROFILE_FUNC();
1440 : /*
1441 : We have to handle elements as follows:
1442 :
1443 : normal edge:
1444 : All nbrs marked: remove
1445 : Some nbrs marked (hnode):
1446 : If has (unconstrained) parent:
1447 : convert to constrained. If parent is normal, then make parent constraining.
1448 : add edge to constrained elements of parent.
1449 : else
1450 : keep edge
1451 :
1452 : constrained edge:
1453 : All nbrs marked: remove
1454 : Some nbrs marked (hnode): keep (important for 3d)
1455 :
1456 : constraining edge (hnode) (only marked, if associated constrained edges are marked, too):
1457 : (All nbrs can't be marked)
1458 : Some nbrs marked:
1459 : If has (unconstrained) parent:
1460 : convert to constrained. If parent is normal, then make parent constraining.
1461 : add edge to constrained elements of parent.
1462 : else
1463 : convert to normal edge
1464 :
1465 : constraining edge (no hnode - no nbrs marked, but children are marked):
1466 : Convert to normal edge
1467 : */
1468 0 : if(!m_pMG)
1469 : return false;
1470 :
1471 : MultiGrid& mg = *m_pMG;
1472 : selector_t& sel = get_refmark_selector();
1473 :
1474 : // We let collect-objects-for-coarsen schedule the message that coarsening begins
1475 0 : collect_objects_for_coarsen(true);
1476 :
1477 : // if a debug file was specified, we'll now save the marks to that file
1478 0 : if(!m_adjustedMarksDebugFilename.empty())
1479 0 : save_coarsen_marks_to_file(sel, m_adjustedMarksDebugFilename.c_str());
1480 :
1481 : // inform derived classes that coarsening begins
1482 0 : pre_coarsen();
1483 :
1484 : ////////////////////
1485 : // VOLUMES
1486 : // erase elements and introduce constraining/constrained elements as required.
1487 0 : mg.erase(sel.begin<Volume>(), sel.end<Volume>());
1488 :
1489 : ////////////////////
1490 : // FACES
1491 0 : for(selector_t::traits<Face>::iterator iter = sel.begin<Face>();
1492 0 : iter != sel.end<Face>();)
1493 : {
1494 : Face* elem = *iter;
1495 : ++iter;
1496 : ISelector::status_t selState = sel.get_selection_status(elem);
1497 0 : sel.deselect(elem);
1498 :
1499 0 : switch(selState){
1500 0 : case RM_COARSEN:
1501 : case HNCM_ALL:{
1502 : // this should only be set on normal or constrained faces
1503 : assert(!elem->is_constraining());
1504 0 : mg.erase(elem);
1505 : }break;
1506 :
1507 0 : case HNCM_PARTIAL:{
1508 0 : if(elem->is_constrained()){
1509 : // a constrained element where not all nbrs are selected
1510 : // will be kept as it is.
1511 0 : continue;
1512 : }
1513 :
1514 : UG_ASSERT(mg.parent_type(elem) == FACE,
1515 : "Only a face with a parent-face may be marked for "
1516 : "partial coarsening!");
1517 :
1518 : // replace elem by a constrained face
1519 : ConstrainedFace* cdf = NULL;
1520 0 : switch(elem->reference_object_id()){
1521 0 : case ROID_TRIANGLE:
1522 0 : cdf = *mg.create_and_replace<ConstrainedTriangle>(elem);
1523 0 : break;
1524 0 : case ROID_QUADRILATERAL:
1525 0 : cdf = *mg.create_and_replace<ConstrainedQuadrilateral>(elem);
1526 0 : break;
1527 0 : default:
1528 0 : UG_THROW("Unknown face reference object type in "
1529 : "HangingNodeRefiner_MultiGrid::coarsen");
1530 : break;
1531 : }
1532 :
1533 : elem = cdf;
1534 :
1535 :
1536 : // the face has to be replaced by a constrained face.
1537 : // Volume-children don't have to be considered here, since
1538 : // they should all be marked with HNCM_ALL_NBRS_SELECTED
1539 0 : Face* parent = dynamic_cast<Face*>(mg.get_parent(elem));
1540 :
1541 0 : if(parent){
1542 : // the parent may not be a constrained object
1543 : UG_ASSERT(!parent->is_constrained(),
1544 : "An element with a constrained parent may not be marked"
1545 : " for partial coarsening!");
1546 :
1547 : // check whether the parent is already constraining.
1548 : // if not, it will be transformed to a constraining element later on
1549 : // and the connection is established then.
1550 0 : ConstrainingFace* cf = dynamic_cast<ConstrainingFace*>(parent);
1551 0 : if(cf){
1552 0 : cf->add_constrained_object(cdf);
1553 : cdf->set_constraining_object(cf);
1554 : }
1555 : }
1556 : else{
1557 : // There is no parent on the local process. We thus only set the
1558 : // type of the constraining element
1559 0 : cdf->set_parent_base_object_id(mg.parent_type(cdf));
1560 : }
1561 :
1562 : }break;
1563 :
1564 0 : case HNCM_REPLACE:{
1565 : assert(!elem->is_constrained());
1566 0 : if(elem->is_constraining()){
1567 : #ifdef UG_DEBUG
1568 : // make sure that the associated constrained vertex will be removed
1569 : ConstrainingFace* cf = dynamic_cast<ConstrainingFace*>(elem);
1570 : if(cf->num_constrained_vertices() > 0){
1571 : Vertex* hv = cf->constrained_vertex(0);
1572 : UG_ASSERT(sel.get_selection_status(hv) == HNCM_ALL,
1573 : "A constraining face may only be replaced by a normal face"
1574 : " during coarsening, if its constrained vertex will be removed."
1575 : " At: " << ElementDebugInfo(mg, elem));
1576 : }
1577 : #endif
1578 :
1579 : // the constraining face has to be transformed to a normal face
1580 0 : switch(elem->reference_object_id()){
1581 0 : case ROID_TRIANGLE:
1582 0 : mg.create_and_replace<Triangle>(elem);
1583 0 : break;
1584 0 : case ROID_QUADRILATERAL:
1585 0 : mg.create_and_replace<Quadrilateral>(elem);
1586 0 : break;
1587 0 : default:
1588 0 : UG_THROW("Unknown face reference object type in "
1589 : "HangingNodeRefiner_MultiGrid::coarsen");
1590 : break;
1591 : }
1592 : }
1593 : else{
1594 : // transform the face to a constraining face
1595 : ConstrainingFace* cf = NULL;
1596 0 : switch(elem->reference_object_id()){
1597 0 : case ROID_TRIANGLE:
1598 0 : cf = *mg.create_and_replace<ConstrainingTriangle>(elem);
1599 0 : break;
1600 0 : case ROID_QUADRILATERAL:
1601 0 : cf = *mg.create_and_replace<ConstrainingQuadrilateral>(elem);
1602 0 : break;
1603 0 : default:
1604 0 : UG_THROW("Unknown face reference object type in "
1605 : "HangingNodeRefiner_MultiGrid::coarsen");
1606 : break;
1607 : }
1608 :
1609 : // make sure that a connection with constrained children is established
1610 : // note - associated constrained vertices and constrained edges can
1611 : // not exist at this point.
1612 0 : for(size_t i = 0; i < mg.num_children<Face>(cf); ++i){
1613 0 : if(ConstrainedFace* cdf =
1614 0 : dynamic_cast<ConstrainedFace*>(mg.get_child<Face>(cf, i)))
1615 : {
1616 0 : cf->add_constrained_object(cdf);
1617 : cdf->set_constraining_object(cf);
1618 : }
1619 : }
1620 : }
1621 :
1622 : }break;
1623 :
1624 0 : default:
1625 0 : UG_THROW("Bad selection mark on face during "
1626 : "HangingNodeRefiner_MultiGrid::coarsen. Internal error.");
1627 : break;
1628 : }
1629 : }
1630 :
1631 : ////////////////////
1632 : // EDGES
1633 0 : for(selector_t::traits<Edge>::iterator iter = sel.begin<Edge>();
1634 0 : iter != sel.end<Edge>();)
1635 : {
1636 : Edge* elem = *iter;
1637 : ++iter;
1638 : ISelector::status_t selState = sel.get_selection_status(elem);
1639 0 : sel.deselect(elem);
1640 :
1641 0 : switch(selState){
1642 0 : case RM_COARSEN:
1643 : case HNCM_ALL:{
1644 : // this should only be set on normal or constrained edges
1645 : assert(!elem->is_constraining());
1646 0 : mg.erase(elem);
1647 : }break;
1648 :
1649 0 : case HNCM_PARTIAL:{
1650 0 : if(elem->is_constrained()){
1651 : // a constrained element where not all nbrs are selected
1652 : // will be kept as it is.
1653 0 : continue;
1654 : }
1655 :
1656 : UG_ASSERT(mg.parent_type(elem) != VOLUME,
1657 : "Children of volume elements may not be marked for"
1658 : " partial coarsening!");
1659 :
1660 : // replace the edge by a constrained edge
1661 0 : ConstrainedEdge* cde = *mg.create_and_replace<ConstrainedEdge>(elem);
1662 : elem = cde;
1663 :
1664 : // establish connection to constraining element
1665 : GridObject* parent = mg.get_parent(cde);
1666 0 : if(parent){
1667 0 : switch(parent->base_object_id()){
1668 : case EDGE:{
1669 0 : Edge* edgeParent = dynamic_cast<Edge*>(parent);
1670 : UG_ASSERT(edgeParent, "Severe type error during coarsing - bad parent type: "
1671 : << ElementDebugInfo(mg, elem));
1672 :
1673 : // check whether the parent is already constraining
1674 0 : ConstrainingEdge* ce = dynamic_cast<ConstrainingEdge*>(edgeParent);
1675 0 : if(ce){
1676 0 : ce->add_constrained_object(cde);
1677 : cde->set_constraining_object(ce);
1678 : }
1679 : }break;
1680 : case FACE:{
1681 : // if the parent is a face, then it the face should have been converted
1682 : // to a constraining face during the face-coarsening step. If it hasn't
1683 : // been converted, then the somethings wrong with the selection states.
1684 0 : ConstrainingFace* cf = dynamic_cast<ConstrainingFace*>(parent);
1685 : UG_ASSERT(cf, "The parent face should have been transformed"
1686 : " to a constraining face already: "
1687 : << ElementDebugInfo(mg, cf));
1688 0 : cf->add_constrained_object(cde);
1689 : cde->set_constraining_object(cf);
1690 0 : continue;
1691 0 : }break;
1692 0 : default:
1693 0 : UG_THROW("Unsupported parent type during partial coarsening: "
1694 : << ElementDebugInfo(mg, parent));
1695 : break;
1696 : }
1697 : }
1698 : else{
1699 : // since no parent is present on the local process, we have to
1700 : // manually set the type of the constraining element.
1701 0 : cde->set_parent_base_object_id(mg.parent_type(cde));
1702 : }
1703 : }break;
1704 :
1705 0 : case HNCM_REPLACE:{
1706 : // this should only not be set on constrained edges
1707 : UG_ASSERT(!elem->is_constrained(), "RegularEdge should not be constrained: "
1708 : << ElementDebugInfo(mg, elem));
1709 0 : if(elem->is_constraining()){
1710 : // the constraining edge has to be transformed to a normal edge
1711 : #ifdef UG_DEBUG
1712 : // make sure that the associated constrained vertex will be removed
1713 : ConstrainingEdge* ce = dynamic_cast<ConstrainingEdge*>(elem);
1714 : if(ce->num_constrained_vertices() > 0){
1715 : Vertex* hv = ce->constrained_vertex(0);
1716 : UG_ASSERT(sel.get_selection_status(hv) == HNCM_ALL,
1717 : "A constraining edge may only be replaced by a normal edge"
1718 : " during coarsening, if its constrained vertex will be removed."
1719 : " At: " << ElementDebugInfo(mg, elem));
1720 : }
1721 : #endif
1722 :
1723 0 : mg.create_and_replace<RegularEdge>(elem);
1724 : }
1725 : else{
1726 : // transform the edge to a constraining edge
1727 0 : ConstrainingEdge* cge = *mg.create_and_replace<ConstrainingEdge>(elem);
1728 : // make sure that a connection with constrained children is established
1729 : // note - associated constrained vertices and constrained edges can
1730 : // not exist at this point.
1731 0 : for(size_t i = 0; i < mg.num_children<Edge>(cge); ++i){
1732 0 : if(ConstrainedEdge* cde =
1733 0 : dynamic_cast<ConstrainedEdge*>(mg.get_child<Edge>(cge, i)))
1734 : {
1735 0 : cge->add_constrained_object(cde);
1736 : cde->set_constraining_object(cge);
1737 : }
1738 : }
1739 : }
1740 : }break;
1741 :
1742 0 : default:
1743 0 : UG_THROW("Bad selection mark on edge during "
1744 : "HangingNodeRefiner_MultiGrid::coarsen. Internal error: "
1745 : << ElementDebugInfo(mg, elem));
1746 : break;
1747 : }
1748 : }
1749 :
1750 : ////////////////////
1751 : // VERTICES
1752 0 : for(selector_t::traits<Vertex>::iterator iter = sel.begin<Vertex>();
1753 0 : iter != sel.end<Vertex>();)
1754 : {
1755 : Vertex* elem = *iter;
1756 : ++iter;
1757 : ISelector::status_t selState = sel.get_selection_status(elem);
1758 0 : sel.deselect(elem);
1759 :
1760 0 : switch(selState){
1761 0 : case RM_COARSEN:
1762 : case HNCM_ALL:{
1763 0 : mg.erase(elem);
1764 : }break;
1765 :
1766 0 : case HNCM_PARTIAL:{
1767 0 : if(elem->is_constrained()){
1768 : // a constrained element where not all nbrs are selected
1769 : // will be kept as it is.
1770 0 : continue;
1771 : }
1772 :
1773 : // a constrained vertex is only created if the parent is an edge or a face
1774 : ConstrainedVertex* hv = NULL;
1775 : char parentType = mg.parent_type(elem);
1776 0 : if((parentType == EDGE) || (parentType == FACE)){
1777 0 : hv = *mg.create_and_replace<ConstrainedVertex>(elem);
1778 : // elem = hv; // never used
1779 : }
1780 : else
1781 : break;
1782 :
1783 : // associated constrained and constraining elements
1784 : GridObject* parent = mg.get_parent(hv);
1785 :
1786 0 : if(parent){
1787 : // the parent may not be a constrained object
1788 : UG_ASSERT(!parent->is_constrained(),
1789 : "An element with a constrained parent may not be marked"
1790 : " for partial coarsening: " << ElementDebugInfo(mg, parent));
1791 :
1792 : // if the parent is an edge or a face, the vertex will be converted
1793 : // to a hanging-vertex
1794 0 : switch(parent->base_object_id()){
1795 : case EDGE:{
1796 : // the parent already has to be a constraining edge, due to
1797 : // the operations performed during edge-coarsening
1798 0 : ConstrainingEdge* ce = dynamic_cast<ConstrainingEdge*>(parent);
1799 : UG_ASSERT(ce, "The parent edge already has to be constraining: "
1800 : << ElementDebugInfo(mg, parent));
1801 : UG_ASSERT(ce->num_constrained_vertices() == 0,
1802 : "The parent edge may not yet have any constrained vertices: "
1803 : << ElementDebugInfo(mg, parent));
1804 :
1805 0 : ce->add_constrained_object(hv);
1806 : hv->set_constraining_object(ce);
1807 : hv->set_local_coordinate_1(0.5);
1808 : }break;
1809 :
1810 : case FACE:{
1811 : // the parent already has to be a constraining face, due to
1812 : // the operations performed during face-coarsening
1813 0 : ConstrainingFace* cf = dynamic_cast<ConstrainingFace*>(parent);
1814 : UG_ASSERT(cf, "The parent face already has to be constraining: "
1815 : << ElementDebugInfo(mg, parent));
1816 : UG_ASSERT(cf->num_constrained_vertices() == 0,
1817 : "The parent face may not yet have any constrained vertices: "
1818 : << ElementDebugInfo(mg, parent));
1819 :
1820 0 : cf->add_constrained_object(hv);
1821 : hv->set_constraining_object(cf);
1822 : hv->set_local_coordinates(0.5, 0.5);
1823 : }break;
1824 :
1825 0 : default:
1826 0 : UG_THROW("Unsupported parent type during partial"
1827 : " vertex coarsening: " << ElementDebugInfo(mg, parent));
1828 : break;
1829 : }
1830 : }
1831 : else{
1832 : // There is no parent on the local process. We thus only set the
1833 : // type of the constraining element
1834 0 : hv->set_parent_base_object_id(mg.parent_type(hv));
1835 : }
1836 :
1837 : }break;
1838 :
1839 0 : default:
1840 0 : UG_THROW("Bad selection mark on vertex during "
1841 : "HangingNodeRefiner_MultiGrid::coarsen. Internal error: "
1842 : << ElementDebugInfo(mg, elem));
1843 : break;
1844 : }
1845 : }
1846 :
1847 0 : debug_save(sel, "coarsen_marks_06_coarsening_done");
1848 0 : if(debugging_enabled()){
1849 0 : ParallelLayoutDebugSave(mg);
1850 : }
1851 :
1852 : //ContinuousDebugSave(sel);
1853 : // inform derived classes that coarsening is done
1854 0 : post_coarsen();
1855 0 : clear_marks();
1856 0 : m_messageHub->post_message(GridMessage_Adaption(GMAT_HNODE_COARSENING_ENDS));
1857 :
1858 0 : return true;
1859 : }
1860 :
1861 : }// end of namespace
1862 :
|
