Line data Source code
1 : /*
2 : * Copyright (c) 2015: G-CSC, Goethe University Frankfurt
3 : * Author: Dmitry Logashenko
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 : /*
34 : * Implementation of the degenerated layer subset manager.
35 : */
36 : // ug4 headers
37 : #include "common/util/string_util.h"
38 :
39 : namespace ug {
40 :
41 : /**
42 : * Class constructor: Parses the subset names, attaches itself to the grid
43 : * to be updated if the grid is refined, and writes itself to the properties
44 : * of the subsets.
45 : */
46 : template <int dim>
47 0 : DegeneratedLayerManager<dim>::DegeneratedLayerManager
48 : (
49 : SmartPtr<MultiGridSubsetHandler> spSH ///< [in] subset handler of the grid
50 : )
51 0 : : m_spSH (spSH), m_bClosed (false)
52 : {
53 : // Check the subset handler:
54 0 : if (m_spSH.invalid ())
55 0 : UG_THROW ("DegeneratedLayerManager: Invalid subset handler specified!");
56 :
57 : // Initialize the data:
58 0 : m_layerSsGrp.set_subset_handler (m_spSH);
59 :
60 : // Initialize the attachment:
61 : MultiGrid * pMG = (MultiGrid *) (m_spSH->multi_grid ());
62 0 : pMG->template attach_to_dv<Vertex> (m_aVertexMarks, D_LAYER_UNKNOWN);
63 : m_aaVertMarks.access (*pMG, m_aVertexMarks);
64 :
65 : // Register the callbacks in the message hub:
66 0 : m_spGridAdaptionCallbackID =
67 : pMG->message_hub()->register_class_callback
68 : (this, & DegeneratedLayerManager<dim>::grid_adaption_callback);
69 0 : m_spGridDistributionCallbackID =
70 : pMG->message_hub()->register_class_callback
71 : (this, & DegeneratedLayerManager<dim>::grid_distribution_callback);
72 0 : }
73 :
74 : /**
75 : * Class destructor: Detaches the attachment.
76 : */
77 : template <int dim>
78 0 : DegeneratedLayerManager<dim>::~DegeneratedLayerManager ()
79 : {
80 : MultiGrid * pMG = (MultiGrid *) (m_spSH->multi_grid ());
81 : m_aaVertMarks.invalidate ();
82 0 : pMG->template detach_from<Vertex> (m_aVertexMarks);
83 0 : m_bClosed = false; // (to be on the very safe side...)
84 0 : }
85 :
86 : /**
87 : * Adds subsets to the manager
88 : */
89 : template <int dim>
90 0 : void DegeneratedLayerManager<dim>::add
91 : (
92 : const char * ss_names ///< [in] subset names of the degenerated layers
93 : )
94 : {
95 0 : if (m_bClosed)
96 : UG_LOG ("DegeneratedLayerManager: Adding subsets to a closed manager.");
97 :
98 : // If extended - then not closed any more
99 0 : m_bClosed = false;
100 :
101 : // Parse the subset names:
102 : std::vector<std::string> vssNames;
103 0 : TokenizeString (ss_names, vssNames);
104 :
105 : // Add the subsets to the group:
106 0 : m_layerSsGrp.add (vssNames);
107 0 : }
108 :
109 : /**
110 : * Removes subsets from the manager
111 : */
112 : template <int dim>
113 0 : void DegeneratedLayerManager<dim>::remove
114 : (
115 : const char * ss_names ///< [in] subset names of the fractures
116 : )
117 : {
118 : // If extended - then not closed any more
119 0 : m_bClosed = false;
120 :
121 : // Parse the subset names:
122 : std::vector<std::string> vssNames;
123 0 : TokenizeString (ss_names, vssNames);
124 :
125 : // Add the subsets to the group:
126 0 : m_layerSsGrp.remove (vssNames);
127 0 : }
128 :
129 : /**
130 : * Closes the manager: Computes the vertex marks
131 : */
132 : template <int dim>
133 0 : void DegeneratedLayerManager<dim>::close ()
134 : {
135 : size_t init_n_subsets = m_layerSsGrp.size ();
136 :
137 : // Check the dimensionality of the subsets:
138 : size_t i = 0;
139 0 : while (i < m_layerSsGrp.size ())
140 0 : if (DimensionOfSubset (*m_spSH, m_layerSsGrp [i]) != dim)
141 : {
142 0 : UG_LOG ("DegeneratedLayerManager: Low-dimensional subset '"
143 : << m_spSH->get_subset_name (m_layerSsGrp [i])
144 : << "' specified as a degenerated layer. Removed from the list.\n");
145 0 : m_layerSsGrp.remove (m_layerSsGrp [i]);
146 : i = 0;
147 : }
148 0 : else i++;
149 :
150 : // Check if all the subsets are removed due to the illegal dimensionality
151 0 : if (init_n_subsets != 0 && m_layerSsGrp.size () == 0)
152 0 : UG_THROW ("DegeneratedLayerManager: All the specified subsets have illegal dimensionality!\n");
153 :
154 : // Mark the vertices:
155 0 : mark_vertices ();
156 :
157 : // Done:
158 0 : m_bClosed = true;
159 0 : }
160 :
161 : /**
162 : * Adds the registered subsets to the refiner
163 : */
164 : template <int dim>
165 0 : void DegeneratedLayerManager<dim>::init_refiner
166 : (
167 : SmartPtr<GlobalFracturedMediaRefiner> refiner,
168 : bool as_low_dim
169 : )
170 : {
171 0 : if (! is_closed ())
172 0 : UG_THROW ("DegeneratedLayerManager: The manager must be closed before use.");
173 0 : for (size_t i = 0; i < m_layerSsGrp.size (); i++)
174 0 : refiner->mark_as_fracture (m_layerSsGrp[i], as_low_dim);
175 0 : }
176 :
177 : /**
178 : * Marks the vertices in such a way that the inner layer vertices are
179 : * marked with D_LAYER_INNER, whereas all other vertices (at the fractures
180 : * and in the bulk medium) are marked with D_LAYER_OUTER. After the call,
181 : * there should be no vertices marked with D_LAYER_UNKNOWN.
182 : */
183 : template <int dim>
184 0 : void DegeneratedLayerManager<dim>::mark_vertices ()
185 : {
186 : typedef typename geometry_traits<element_type>::const_iterator t_elem_iter;
187 :
188 : MultiGrid * pMG = (MultiGrid *) (m_spSH->multi_grid ());
189 :
190 : // Mark all vertices of the fracture elements:
191 0 : for (size_t i = 0; i < m_layerSsGrp.size (); i++)
192 : {
193 : int si = m_layerSsGrp [i];
194 0 : for (size_t lev = 0; lev < pMG->num_levels (); lev++)
195 : {
196 0 : t_elem_iter list_end = m_spSH->template end<element_type> (si, lev);
197 0 : for (t_elem_iter iElem = m_spSH->template begin<element_type> (si, lev);
198 0 : iElem != list_end; ++iElem)
199 : {
200 : element_type * elem = *iElem;
201 0 : for (size_t i = 0; i < elem->num_vertices (); i++)
202 : {
203 0 : Vertex * vert = elem->vertex(i);
204 0 : m_aaVertMarks [vert] = (! vert->is_constrained ())? D_LAYER_INNER : D_LAYER_OUTER;
205 : // REMARK: Constrained vertices cannot be inner for degenerated layers!
206 : // We cannot catch this situation somewhere else, so we do it here.
207 : }
208 : }
209 : }
210 : }
211 :
212 : // Unmark all vertices in the other (not-degenerated-layer) elements:
213 0 : for (int si = 0; si < m_spSH->num_subsets (); si++)
214 0 : if (DimensionOfSubset (*m_spSH, si) == dim && ! m_layerSsGrp.contains (si))
215 0 : for (size_t lev = 0; lev < pMG->num_levels (); lev++)
216 : {
217 0 : t_elem_iter list_end = m_spSH->template end<element_type> (si, lev);
218 0 : for (t_elem_iter iElem = m_spSH->template begin<element_type> (si, lev);
219 0 : iElem != list_end; ++iElem)
220 : {
221 : element_type * elem = *iElem;
222 0 : for (size_t i = 0; i < elem->num_vertices (); i++)
223 0 : m_aaVertMarks [elem->vertex(i)] = D_LAYER_OUTER;
224 : }
225 : }
226 :
227 : #ifdef UG_PARALLEL
228 : // Make the marking consistent
229 : AttachmentAllReduce<Vertex> (*pMG, m_aVertexMarks, PCL_RO_MIN);
230 : #endif // UG_PARALLEL
231 0 : }
232 :
233 : /**
234 : * The grid adaption callback: Catches the messages about grid refinements, ...
235 : */
236 : template <int dim>
237 0 : void DegeneratedLayerManager<dim>::grid_adaption_callback
238 : (
239 : const GridMessage_Adaption & msg ///< the message from the message hub
240 : )
241 : {
242 0 : if (m_bClosed && msg.adaption_ends ())
243 0 : mark_vertices ();
244 0 : }
245 :
246 : /**
247 : * The grid adaption callback: Catches the messages about the distribution of the grid
248 : */
249 : template <int dim>
250 0 : void DegeneratedLayerManager<dim>::grid_distribution_callback
251 : (
252 : const GridMessage_Distribution & msg ///< the message from the message hub
253 : )
254 : {
255 0 : if (m_bClosed && msg.msg () == GMDT_DISTRIBUTION_STOPS)
256 0 : mark_vertices ();
257 0 : }
258 :
259 : /**
260 : * For a given element, finds its inner and outer layer sides. For these
261 : * sides, gets the correspondence of the vertices in the degenerated element
262 : * (i.e. considering all other sides as degenerated).
263 : */
264 : template <int dim>
265 0 : void DegeneratedLayerManager<dim>::get_layer_sides
266 : (
267 : element_type * elem, ///< [in] the element
268 : size_t & num_side_co, ///< [out] number of corners of the inner/outer sides
269 : side_type * & inner_side, ///< [out] its fracture inner side
270 : size_t & inner_side_idx, ///< [out] index of the inner side in the reference element
271 : size_t inner_side_corners [], ///< [out] inner side corner idx -> elem. corner idx (maxLayerSideCorners elements)
272 : side_type * & outer_side, ///< [out] its fracture outer side
273 : size_t & outer_side_idx, ///< [out] index of the outer side in the reference element
274 : size_t outer_side_corners [], ///< [out] outer side corner idx -> elem. corner idx (maxLayerSideCorners elements)
275 : size_t ass_co [] ///< [out] correspondence of the corners of the sides (2 * maxLayerSideCorners elements or NULL)
276 : )
277 : {
278 : size_t n_inner, n_outer;
279 0 : size_t n_co = elem->num_vertices ();
280 : MultiGrid * pMG = (MultiGrid *) (m_spSH->multi_grid ());
281 : typename Grid::traits<Edge>::secure_container edge_list;
282 :
283 : // First, we check, whether the inner and outer sides are well defined.
284 : // Loop over the sides: We look for sides with only inner or only outer corners
285 0 : inner_side = NULL; n_inner = 0;
286 0 : outer_side = NULL; n_outer = 0;
287 0 : for (size_t i = 0; i < elem->num_sides (); i++)
288 : {
289 0 : side_type * side = pMG->get_side (elem, i);
290 : bool has_inner = false, has_outer = false;
291 0 : size_t side_n_co = side->num_vertices ();
292 0 : for (size_t co = 0; co < side_n_co; co++)
293 : {
294 0 : int mark = vert_mark (side->vertex (co));
295 0 : if (mark == D_LAYER_OUTER)
296 : has_outer = true;
297 0 : else if (mark == D_LAYER_INNER)
298 : has_inner = true;
299 : else
300 0 : UG_THROW ("DegeneratedLayerManager: Some vertices are not marked!");
301 : }
302 :
303 0 : if (has_inner && has_outer) continue; // this is a 'degenerated' side
304 :
305 0 : if (has_inner)
306 : { // this is the inner side
307 0 : if (inner_side != NULL)
308 0 : UG_THROW ("DegeneratedLayerManager: Two inner sides of a degenerated layer element found!");
309 0 : inner_side = side;
310 0 : inner_side_idx = i;
311 : n_inner = side_n_co;
312 : }
313 : else
314 : { // this is the outer side
315 0 : if (outer_side != NULL)
316 0 : UG_THROW ("DegeneratedLayerManager: Two outer sides of a degenerated layer element found!");
317 0 : outer_side = side;
318 0 : outer_side_idx = i;
319 : n_outer = side_n_co;
320 : }
321 : }
322 0 : if (n_inner != n_outer || 2 * n_inner != n_co)
323 0 : UG_THROW ("DegeneratedLayerManager: Unrecognized degenerated layer element!");
324 0 : num_side_co = n_inner;
325 :
326 : // Now: Corner of inner side are exactly those marked as inner, and the same for the outer side
327 : Vertex * inner_vert_ptr [maxLayerSideCorners], * outer_vert_ptr [maxLayerSideCorners];
328 : size_t inner_corners [maxLayerSideCorners], outer_corners [maxLayerSideCorners];
329 : size_t inner_i = 0, outer_i = 0;
330 0 : for (size_t co = 0; co < n_co; co++)
331 : {
332 0 : Vertex * vrt = elem->vertex (co);
333 0 : if (vert_mark (vrt) == D_LAYER_INNER)
334 : {
335 0 : inner_vert_ptr [inner_i] = vrt;
336 0 : inner_corners [inner_i] = co;
337 0 : inner_i++;
338 : }
339 : else
340 : {
341 0 : outer_vert_ptr [outer_i] = vrt;
342 0 : outer_corners [outer_i] = co;
343 0 : outer_i++;
344 : }
345 : }
346 0 : if (inner_i != n_inner || outer_i != n_outer)
347 0 : UG_THROW ("DegeneratedLayerManager: Failed to get the vertex-side correspondence!");
348 :
349 : // Order the corner indices according to the corners of the sides
350 0 : for (size_t i = 0; i < n_inner; i++)
351 : {
352 : size_t co;
353 :
354 : co = 0;
355 0 : while (inner_side->vertex (co) != inner_vert_ptr [i])
356 0 : if (++co == n_inner)
357 0 : UG_THROW ("DegeneratedLayerManager: Failed! Missing inner vertex?");
358 0 : inner_side_corners [co] = inner_corners [i];
359 :
360 : co = 0;
361 0 : while (outer_side->vertex (co) != outer_vert_ptr [i])
362 0 : if (++co == n_outer)
363 0 : UG_THROW ("DegeneratedLayerManager: Failed! Missing outer vertex?");
364 0 : outer_side_corners [co] = outer_corners [i];
365 : }
366 :
367 : // Loop over the edges: Find out the associated corners
368 0 : if (ass_co != NULL)
369 : {
370 0 : pMG->associated_elements (edge_list, elem);
371 0 : for (size_t i = 0; i < edge_list.size (); i++)
372 : {
373 : Edge * e = edge_list [i];
374 0 : Vertex * v_1 = e->vertex (0); int mark_1 = vert_mark (v_1);
375 0 : Vertex * v_2 = e->vertex (1); int mark_2 = vert_mark (v_2);
376 : Vertex * t;
377 : size_t co_1, co_2;
378 :
379 0 : if (mark_1 == mark_2)
380 0 : continue; // an edge of the inner or the outer size
381 :
382 0 : if (mark_1 == D_LAYER_OUTER)
383 : { // Vertex # 1 must be inner, vertex # 2 must be outer
384 : t = v_1; v_1 = v_2; v_2 = t;
385 : }
386 :
387 : co_1 = 0;
388 0 : while (inner_vert_ptr[co_1] != v_1)
389 0 : if (++co_1 == n_inner)
390 0 : UG_THROW ("DegeneratedLayerManager: Cannot find an inner node by vertex!");
391 :
392 : co_2 = 0;
393 0 : while (outer_vert_ptr[co_2] != v_2)
394 0 : if (++co_2 == n_outer)
395 0 : UG_THROW ("DegeneratedLayerManager: Cannot find an outer node by vertex!");
396 :
397 0 : ass_co [inner_corners [co_1]] = outer_corners [co_2];
398 0 : ass_co [outer_corners [co_2]] = inner_corners [co_1];
399 : }
400 : }
401 0 : }
402 :
403 : /**
404 : * Assigns a different subset index to the inner sides of a layer. (Only to
405 : * the sides, not vertices, or edges in 3d!) If -1 passed instead of the index
406 : * of the subset to assign, new subset is created. The function returns the
407 : * index of the assigned subset.
408 : */
409 : template <int dim>
410 0 : int DegeneratedLayerManager<dim>::assign_middle_subset
411 : (
412 : int layer_si, ///< subset index of the layer
413 : int middle_si ///< the subset index to assign (or -1 to create a new subset)
414 : )
415 : {
416 : typedef typename geometry_traits<element_type>::const_iterator t_elem_iter;
417 :
418 : MultiGrid * pMG = (MultiGrid *) (m_spSH->multi_grid ());
419 :
420 : // Assign the default value of the subset (if necessary)
421 0 : if (middle_si < 0) middle_si = m_spSH->num_subsets ();
422 :
423 : // Loop the elements of the layer in all the grid levels
424 0 : pMG->message_hub()->post_message (GridMessage_Creation (GMCT_CREATION_STARTS, 0));
425 0 : for (size_t lev = 0; lev < pMG->num_levels (); lev++)
426 : {
427 0 : t_elem_iter list_end = m_spSH->template end<element_type> (layer_si, lev);
428 0 : for (t_elem_iter iElem = m_spSH->template begin<element_type> (layer_si, lev);
429 0 : iElem != list_end; ++iElem)
430 : {
431 : element_type * elem = *iElem;
432 : size_t num_side_co;
433 : side_type * inner_side, * outer_side;
434 : size_t inner_side_idx, outer_side_idx;
435 : size_t inner_side_corners [maxLayerSideCorners], outer_side_corners [maxLayerSideCorners];
436 :
437 : // get the inner side
438 0 : get_layer_sides (elem, num_side_co,
439 : inner_side, inner_side_idx, inner_side_corners,
440 : outer_side, outer_side_idx, outer_side_corners);
441 :
442 : // assign the subset
443 0 : m_spSH->assign_subset (inner_side, middle_si);
444 : }
445 : }
446 0 : pMG->message_hub()->post_message (GridMessage_Creation (GMCT_CREATION_STOPS, 0));
447 :
448 0 : return middle_si;
449 : }
450 :
451 : /**
452 : * Assigns a different subset to the inner sides of a layer. (Only to
453 : * the sides, not vertices, or edges in 3d!) If the given subset does not exist,
454 : * it is created. The function returns the index of the assigned subset.
455 : */
456 : template <int dim>
457 0 : int DegeneratedLayerManager<dim>::assign_middle_subset
458 : (
459 : int layer_si, ///< subset index of the layer
460 : const char* middle_ss_name ///< name of the subset to assign
461 : )
462 : {
463 : // Look for the subset to assign
464 : int middle_si = -1;
465 : bool new_ss = true;
466 0 : for (int si = 0; si < m_spSH->num_subsets (); si++)
467 0 : if (strcmp (middle_ss_name, m_spSH->get_subset_name (si)) == 0)
468 : {
469 : middle_si = si;
470 : new_ss = false;
471 : break;
472 : }
473 :
474 : // Assign the subset
475 0 : middle_si = assign_middle_subset (layer_si, middle_si);
476 :
477 : // Set the name
478 0 : if (new_ss)
479 0 : m_spSH->set_subset_name (middle_ss_name, middle_si);
480 :
481 0 : return middle_si;
482 : }
483 :
484 : /**
485 : * Assigns a different subset to the inner sides of a layer. (Only to
486 : * the sides, not vertices, or edges in 3d!) If the given subset does not exist,
487 : * it is created. The function returns the index of the assigned subset.
488 : */
489 : template <int dim>
490 0 : int DegeneratedLayerManager<dim>::assign_middle_subset
491 : (
492 : const char* layer_ss_name, ///< subset name of the layer
493 : const char* middle_ss_name ///< name of the subset to assign
494 : )
495 : {
496 : // Look for the subset of the layer
497 0 : for (int si = 0; si < m_spSH->num_subsets (); si++)
498 0 : if (strcmp (layer_ss_name, m_spSH->get_subset_name (si)) == 0)
499 0 : return assign_middle_subset (si, middle_ss_name);
500 0 : UG_THROW ("DegeneratedLayerManager: Subset '" << layer_ss_name << "' not found");
501 : }
502 :
503 : } // end namespace ug
504 :
505 : /* End of File */
|
