Line data Source code
1 : /*
2 : * Copyright (c) 2012-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #include "lib_grid/algorithms/selection_util.h"
34 : #include "adaptive_regular_mg_refiner.h"
35 :
36 : using namespace std;
37 :
38 : namespace ug{
39 :
40 0 : AdaptiveRegularRefiner_MultiGrid::
41 0 : AdaptiveRegularRefiner_MultiGrid(SPRefinementProjector projector) :
42 0 : HangingNodeRefiner_MultiGrid(projector)
43 : {
44 0 : }
45 :
46 0 : AdaptiveRegularRefiner_MultiGrid::
47 : AdaptiveRegularRefiner_MultiGrid(MultiGrid& mg,
48 0 : SPRefinementProjector projector) :
49 0 : HangingNodeRefiner_MultiGrid(projector)
50 : {
51 0 : set_grid(&mg);
52 0 : }
53 :
54 0 : AdaptiveRegularRefiner_MultiGrid::
55 0 : ~AdaptiveRegularRefiner_MultiGrid()
56 : {
57 0 : }
58 :
59 0 : void AdaptiveRegularRefiner_MultiGrid::
60 : assign_grid(MultiGrid& mg)
61 : {
62 0 : set_grid(&mg);
63 0 : }
64 :
65 0 : void AdaptiveRegularRefiner_MultiGrid::
66 : set_grid(MultiGrid* mg)
67 : {
68 0 : HangingNodeRefiner_MultiGrid::set_grid(mg);
69 :
70 0 : m_closureElems.enable_autoselection(false);
71 0 : m_closureElems.enable_selection_inheritance(false);
72 0 : m_closureElems.enable_strict_inheritance(false);
73 :
74 0 : m_closureElems.assign_grid(mg);
75 0 : }
76 :
77 :
78 0 : void AdaptiveRegularRefiner_MultiGrid::
79 : remove_closure_elements()
80 : {
81 : // remove all closure elements from the current grid
82 :
83 0 : if(!multi_grid()){
84 0 : UG_THROW("AdaptiveRegularRefiner_MultiGrid has to be associated with a grid!");
85 : }
86 :
87 0 : EraseSelectedObjects(m_closureElems);
88 0 : }
89 :
90 :
91 0 : void AdaptiveRegularRefiner_MultiGrid::
92 : create_closure_elements()
93 : {
94 0 : if(!multi_grid()){
95 0 : UG_THROW("AdaptiveRegularRefiner_MultiGrid has to be associated with a grid!");
96 : }
97 :
98 0 : MultiGrid& mg = *multi_grid();
99 :
100 0 : if(mg.num<Volume>() > 0)
101 0 : create_closure_elements_3d();
102 0 : else if(mg.num<Face>() > 0)
103 0 : create_closure_elements_2d();
104 0 : }
105 :
106 :
107 0 : void AdaptiveRegularRefiner_MultiGrid::
108 : create_closure_elements_2d()
109 : {
110 : //todo: This method currently only works for one type of elements. No manifolds
111 : // are currently supported.
112 :
113 : // for each surface element (faces in 2d, volumes in 3d) adjacent to a constraining
114 : // element, we generate closure elements in the level above.
115 :
116 0 : if(!multi_grid()){
117 0 : UG_THROW("AdaptiveRegularRefiner_MultiGrid has to be associated with a grid!");
118 : }
119 :
120 : // remove all closure elements. This is currently required, since the selector
121 : // m_closureElems is currently also used to store non-closer elements, which are
122 : // to be refined.
123 0 : remove_closure_elements();
124 :
125 0 : MultiGrid& mg = *multi_grid();
126 :
127 : // iterate over all constraining edges and collect associated surface faces
128 : // Those then have to be refined to generate a closure.
129 :
130 : Grid::face_traits::secure_container assElems;
131 : Grid::edge_traits::secure_container assEdges;
132 : Face::ConstVertexArray vrts;
133 : std::vector<Vertex*> newVrtVrts;
134 : std::vector<Vertex*> newEdgeVrts;
135 : std::vector<Face*> newFaces;
136 0 : EdgeDescriptor ed;
137 :
138 : // we'll select all new elements on the fly
139 0 : m_closureElems.enable_autoselection(true);
140 :
141 : for(Grid::traits<ConstrainingEdge>::iterator i_edge = mg.begin<ConstrainingEdge>();
142 0 : i_edge != mg.end<ConstrainingEdge>(); ++i_edge)
143 : {
144 : // check all associated elements of i_edge, whether one is a surface element.
145 : // If so, create the closure.
146 0 : mg.associated_elements(assElems, *i_edge);
147 :
148 0 : for(size_t i_elem = 0; i_elem < assElems.size(); ++i_elem){
149 : Face* elem = assElems[i_elem];
150 0 : if(mg.has_children(elem))
151 0 : continue;
152 :
153 : newVrtVrts.clear();
154 : newEdgeVrts.clear();
155 :
156 : // copy associated vertices and edges to the next level
157 0 : vrts = elem->vertices();
158 0 : size_t numVrts = elem->num_vertices();
159 0 : for(size_t i_vrt = 0; i_vrt < numVrts; ++i_vrt){
160 0 : Vertex* vrt = vrts[i_vrt];
161 0 : if(!mg.has_children(vrt)){
162 0 : newVrtVrts.push_back(*mg.create<RegularVertex>(vrt));
163 0 : if(m_projector.valid())
164 0 : m_projector->new_vertex(newVrtVrts.back(), vrt);
165 : }
166 : else
167 0 : newVrtVrts.push_back(mg.get_child_vertex(vrt));
168 : }
169 :
170 : mg.associated_elements_sorted(assEdges, elem);
171 0 : for(size_t i = 0; i < assEdges.size(); ++i){
172 : Edge* e = assEdges[i];
173 : if(!mg.has_children(e)){
174 0 : ed.set_vertices(mg.get_child_vertex(e->vertex(0)),
175 0 : mg.get_child_vertex(e->vertex(1)));
176 0 : mg.create<RegularEdge>(ed, e);
177 0 : newEdgeVrts.push_back(NULL);
178 : }
179 : else
180 0 : newEdgeVrts.push_back(mg.get_child_vertex(e));
181 : }
182 :
183 : // refine the element
184 : Vertex* newFaceVrt;
185 0 : if(elem->refine(newFaces, &newFaceVrt, &newEdgeVrts.front(),
186 : NULL, &newVrtVrts.front()))
187 : {
188 0 : if(newFaceVrt){
189 : mg.register_element(newFaceVrt, elem);
190 0 : if(m_projector.valid())
191 0 : m_projector->new_vertex(newFaceVrt, elem);
192 : }
193 0 : for(size_t i = 0; i < newFaces.size(); ++i)
194 0 : mg.register_element(newFaces[i], elem);
195 : }
196 : }
197 : }
198 :
199 : // stop selecting new elements
200 0 : m_closureElems.enable_autoselection(false);
201 0 : }
202 :
203 0 : void AdaptiveRegularRefiner_MultiGrid::
204 : create_closure_elements_3d()
205 : {
206 : //todo: This method currently only works for one type of elements. No manifolds
207 : // are currently supported.
208 :
209 : // for each surface element (faces in 2d, volumes in 3d) adjacent to a constraining
210 : // element, we generate closure elements in the level above.
211 :
212 0 : if(!multi_grid()){
213 0 : UG_THROW("AdaptiveRegularRefiner_MultiGrid has to be associated with a grid!");
214 : }
215 :
216 : // remove all closure elements. This is currently required, since the selector
217 : // m_closureElems is currently also used to store non-closer elements, which are
218 : // to be refined.
219 0 : remove_closure_elements();
220 :
221 0 : MultiGrid& mg = *multi_grid();
222 :
223 : // iterate over all constraining edges and collect associated surface faces / volumes.
224 : // Those then have to be refined to generate a closure.
225 :
226 : Grid::volume_traits::secure_container assElems;
227 : Grid::edge_traits::secure_container assEdges;
228 : Grid::face_traits::secure_container assFaces;
229 : Volume::ConstVertexArray vrts;
230 : std::vector<Vertex*> newVolVrtVrts;
231 : std::vector<Vertex*> newVolEdgeVrts;
232 : std::vector<Vertex*> newVolFaceVrts;
233 : std::vector<Volume*> newVols;
234 0 : EdgeDescriptor ed;
235 0 : FaceDescriptor fd;
236 :
237 : // when refining the associated faces, we need some structs, too
238 : Grid::edge_traits::secure_container assFaceEdges;
239 : std::vector<Vertex*> newFaceVrtVrts;
240 : std::vector<Vertex*> newFaceEdgeVrts;
241 : std::vector<Face*> newFaces;
242 :
243 : // we'll select all new elements on the fly
244 0 : m_closureElems.enable_autoselection(true);
245 :
246 : for(Grid::traits<ConstrainingEdge>::iterator i_edge = mg.begin<ConstrainingEdge>();
247 0 : i_edge != mg.end<ConstrainingEdge>(); ++i_edge)
248 : {
249 : // check all associated elements of i_edge, whether one is a surface element.
250 : // If so, select it into m_closureElems. This is only temporary, since it
251 : // isn't a real closure element.
252 0 : mg.associated_elements(assElems, *i_edge);
253 :
254 0 : for(size_t i_elem = 0; i_elem < assElems.size(); ++i_elem){
255 : Volume* elem = assElems[i_elem];
256 0 : if(mg.has_children(elem))
257 0 : continue;
258 :
259 : newVolVrtVrts.clear();
260 : newVolEdgeVrts.clear();
261 : newVolFaceVrts.clear();
262 :
263 : // copy associated vertices and edges to the next level
264 0 : vrts = elem->vertices();
265 0 : size_t numVrts = elem->num_vertices();
266 0 : for(size_t i_vrt = 0; i_vrt < numVrts; ++i_vrt){
267 0 : Vertex* vrt = vrts[i_vrt];
268 0 : if(!mg.has_children(vrt)){
269 0 : newVolVrtVrts.push_back(*mg.create<RegularVertex>(vrt));
270 0 : if(m_projector.valid())
271 0 : m_projector->new_vertex(newVolVrtVrts.back(), vrt);
272 : }
273 : else
274 0 : newVolVrtVrts.push_back(mg.get_child_vertex(vrt));
275 : }
276 :
277 : mg.associated_elements_sorted(assEdges, elem);
278 0 : for(size_t i = 0; i < assEdges.size(); ++i){
279 : Edge* e = assEdges[i];
280 : if(!mg.has_children(e)){
281 0 : ed.set_vertices(mg.get_child_vertex(e->vertex(0)),
282 0 : mg.get_child_vertex(e->vertex(1)));
283 0 : mg.create<RegularEdge>(ed, e);
284 0 : newVolEdgeVrts.push_back(NULL);
285 : }
286 : else
287 0 : newVolEdgeVrts.push_back(mg.get_child_vertex(e));
288 : }
289 :
290 : // we have to either refine or copy associated faces
291 : mg.associated_elements_sorted(assFaces, elem);
292 0 : for(size_t i_face = 0; i_face < assFaces.size(); ++i_face){
293 :
294 : Face* f = assFaces[i_face];
295 0 : if(mg.has_children(f)){
296 0 : newVolFaceVrts.push_back(mg.get_child_vertex(f));
297 0 : continue;
298 : }
299 :
300 : // collect child vertices of associated edges
301 : mg.associated_elements_sorted(assFaceEdges, f);
302 : newFaceEdgeVrts.clear();
303 :
304 : bool faceRefinement = false;
305 0 : for(size_t i = 0; i < assFaceEdges.size(); ++i){
306 0 : Vertex* child = mg.get_child_vertex(assFaceEdges[i]);
307 0 : newFaceEdgeVrts.push_back(child);
308 0 : faceRefinement |= (child != NULL);
309 : }
310 :
311 0 : if(faceRefinement){
312 : newFaceVrtVrts.clear();
313 0 : for(size_t i = 0; i < f->num_vertices(); ++i)
314 0 : newFaceVrtVrts.push_back(mg.get_child_vertex(f->vertex(i)));
315 :
316 0 : Vertex* newFaceVrt = NULL;
317 0 : if(f->refine(newFaces, &newFaceVrt, &newFaceEdgeVrts.front(),
318 : NULL, &newFaceVrtVrts.front()))
319 : {
320 0 : if(newFaceVrt){
321 : mg.register_element(newFaceVrt, f);
322 0 : if(m_projector.valid())
323 0 : m_projector->new_vertex(newFaceVrt, f);
324 : }
325 0 : for(size_t i = 0; i < newFaces.size(); ++i)
326 0 : mg.register_element(newFaces[i], f);
327 : }
328 0 : newVolFaceVrts.push_back(newFaceVrt);
329 : }
330 : else{
331 0 : Face::ConstVertexArray fvrts = f->vertices();
332 0 : if(f->num_vertices() == 3)
333 0 : mg.create<Triangle>(TriangleDescriptor(
334 : mg.get_child_vertex(fvrts[0]),
335 : mg.get_child_vertex(fvrts[1]),
336 : mg.get_child_vertex(fvrts[2])),
337 : f);
338 0 : else if(f->num_vertices() == 4)
339 0 : mg.create<Quadrilateral>(QuadrilateralDescriptor(
340 : mg.get_child_vertex(fvrts[0]),
341 : mg.get_child_vertex(fvrts[1]),
342 : mg.get_child_vertex(fvrts[2]),
343 : mg.get_child_vertex(fvrts[3])),
344 : f);
345 0 : newVolFaceVrts.push_back(NULL);
346 : }
347 : }
348 :
349 : // refine the element
350 : // if we're performing tetrahedral refinement, we have to collect
351 : // the corner coordinates, so that the refinement algorithm may choose
352 : // the best interior diagonal.
353 0 : vector3 corners[4];
354 : vector3* pCorners = NULL;
355 0 : if((elem->num_vertices() == 4) && m_projector.valid()){
356 0 : for(size_t i = 0; i < 4; ++i){
357 0 : corners[i] = m_projector->geometry()->pos(vrts[i]);
358 : }
359 : pCorners = corners;
360 : }
361 :
362 : Vertex* newVolVrt;
363 0 : if(elem->refine(newVols, &newVolVrt, &newVolEdgeVrts.front(),
364 0 : &newVolFaceVrts.front(), NULL, RegularVertex(),
365 : &newVolVrtVrts.front(), pCorners))
366 : {
367 0 : if(newVolVrt){
368 : mg.register_element(newVolVrt, elem);
369 0 : if(m_projector.valid())
370 0 : m_projector->new_vertex(newVolVrt, elem);
371 : }
372 :
373 0 : for(size_t i = 0; i < newVols.size(); ++i)
374 0 : mg.register_element(newVols[i], elem);
375 : }
376 : }
377 : }
378 :
379 : // stop selecting new elements
380 0 : m_closureElems.enable_autoselection(false);
381 0 : }
382 :
383 : template <class TElem>
384 0 : void AdaptiveRegularRefiner_MultiGrid::
385 : get_parents_of_marked_closure_elements(std::vector<GridObject*>& parents,
386 : Selector::status_t mark)
387 : {
388 : UG_ASSERT(multi_grid(), "A multi grid has to be assigned to the refiner.");
389 0 : MultiGrid& mg = *multi_grid();
390 :
391 : typedef typename BaseClass::selector_t::template traits<TElem>::iterator TIter;
392 0 : for(TIter iter = m_selMarkedElements.begin<TElem>();
393 0 : iter != m_selMarkedElements.end<TElem>(); ++iter)
394 : {
395 : TElem* e = *iter;
396 0 : if(!m_closureElems.is_selected(e))
397 0 : continue;
398 :
399 0 : if(get_mark(e) & mark){
400 0 : GridObject* parent = mg.get_parent(e);
401 0 : if(parent && !m_closureElems.is_selected(parent))
402 0 : parents.push_back(parent);
403 : }
404 : }
405 0 : }
406 :
407 0 : void AdaptiveRegularRefiner_MultiGrid::
408 : perform_refinement()
409 : {
410 : // todo: copy refinement marks from closure elements to their parents
411 : vector<GridObject*> parents;
412 : Selector::status_t refMark = RM_REFINE | RM_ANISOTROPIC;
413 0 : get_parents_of_marked_closure_elements<Vertex>(parents, refMark);
414 0 : get_parents_of_marked_closure_elements<Edge>(parents, refMark);
415 0 : get_parents_of_marked_closure_elements<Face>(parents, refMark);
416 0 : get_parents_of_marked_closure_elements<Volume>(parents, refMark);
417 :
418 0 : remove_closure_elements();
419 :
420 : // mark parents of formerly marked closure elements for refinement
421 0 : mark(parents.begin(), parents.end(), RM_REFINE);
422 :
423 0 : HangingNodeRefiner_MultiGrid::perform_refinement();
424 0 : create_closure_elements();
425 0 : }
426 :
427 0 : bool AdaptiveRegularRefiner_MultiGrid::
428 : perform_coarsening()
429 : {
430 : // todo: copy coarsen marks from closure elements to their parents
431 0 : remove_closure_elements();
432 0 : bool bSuccess = HangingNodeRefiner_MultiGrid::perform_coarsening();
433 0 : create_closure_elements();
434 0 : return bSuccess;
435 : }
436 :
437 : }// end of namespace
|
