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 <vector>
34 : #include "hanging_node_refiner_base.h"
35 : #include "lib_grid/algorithms/geom_obj_util/geom_obj_util.h"
36 : #include "lib_grid/algorithms/debug_util.h"
37 : #include "lib_grid/algorithms/subset_util.h"
38 : #include "lib_grid/file_io/file_io.h"
39 : #include "ref_mark_adjusters/std_hnode_adjuster.h"
40 : #include "lib_grid/tools/selector_multi_grid.h"
41 : #include "lib_grid/tools/periodic_boundary_manager.h"
42 :
43 : #ifdef UG_PARALLEL
44 : #include "pcl/pcl_util.h"
45 : #endif
46 :
47 : //define PROFILE_HANGING_NODE_REFINER if you want to profile
48 : //the refinement code.
49 : #define PROFILE_HANGING_NODE_REFINER
50 : #ifdef PROFILE_HANGING_NODE_REFINER
51 : #define HNODE_PROFILE_FUNC() PROFILE_FUNC()
52 : #define HNODE_PROFILE_BEGIN(name) PROFILE_BEGIN(name)
53 : #define HNODE_PROFILE_END() PROFILE_END()
54 : #else
55 : #define HNODE_PROFILE_FUNC()
56 : #define HNODE_PROFILE_BEGIN(name)
57 : #define HNODE_PROFILE_END()
58 : #endif
59 :
60 :
61 : using namespace std;
62 :
63 : namespace ug{
64 : template <class TSelector>
65 0 : HangingNodeRefinerBase<TSelector>::
66 : HangingNodeRefinerBase(SPRefinementProjector projector) :
67 : IRefiner(projector),
68 0 : m_pGrid(NULL),
69 0 : m_nodeDependencyOrder1(true),
70 0 : m_adjustingRefMarks(false)
71 : //,m_automarkHigherDimensionalObjects(false)
72 : {
73 0 : add_ref_mark_adjuster(StdHNodeAdjuster::create());
74 0 : }
75 :
76 : template <class TSelector>
77 0 : HangingNodeRefinerBase<TSelector>::
78 0 : HangingNodeRefinerBase(const HangingNodeRefinerBase&)
79 : {
80 0 : throw(UGError("no copy construction allowed."));
81 0 : }
82 :
83 : template <class TSelector>
84 0 : HangingNodeRefinerBase<TSelector>::
85 : ~HangingNodeRefinerBase()
86 : {
87 0 : if(m_pGrid)
88 : {
89 0 : m_pGrid->unregister_observer(this);
90 : }
91 0 : }
92 :
93 : template <class TSelector>
94 0 : void HangingNodeRefinerBase<TSelector>::
95 : set_grid(typename TSelector::grid_type* grid)
96 : {
97 0 : if(m_pGrid)
98 : {
99 0 : m_pGrid->unregister_observer(this);
100 0 : m_selMarkedElements.assign_grid(NULL);
101 0 : m_pGrid = NULL;
102 : }
103 :
104 0 : if(grid){
105 0 : m_pGrid = grid;
106 0 : grid->register_observer(this, OT_GRID_OBSERVER);
107 0 : m_selMarkedElements.assign_grid(*grid);
108 0 : m_selMarkedElements.enable_autoselection(false);
109 0 : m_selMarkedElements.enable_selection_inheritance(false);
110 0 : set_message_hub(grid->message_hub());
111 : }
112 0 : }
113 :
114 : template <class TSelector> void
115 0 : HangingNodeRefinerBase<TSelector>::grid_to_be_destroyed(Grid* grid)
116 : {
117 0 : m_pGrid = NULL;
118 0 : }
119 :
120 : template <class TSelector> void
121 0 : HangingNodeRefinerBase<TSelector>::clear_marks()
122 : {
123 0 : m_selMarkedElements.clear();
124 : m_newlyMarkedRefVrts.clear();
125 : m_newlyMarkedRefEdges.clear();
126 : m_newlyMarkedRefFaces.clear();
127 : m_newlyMarkedRefVols.clear();
128 0 : }
129 :
130 : template <class TSelector>
131 0 : bool HangingNodeRefinerBase<TSelector>::
132 : mark(Vertex* v, RefinementMark refMark)
133 : {
134 : assert(m_pGrid && "ERROR in HangingNodeRefinerBase::mark_for_refinement(...): No grid assigned.");
135 0 : if(get_mark(v) != refMark){
136 : // if ref-mark-adjustment is taking place, we'll also allow for the selection
137 : // of non-surface vertices
138 0 : if(refinement_is_allowed(v) || m_adjustingRefMarks){
139 0 : m_selMarkedElements.select(v, refMark);
140 0 : if(m_adjustingRefMarks && (refMark & (RM_FULL | RM_LOCAL | RM_CLOSURE | RM_DUMMY)))
141 0 : m_newlyMarkedRefVrts.push_back(v);
142 0 : return true;
143 : }
144 : }
145 : return false;
146 : }
147 :
148 : template <class TSelector>
149 0 : bool HangingNodeRefinerBase<TSelector>::mark(Edge* e, RefinementMark refMark)
150 : {
151 : assert(m_pGrid && "ERROR in HangingNodeRefinerBase::mark_for_refinement(...): No grid assigned.");
152 0 : if(get_mark(e) != refMark){
153 0 : if(refinement_is_allowed(e)){
154 0 : m_selMarkedElements.select(e, refMark);
155 0 : if(m_adjustingRefMarks && (refMark & (RM_REFINE | RM_LOCAL | RM_CLOSURE | RM_DUMMY)))
156 0 : m_newlyMarkedRefEdges.push_back(e);
157 0 : return true;
158 : }
159 : }
160 : return false;
161 : }
162 :
163 : template <class TSelector>
164 0 : bool HangingNodeRefinerBase<TSelector>::mark(Face* f, RefinementMark refMark)
165 : {
166 : assert(m_pGrid && "ERROR in HangingNodeRefinerBase::mark_for_refinement(...): No grid assigned.");
167 :
168 0 : if(get_mark(f) != refMark){
169 0 : if(refinement_is_allowed(f)){
170 0 : m_selMarkedElements.select(f, refMark);
171 0 : if(m_adjustingRefMarks && (refMark & (RM_REFINE | RM_LOCAL | RM_CLOSURE | RM_DUMMY)))
172 0 : m_newlyMarkedRefFaces.push_back(f);
173 0 : return true;
174 : }
175 : }
176 : return false;
177 : }
178 :
179 : template <class TSelector>
180 0 : bool HangingNodeRefinerBase<TSelector>::mark(Volume* v, RefinementMark refMark)
181 : {
182 : assert(m_pGrid && "ERROR in HangingNodeRefinerBase::mark_for_refinement(...): No grid assigned.");
183 0 : if(get_mark(v) != refMark){
184 0 : if(refinement_is_allowed(v)){
185 0 : m_selMarkedElements.select(v, refMark);
186 0 : if(m_adjustingRefMarks && (refMark & (RM_REFINE | RM_LOCAL | RM_CLOSURE | RM_DUMMY)))
187 0 : m_newlyMarkedRefVols.push_back(v);
188 0 : return true;
189 : }
190 : }
191 : return false;
192 : }
193 :
194 : template <class TSelector>
195 0 : void HangingNodeRefinerBase<TSelector>::
196 : mark_neighborhood(size_t numIterations, RefinementMark refMark, bool sideNbrsOnly)
197 : {
198 0 : if(!m_pGrid)
199 0 : return;
200 :
201 : typedef typename TSelector::template traits<Vertex>::iterator SelVrtIter;
202 : typedef typename TSelector::template traits<Edge>::iterator SelEdgeIter;
203 : typedef typename TSelector::template traits<Face>::iterator SelFaceIter;
204 : typedef typename TSelector::template traits<Volume>::iterator SelVolIter;
205 :
206 : Grid::edge_traits::secure_container edges;
207 : Grid::face_traits::secure_container faces;
208 : Grid::volume_traits::secure_container vols;
209 :
210 0 : EdgeDescriptor edesc;
211 :
212 0 : TSelector& sel = m_selMarkedElements;
213 0 : Grid& g = *m_pGrid;
214 :
215 : #ifdef UG_PARALLEL
216 : bool hasVolumes = pcl::OneProcTrue(g.num_volumes() > 0);
217 : bool hasFaces;
218 : if(hasVolumes)
219 : hasFaces = true;
220 : else
221 : hasFaces = pcl::OneProcTrue(g.num_faces() > 0);
222 : bool hasEdges;
223 : if(hasFaces)
224 : hasEdges = true;
225 : else
226 : hasEdges = pcl::OneProcTrue(g.num_edges() > 0);
227 : #else
228 : bool hasEdges = g.num_edges() > 0;
229 : bool hasFaces = g.num_faces() > 0;
230 : bool hasVolumes = g.num_volumes() > 0;
231 : #endif
232 :
233 : //todo: Speed could be improved by only considering newly marked elements after each iteration.
234 0 : for(size_t iteration = 0; iteration < numIterations; ++iteration){
235 0 : if(sideNbrsOnly){
236 : // first we'll select all unselected sides which are connected to marked elements
237 0 : for(SelVolIter iter = sel.template begin<Volume>();
238 0 : iter != sel.template end<Volume>(); ++iter)
239 : {
240 : Volume* vol = *iter;
241 0 : RefinementMark s = get_mark(vol);
242 : g.associated_elements(faces, vol);
243 0 : for_each_in_vec(Face* f, faces){
244 0 : if(!sel.is_selected(f) && this->refinement_is_allowed(f))
245 0 : this->mark(f, s);
246 : }end_for;
247 : }
248 :
249 0 : for(SelFaceIter iter = sel.template begin<Face>();
250 0 : iter != sel.template end<Face>(); ++iter)
251 : {
252 : Face* f = *iter;
253 0 : RefinementMark s = get_mark(f);
254 : g.associated_elements(edges, f);
255 0 : for_each_in_vec(Edge* e, edges){
256 0 : if(!sel.is_selected(e) && this->refinement_is_allowed(e))
257 0 : this->mark(e, s);
258 : }end_for;
259 : }
260 :
261 0 : for(SelEdgeIter iter = sel.template begin<Edge>();
262 0 : iter != sel.template end<Edge>(); ++iter)
263 : {
264 : Edge* e = *iter;
265 0 : RefinementMark s = get_mark(e);
266 0 : for(size_t i = 0; i < 2; ++i){
267 0 : Vertex* vrt = e->vertex(i);
268 0 : if(!sel.is_selected(vrt) && this->refinement_is_allowed(vrt))
269 0 : mark(vrt, s);
270 : }
271 : }
272 : }
273 : else{
274 : // first we'll select all vertices which are connected to marked elements
275 0 : for(SelEdgeIter iter = sel.template begin<Edge>();
276 0 : iter != sel.template end<Edge>(); ++iter)
277 : {
278 : Edge* e = *iter;
279 0 : sel.select(e->vertex(0), sel.get_selection_status(e));
280 0 : sel.select(e->vertex(1), sel.get_selection_status(e));
281 : }
282 :
283 0 : for(SelFaceIter iter = sel.template begin<Face>();
284 0 : iter != sel.template end<Face>(); ++iter)
285 : {
286 : Face* f = *iter;
287 : ISelector::status_t s = sel.get_selection_status(f);
288 0 : Face::ConstVertexArray faceVrts = f->vertices();
289 0 : for(size_t i = 0; i < f->num_vertices(); ++i)
290 0 : sel.select(faceVrts[i], s);
291 : }
292 :
293 0 : for(SelVolIter iter = sel.template begin<Volume>();
294 0 : iter != sel.template end<Volume>(); ++iter)
295 : {
296 : Volume* vol = *iter;
297 : ISelector::status_t s = sel.get_selection_status(vol);
298 0 : Volume::ConstVertexArray volVrts = vol->vertices();
299 0 : for(size_t i = 0; i < vol->num_vertices(); ++i)
300 0 : sel.select(volVrts[i], s);
301 : }
302 : }
303 :
304 : // if we're in a parallel environment, we have to broadcast the current selection
305 0 : sel.broadcast_selection_states();
306 :
307 0 : if(sideNbrsOnly){
308 : // mark those elements which have a marked side
309 0 : if(hasVolumes){
310 : vector<Edge*> markedDummyEdges;
311 0 : lg_for_each(Volume, vol, g){
312 0 : if(this->refinement_is_allowed(vol)){
313 0 : RefinementMark s = get_mark(vol);
314 0 : if(s < refMark){
315 : g.associated_elements(faces, vol);
316 : RefinementMark rm = RM_NONE;
317 0 : for_each_in_vec(Face* f, faces){
318 0 : RefinementMark sSide = get_mark(f);
319 0 : if(sSide){
320 : rm = refMark;
321 0 : if(rm == RM_NONE){
322 : switch(sSide){
323 : case RM_REFINE: rm = RM_REFINE; break;
324 : case RM_CLOSURE: rm = RM_CLOSURE; break;
325 : case RM_LOCAL: rm = RM_LOCAL; break;
326 : case RM_COARSEN: rm = RM_COARSEN; break;
327 : default: break;
328 : }
329 : }
330 0 : mark(vol, rm);
331 : break;
332 : }
333 : }end_for;
334 :
335 0 : if(rm == RM_ANISOTROPIC){
336 : // we'll also copy edge-marks to guarantee an anisotropic refinement
337 : bool copying = true;
338 0 : while(copying){
339 : copying = false;
340 0 : for_each_in_vec(Face* f, faces){
341 : g.associated_elements(edges, f);
342 0 : for_each_in_vec(Edge* e, edges){
343 0 : RefinementMark erm = get_mark(e);
344 0 : if(erm == RM_REFINE){
345 0 : if(f->get_opposing_side(e, edesc)){
346 0 : Edge* oe = g.get_edge(edesc);
347 0 : if(oe && get_mark(oe) < erm){
348 : copying = true;
349 0 : mark(oe, RM_DUMMY);
350 0 : markedDummyEdges.push_back(oe);
351 : }
352 : }
353 : }
354 : }end_for;
355 : }end_for;
356 : }
357 : }
358 : }
359 : }
360 : }lg_end_for;
361 :
362 0 : for_each_in_vec(Edge* e, markedDummyEdges){
363 0 : mark(e, RM_REFINE);
364 : }end_for;
365 0 : }
366 0 : else if(hasFaces){
367 0 : lg_for_each(Face, f, g){
368 0 : if(this->refinement_is_allowed(f)){
369 0 : RefinementMark s = get_mark(f);
370 0 : if(s < refMark){
371 : g.associated_elements(edges, f);
372 : RefinementMark rm = RM_NONE;
373 0 : for_each_in_vec(Edge* e, edges){
374 0 : RefinementMark sSide = get_mark(e);
375 0 : if(sSide){
376 : rm = refMark;
377 0 : if(rm == RM_NONE){
378 : switch(sSide){
379 : case RM_REFINE: rm = RM_REFINE; break;
380 : case RM_CLOSURE: rm = RM_CLOSURE; break;
381 : case RM_LOCAL: rm = RM_LOCAL; break;
382 : case RM_COARSEN: rm = RM_COARSEN; break;
383 : default: break;
384 : }
385 : }
386 0 : mark(f, rm);
387 : break;
388 : }
389 : }end_for;
390 :
391 : // todo: if rm == RM_ANISOTROPIC: mark opposite edges (cf hasVolumes)
392 : }
393 : }
394 : }lg_end_for;
395 : }
396 0 : else if(hasEdges){
397 0 : lg_for_each(Edge, e, g){
398 0 : if(this->refinement_is_allowed(e)){
399 0 : RefinementMark s = get_mark(e);
400 0 : if(s < refMark){
401 0 : Edge::ConstVertexArray vrts = e->vertices();
402 0 : for(size_t i = 0; i < 2; ++i){
403 0 : ISelector::status_t sSide = sel.get_selection_status(vrts[i]);
404 0 : if(sSide){
405 : RefinementMark rm = refMark;
406 0 : if(rm == RM_NONE){
407 : switch(sSide){
408 : case RM_REFINE: rm = RM_REFINE; break;
409 : case RM_CLOSURE: rm = RM_CLOSURE; break;
410 : case RM_LOCAL: rm = RM_LOCAL; break;
411 : case RM_COARSEN: rm = RM_COARSEN; break;
412 : default: break;
413 : }
414 : }
415 0 : mark(e, rm);
416 : }
417 : }
418 : }
419 : }
420 : }lg_end_for;
421 : }
422 : }
423 : else{
424 : // mark all associated elements of marked vertices
425 0 : for(SelVrtIter iter = sel.template begin<Vertex>();
426 0 : iter != sel.template end<Vertex>(); ++iter)
427 : {
428 : ISelector::status_t s = sel.get_selection_status(*iter);
429 : RefinementMark rm = refMark;
430 0 : if(rm == RM_NONE){
431 : switch(s){
432 : case RM_REFINE: rm = RM_REFINE; break;
433 : case RM_CLOSURE: rm = RM_CLOSURE; break;
434 : case RM_LOCAL: rm = RM_LOCAL; break;
435 : case RM_COARSEN: rm = RM_COARSEN; break;
436 : default: break;
437 : }
438 : }
439 :
440 : g.associated_elements(edges, *iter);
441 0 : for(size_t i = 0; i < edges.size(); ++i)
442 0 : if(!this->is_marked(edges[i]))
443 0 : this->mark(edges[i], rm);
444 :
445 : g.associated_elements(faces, *iter);
446 0 : for(size_t i = 0; i < faces.size(); ++i)
447 0 : if(!this->is_marked(faces[i]))
448 0 : this->mark(faces[i], rm);
449 :
450 : g.associated_elements(vols, *iter);
451 0 : for(size_t i = 0; i < vols.size(); ++i)
452 0 : if(!this->is_marked(vols[i]))
453 0 : this->mark(vols[i], rm);
454 : }
455 :
456 : // since we selected vertices which possibly may not be refined, we have to
457 : // deselect those now.
458 0 : for(SelVrtIter iter = sel.template begin<Vertex>();
459 0 : iter != sel.template end<Vertex>();)
460 : {
461 : Vertex* vrt = *iter;
462 : ++iter;
463 :
464 0 : if(!refinement_is_allowed(vrt))
465 0 : sel.deselect(vrt);
466 : }
467 : }
468 :
469 : }
470 : }
471 :
472 : template <class TSelector>
473 0 : RefinementMark HangingNodeRefinerBase<TSelector>::
474 : get_mark(Vertex* v) const
475 : {
476 : return (RefinementMark)(m_selMarkedElements.get_selection_status(v)
477 0 : & (RM_REFINE | RM_LOCAL | RM_CLOSURE | RM_COARSEN | RM_DUMMY));
478 : }
479 :
480 : template <class TSelector>
481 0 : RefinementMark HangingNodeRefinerBase<TSelector>::
482 : get_mark(Edge* e) const
483 : {
484 : return (RefinementMark)(m_selMarkedElements.get_selection_status(e)
485 0 : & (RM_REFINE | RM_LOCAL | RM_CLOSURE | RM_COARSEN | RM_DUMMY));
486 : }
487 :
488 : template <class TSelector>
489 0 : RefinementMark HangingNodeRefinerBase<TSelector>::
490 : get_mark(Face* f) const
491 : {
492 : return (RefinementMark)(m_selMarkedElements.get_selection_status(f)
493 0 : & (RM_REFINE | RM_LOCAL | RM_CLOSURE | RM_COARSEN | RM_DUMMY));
494 : }
495 :
496 : template <class TSelector>
497 0 : RefinementMark HangingNodeRefinerBase<TSelector>::
498 : get_mark(Volume* v) const
499 : {
500 : return (RefinementMark)(m_selMarkedElements.get_selection_status(v)
501 0 : & (RM_REFINE | RM_LOCAL | RM_CLOSURE | RM_COARSEN | RM_DUMMY));
502 : }
503 :
504 :
505 : template <class TSelector>
506 0 : bool HangingNodeRefinerBase<TSelector>::
507 : save_marks_to_file(const char* filename)
508 : {
509 : UG_DLOG(LIB_GRID, 1, " saving marks to file...\n");
510 0 : if(!m_pGrid){
511 0 : UG_THROW("ERROR in HangingNodeRefinerBase::save_marks_to_file: No grid assigned!");
512 : }
513 :
514 : Grid& g = *m_pGrid;
515 0 : SubsetHandler sh(g);
516 :
517 0 : AssignGridToSubset(g, sh, 0);
518 : selector_t& sel = get_refmark_selector();
519 :
520 0 : for(VertexIterator iter = g.vertices_begin(); iter != g.vertices_end(); ++iter){
521 : typename selector_t::status_t status = sel.get_selection_status(*iter);
522 0 : sh.assign_subset(*iter, status);
523 : }
524 0 : for(EdgeIterator iter = g.edges_begin(); iter != g.edges_end(); ++iter){
525 : typename selector_t::status_t status = sel.get_selection_status(*iter);
526 0 : sh.assign_subset(*iter, status);
527 : }
528 0 : for(FaceIterator iter = g.faces_begin(); iter != g.faces_end(); ++iter){
529 : typename selector_t::status_t status = sel.get_selection_status(*iter);
530 0 : sh.assign_subset(*iter, status);
531 : }
532 0 : for(VolumeIterator iter = g.volumes_begin(); iter != g.volumes_end(); ++iter){
533 : typename selector_t::status_t status = sel.get_selection_status(*iter);
534 0 : sh.assign_subset(*iter, status);
535 : }
536 :
537 0 : sh.subset_info(0).name = "_NONE_";
538 0 : for (byte si = 1; si > 0; ++si)
539 : {
540 0 : sh.subset_info(si).name = "_";
541 0 : for (byte b = 1; b != 0; b = b << 1)
542 0 : if (si & b)
543 0 : switch (b)
544 : {
545 0 : case RM_CLOSURE: sh.subset_info(si).name.append("CLOSURE_"); break;
546 0 : case RM_LOCAL: sh.subset_info(si).name.append("LOCAL_"); break;
547 0 : case RM_FULL: sh.subset_info(si).name.append("FULL_"); break;
548 0 : case RM_COARSEN: sh.subset_info(si).name.append("COARSEN_"); break;
549 0 : case RM_DUMMY: sh.subset_info(si).name.append("DUMMY_"); break;
550 0 : case HNRM_TO_NORMAL: sh.subset_info(si).name.append("TONORMAL_"); break;
551 0 : case HNRM_TO_CONSTRAINED: sh.subset_info(si).name.append("TOCONSTRAINED_"); break;
552 0 : case HNRM_TO_CONSTRAINING: sh.subset_info(si).name.append("TOCONSTRAINING_"); break;
553 0 : default: UG_THROW("Invalid status.")
554 : }
555 : }
556 : #if 0
557 : AssignGridToSubset(g, sh, 6);
558 :
559 : selector_t& sel = get_refmark_selector();
560 :
561 : for(VertexIterator iter = g.vertices_begin(); iter != g.vertices_end(); ++iter){
562 : typename selector_t::status_t status = sel.get_selection_status(*iter);
563 : switch(status){
564 : case RM_NONE: break;
565 : case RM_REFINE: sh.assign_subset(*iter, 0); break;
566 : case RM_ANISOTROPIC: sh.assign_subset(*iter, 1); break;
567 : case RM_COPY: sh.assign_subset(*iter, 2); break;
568 : case RM_CLOSURE: sh.assign_subset(*iter, 3); break;
569 : case RM_COARSEN: sh.assign_subset(*iter, 4); break;
570 : default: sh.assign_subset(*iter, 5); break;
571 : }
572 : }
573 :
574 : for(EdgeIterator iter = g.edges_begin(); iter != g.edges_end(); ++iter){
575 : typename selector_t::status_t status = sel.get_selection_status(*iter);
576 : switch(status){
577 : case RM_NONE: break;
578 : case RM_REFINE: sh.assign_subset(*iter, 0); break;
579 : case RM_ANISOTROPIC: sh.assign_subset(*iter, 1); break;
580 : case RM_COPY: sh.assign_subset(*iter, 2); break;
581 : case RM_CLOSURE: sh.assign_subset(*iter, 3); break;
582 : case RM_COARSEN: sh.assign_subset(*iter, 4); break;
583 : default: sh.assign_subset(*iter, 5); break;
584 : }
585 : }
586 :
587 : for(FaceIterator iter = g.faces_begin(); iter != g.faces_end(); ++iter){
588 : typename selector_t::status_t status = sel.get_selection_status(*iter);
589 : switch(status){
590 : case RM_NONE: break;
591 : case RM_REFINE: sh.assign_subset(*iter, 0); break;
592 : case RM_ANISOTROPIC: sh.assign_subset(*iter, 1); break;
593 : case RM_COPY: sh.assign_subset(*iter, 2); break;
594 : case RM_CLOSURE: sh.assign_subset(*iter, 3); break;
595 : case RM_COARSEN: sh.assign_subset(*iter, 4); break;
596 : default: sh.assign_subset(*iter, 5); break;
597 : }
598 : }
599 :
600 : for(VolumeIterator iter = g.volumes_begin(); iter != g.volumes_end(); ++iter){
601 : typename selector_t::status_t status = sel.get_selection_status(*iter);
602 : switch(status){
603 : case RM_NONE: break;
604 : case RM_REFINE: sh.assign_subset(*iter, 0); break;
605 : case RM_ANISOTROPIC: sh.assign_subset(*iter, 1); break;
606 : case RM_COPY: sh.assign_subset(*iter, 2); break;
607 : case RM_CLOSURE: sh.assign_subset(*iter, 3); break;
608 : case RM_COARSEN: sh.assign_subset(*iter, 4); break;
609 : default: sh.assign_subset(*iter, 5); break;
610 : }
611 : }
612 :
613 : sh.subset_info(0).name = "refine regular";
614 : sh.subset_info(1).name = "refine anisotropic";
615 : sh.subset_info(2).name = "refine copy";
616 : sh.subset_info(3).name = "refine closure";
617 : sh.subset_info(4).name = "coarsen";
618 : sh.subset_info(5).name = "combi-mark";
619 : sh.subset_info(6).name = "no marks";
620 :
621 : #endif
622 :
623 0 : AssignSubsetColors(sh);
624 0 : EraseEmptySubsets(sh);
625 :
626 0 : if(MultiGrid* pmg = dynamic_cast<MultiGrid*>(&g))
627 0 : return SaveGridHierarchyTransformed(*pmg, sh, filename, 0.1);
628 : else
629 0 : return SaveGridToFile(g, sh, filename);
630 0 : }
631 :
632 : template <class TSelector>
633 0 : void HangingNodeRefinerBase<TSelector>::
634 : enable_node_dependency_order_1(bool bEnable)
635 : {
636 0 : m_nodeDependencyOrder1 = bEnable;
637 0 : for(size_t i = 0; i < m_refMarkAdjusters.size(); ++i){
638 : m_refMarkAdjusters[i]->enable_node_dependency_order_1(bEnable);
639 : }
640 0 : }
641 :
642 : template <class TSelector>
643 0 : void HangingNodeRefinerBase<TSelector>::perform_refinement()
644 : {
645 : HNODE_PROFILE_BEGIN(perform_hnode_refinement);
646 : UG_DLOG(LIB_GRID, 1, "performing hanging-node-refine:\n");
647 :
648 0 : if(!m_pGrid)
649 0 : throw(UGError("ERROR in HangingNodeRefinerBase::refine(...): No grid assigned."));
650 :
651 0 : if(m_selMarkedElements.grid() != m_pGrid)
652 0 : throw(UGError("selector not initialized properly. Use HangingNodeRefinerBase::set_grid."));
653 :
654 : Grid& grid = *m_pGrid;
655 :
656 : // check grid options.
657 0 : if(!grid.option_is_enabled(GRIDOPT_AUTOGENERATE_SIDES))
658 : {
659 : LOG("WARNING in HangingNodeRefiner::refine(): grid option GRIDOPT_AUTOGENERATE_SIDES auto-enabled." << endl);
660 0 : grid.enable_options(GRIDOPT_AUTOGENERATE_SIDES);
661 : }
662 :
663 : // containers used for temporary results
664 : vector<Edge*> vEdges;
665 : vector<Face*> vFaces;
666 : vector<Volume*> vVols;
667 :
668 : HNODE_PROFILE_BEGIN(href_AdjustingMarks);
669 : // fills the queues with the elements that have to be refined.
670 0 : collect_objects_for_refine();
671 :
672 : // assigns hnode marks
673 0 : assign_hnode_marks();
674 : HNODE_PROFILE_END();
675 :
676 : // {
677 : // UG_LOG("DEBUG SAVE...\n");
678 : // static int refselCount = 0;
679 : // stringstream ss;
680 : // ss << "refselbefore_" << refselCount << ".ugx";
681 : // save_marks_to_file(ss.str().c_str());
682 : // ++refselCount;
683 : // }
684 :
685 : // if a debug file was specified, we'll now save the marks to that file
686 0 : if(!m_adjustedMarksDebugFilename.empty())
687 0 : save_marks_to_file(m_adjustedMarksDebugFilename.c_str());
688 :
689 0 : m_messageHub->post_message(GridMessage_Adaption(GMAT_HNODE_REFINEMENT_BEGINS,
690 : m_selMarkedElements.get_grid_objects()));
691 0 : if(projector()->refinement_begins_requires_subgrid()){
692 0 : SubGrid<IsSelected> sg(m_selMarkedElements.get_grid_objects(),
693 0 : IsSelected(m_selMarkedElements));
694 0 : projector()->refinement_begins(&sg);
695 : }
696 : else
697 0 : projector()->refinement_begins(NULL);
698 :
699 : // call pre_refine to allow derived classes to perform some actions
700 : HNODE_PROFILE_BEGIN(href_PreRefine);
701 0 : pre_refine();
702 : HNODE_PROFILE_END();
703 :
704 : ////////////////////////////////
705 : // ConstrainedVertices
706 : UG_DLOG(LIB_GRID, 1, " constrained vertices.\n");
707 : HNODE_PROFILE_BEGIN(href_ConstrainedVertices);
708 : for(typename selector_t::template traits<ConstrainedVertex>::iterator
709 0 : iter = m_selMarkedElements.template begin<ConstrainedVertex>();
710 0 : iter != m_selMarkedElements.template end<ConstrainedVertex>();)
711 : {
712 : ConstrainedVertex* cdv = *iter;
713 : ++iter;
714 0 : process_constrained_vertex(cdv);
715 : }
716 : HNODE_PROFILE_END();
717 :
718 : ////////////////////////////////
719 : // ConstrainedEdges
720 : UG_DLOG(LIB_GRID, 1, " constrained edges.\n");
721 : HNODE_PROFILE_BEGIN(href_ConstrainedEdges);
722 0 : for(typename selector_t::template traits<ConstrainedEdge>::iterator iter =
723 : m_selMarkedElements.template begin<ConstrainedEdge>();
724 0 : iter != m_selMarkedElements.template end<ConstrainedEdge>();)
725 : {
726 : ConstrainedEdge* cde = *iter;
727 : ++iter;
728 0 : if(marked_to_normal(cde))
729 0 : remove_hmark(cde, RM_REFINE);
730 0 : process_constrained_edge(cde);
731 : }
732 : HNODE_PROFILE_END();
733 :
734 : ////////////////////////////////
735 : // ConstrainingEdges
736 : // iterate through scheduled cg-edges
737 : UG_DLOG(LIB_GRID, 1, " constraining edges.\n");
738 : HNODE_PROFILE_BEGIN(href_ConstrainingEdges);
739 : {
740 0 : typename selector_t::template traits<ConstrainingEdge>::iterator iter =
741 : m_selMarkedElements.template begin<ConstrainingEdge>();
742 :
743 0 : while(iter != m_selMarkedElements.template end<ConstrainingEdge>()){
744 : ConstrainingEdge* cge = *iter;
745 : ++iter;
746 0 : if(marked_to_normal(cge))
747 0 : remove_hmark(cge, RM_REFINE);
748 : // if a constraining edge was previously created through replacement of
749 : // a constrained edge, we won't process it here
750 0 : if(!marked_to_constraining(cge))
751 0 : process_constraining_edge(cge);
752 : }
753 : }
754 : HNODE_PROFILE_END();
755 :
756 : ////////////////////////////////
757 : // Normal Edges
758 : // iterate through scheduled edges
759 : UG_DLOG(LIB_GRID, 1, " normal edges.\n");
760 : HNODE_PROFILE_BEGIN(href_NormalEdges);
761 : {
762 : typename selector_t::template traits<RegularEdge>::iterator
763 0 : iter = m_selMarkedElements.template begin<RegularEdge>();
764 0 : while(iter != m_selMarkedElements.template end<RegularEdge>())
765 : {
766 : RegularEdge* e = *iter;
767 : ++iter;
768 :
769 : // a normal edge may have previously been created by replacing a
770 : // constrained or constraining edge. Those edges won't be considered here
771 0 : if(!refinement_is_allowed(e)){
772 0 : continue;
773 : }
774 :
775 0 : if(marked_to_constraining(e)){
776 0 : refine_edge_with_hanging_vertex(e);
777 : }
778 0 : else if(marked_refine(e)){
779 0 : refine_edge_with_normal_vertex(e);
780 : }
781 : }
782 : }
783 : HNODE_PROFILE_END();
784 :
785 : ////////////////////////////////
786 : // Faces
787 : HNODE_PROFILE_BEGIN(href_ConstrainedFaces);
788 : UG_DLOG(LIB_GRID, 1, " constrained triangles.\n");
789 : for(typename selector_t::template traits<ConstrainedTriangle>::iterator
790 0 : iter = m_selMarkedElements.template begin<ConstrainedTriangle>();
791 0 : iter != m_selMarkedElements.template end<ConstrainedTriangle>();)
792 : {
793 : ConstrainedTriangle* cdf = *iter;
794 : ++iter;
795 0 : if(marked_to_normal(cdf))
796 0 : remove_hmark(cdf, RM_REFINE);
797 0 : process_constrained_face(cdf);
798 : }
799 :
800 : UG_DLOG(LIB_GRID, 1, " constrained quadrilaterals.\n");
801 : for(typename selector_t::template traits<ConstrainedQuadrilateral>::iterator
802 0 : iter = m_selMarkedElements.template begin<ConstrainedQuadrilateral>();
803 0 : iter != m_selMarkedElements.template end<ConstrainedQuadrilateral>();)
804 : {
805 : ConstrainedQuadrilateral* cdf = *iter;
806 : ++iter;
807 0 : if(marked_to_normal(cdf))
808 0 : remove_hmark(cdf, RM_REFINE);
809 0 : process_constrained_face(cdf);
810 : }
811 : HNODE_PROFILE_END();
812 :
813 : HNODE_PROFILE_BEGIN(href_ConstrainingFaces);
814 : UG_DLOG(LIB_GRID, 1, " constraining triangles.\n");
815 : // constraining triangles
816 : {
817 : typename selector_t::template traits<ConstrainingTriangle>::iterator
818 0 : iter = m_selMarkedElements.template begin<ConstrainingTriangle>();
819 0 : while(iter != m_selMarkedElements.template end<ConstrainingTriangle>()){
820 : ConstrainingTriangle* cgf = *iter;
821 : ++iter;
822 0 : if(marked_to_normal(cgf))
823 0 : remove_hmark(cgf, RM_REFINE);
824 : // if a constraining face was previously created through replacement of
825 : // a constrained face, we won't process it here
826 0 : if(!marked_to_constraining(cgf))
827 0 : process_constraining_face(cgf);
828 : }
829 : }
830 :
831 : UG_DLOG(LIB_GRID, 1, " constraining quadrilaterals.\n");
832 : // constraining quadrilaterals
833 : {
834 : typename selector_t::template traits<ConstrainingQuadrilateral>::iterator
835 0 : iter = m_selMarkedElements.template begin<ConstrainingQuadrilateral>();
836 0 : while(iter != m_selMarkedElements.template end<ConstrainingQuadrilateral>()){
837 : ConstrainingQuadrilateral* cgf = *iter;
838 : ++iter;
839 0 : if(marked_to_normal(cgf))
840 0 : remove_hmark(cgf, RM_REFINE);
841 : // if a constraining face was previously created through replacement of
842 : // a constrained face, we won't process it here
843 0 : if(!marked_to_constraining(cgf))
844 0 : process_constraining_face(cgf);
845 : }
846 : }
847 : HNODE_PROFILE_END();
848 :
849 : HNODE_PROFILE_BEGIN(href_NormalFaces);
850 : UG_DLOG(LIB_GRID, 1, " normal triangles.\n");
851 : // normal triangles
852 : {
853 : typename selector_t::template traits<Triangle>::iterator
854 0 : iter = m_selMarkedElements.template begin<Triangle>();
855 0 : while(iter != m_selMarkedElements.template end<Triangle>()){
856 : Face* f = *iter;
857 : ++iter;
858 :
859 : // a normal face may have previously been created by replacing a
860 : // constrained or constraining face. Those faces won't be considered here
861 0 : if((!refinement_is_allowed(f)) || marked_to_normal(f)){
862 0 : continue;
863 : }
864 :
865 0 : if(marked_to_constraining(f))
866 0 : refine_face_with_hanging_vertex(f);
867 0 : else if(marked_refine(f))
868 0 : refine_face_with_normal_vertex(f);
869 : }
870 : }
871 :
872 : UG_DLOG(LIB_GRID, 1, " normal quadrilaterals.\n");
873 : // normal quadrilaterals
874 : {
875 : typename selector_t::template traits<Quadrilateral>::iterator
876 0 : iter = m_selMarkedElements.template begin<Quadrilateral>();
877 0 : while(iter != m_selMarkedElements.template end<Quadrilateral>()){
878 : Face* f = *iter;
879 : ++iter;
880 :
881 : // a normal face may have previously been created by replacing a
882 : // constrained or constraining face. Those faces won't be considered here
883 0 : if((!refinement_is_allowed(f)) || marked_to_normal(f)){
884 0 : continue;
885 : }
886 :
887 0 : if(marked_to_constraining(f))
888 0 : refine_face_with_hanging_vertex(f);
889 0 : else if(marked_refine(f))
890 0 : refine_face_with_normal_vertex(f);
891 : }
892 : }
893 : HNODE_PROFILE_END();
894 :
895 : ////////////////////////////////
896 : // Volumes
897 : UG_DLOG(LIB_GRID, 1, " volumes.\n");
898 : HNODE_PROFILE_BEGIN(href_Volumes);
899 : {
900 : typename selector_t::template traits<Volume>::iterator
901 0 : iter = m_selMarkedElements.template begin<Volume>();
902 0 : while(iter != m_selMarkedElements.template end<Volume>()){
903 : Volume* vol = *iter;
904 : ++iter;
905 0 : if(refinement_is_allowed(vol)){
906 0 : refine_volume_with_normal_vertex(vol);
907 : }
908 : }
909 : }
910 : HNODE_PROFILE_END();
911 :
912 : UG_DLOG(LIB_GRID, 1, " refinement done.\n");
913 :
914 : ////////////////////////////////
915 : // call post_refine to allow derived classes to perform some actions
916 : HNODE_PROFILE_BEGIN(href_PostRefine);
917 0 : post_refine();
918 : HNODE_PROFILE_END();
919 :
920 :
921 :
922 : ////////////////////////////////
923 : // before calling the refinement callback, we make sure that only elements are
924 : // marked, which actually received new children
925 : for(typename selector_t::template traits<Vertex>::iterator
926 0 : iter = m_selMarkedElements.template begin<Vertex>();
927 0 : iter != m_selMarkedElements.template end<Vertex>();)
928 : {
929 : Vertex* e = *iter;
930 : ++iter;
931 0 : if(!(marked_refine(e)))
932 0 : m_selMarkedElements.deselect(e);
933 : }
934 : for(typename selector_t::template traits<Edge>::iterator
935 0 : iter = m_selMarkedElements.template begin<Edge>();
936 0 : iter != m_selMarkedElements.template end<Edge>();)
937 : {
938 : Edge* e = *iter;
939 : ++iter;
940 0 : if(!(marked_refine(e)))
941 0 : m_selMarkedElements.deselect(e);
942 : }
943 : for(typename selector_t::template traits<Face>::iterator
944 0 : iter = m_selMarkedElements.template begin<Face>();
945 0 : iter != m_selMarkedElements.template end<Face>();)
946 : {
947 : Face* e = *iter;
948 : ++iter;
949 0 : if(!(marked_refine(e)))
950 0 : m_selMarkedElements.deselect(e);
951 : }
952 :
953 : // {
954 : // UG_LOG("DEBUG SAVE...\n");
955 : // static int refselCount = 0;
956 : // stringstream ss;
957 : // ss << "refselafter_" << refselCount << ".ugx";
958 : // save_marks_to_file(ss.str().c_str());
959 : // ++refselCount;
960 : // }
961 :
962 0 : projector()->refinement_ends();
963 :
964 : // notify the grid's message hub that refinement ends
965 : HNODE_PROFILE_BEGIN(href_AdaptionEndsMessage);
966 0 : m_messageHub->post_message(GridMessage_Adaption(GMAT_HNODE_REFINEMENT_ENDS,
967 : m_selMarkedElements.get_grid_objects()));
968 : HNODE_PROFILE_END();
969 :
970 : ////////////////////////////////
971 : // Clean up
972 : // clear the refinement-callback if necessary
973 : HNODE_PROFILE_BEGIN(href_CleanUp);
974 0 : clear_marks();
975 : HNODE_PROFILE_END();
976 :
977 : // debugging utilities for the periodic boundary manager
978 : // if(m_pGrid->periodic_boundary_manager()){
979 : // m_pGrid->periodic_boundary_manager()->print_identification<Vertex>();
980 : // m_pGrid->periodic_boundary_manager()->print_identification<Edge>();
981 : // m_pGrid->periodic_boundary_manager()->print_identification<Face>();
982 :
983 : // UG_LOG("DEBUGGING: Checking validity of PeriodicBoundaryManager:\n");
984 : // if(m_pGrid->periodic_boundary_manager())
985 : // m_pGrid->periodic_boundary_manager()->validity_check();
986 : // }
987 :
988 : UG_DLOG(LIB_GRID, 1, " done.\n");
989 0 : }
990 :
991 : template <class TSelector>
992 : template <class TElem>
993 0 : bool HangingNodeRefinerBase<TSelector>::remove_coarsen_marks()
994 : {
995 : typedef typename selector_t::template traits<TElem>::iterator ElemIter;
996 : bool removedCoarsenMark = false;
997 0 : for(ElemIter iter = m_selMarkedElements.template begin<TElem>();
998 0 : iter != m_selMarkedElements.template end<TElem>();)
999 : {
1000 : TElem* e = *iter;
1001 : ++iter;
1002 0 : if(get_mark(e) == RM_COARSEN){
1003 0 : m_selMarkedElements.deselect(e);
1004 : removedCoarsenMark = true;
1005 : }
1006 : }
1007 :
1008 0 : return removedCoarsenMark;
1009 : }
1010 :
1011 : template <class TSelector>
1012 0 : void HangingNodeRefinerBase<TSelector>::collect_objects_for_refine()
1013 : {
1014 : HNODE_PROFILE_FUNC();
1015 : UG_DLOG(LIB_GRID, 1, "hnode_ref-start: collect_objects_for_refine\n");
1016 :
1017 0 : m_adjustingRefMarks = true;
1018 :
1019 : // build correct selection. see HangingVertexRefiner description.
1020 : // unmark all elements which are marked for coarsening
1021 :
1022 0 : bool removedCoarseMarks = remove_coarsen_marks<Vertex>();
1023 0 : removedCoarseMarks |= remove_coarsen_marks<Edge>();
1024 0 : removedCoarseMarks |= remove_coarsen_marks<Face>();
1025 0 : removedCoarseMarks |= remove_coarsen_marks<Volume>();
1026 :
1027 0 : if(removedCoarseMarks){
1028 : UG_LOG("WARNING in HangingNodeRefinerBase::collect_objects_for_refine: "
1029 : "Removed coarsen marks.\n");
1030 : }
1031 :
1032 : std::vector<Vertex*> newlyMarkedVrts;
1033 : std::vector<Edge*> newlyMarkedEdges;
1034 : std::vector<Face*> newlyMarkedFaces;
1035 : std::vector<Volume*> newlyMarkedVols;
1036 :
1037 0 : newlyMarkedVrts.assign(m_selMarkedElements.template begin<Vertex>(),
1038 : m_selMarkedElements.template end<Vertex>());
1039 0 : newlyMarkedEdges.assign(m_selMarkedElements.template begin<Edge>(),
1040 : m_selMarkedElements.template end<Edge>());
1041 0 : newlyMarkedFaces.assign(m_selMarkedElements.template begin<Face>(),
1042 : m_selMarkedElements.template end<Face>());
1043 0 : newlyMarkedVols.assign(m_selMarkedElements.template begin<Volume>(),
1044 : m_selMarkedElements.template end<Volume>());
1045 :
1046 : bool continueAdjustment = true;
1047 : bool firstAdjustment = true;
1048 :
1049 0 : while(continueAdjustment){
1050 0 : if(!firstAdjustment){
1051 : // we don't simply pass m_newlyMarkedXXX to the adjusters, since we want to
1052 : // record newly marked elements during adjustment. Those newly marked elems
1053 : // are then used for the next calls, etc.
1054 : // This is only necessary if we're not in the first adjustment iteration.
1055 : newlyMarkedVrts.swap(m_newlyMarkedRefVrts);
1056 : newlyMarkedEdges.swap(m_newlyMarkedRefEdges);
1057 : newlyMarkedFaces.swap(m_newlyMarkedRefFaces);
1058 : newlyMarkedVols.swap(m_newlyMarkedRefVols);
1059 : }
1060 :
1061 : firstAdjustment = false;
1062 :
1063 : m_newlyMarkedRefVrts.clear();
1064 : m_newlyMarkedRefEdges.clear();
1065 : m_newlyMarkedRefFaces.clear();
1066 : m_newlyMarkedRefVols.clear();
1067 :
1068 : // call the adjusters
1069 0 : for(size_t i = 0; i < m_refMarkAdjusters.size(); ++i){
1070 0 : if(m_refMarkAdjusters[i]->enabled()){
1071 0 : m_refMarkAdjusters[i]->ref_marks_changed(*this, newlyMarkedVrts,
1072 : newlyMarkedEdges,
1073 : newlyMarkedFaces,
1074 : newlyMarkedVols);
1075 : }
1076 : }
1077 :
1078 : bool continueRequired = (!m_newlyMarkedRefVrts.empty())
1079 0 : || (!m_newlyMarkedRefEdges.empty())
1080 0 : || (!m_newlyMarkedRefFaces.empty())
1081 0 : || (!m_newlyMarkedRefVols.empty());
1082 :
1083 0 : continueAdjustment = continue_collect_objects_for_refine(continueRequired);
1084 : }
1085 :
1086 0 : m_adjustingRefMarks = false;
1087 : UG_DLOG(LIB_GRID, 1, "hnode_ref-stop: collect_objects_for_refine\n");
1088 0 : }
1089 :
1090 : template <class TSelector>
1091 0 : void HangingNodeRefinerBase<TSelector>::
1092 : assign_hnode_marks()
1093 : {
1094 : UG_DLOG(LIB_GRID, 1, " assigning hnode marks...\n");
1095 : // iterate over all faces and volumes. If the element is not marked, but
1096 : // a side is marked, the side has to be marked for hnode refinement.
1097 : // Note that we won't mark any new elements for refinement here - we only adjust the marks
1098 : // or add conversion marks (HNRM_TO_NORMAL etc).
1099 : // Note also that we won't remove any marks during this algorithm (neither normal
1100 : // nor hnode marks).
1101 : Grid::edge_traits::secure_container edges;
1102 : vector<Face*> faces;
1103 : vector<Volume*> vols;
1104 :
1105 : std::vector<int> vinds;
1106 0 : vinds.reserve(4);
1107 :
1108 : // the grid
1109 : UG_ASSERT(m_pGrid, "A grid is required to perform this operation!");
1110 0 : Grid& grid = *m_pGrid;
1111 :
1112 0 : if(grid.num<Volume>() > 0){
1113 0 : for(sel_face_iter iter = m_selMarkedElements.template begin<Face>();
1114 0 : iter != m_selMarkedElements.template end<Face>(); ++iter)
1115 : {
1116 : Face* f = *iter;
1117 : CollectAssociated(vols, grid, f);
1118 : // if one of the volumes is marked with RM_LOCAL, we'll have to
1119 : // set up a local mark for the face and compare it.
1120 : bool gotLocal = false;
1121 0 : for(size_t i = 0; i < vols.size(); ++i){
1122 0 : if(marked_local(vols[i])){
1123 : gotLocal = true;
1124 : break;
1125 : }
1126 : }
1127 :
1128 0 : if(gotLocal){
1129 : int localFaceMark = 0;
1130 0 : if(marked_full(f)){
1131 0 : if(f->num_vertices() == 3)
1132 : localFaceMark = 7;
1133 0 : else if(f->num_vertices() == 4)
1134 : localFaceMark = 15;
1135 : else{
1136 0 : UG_THROW("Unsupported face type with too many vertices: " << f->num_vertices());
1137 : }
1138 : }
1139 0 : else if(marked_local(f))
1140 0 : localFaceMark = get_local_mark(f);
1141 0 : else if(marked_closure(f)){
1142 : grid.associated_elements_sorted(edges, f);
1143 0 : for(size_t i = 0; i < edges.size(); ++i){
1144 0 : if(marked_full(edges[i]))
1145 0 : localFaceMark |= (1<<i);
1146 : }
1147 : }
1148 :
1149 0 : for(size_t i_vol = 0; i_vol < vols.size(); ++i_vol){
1150 0 : Volume* v = vols[i_vol];
1151 0 : if(!refinement_is_allowed(v))
1152 0 : continue;
1153 :
1154 0 : if(!m_selMarkedElements.is_selected(v))
1155 : {
1156 0 : add_hmark(f, HNRM_TO_CONSTRAINING);
1157 : break;
1158 : }
1159 : else{
1160 0 : if(!marked_local(v))
1161 0 : continue;
1162 :
1163 0 : int volLocalMark = get_local_mark(v);
1164 : int sideMark = 0;
1165 :
1166 0 : Face::ConstVertexArray vrts = f->vertices();
1167 0 : const size_t numVrts = f->num_vertices();
1168 :
1169 0 : vinds.resize(numVrts);
1170 0 : for(size_t j = 0; j < numVrts; ++j){
1171 0 : vinds[j] = GetVertexIndex(v, vrts[j]);
1172 : }
1173 :
1174 :
1175 0 : for(size_t j = 0; j < numVrts; ++j){
1176 : const int edgeInd =
1177 0 : v->get_edge_index_from_vertices(
1178 0 : vinds[j], vinds[(j+1)%numVrts]);
1179 :
1180 0 : sideMark |= ((volLocalMark >> edgeInd) & 1) << j;
1181 : }
1182 :
1183 :
1184 0 : if(sideMark != localFaceMark){
1185 0 : add_hmark(f, HNRM_TO_CONSTRAINING);
1186 : break;
1187 : }
1188 : }
1189 : }
1190 : }
1191 : else{
1192 0 : for(size_t i = 0; i < vols.size(); ++i){
1193 0 : if(refinement_is_allowed(vols[i])
1194 0 : && (!m_selMarkedElements.is_selected(vols[i])))
1195 : {
1196 0 : add_hmark(f, HNRM_TO_CONSTRAINING);
1197 : break;
1198 : }
1199 : }
1200 : }
1201 :
1202 : }
1203 :
1204 : // make sure, that constrained faces, edges and vertices of selected
1205 : // constraining faces will be converted to normal edges, if they are not
1206 : // yet marked and to constraining edges, if they are marked.
1207 : for(typename selector_t::template traits<ConstrainingTriangle>::iterator
1208 0 : iter = m_selMarkedElements.template begin<ConstrainingTriangle>();
1209 0 : iter != m_selMarkedElements.template end<ConstrainingTriangle>(); ++iter)
1210 : {
1211 : ConstrainingTriangle* cgf = *iter;
1212 0 : if(marked_refine(cgf)){
1213 0 : add_hmark(cgf, HNRM_TO_NORMAL);
1214 0 : for(size_t i = 0; i < cgf->num_constrained_vertices(); ++i)
1215 0 : add_hmark(cgf->constrained_vertex(i), HNRM_TO_NORMAL);
1216 0 : for(size_t i = 0; i < cgf->num_constrained_edges(); ++i){
1217 : Edge* cde = cgf->constrained_edge(i);
1218 0 : if(marked_refine(cde))
1219 0 : add_hmark(cde, HNRM_TO_CONSTRAINING);
1220 : else
1221 0 : add_hmark(cde, HNRM_TO_NORMAL);
1222 : }
1223 0 : for(size_t i = 0; i < cgf->num_constrained_faces(); ++i){
1224 : Face* cdf = cgf->constrained_face(i);
1225 0 : if(marked_refine(cdf))
1226 0 : add_hmark(cdf, HNRM_TO_CONSTRAINING);
1227 : else
1228 0 : add_hmark(cdf, HNRM_TO_NORMAL);
1229 : }
1230 : }
1231 : }
1232 :
1233 : for(typename selector_t::template traits<ConstrainingQuadrilateral>::iterator
1234 0 : iter = m_selMarkedElements.template begin<ConstrainingQuadrilateral>();
1235 0 : iter != m_selMarkedElements.template end<ConstrainingQuadrilateral>(); ++iter)
1236 : {
1237 : ConstrainingQuadrilateral* cgf = *iter;
1238 0 : if(marked_refine(cgf)){
1239 0 : add_hmark(cgf, HNRM_TO_NORMAL);
1240 0 : for(size_t i = 0; i < cgf->num_constrained_vertices(); ++i)
1241 0 : add_hmark(cgf->constrained_vertex(i), HNRM_TO_NORMAL);
1242 0 : for(size_t i = 0; i < cgf->num_constrained_edges(); ++i){
1243 : Edge* cde = cgf->constrained_edge(i);
1244 0 : if(marked_refine(cde))
1245 0 : add_hmark(cde, HNRM_TO_CONSTRAINING);
1246 : else
1247 0 : add_hmark(cde, HNRM_TO_NORMAL);
1248 : }
1249 0 : for(size_t i = 0; i < cgf->num_constrained_faces(); ++i){
1250 : Face* cdf = cgf->constrained_face(i);
1251 0 : if(marked_refine(cdf))
1252 0 : add_hmark(cdf, HNRM_TO_CONSTRAINING);
1253 : else
1254 0 : add_hmark(cdf, HNRM_TO_NORMAL);
1255 : }
1256 : }
1257 : }
1258 : }
1259 :
1260 0 : if(grid.num<Face>() > 0){
1261 0 : for(sel_edge_iter iter = m_selMarkedElements.template begin<Edge>();
1262 0 : iter != m_selMarkedElements.template end<Edge>(); ++iter)
1263 : {
1264 : Edge* e = *iter;
1265 0 : if(!marked_refine(e))
1266 0 : continue;
1267 :
1268 : CollectAssociated(faces, grid, e);
1269 0 : for(size_t i = 0; i < faces.size(); ++i){
1270 0 : Face* f = faces[i];
1271 0 : if(marked_to_constraining(f))
1272 : {
1273 0 : add_hmark(e, HNRM_TO_CONSTRAINING);
1274 : break;
1275 : }
1276 :
1277 0 : if(!refinement_is_allowed(f))
1278 0 : continue;
1279 :
1280 0 : if(!m_selMarkedElements.is_selected(f))
1281 : {
1282 0 : add_hmark(e, HNRM_TO_CONSTRAINING);
1283 : break;
1284 : }
1285 0 : else if(marked_local(f)){
1286 0 : int localMark = get_local_mark(f);
1287 0 : int ei = GetEdgeIndex(f, e);
1288 0 : if(!(localMark & 1<<ei)){
1289 0 : add_hmark(e, HNRM_TO_CONSTRAINING);
1290 : break;
1291 : }
1292 : }
1293 : }
1294 : }
1295 :
1296 : for(typename selector_t::template traits<ConstrainingEdge>::iterator
1297 0 : iter = m_selMarkedElements.template begin<ConstrainingEdge>();
1298 0 : iter != m_selMarkedElements.template end<ConstrainingEdge>(); ++iter)
1299 : {
1300 : ConstrainingEdge* cge = *iter;
1301 0 : if(marked_refine(cge)){
1302 0 : add_hmark(cge, HNRM_TO_NORMAL);
1303 0 : if(!marked_to_constraining(cge)){
1304 : // this happens in volume geometries, where only some of the
1305 : // associated volumes are refined
1306 0 : for(size_t i = 0; i < cge->num_constrained_vertices(); ++i)
1307 0 : add_hmark(cge->constrained_vertex(i), HNRM_TO_NORMAL);
1308 :
1309 0 : for(size_t i = 0; i < cge->num_constrained_edges(); ++i){
1310 : Edge* cde = cge->constrained_edge(i);
1311 0 : if(marked_refine(cde))
1312 0 : add_hmark(cde, HNRM_TO_CONSTRAINING);
1313 : else
1314 0 : add_hmark(cde, HNRM_TO_NORMAL);
1315 : }
1316 : }
1317 : }
1318 : }
1319 :
1320 : // FIXME: Also take care of (unmarkable) SHADOW_COPY edges
1321 : // whose children are marked for refinement:
1322 : // These children need to be marked HNRM_TO_CONSTRAINING!
1323 : // How can these be identified!?
1324 : // This is now only implemented in HaningNodeRefiner_MultiGrid::assign_hnode_marks
1325 :
1326 : }
1327 0 : }
1328 :
1329 : template <class TSelector>
1330 : template <class TElem>
1331 0 : void HangingNodeRefinerBase<TSelector>::
1332 : add_hmark(TElem* elem, HNodeRefMarks mark)
1333 : {
1334 : // we have to consider periodic boundaries
1335 0 : PeriodicBoundaryManager* pbm = m_pGrid->periodic_boundary_manager();
1336 : typedef typename TElem::grid_base_object base_elem_t;
1337 : base_elem_t* e = elem;
1338 0 : if(pbm && pbm->is_periodic(e)){
1339 0 : base_elem_t* master = pbm->master(e);
1340 0 : m_selMarkedElements.select(master,
1341 0 : m_selMarkedElements.get_selection_status(master) | mark);
1342 :
1343 : typedef typename PeriodicBoundaryManager::Group<base_elem_t>::SlaveContainer
1344 : slave_container_t;
1345 0 : slave_container_t* slaves = pbm->slaves(master);
1346 0 : for(typename slave_container_t::iterator i = slaves->begin();
1347 0 : i != slaves->end(); ++i)
1348 : {
1349 0 : m_selMarkedElements.select(*i,
1350 0 : m_selMarkedElements.get_selection_status(*i) | mark);
1351 : }
1352 : }
1353 : else{
1354 0 : m_selMarkedElements.select(elem,
1355 0 : m_selMarkedElements.get_selection_status(elem) | mark);
1356 : }
1357 0 : }
1358 :
1359 : template <class TSelector>
1360 : template <class TElem>
1361 0 : void HangingNodeRefinerBase<TSelector>::
1362 : remove_hmark(TElem* elem, uint mark)
1363 : {
1364 : // we have to consider periodic boundaries
1365 0 : PeriodicBoundaryManager* pbm = m_pGrid->periodic_boundary_manager();
1366 : typedef typename TElem::grid_base_object base_elem_t;
1367 : base_elem_t* e = elem;
1368 0 : if(pbm && pbm->is_periodic(e)){
1369 0 : base_elem_t* master = pbm->master(e);
1370 0 : m_selMarkedElements.select(master,
1371 0 : m_selMarkedElements.get_selection_status(master) & (~mark));
1372 :
1373 : typedef typename PeriodicBoundaryManager::Group<base_elem_t>::SlaveContainer
1374 : slave_container_t;
1375 0 : slave_container_t* slaves = pbm->slaves(master);
1376 0 : for(typename slave_container_t::iterator i = slaves->begin();
1377 0 : i != slaves->end(); ++i)
1378 : {
1379 0 : m_selMarkedElements.select(*i,
1380 0 : m_selMarkedElements.get_selection_status(*i) & (~mark));
1381 : }
1382 : }
1383 : else{
1384 0 : m_selMarkedElements.select(elem,
1385 0 : m_selMarkedElements.get_selection_status(elem) & (~mark));
1386 : }
1387 0 : }
1388 :
1389 : ////////////////////////////////////////////////////////////////////////
1390 : // implementation of refine-methods.
1391 : template <class TSelector>
1392 0 : void HangingNodeRefinerBase<TSelector>::
1393 : process_constrained_vertex(ConstrainedVertex* cdv)
1394 : {
1395 :
1396 0 : if(marked_to_normal(cdv)){
1397 : Edge* parentEdge = NULL;
1398 : Face* parentFace = NULL;
1399 : GridObject* constrObj = cdv->get_constraining_object();
1400 :
1401 0 : if(constrObj){
1402 0 : switch(cdv->get_parent_base_object_id()){
1403 : case EDGE:{
1404 0 : ConstrainingEdge* cge = dynamic_cast<ConstrainingEdge*>(constrObj);
1405 : UG_ASSERT(cge, "Constraining edge has invalid type!");
1406 : UG_ASSERT(marked_to_normal(cge),
1407 : "Constraining edge has to be converted to a normal edge!");
1408 : cge->clear_constrained_vertices();
1409 : parentEdge = cge;
1410 : UG_ASSERT(get_center_vertex(cge) == cdv,
1411 : "Center vertex and constrained vertex do not match!");
1412 : }break;
1413 : case FACE:{
1414 0 : ConstrainingFace* cgf = dynamic_cast<ConstrainingFace*>(constrObj);
1415 : UG_ASSERT(cgf, "Constraining face has invalid type!");
1416 : UG_ASSERT(marked_to_normal(cgf),
1417 : "Constraining edge has to be converted to a normal edge!");
1418 : cgf->clear_constrained_vertices();
1419 : parentFace = cgf;
1420 : UG_ASSERT((get_center_vertex(cgf) == NULL) || (get_center_vertex(cgf) == cdv),
1421 : "Center vertex and constrained vertex do not match!");
1422 : }break;
1423 : default:
1424 : break;
1425 : }
1426 : }
1427 :
1428 0 : Grid& grid = *m_pGrid;
1429 0 : Vertex* nVrt = *grid.create_and_replace<RegularVertex>(cdv);
1430 0 : if(parentEdge)
1431 0 : set_center_vertex(parentEdge, nVrt);
1432 0 : else if(parentFace)
1433 0 : set_center_vertex(parentFace, nVrt);
1434 : }
1435 0 : }
1436 :
1437 : template <class TSelector>
1438 0 : void HangingNodeRefinerBase<TSelector>::
1439 : process_constrained_edge(ConstrainedEdge* cde)
1440 : {
1441 : GridObject* constrObj = cde->get_constraining_object();
1442 0 : if(constrObj && (marked_to_normal(cde) || marked_to_constraining(cde))){
1443 0 : switch(cde->get_parent_base_object_id()){
1444 : case EDGE:{
1445 0 : ConstrainingEdge* cge = dynamic_cast<ConstrainingEdge*>(constrObj);
1446 : UG_ASSERT(cge, "Constraining edge has invalid type!");
1447 : UG_ASSERT(marked_to_normal(cge),
1448 : "Constraining edge has to be converted to a normal edge!");
1449 : cge->clear_constrained_edges();
1450 : }break;
1451 : case FACE:{
1452 0 : ConstrainingFace* cgf = dynamic_cast<ConstrainingFace*>(constrObj);
1453 : UG_ASSERT(cgf, "Constraining face has invalid type!");
1454 : UG_ASSERT(marked_to_normal(cgf),
1455 : "Constraining face has to be converted to a normal face!");
1456 : cgf->clear_constrained_edges();
1457 : }break;
1458 : default:
1459 : break;
1460 : }
1461 : }
1462 :
1463 0 : Grid& grid = *m_pGrid;
1464 0 : if(marked_to_normal(cde)){
1465 0 : grid.create_and_replace<RegularEdge>(cde);
1466 : }
1467 0 : else if(marked_to_constraining(cde)){
1468 0 : refine_edge_with_hanging_vertex(cde);
1469 : }
1470 0 : }
1471 :
1472 : template <class TSelector>
1473 0 : void HangingNodeRefinerBase<TSelector>::
1474 : process_constraining_edge(ConstrainingEdge* cge)
1475 : {
1476 : /* HNODE_PROFILE_FUNC();
1477 :
1478 : // make sure that there is one hanging vertex and two constrained edges.
1479 : assert(cge->num_constrained_vertices() == 1 && "bad number of constrained vertices. There has to be exactly 1.");
1480 : assert(cge->num_constrained_edges() == 2 && "bad number of constrained edges. There have to be exactly 2.");
1481 :
1482 : // the grid
1483 : Grid& grid = *m_pGrid;
1484 :
1485 : // the central hanging vertex has to be transformed into a normal vertex
1486 : ConstrainedVertex* centralHV = NULL;
1487 : if(cge->num_constrained_vertices() > 0)
1488 : centralHV = dynamic_cast<ConstrainedVertex*>(cge->constrained_vertex(0));
1489 :
1490 : if(!centralHV){
1491 : UG_LOG("The central hanging vertex of a constraining edge is missing. ignoring edge.\n");
1492 : return;
1493 : }
1494 :
1495 : // replace the central vertex with a normal vertex
1496 : RegularVertex* centerVrt = *grid.create_and_replace<RegularVertex>(centralHV);
1497 :
1498 : // Iterate over the constrained edges.
1499 : // Unmarked constrained edges will be replaced by a normal edge.
1500 : // Marked ones will be replaced by a ConstrainingEdge. Additionally
1501 : // associated constrained edges will be created together with the
1502 : // new central vertex.
1503 : for(size_t i = 0; i < cge->num_constrained_edges(); ++i){
1504 : Edge* cde = cge->constrained_edge(i);
1505 : if(is_marked(cde)){
1506 : refine_edge_with_hanging_vertex(cde);
1507 : }
1508 : else{
1509 : // the constrained-edge can be transformed to a normal edge
1510 : grid.create_and_replace<RegularEdge>(cde);
1511 : }
1512 : }
1513 :
1514 : cge->clear_constrained_objects();
1515 : set_center_vertex(cge, centerVrt);*/
1516 0 : }
1517 :
1518 : template <class TSelector>
1519 0 : void HangingNodeRefinerBase<TSelector>::
1520 : refine_edge_with_normal_vertex(Edge* e, Vertex** newCornerVrts)
1521 : {
1522 : UG_ASSERT(refinement_is_allowed(e), "Refinement of given edge not allowed!");
1523 :
1524 : // the grid
1525 0 : Grid& grid = *m_pGrid;
1526 :
1527 0 : if(marked_adaptive(e)){
1528 0 : if(newCornerVrts){
1529 0 : grid.create_by_cloning(e, EdgeDescriptor(newCornerVrts[0], newCornerVrts[1]), e);
1530 : }
1531 0 : return;
1532 : }
1533 :
1534 0 : RegularVertex* nVrt = *grid.create<RegularVertex>(e);
1535 0 : set_center_vertex(e, nVrt);
1536 :
1537 : // allow projector to calculate a new position
1538 0 : if(m_projector.valid())
1539 0 : m_projector->new_vertex(nVrt, e);
1540 :
1541 : // split the edge
1542 0 : vector<Edge*> vEdges(2);
1543 0 : e->refine(vEdges, nVrt, newCornerVrts);
1544 : assert((vEdges.size() == 2) && "Edge::refine - produced wrong number of edges.");
1545 0 : grid.register_element(vEdges[0], e);
1546 0 : grid.register_element(vEdges[1], e);
1547 0 : }
1548 :
1549 : template <class TSelector>
1550 0 : void HangingNodeRefinerBase<TSelector>::
1551 : refine_edge_with_hanging_vertex(Edge* e, Vertex** newCornerVrts)
1552 : {
1553 : UG_ASSERT(refinement_is_allowed(e), "Refinement of given edge not allowed!");
1554 :
1555 0 : Grid& grid = *m_pGrid;
1556 : // we have to insert a hanging node.
1557 : // e has to be transformed to a constraining edge at the same time.
1558 : assert(!ConstrainingEdge::type_match(e) && "invalid operation. e is a constraining edge.");
1559 : //assert(!ConstrainedEdge::type_match(e) && "invalid operation. e is a constrained edge.");
1560 :
1561 0 : ConstrainingEdge* ce = *grid.create_and_replace<ConstrainingEdge>(e);
1562 :
1563 0 : ConstrainedVertex* hv = *grid.create<ConstrainedVertex>(ce);
1564 :
1565 : // allow projector to calculate a new position
1566 0 : if(m_projector.valid())
1567 0 : m_projector->new_vertex(hv, ce);
1568 :
1569 0 : set_center_vertex(ce, hv);
1570 : hv->set_constraining_object(ce);
1571 0 : ce->add_constrained_object(hv);
1572 :
1573 : hv->set_local_coordinate_1(0.5);
1574 :
1575 : // two constrained edges have to be created.
1576 : ConstrainedEdge* nEdge[2];
1577 0 : if(newCornerVrts){
1578 0 : nEdge[0] = *grid.create<ConstrainedEdge>(EdgeDescriptor(newCornerVrts[0], hv), ce);
1579 0 : nEdge[1] = *grid.create<ConstrainedEdge>(EdgeDescriptor(hv, newCornerVrts[1]), ce);
1580 : }
1581 : else{
1582 0 : nEdge[0] = *grid.create<ConstrainedEdge>(EdgeDescriptor(ce->vertex(0), hv), ce);
1583 0 : nEdge[1] = *grid.create<ConstrainedEdge>(EdgeDescriptor(hv, ce->vertex(1)), ce);
1584 : }
1585 :
1586 0 : for(uint i = 0; i < 2; ++i)
1587 : {
1588 0 : ce->add_constrained_object(nEdge[i]);
1589 : nEdge[i]->set_constraining_object(ce);
1590 : }
1591 0 : }
1592 :
1593 : template <class TSelector>
1594 0 : void HangingNodeRefinerBase<TSelector>::
1595 : refine_face_with_normal_vertex(Face* f, Vertex** newCornerVrts)
1596 : {
1597 : UG_ASSERT(refinement_is_allowed(f), "Refinement of given face not allowed!");
1598 :
1599 : //UG_LOG("refine_face_with_normal_vertex\n");
1600 0 : Grid& grid = *m_pGrid;
1601 :
1602 : Vertex* vNewEdgeVertices[MAX_FACE_VERTICES];
1603 0 : vector<Face*> vFaces(f->num_vertices());// heuristic
1604 :
1605 0 : const size_t numVrts = f->num_vertices();
1606 : const size_t numEdges = f->num_edges();
1607 : bool noEdgeVrts = true;
1608 :
1609 0 : const int localMark = get_local_mark(f);
1610 0 : if(localMark && marked_local(f)){
1611 0 : for(size_t i = 0; i < numEdges; ++i){
1612 0 : if(localMark & (1<<i)){
1613 0 : vNewEdgeVertices[i] = get_center_vertex(grid.get_edge(f, i));
1614 0 : if(vNewEdgeVertices[i])
1615 : noEdgeVrts = false;
1616 : }
1617 : else
1618 0 : vNewEdgeVertices[i] = NULL;
1619 : }
1620 : }
1621 : else{
1622 0 : for(size_t i = 0; i < numEdges; ++i){
1623 0 : Edge* e = grid.get_edge(f, i);
1624 :
1625 : // if the face is refined with a regular rule, then every edge has to have
1626 : // an associated center vertex
1627 : UG_ASSERT(marked_adaptive(f) ||
1628 : (get_mark(f) == RM_REFINE && get_center_vertex(e) != NULL),
1629 : "violated for " << ElementDebugInfo(grid, f));
1630 :
1631 : // assign the center vertex
1632 0 : vNewEdgeVertices[i] = get_center_vertex(e);
1633 0 : if(vNewEdgeVertices[i])
1634 : noEdgeVrts = false;
1635 : }
1636 : }
1637 :
1638 0 : if(marked_adaptive(f) && noEdgeVrts){
1639 0 : if(newCornerVrts){
1640 0 : FaceDescriptor desc(numVrts);
1641 0 : for(size_t i = 0; i < numVrts; ++i)
1642 0 : desc.set_vertex(i, newCornerVrts[i]);
1643 0 : grid.create_by_cloning(f, desc, f);
1644 : }
1645 : return;
1646 : }
1647 :
1648 : // we'll perform a regular refine
1649 0 : Vertex* nVrt = NULL;
1650 : /*f->refine_regular(vFaces, &nVrt, vNewEdgeVertices, NULL,
1651 : RegularVertex(), newCornerVrts);*/
1652 0 : f->refine(vFaces, &nVrt, vNewEdgeVertices, NULL, newCornerVrts);
1653 :
1654 : // if a new vertex has been created during refine, then register it at the grid.
1655 0 : if(nVrt)
1656 : {
1657 : grid.register_element(nVrt, f);
1658 :
1659 : // allow projector to calculate a new position
1660 0 : if(m_projector.valid())
1661 0 : m_projector->new_vertex(nVrt, f);
1662 : }
1663 :
1664 0 : for(uint i = 0; i < vFaces.size(); ++i)
1665 : {
1666 0 : grid.register_element(vFaces[i], f);
1667 : }
1668 :
1669 0 : set_center_vertex(f, nVrt);
1670 0 : }
1671 :
1672 : template <class TSelector>
1673 0 : void HangingNodeRefinerBase<TSelector>::
1674 : refine_face_with_hanging_vertex(Face* f, Vertex** newCornerVrts)
1675 : {
1676 : UG_ASSERT(refinement_is_allowed(f), "Refinement of given face not allowed!");
1677 :
1678 : //UG_LOG("refine_face_with_hanging_vertex\n");
1679 0 : Grid& grid = *m_pGrid;
1680 :
1681 0 : size_t numVrts = f->num_vertices();
1682 : /*
1683 : vector<Edge*> vEdges(f->num_edges());
1684 : vector<Vertex*> vNewEdgeVertices(f->num_edges());
1685 : vector<Face*> vFaces(numVrts);// heuristic
1686 :
1687 : //todo: iterate over edges directly
1688 : // collect all associated edges.
1689 : CollectEdges(vEdges, grid, f);
1690 : size_t numEdges = vEdges.size();
1691 :
1692 : assert(numEdges == f->num_edges() && "ERROR in RefineFaceWithNormalVertex(...): associated edges missing.");
1693 :
1694 : // each should have an associated vertex. sort them into vNewEdgeVertices.
1695 : for(size_t i = 0; i < numEdges; ++i)
1696 : {
1697 : Edge* e = vEdges[i];
1698 : int edgeIndex = GetEdgeIndex(f, e);
1699 :
1700 : assert((edgeIndex >= 0) && (edgeIndex < (int)vEdges.size()) && "ERROR in RefineFaceWithNormalVertex(...): unknown problem in CollectEdges / GetEdgeIndex.");
1701 : //assert((get_center_vertex(e) != NULL) && "ERROR in RefineFaceWithNormalVertex(...): no new vertex on refined edge.");
1702 : vNewEdgeVertices[edgeIndex] = get_center_vertex(e);
1703 : }
1704 : */
1705 : Vertex* vNewEdgeVertices[MAX_FACE_VERTICES];
1706 0 : vector<Face*> vFaces(f->num_vertices());// heuristic
1707 :
1708 : size_t numEdges = f->num_edges();
1709 : size_t numMarkedEdges = 0;
1710 :
1711 0 : const int localMark = get_local_mark(f);
1712 0 : if(localMark && marked_local(f)){
1713 0 : for(size_t i = 0; i < numEdges; ++i){
1714 0 : if(localMark & (1<<i)){
1715 0 : vNewEdgeVertices[i] = get_center_vertex(grid.get_edge(f, i));
1716 0 : if(vNewEdgeVertices[i])
1717 0 : ++numMarkedEdges;
1718 : }
1719 : else
1720 0 : vNewEdgeVertices[i] = NULL;
1721 : }
1722 : }
1723 : else{
1724 0 : for(size_t i = 0; i < numEdges; ++i){
1725 0 : Edge* e = grid.get_edge(f, i);
1726 :
1727 : // if the face is refined with a regular rule, then every edge has to have
1728 : // an associated center vertex
1729 : assert(marked_adaptive(f) ||
1730 : ((get_mark(f) == RM_REFINE)));
1731 :
1732 : // assign the center vertex
1733 0 : vNewEdgeVertices[i] = get_center_vertex(e);
1734 0 : if(vNewEdgeVertices[i])
1735 0 : ++numMarkedEdges;
1736 : }
1737 : }
1738 :
1739 : ConstrainingFace* cgf = NULL;
1740 : ConstrainedVertex* hv = NULL;
1741 :
1742 : // the face has to be replaced by a constraining face.
1743 : // we'll perform a switch here depending on the number of vertices
1744 0 : switch(numVrts)
1745 : {
1746 0 : case 3:
1747 : {
1748 : // create the constraining triangle and replace the old face.
1749 0 : cgf = *grid.create_and_replace<ConstrainingTriangle>(f);
1750 :
1751 : // create the constrained faces.
1752 : // the following triangle will not be registered at the grid. Just a temporary one.
1753 0 : ConstrainedTriangle constrainedTri(cgf->vertex(0),
1754 0 : cgf->vertex(1),
1755 0 : cgf->vertex(2));
1756 :
1757 : // refine the constrained tri
1758 : Vertex* tmpVrt;
1759 0 : constrainedTri.refine(vFaces, &tmpVrt, vNewEdgeVertices,
1760 : NULL, newCornerVrts);
1761 : }
1762 0 : break;
1763 0 : case 4:
1764 : {
1765 0 : cgf = *grid.create_and_replace<ConstrainingQuadrilateral>(f);
1766 :
1767 : // a central hanging vertex is required
1768 0 : if(numMarkedEdges == 4){
1769 0 : hv = *grid.create<ConstrainedVertex>(cgf);
1770 :
1771 : // allow projector to calculate a new position
1772 0 : if(m_projector.valid())
1773 0 : m_projector->new_vertex(hv, cgf);
1774 :
1775 0 : set_center_vertex(cgf, hv);
1776 :
1777 : hv->set_constraining_object(cgf);
1778 0 : cgf->add_constrained_object(hv);
1779 : hv->set_local_coordinates(0.5, 0.5);
1780 : }
1781 :
1782 : // create the constrained faces.
1783 : // the following quadrilateral will not be registered at the grid. Just a temporary one.
1784 0 : ConstrainedQuadrilateral cdf(cgf->vertex(0), cgf->vertex(1),
1785 0 : cgf->vertex(2), cgf->vertex(3));
1786 :
1787 : // refine the constrained quad
1788 : Vertex* tmpVrt;
1789 0 : cdf.refine(vFaces, &tmpVrt, vNewEdgeVertices, hv, newCornerVrts);
1790 :
1791 0 : UG_COND_THROW(tmpVrt != hv, "Internal refinement error.");
1792 : }
1793 0 : break;
1794 : default:
1795 : assert(!"unsupported element type.");
1796 : break;
1797 : }
1798 :
1799 : if(cgf)
1800 : {
1801 : // register the new faces
1802 0 : for(size_t i = 0; i < vFaces.size(); ++i)
1803 : {
1804 0 : ConstrainedFace* cdf = dynamic_cast<ConstrainedFace*>(vFaces[i]);
1805 : assert(cdf && "constrained face refine did produce faces which are not constrained.");
1806 0 : if(cdf)
1807 : {
1808 : grid.register_element(cdf, cgf);
1809 : cdf->set_constraining_object(cgf);
1810 0 : cgf->add_constrained_object(cdf);
1811 : }
1812 : }
1813 :
1814 : // we have to link the new constrained edges which have been auto-generated between the constrained faces.
1815 : // Since only edges that lie inside of the constraining face are newly created, and since only those
1816 : // have to be linked with the constraining face, the following algorithm will be ok for
1817 : // triangles and quadrilaterals.
1818 : // Check for each new edge-vertex, if an edge exists with the new center vertex or with it's next neighbor.
1819 0 : if(hv)
1820 : {
1821 0 : for(size_t i = 0; i < numEdges; ++i) {
1822 0 : ConstrainedEdge* e = dynamic_cast<ConstrainedEdge*>(grid.get_edge(vNewEdgeVertices[i], hv));
1823 0 : if(e)
1824 : {
1825 : // link e with the constraining face
1826 : e->set_constraining_object(cgf);
1827 0 : cgf->add_constrained_object(e);
1828 : }
1829 : }
1830 : }
1831 0 : else if(numVrts == 3 && numMarkedEdges == 3){
1832 0 : for(size_t i = 0; i < numEdges; ++i) {
1833 : // check if a constrained edge exists between the vertex and its next neighbor
1834 0 : Vertex* vNext = vNewEdgeVertices[(i + 1) % numEdges];
1835 0 : ConstrainedEdge* e = dynamic_cast<ConstrainedEdge*>(grid.get_edge(vNewEdgeVertices[i], vNext));
1836 0 : if(e)
1837 : {
1838 : // link e with the constraining face
1839 : e->set_constraining_object(cgf);
1840 0 : cgf->add_constrained_object(e);
1841 : }
1842 : }
1843 : }
1844 : else{
1845 0 : for(size_t i = 0; i < numEdges; ++i) {
1846 0 : Vertex* vCur = vNewEdgeVertices[i];
1847 0 : if(!vCur)
1848 0 : continue;
1849 :
1850 : // check if constrained edges exist between pairs of new edge vertices
1851 : // and pairs of new edge vertices and normal vertices
1852 0 : for(size_t j = i + 1; j < numEdges; ++j){
1853 0 : Vertex* vNext = vNewEdgeVertices[j];
1854 0 : if(vNext){
1855 0 : ConstrainedEdge* e = dynamic_cast<ConstrainedEdge*>(grid.get_edge(vCur, vNext));
1856 0 : if(e)
1857 : {
1858 : // link e with the constraining face
1859 : e->set_constraining_object(cgf);
1860 0 : cgf->add_constrained_object(e);
1861 : }
1862 : }
1863 : }
1864 :
1865 0 : for(size_t j = i + 2; j < i + 2 + numVrts - 2; ++j){
1866 : Vertex* vNext;
1867 0 : if(newCornerVrts)
1868 0 : vNext = newCornerVrts[j%numVrts];
1869 : else
1870 0 : vNext = cgf->vertex(j%numVrts);
1871 :
1872 0 : ConstrainedEdge* e = dynamic_cast<ConstrainedEdge*>(grid.get_edge(vCur, vNext));
1873 0 : if(e)
1874 : {
1875 : // link e with the constraining face
1876 : e->set_constraining_object(cgf);
1877 0 : cgf->add_constrained_object(e);
1878 : }
1879 : }
1880 : }
1881 : }
1882 : }
1883 0 : }
1884 :
1885 : template <class TSelector>
1886 0 : void HangingNodeRefinerBase<TSelector>::
1887 : process_constrained_face(ConstrainedFace* cdf)
1888 : {
1889 : GridObject* constrObj = cdf->get_constraining_object();
1890 :
1891 0 : if(constrObj && (marked_to_normal(cdf) || marked_to_constraining(cdf))){
1892 0 : switch(cdf->get_parent_base_object_id()){
1893 : case FACE:{
1894 0 : ConstrainingFace* cgf = dynamic_cast<ConstrainingFace*>(constrObj);
1895 : UG_ASSERT(cgf, "Constraining face has invalid type!");
1896 : UG_ASSERT(marked_to_normal(cgf),
1897 : "Constraining face has to be converted to a normal face!");
1898 : cgf->clear_constrained_faces();
1899 : }break;
1900 : default:
1901 : UG_ASSERT(cdf->get_constraining_object() == NULL,
1902 : "Invalid constraining object encountered!");
1903 : break;
1904 : }
1905 : }
1906 :
1907 0 : Grid& grid = *m_pGrid;
1908 0 : if(marked_to_normal(cdf)){
1909 0 : if(cdf->num_vertices() == 3)
1910 0 : grid.create_and_replace<Triangle>(cdf);
1911 : else
1912 0 : grid.create_and_replace<Quadrilateral>(cdf);
1913 : }
1914 0 : else if(marked_to_constraining(cdf)){
1915 0 : refine_face_with_hanging_vertex(cdf);
1916 : }
1917 0 : }
1918 :
1919 : template <class TSelector>
1920 0 : void HangingNodeRefinerBase<TSelector>::
1921 : process_constraining_face(ConstrainingFace* cgf)
1922 : {
1923 0 : Grid& grid = *m_pGrid;
1924 : /*
1925 : size_t numVrts = cgf->num_vertices();
1926 :
1927 : // make sure that there is one hanging vertex and two constrained edges.
1928 : UG_ASSERT(cgf->num_constrained_edges() == numVrts,
1929 : "bad number of constrained edges: " << cgf->num_constrained_edges()
1930 : << ". There have to be as many as vertices: " << numVrts << "."
1931 : << "At face with center " << GetGridObjectCenter(grid, cgf));
1932 : UG_ASSERT(cgf->num_constrained_faces() == 4,
1933 : "bad number of constrained faces. There have to be exactly 4. "
1934 : << "At face with center " << GetGridObjectCenter(grid, cgf));
1935 :
1936 : ConstrainedVertex* centralHV = NULL;
1937 : RegularVertex* centerVrt = NULL;
1938 :
1939 : if(numVrts == 4){
1940 : // the central hanging vertex has to be transformed into a normal vertex
1941 : centralHV = NULL;
1942 : if(cgf->num_constrained_vertices() > 0)
1943 : centralHV = dynamic_cast<ConstrainedVertex*>(cgf->constrained_vertex(0));
1944 :
1945 : // replace the central vertex with a normal vertex
1946 : if(centralHV)
1947 : centerVrt = *grid.create_and_replace<RegularVertex>(centralHV);
1948 : }
1949 :
1950 : // Iterate over the constrained edges.
1951 : // Unmarked constrained edges will be replaced by a normal edge.
1952 : // Marked ones will be replaced by a ConstrainingEdge. Additionally
1953 : // associated constrained edges will be created together with the
1954 : // new central vertex.
1955 : for(size_t i = 0; i < cgf->num_constrained_edges(); ++i){
1956 : Edge* cde = cgf->constrained_edge(i);
1957 : if(is_marked(cde)){
1958 : refine_edge_with_hanging_vertex(cde);
1959 : }
1960 : else{
1961 : // the constrained-edge can be transformed to a normal edge
1962 : grid.create_and_replace<RegularEdge>(cde);
1963 : }
1964 : }
1965 :
1966 : // iterate over the constrained faces.
1967 : // If it is marked, we'll replace it by a constraining face and create
1968 : // associated constrained faces.
1969 : // if not, it will simply be transformed to a normal face.
1970 : // To ease implementation we will transform it anyway and if required we will
1971 : // call refine_face_with_hanging_vertex(...).
1972 : for(size_t i = 0; i < cgf->num_constrained_faces(); ++i){
1973 : Face* f = cgf->constrained_face(i);
1974 : if(is_marked(f)){
1975 : // refine it using hanging_node_refinement.
1976 : refine_face_with_hanging_vertex(f);
1977 : }
1978 : else{
1979 : // replace it by a normal face
1980 : if(f->num_vertices() == 3)
1981 : f = *grid.create_and_replace<Triangle>(f);
1982 : else
1983 : f = *grid.create_and_replace<Quadrilateral>(f);
1984 : }
1985 : }
1986 : */
1987 : // cgf->clear_constrained_objects();
1988 : // cgf itself now has to be transformed to a normal face
1989 : UG_ASSERT(marked_to_normal(cgf), "A constraining face has to be converted to"
1990 : " a normal face when it is refined.");
1991 0 : Vertex* centerVrt = get_center_vertex(cgf);
1992 : Face* nFace;
1993 0 : if(cgf->num_vertices() == 3)
1994 0 : nFace = *grid.create_and_replace<Triangle>(cgf);
1995 : else
1996 0 : nFace = *grid.create_and_replace<Quadrilateral>(cgf);
1997 :
1998 0 : if(centerVrt)
1999 0 : set_center_vertex(nFace, centerVrt);
2000 0 : }
2001 :
2002 : template <class TSelector>
2003 0 : void HangingNodeRefinerBase<TSelector>::
2004 : refine_volume_with_normal_vertex(Volume* v, Vertex** newCornerVrts)
2005 : {
2006 : UG_ASSERT(refinement_is_allowed(v), "Refinement of given volume not allowed!");
2007 :
2008 0 : Grid& grid = *m_pGrid;
2009 :
2010 : //vector<Edge*> vEdges(v->num_edges());
2011 0 : vector<Vertex*> vNewEdgeVertices(v->num_edges());
2012 : //vector<Face*> vFaces(v->num_faces());
2013 0 : vector<Vertex*> vNewFaceVertices(v->num_faces());
2014 0 : vector<Volume*> vVolumes(8);// heuristic
2015 : // collect all associated edges.
2016 :
2017 0 : const size_t numEdges = v->num_edges();
2018 : bool noEdgeVrts = true;
2019 :
2020 :
2021 0 : const int localMark = get_local_mark(v);
2022 0 : if(localMark && marked_local(v)){
2023 0 : for(size_t i = 0; i < numEdges; ++i){
2024 0 : if(localMark & (1<<i)){
2025 0 : vNewEdgeVertices[i] = get_center_vertex(grid.get_edge(v, i));
2026 0 : if(vNewEdgeVertices[i])
2027 : noEdgeVrts = false;
2028 : }
2029 : else
2030 0 : vNewEdgeVertices[i] = NULL;
2031 : }
2032 : }
2033 : else{
2034 0 : for(size_t i = 0; i < numEdges; ++i){
2035 0 : Edge* e = grid.get_edge(v, i);
2036 0 : vNewEdgeVertices[i] = get_center_vertex(e);
2037 0 : if(vNewEdgeVertices[i])
2038 : noEdgeVrts = false;
2039 0 : UG_COND_THROW(!(marked_adaptive(v) || vNewEdgeVertices[i]),
2040 : "In order to fully refine a volume, all edges have "
2041 : "to contain a new vertex!");
2042 : }
2043 : }
2044 :
2045 0 : if(marked_adaptive(v) && noEdgeVrts){
2046 0 : if(newCornerVrts){
2047 0 : size_t numVrts = v->num_vertices();
2048 0 : VolumeDescriptor desc(numVrts);
2049 0 : for(size_t i = 0; i < numVrts; ++i)
2050 0 : desc.set_vertex(i, newCornerVrts[i]);
2051 0 : grid.create_by_cloning(v, desc, v);
2052 : }
2053 : return;
2054 : }
2055 :
2056 0 : size_t numFaces = v->num_faces();
2057 0 : for(size_t i = 0; i < numFaces; ++i){
2058 0 : Face* f = grid.get_face(v, i);
2059 :
2060 : /*if(!VolumeContains(v, f))
2061 : {
2062 : UG_LOG("Grid::get_face(vol, ind) returned bad face.");
2063 : MultiGrid* pmg = dynamic_cast<MultiGrid*>(m_pGrid);
2064 : UG_LOG("Vol in level " << pmg->get_level(v));
2065 : UG_LOG(", face in level " << pmg->get_level(f) << endl);
2066 : UG_LOG("positions of volume vertices:");
2067 : for(size_t i_c = 0; i_c < v->num_vertices(); ++i_c){
2068 : UG_LOG(" " << GetGridObjectCenter(grid, v->vertex(i_c)));
2069 : }
2070 : UG_LOG("\npositions of face vertices:");
2071 : for(size_t i_c = 0; i_c < f->num_vertices(); ++i_c){
2072 : UG_LOG(" " << GetGridObjectCenter(grid, f->vertex(i_c)));
2073 : }
2074 : UG_LOG(endl);
2075 : }*/
2076 :
2077 0 : if(f->num_vertices() == 3)
2078 0 : vNewFaceVertices[i] = NULL;
2079 : else{
2080 0 : vNewFaceVertices[i] = get_center_vertex(f);
2081 0 : UG_COND_THROW(!(marked_adaptive(v) || vNewFaceVertices[i]),
2082 : "In order to fully refine a volume, all quadrilateral sides have"
2083 : " to contain a new vertex!"
2084 : << ", vol-center " << GetGridObjectCenter(grid, v)
2085 : << ", face-center " << GetGridObjectCenter(grid, f));
2086 : //todo:remove this!!!
2087 : /*if(!vNewFaceVertices[i]){
2088 : UG_LOG("missing face-vertex: " << GetGridObjectCenter(grid, f) << endl);
2089 : UG_LOG("corners of face:");
2090 : for(size_t i_c = 0; i_c < f->num_vertices(); ++i_c){
2091 : UG_LOG(" " << GetGridObjectCenter(grid, f->vertex(i_c)));
2092 : }
2093 : UG_LOG(endl);
2094 : UG_LOG("during refinement of volume: " << GetGridObjectCenter(grid, v) << endl);
2095 : }*/
2096 : }
2097 : }
2098 :
2099 : // if we're performing tetrahedral refinement, we have to collect
2100 : // the corner coordinates, so that the refinement algorithm may choose
2101 : // the best interior diagonal.
2102 0 : vector3 corners[4];
2103 : vector3* pCorners = NULL;
2104 0 : if((v->num_vertices() == 4) && m_projector.valid()){
2105 0 : for(size_t i = 0; i < 4; ++i){
2106 0 : corners[i] = m_projector->geometry()->pos(v->vertex(i));
2107 : }
2108 : pCorners = corners;
2109 : }
2110 :
2111 : // refine the volume and register new volumes at the grid.
2112 0 : Vertex* createdVrt = NULL;
2113 0 : v->refine(vVolumes, &createdVrt, &vNewEdgeVertices.front(),
2114 0 : &vNewFaceVertices.front(), NULL, RegularVertex(), newCornerVrts, pCorners);
2115 :
2116 0 : if(createdVrt){
2117 : // register the new vertex
2118 : grid.register_element(createdVrt, v);
2119 :
2120 : // allow projector to calculate a new position
2121 0 : if(m_projector.valid())
2122 0 : m_projector->new_vertex(createdVrt, v);
2123 : }
2124 :
2125 0 : for(uint i = 0; i < vVolumes.size(); ++i)
2126 0 : grid.register_element(vVolumes[i], v);
2127 0 : }
2128 :
2129 : template class HangingNodeRefinerBase<Selector>;
2130 : template void HangingNodeRefinerBase<Selector>::add_hmark(Vertex*, HNodeRefMarks);
2131 : template void HangingNodeRefinerBase<Selector>::add_hmark(Edge*, HNodeRefMarks);
2132 : template void HangingNodeRefinerBase<Selector>::add_hmark(Face*, HNodeRefMarks);
2133 : template void HangingNodeRefinerBase<Selector>::add_hmark(Volume*, HNodeRefMarks);
2134 :
2135 : template class HangingNodeRefinerBase<MGSelector>;
2136 : template void HangingNodeRefinerBase<MGSelector>::add_hmark(Vertex*, HNodeRefMarks);
2137 : template void HangingNodeRefinerBase<MGSelector>::add_hmark(Edge*, HNodeRefMarks);
2138 : template void HangingNodeRefinerBase<MGSelector>::add_hmark(Face*, HNodeRefMarks);
2139 : template void HangingNodeRefinerBase<MGSelector>::add_hmark(Volume*, HNodeRefMarks);
2140 :
2141 : }// end of namespace
|
