Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #include <cassert>
34 : #include "common/profiler/profiler.h"
35 : #include "global_multi_grid_refiner.h"
36 : #include "lib_grid/algorithms/algorithms.h"
37 : #include "lib_grid/file_io/file_io.h"
38 :
39 : //define PROFILE_GLOBAL_MULTI_GRID_REFINER if you want to profile
40 : //the refinement code.
41 : #define PROFILE_GLOBAL_MULTI_GRID_REFINER
42 : #ifdef PROFILE_GLOBAL_MULTI_GRID_REFINER
43 : #define GMGR_PROFILE_FUNC() PROFILE_FUNC_GROUP("grid")
44 : #define GMGR_PROFILE(name) PROFILE_BEGIN_GROUP(name, "grid")
45 : #define GMGR_PROFILE_END() PROFILE_END()
46 : #else
47 : #define GMGR_PROFILE_FUNC()
48 : #define GMGR_PROFILE(name)
49 : #define GMGR_PROFILE_END()
50 : #endif
51 :
52 : using namespace std;
53 :
54 : namespace ug
55 : {
56 :
57 0 : GlobalMultiGridRefiner::
58 0 : GlobalMultiGridRefiner(SPRefinementProjector projector) :
59 : IRefiner(projector),
60 0 : m_pMG(NULL)
61 : {
62 0 : }
63 :
64 0 : GlobalMultiGridRefiner::
65 0 : GlobalMultiGridRefiner(MultiGrid& mg, SPRefinementProjector projector) :
66 0 : IRefiner(projector)
67 : {
68 0 : m_pMG = NULL;
69 0 : assign_grid(mg);
70 0 : }
71 :
72 0 : GlobalMultiGridRefiner::~GlobalMultiGridRefiner()
73 : {
74 0 : if(m_pMG)
75 0 : m_pMG->unregister_observer(this);
76 0 : }
77 :
78 0 : void GlobalMultiGridRefiner::grid_to_be_destroyed(Grid* grid)
79 : {
80 0 : m_pMG = NULL;
81 0 : }
82 :
83 0 : void GlobalMultiGridRefiner::assign_grid(MultiGrid& mg)
84 : {
85 0 : assign_grid(&mg);
86 0 : }
87 :
88 0 : void GlobalMultiGridRefiner::assign_grid(MultiGrid* mg)
89 : {
90 0 : if(m_pMG){
91 0 : m_pMG->unregister_observer(this);
92 0 : m_pMG = NULL;
93 : }
94 :
95 0 : if(mg){
96 0 : m_pMG = mg;
97 0 : set_message_hub(mg->message_hub());
98 0 : m_pMG->register_observer(this, OT_GRID_OBSERVER);
99 : }
100 0 : }
101 :
102 :
103 0 : void GlobalMultiGridRefiner::
104 : num_marked_edges_local(std::vector<int>& numMarkedEdgesOut)
105 : {
106 0 : num_marked_elems<Edge>(numMarkedEdgesOut);
107 0 : }
108 :
109 0 : void GlobalMultiGridRefiner::
110 : num_marked_faces_local(std::vector<int>& numMarkedFacesOut)
111 : {
112 0 : num_marked_elems<Face>(numMarkedFacesOut);
113 0 : }
114 :
115 0 : void GlobalMultiGridRefiner::
116 : num_marked_volumes_local(std::vector<int>& numMarkedVolsOut)
117 : {
118 0 : num_marked_elems<Volume>(numMarkedVolsOut);
119 0 : }
120 :
121 :
122 : template <class TElem>
123 0 : void GlobalMultiGridRefiner::
124 : num_marked_elems(std::vector<int>& numMarkedElemsOut)
125 : {
126 : numMarkedElemsOut.clear();
127 0 : if(!m_pMG)
128 : return;
129 0 : numMarkedElemsOut.resize(m_pMG->num_levels(), 0);
130 0 : if(m_pMG->num_levels() > 0)
131 0 : numMarkedElemsOut.back() = m_pMG->num<TElem>(m_pMG->top_level());
132 : }
133 :
134 : ////////////////////////////////////////////////////////////////////////
135 0 : void GlobalMultiGridRefiner::perform_refinement()
136 : {
137 : UG_DLOG(LIB_GRID, 1, "GlobalMultiGridRefiner\n");
138 :
139 : GMGR_PROFILE_FUNC();
140 :
141 : assert(m_pMG && "refiner has to be assigned to a multi-grid!");
142 0 : if(!m_pMG)
143 0 : return;
144 :
145 : // the multi-grid
146 : MultiGrid& mg = *m_pMG;
147 :
148 : // make sure that the required options are enabled.
149 0 : if(mg.num_volumes() > 0){
150 0 : if(!mg.option_is_enabled(VOLOPT_AUTOGENERATE_FACES))
151 : {
152 : LOG("WARNING in GlobalMultiGridRefiner::refine(): auto-enabling VOLOPT_AUTOGENERATE_FACES.\n");
153 0 : mg.enable_options(VOLOPT_AUTOGENERATE_FACES);
154 : }
155 : }
156 :
157 0 : if(mg.num_faces() > 0){
158 0 : if(!mg.option_is_enabled(FACEOPT_AUTOGENERATE_EDGES))
159 : {
160 : LOG("WARNING in GlobalMultiGridRefiner::refine(): auto-enabling FACEOPT_AUTOGENERATE_EDGES.\n");
161 0 : mg.enable_options(FACEOPT_AUTOGENERATE_EDGES);
162 : }
163 : }
164 :
165 0 : if(mg.num_levels() == 0)
166 : return;
167 :
168 : // the old top level
169 0 : int oldTopLevel = mg.num_levels() - 1;
170 : //todo: Ghosts have to be ignored in the specified grid-object-collection - otherwise problems may occur
171 : // e.g. in the AdaptionSurfaceGridFunction during prolongation of piecewise const functions... (M. Breit)
172 0 : m_messageHub->post_message(GridMessage_Adaption(GMAT_GLOBAL_REFINEMENT_BEGINS,
173 0 : mg.get_grid_objects(oldTopLevel)));
174 :
175 0 : if(projector()->refinement_begins_requires_subgrid()){
176 0 : SubGrid<ConsiderAll> sg(mg.get_grid_objects(), ConsiderAll());
177 0 : projector()->refinement_begins(&sg);
178 : }
179 : else
180 0 : projector()->refinement_begins(NULL);
181 :
182 : UG_DLOG(LIB_GRID, 1, "REFINER: reserving memory...");
183 :
184 : //todo: Adjust for parallel computations with ghost elements!
185 : // In that case the algorithm reserves too much memory.
186 : // Use e.g. a virtual method GlobalMultiGridRefiner::reserve_memory
187 : GMGR_PROFILE(GMGR_Reserve);
188 : {
189 : int l = oldTopLevel;
190 :
191 : GMGR_PROFILE(GMGR_ReserveVrtData);
192 0 : mg.reserve<Vertex>(mg.num<Vertex>() +
193 0 : + mg.num<Vertex>(l) + mg.num<Edge>(l)
194 0 : + mg.num<Quadrilateral>(l) + mg.num<Hexahedron>(l));
195 : GMGR_PROFILE_END();
196 :
197 : GMGR_PROFILE(GMGR_ReserveEdgeData);
198 0 : mg.reserve<Edge>(mg.num<Edge>()
199 0 : + 2 * mg.num<Edge>(l) + 3 * mg.num<Triangle>(l)
200 0 : + 4 * mg.num<Quadrilateral>(l) + 3 * mg.num<Prism>(l)
201 0 : + mg.num<Tetrahedron>(l)
202 0 : + 4 * mg.num<Pyramid>(l) + 6 * mg.num<Hexahedron>(l));
203 : GMGR_PROFILE_END();
204 :
205 : GMGR_PROFILE(GMGR_ReserveFaceData);
206 0 : mg.reserve<Face>(mg.num<Face>()
207 0 : + 4 * mg.num<Face>(l) + 10 * mg.num<Prism>(l)
208 0 : + 8 * mg.num<Tetrahedron>(l)
209 0 : + 9 * mg.num<Pyramid>(l) + 12 * mg.num<Hexahedron>(l));
210 : GMGR_PROFILE_END();
211 :
212 : GMGR_PROFILE(GMGR_ReserveVolData);
213 0 : mg.reserve<Volume>(mg.num<Volume>()
214 0 : + 8 * mg.num<Tetrahedron>(l) + 8 * mg.num<Prism>(l)
215 0 : + 6 * mg.num<Pyramid>(l) + 8 * mg.num<Hexahedron>(l));
216 : GMGR_PROFILE_END();
217 : }
218 : GMGR_PROFILE_END();
219 : UG_DLOG(LIB_GRID, 1, " done.\n");
220 :
221 : UG_DLOG(LIB_GRID, 1, " refinement begins.\n");
222 : // notify derivates that refinement begins
223 0 : refinement_step_begins();
224 :
225 : // cout << "num marked edges: " << m_selMarks.num<Edge>() << endl;
226 : // cout << "num marked faces: " << m_selMarks.num<Face>() << endl;
227 :
228 : // we want to add new elements in a new layer.
229 : bool bHierarchicalInsertionWasEnabled = mg.hierarchical_insertion_enabled();
230 0 : if(!bHierarchicalInsertionWasEnabled)
231 0 : mg.enable_hierarchical_insertion(true);
232 :
233 :
234 : // some buffers
235 : vector<Vertex*> vVrts;
236 : vector<Vertex*> vEdgeVrts;
237 : vector<Vertex*> vFaceVrts;
238 : vector<Edge*> vEdges;
239 : vector<Face*> vFaces;
240 : vector<Volume*> vVols;
241 :
242 : // some repeatedly used objects
243 0 : EdgeDescriptor ed;
244 0 : FaceDescriptor fd;
245 0 : VolumeDescriptor vd;
246 :
247 : UG_DLOG(LIB_GRID, 1, " creating new vertices\n");
248 :
249 : // create new vertices from marked vertices
250 0 : for(VertexIterator iter = mg.begin<Vertex>(oldTopLevel);
251 0 : iter != mg.end<Vertex>(oldTopLevel); ++iter)
252 : {
253 0 : if(!refinement_is_allowed(*iter))
254 0 : continue;
255 :
256 : Vertex* v = *iter;
257 :
258 : // create a new vertex in the next layer.
259 : //GMGR_PROFILE(GMGR_Refine_CreatingVertices);
260 0 : Vertex* nVrt = *mg.create_by_cloning(v, v);
261 :
262 : // allow refCallback to calculate a new position
263 0 : if(m_projector.valid())
264 0 : m_projector->new_vertex(nVrt, v);
265 : //GMGR_PROFILE_END();
266 : }
267 :
268 :
269 : UG_DLOG(LIB_GRID, 1, " creating new edges\n");
270 :
271 : // create new vertices and edges from marked edges
272 0 : for(EdgeIterator iter = mg.begin<Edge>(oldTopLevel);
273 0 : iter != mg.end<Edge>(oldTopLevel); ++iter)
274 : {
275 0 : if(!refinement_is_allowed(*iter))
276 0 : continue;
277 :
278 : // collect_objects_for_refine removed all edges that already were
279 : // refined. No need to check that again.
280 : Edge* e = *iter;
281 :
282 : // debug: make sure that both vertices may be refined
283 : /* #ifdef UG_DEBUG
284 : if(!refinement_is_allowed(e->vertex(0))
285 : || !refinement_is_allowed(e->vertex(1)))
286 : {
287 : UG_LOG("Can't refine edge between vertices ");
288 : if(mg.has_vertex_attachment(aPosition)){
289 : Grid::VertexAttachmentAccessor<APosition> aaPos(mg, aPosition);
290 : UG_LOG(aaPos[e->vertex(0)] << " and " << aaPos[e->vertex(1)] << endl);
291 : }
292 : else if(mg.has_vertex_attachment(aPosition2)){
293 : Grid::VertexAttachmentAccessor<APosition2> aaPos(mg, aPosition2);
294 : UG_LOG(aaPos[e->vertex(0)] << " and " << aaPos[e->vertex(1)] << endl);
295 : }
296 : else if(mg.has_vertex_attachment(aPosition1)){
297 : Grid::VertexAttachmentAccessor<APosition1> aaPos(mg, aPosition1);
298 : UG_LOG(aaPos[e->vertex(0)] << " and " << aaPos[e->vertex(1)] << endl);
299 : }
300 : }
301 : #endif // UG_DEBUG
302 : */
303 :
304 : assert(refinement_is_allowed(e->vertex(0))
305 : && refinement_is_allowed(e->vertex(1)));
306 :
307 : //GMGR_PROFILE(GMGR_Refine_CreatingEdgeVertices);
308 : // create two new edges by edge-split
309 0 : RegularVertex* nVrt = *mg.create<RegularVertex>(e);
310 :
311 : // allow refCallback to calculate a new position
312 0 : if(m_projector.valid())
313 0 : m_projector->new_vertex(nVrt, e);
314 : //GMGR_PROFILE_END();
315 :
316 : // split the edge
317 : //GMGR_PROFILE(GMGR_Refine_CreatingEdges);
318 : Vertex* substituteVrts[2];
319 0 : substituteVrts[0] = mg.get_child_vertex(e->vertex(0));
320 0 : substituteVrts[1] = mg.get_child_vertex(e->vertex(1));
321 :
322 0 : e->refine(vEdges, nVrt, substituteVrts);
323 : assert((vEdges.size() == 2) && "RegularEdge refine produced wrong number of edges.");
324 0 : mg.register_element(vEdges[0], e);
325 0 : mg.register_element(vEdges[1], e);
326 : //GMGR_PROFILE_END();
327 : }
328 :
329 :
330 : UG_DLOG(LIB_GRID, 1, " creating new faces\n");
331 :
332 : // create new vertices and faces from marked faces
333 0 : for(FaceIterator iter = mg.begin<Face>(oldTopLevel);
334 0 : iter != mg.end<Face>(oldTopLevel); ++iter)
335 : {
336 0 : if(!refinement_is_allowed(*iter))
337 0 : continue;
338 :
339 : Face* f = *iter;
340 : // collect child-vertices
341 : vVrts.clear();
342 0 : for(uint j = 0; j < f->num_vertices(); ++j)
343 0 : vVrts.push_back(mg.get_child_vertex(f->vertex(j)));
344 :
345 : // collect the associated edges
346 : vEdgeVrts.clear();
347 : //bool bIrregular = false;
348 0 : for(uint j = 0; j < f->num_edges(); ++j)
349 0 : vEdgeVrts.push_back(mg.get_child_vertex(mg.get_edge(f, j)));
350 :
351 : //GMGR_PROFILE(GMGR_Refine_CreatingFaces);
352 : Vertex* newVrt;
353 0 : if(f->refine(vFaces, &newVrt, &vEdgeVrts.front(), NULL, &vVrts.front())){
354 : // if a new vertex was generated, we have to register it
355 0 : if(newVrt){
356 : //GMGR_PROFILE(GMGR_Refine_CreatingVertices);
357 0 : mg.register_element(newVrt, f);
358 : // allow refCallback to calculate a new position
359 0 : if(m_projector.valid())
360 0 : m_projector->new_vertex(newVrt, f);
361 : //GMGR_PROFILE_END();
362 : }
363 :
364 : // register the new faces and assign status
365 0 : for(size_t j = 0; j < vFaces.size(); ++j)
366 0 : mg.register_element(vFaces[j], f);
367 : }
368 : else{
369 : LOG(" WARNING in Refine: could not refine face.\n");
370 : }
371 : //GMGR_PROFILE_END();
372 : }
373 :
374 :
375 : UG_DLOG(LIB_GRID, 1, " creating new volumes\n");
376 :
377 : // only used for tetrahedron or octahedron refinement
378 0 : vector<vector3> corners(6, vector3(0, 0, 0));
379 :
380 : // create new vertices and volumes from marked volumes
381 0 : for(VolumeIterator iter = mg.begin<Volume>(oldTopLevel);
382 0 : iter != mg.end<Volume>(oldTopLevel); ++iter)
383 : {
384 0 : if(!refinement_is_allowed(*iter))
385 0 : continue;
386 :
387 : Volume* v = *iter;
388 : //GMGR_PROFILE(GMGR_Refining_Volume);
389 :
390 : // collect child-vertices
391 : //GMGR_PROFILE(GMGR_CollectingVolumeVertices);
392 : vVrts.clear();
393 0 : for(uint j = 0; j < v->num_vertices(); ++j)
394 0 : vVrts.push_back(mg.get_child_vertex(v->vertex(j)));
395 : //GMGR_PROFILE_END();
396 :
397 : // collect the associated edges
398 : vEdgeVrts.clear();
399 : //GMGR_PROFILE(GMGR_CollectingVolumeEdgeVertices);
400 : //bool bIrregular = false;
401 0 : for(uint j = 0; j < v->num_edges(); ++j)
402 0 : vEdgeVrts.push_back(mg.get_child_vertex(mg.get_edge(v, j)));
403 : //GMGR_PROFILE_END();
404 :
405 : // collect associated face-vertices
406 : vFaceVrts.clear();
407 : //GMGR_PROFILE(GMGR_CollectingVolumeFaceVertices);
408 0 : for(uint j = 0; j < v->num_faces(); ++j)
409 0 : vFaceVrts.push_back(mg.get_child_vertex(mg.get_face(v, j)));
410 : //GMGR_PROFILE_END();
411 :
412 : // if we're performing tetrahedral or octahedral refinement, we have to collect
413 : // the corner coordinates, so that the refinement algorithm may choose
414 : // the best interior diagonal.
415 : vector3* pCorners = NULL;
416 0 : if((v->num_vertices() == 4) && m_projector.valid()){
417 0 : for(size_t i = 0; i < 4; ++i){
418 0 : corners[i] = m_projector->geometry()->pos(v->vertex(i));
419 : }
420 : pCorners = &corners.front();
421 : }
422 0 : if((v->reference_object_id() == ROID_OCTAHEDRON) && m_projector.valid()){
423 0 : for(size_t i = 0; i < 6; ++i){
424 0 : corners[i] = m_projector->geometry()->pos(v->vertex(i));
425 : }
426 : pCorners = &corners.front();
427 : }
428 :
429 : Vertex* newVrt;
430 0 : if(v->refine(vVols, &newVrt, &vEdgeVrts.front(), &vFaceVrts.front(),
431 0 : NULL, RegularVertex(), &vVrts.front(), pCorners)){
432 : // if a new vertex was generated, we have to register it
433 0 : if(newVrt){
434 0 : mg.register_element(newVrt, v);
435 : // allow refCallback to calculate a new position
436 0 : if(m_projector.valid())
437 0 : m_projector->new_vertex(newVrt, v);
438 : }
439 :
440 : // register the new faces and assign status
441 0 : for(size_t j = 0; j < vVols.size(); ++j)
442 0 : mg.register_element(vVols[j], v);
443 : }
444 : else{
445 : LOG(" WARNING in Refine: could not refine volume.\n");
446 : }
447 : //GMGR_PROFILE_END();
448 : }
449 :
450 : // done - clean up
451 0 : if(!bHierarchicalInsertionWasEnabled)
452 0 : mg.enable_hierarchical_insertion(false);
453 :
454 : // notify derivates that refinement ends
455 0 : refinement_step_ends();
456 :
457 0 : projector()->refinement_ends();
458 0 : m_messageHub->post_message(GridMessage_Adaption(GMAT_GLOBAL_REFINEMENT_ENDS,
459 0 : mg.get_grid_objects(oldTopLevel)));
460 :
461 : UG_DLOG(LIB_GRID, 1, " refinement done.");
462 0 : }
463 :
464 0 : bool GlobalMultiGridRefiner::save_marks_to_file(const char* filename)
465 : {
466 : GMGR_PROFILE(GlobalMultiGridRefiner_save_marks_to_file);
467 0 : if(!m_pMG){
468 0 : UG_THROW("ERROR in GlobalMultiGridRefiner::save_marks_to_file: No grid assigned!");
469 : }
470 :
471 : MultiGrid& mg = *m_pMG;
472 0 : SubsetHandler sh(mg);
473 :
474 0 : AssignGridToSubset(mg, sh, 2);
475 0 : int lvl = mg.num_levels() - 1;
476 0 : sh.assign_subset(mg.begin<Vertex>(lvl), mg.end<Vertex>(lvl), 0);
477 0 : sh.assign_subset(mg.begin<Edge>(lvl), mg.end<Edge>(lvl), 0);
478 0 : sh.assign_subset(mg.begin<Face>(lvl), mg.end<Face>(lvl), 0);
479 0 : sh.assign_subset(mg.begin<Volume>(lvl), mg.end<Volume>(lvl), 0);
480 :
481 0 : sh.subset_info(0).name = "refine";
482 0 : sh.subset_info(1).name = "no marks";
483 :
484 0 : AssignSubsetColors(sh);
485 :
486 0 : return SaveGridToFile(mg, sh, filename);
487 0 : }
488 :
489 : }// end of namespace
|
