Line data Source code
1 : /*
2 : * Copyright (c) 2009-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 : #ifndef __H__LIB_GRID__SUBSET_UTIL_IMPL__
34 : #define __H__LIB_GRID__SUBSET_UTIL_IMPL__
35 :
36 : #include "subset_util.h"
37 : #include "selection_util.h"
38 : #include "lib_grid/callbacks/basic_callbacks.h"
39 : #include "lib_grid/callbacks/subset_callbacks.h"
40 : #include "lib_grid/callbacks/selection_callbacks.h"
41 :
42 : namespace ug
43 : {
44 : ////////////////////////////////////////////////////////////////////////
45 : // FindFirstFreeSubset
46 : template <class TElem>
47 0 : int GetMaxSubsetIndex(SubsetHandler& sh)
48 : {
49 : // go from back to front
50 0 : for(int i = (int)sh.num_subsets() - 1; i >= 0; --i)
51 : {
52 0 : if(sh.num_elements<TElem>(i) > 0)
53 : {
54 : // this is the highest subset that contains elements of type TElem
55 0 : return i;
56 : }
57 : }
58 :
59 : // no subset contains elements of type TElem.
60 : return -1;
61 : }
62 :
63 : ////////////////////////////////////////////////////////////////////////
64 : // MakeSubsetsConsecutive
65 : template <class TElem>
66 0 : void MakeSubsetsConsecutive(SubsetHandler& sh)
67 : {
68 : // TODO: this algo could be slightly improved regarding runtime.
69 :
70 : // iterate through all subsets.
71 0 : for(int i = 0; i < sh.num_subsets(); ++i)
72 : {
73 : // check whether the subset is empty
74 0 : if(sh.num_elements<TElem>(i) == 0)
75 : {
76 : // it is. find the next filled one.
77 0 : for(int j = i + 1; j < sh.num_subsets(); ++j)
78 : {
79 0 : if(sh.num_elements<TElem>(j) > 0)
80 : {
81 : // this subset is filled. move it to position i.
82 0 : sh.move_subset(j, i);
83 0 : break;
84 : }
85 : }
86 : }
87 : }
88 0 : }
89 :
90 : ////////////////////////////////////////////////////////////////////////
91 : // EraseEmptySubsets
92 : /// Erases all subsets which do not contain any geometric objects
93 : void EraseEmptySubsets(ISubsetHandler& sh);
94 :
95 :
96 : template <class TElem>
97 : void SeparateSubsetsByLowerDimSubsets(Grid& grid, SubsetHandler& sh,
98 : bool appendAtEnd)
99 : {
100 : SeparateSubsetsByLowerDimSeparators<TElem>(grid, sh, appendAtEnd,
101 : IsNotInSubset(sh, -1));
102 : }
103 :
104 : template <class TElem>
105 : void SeparateSubsetsByLowerDimSelection(Grid& grid, SubsetHandler& sh,
106 : Selector& sel, bool appendAtEnd)
107 : {
108 : SeparateSubsetsByLowerDimSeparators<TElem>(grid, sh, appendAtEnd,
109 : IsSelected(sel));
110 : }
111 :
112 : template <class TElem>
113 : void SeparateSubsetsByLowerDimSeparators(Grid& grid, SubsetHandler& sh,
114 : bool appendAtEnd,
115 : boost::function<bool (typename TElem::lower_dim_base_object*)>
116 : cbIsSeparator)
117 :
118 : {
119 : using namespace std;
120 :
121 : // the element type of separating elements
122 : typedef typename TElem::lower_dim_base_object TSide;
123 :
124 : // assign all elements to subset -1
125 : sh.assign_subset(grid.begin<TElem>(), grid.end<TElem>(), -1);
126 :
127 : // we'll keep all unassigned volumes in a selector.
128 : Selector sel(grid);
129 : sel.select(grid.begin<TElem>(), grid.end<TElem>());
130 :
131 : // those vectors will be used to gather element neighbours.
132 : vector<TSide*> vSides;
133 : vector<TElem*> vElems;
134 :
135 : // this stack contains all volumes that we still have to check for neighbours.
136 : stack<TElem*> stkElems;
137 :
138 : // now - while there are unassigned elements.
139 : int subsetIndex = 0;
140 : if(appendAtEnd)
141 : subsetIndex = sh.num_subsets();
142 :
143 : while(!sel.empty())
144 : {
145 : // choose the element with which we want to start
146 : // TODO: if material-points are supplied, this should be the
147 : // the element that contains the i-th material point.
148 : stkElems.push(*sel.begin<TElem>());
149 : while(!stkElems.empty())
150 : {
151 : TElem* elem = stkElems.top();
152 : stkElems.pop();
153 : // if the volume is unselected it has already been processed.
154 : if(!sel.is_selected(elem))
155 : continue;
156 : sel.deselect(elem);
157 :
158 : // assign elem to its new subset
159 : sh.assign_subset(elem, subsetIndex);
160 :
161 : // check neighbour-elements, whether they belong to the same subset.
162 : // iterate through the sides of the element
163 : for(uint i = 0; i < elem->num_sides(); ++i)
164 : {
165 : // get the i-th side
166 : TSide* side = grid.get_side(elem, i);
167 :
168 : // check whether the side is regarded as a separator.
169 : // If not, we'll add all associated elements.
170 : if(!cbIsSeparator(side))
171 : {
172 : CollectAssociated(vElems, grid, side);
173 :
174 : // add all elements that are still selected (elem is not selected anymore).
175 : for(uint j = 0; j < vElems.size(); ++j)
176 : {
177 : if(sel.is_selected(vElems[j]))
178 : stkElems.push(vElems[j]);
179 : }
180 : }
181 : }
182 : }
183 : // the stack is empty. increase subset index.
184 : subsetIndex++;
185 : }
186 : }
187 :
188 :
189 : ////////////////////////////////////////////////////////////////////////
190 : template <class TIterator>
191 0 : void CopySubsetIndicesToSides(ISubsetHandler& sh, TIterator elemsBegin,
192 : TIterator elemsEnd, bool toUnassignedOnly)
193 : {
194 : typedef typename PtrToValueType<typename TIterator::value_type>::base_type TElem;
195 : typedef typename TElem::side TSide;
196 :
197 : if(!TElem::HAS_SIDES)
198 : return;
199 :
200 : UG_ASSERT(sh.grid(), "No grid assigned to subset-handler");
201 :
202 0 : Grid& grid = *sh.grid();
203 :
204 : typename Grid::traits<TSide>::secure_container sides;
205 0 : for(TIterator iter = elemsBegin; iter != elemsEnd; ++iter){
206 : TElem* e = *iter;
207 :
208 : int si = sh.get_subset_index(e);
209 :
210 : grid.associated_elements(sides, e);
211 0 : for(size_t i = 0; i < sides.size(); ++i){
212 : TSide* s = sides[i];
213 0 : if(toUnassignedOnly){
214 0 : if(sh.get_subset_index(s) == -1)
215 0 : sh.assign_subset(s, si);
216 : }
217 : else
218 0 : sh.assign_subset(s, si);
219 : }
220 : }
221 : }
222 :
223 :
224 : ////////////////////////////////////////////////////////////////////////
225 : template <class TElem, class TSubsetHandler>
226 0 : void AssignUnassignedElemsToSubset(TSubsetHandler& sh, int si)
227 : {
228 : typedef typename geometry_traits<TElem>::iterator ElemIter;
229 :
230 : // access the grid on which sh operates.
231 0 : if(!sh.grid())
232 : return;
233 :
234 0 : Grid& grid = *sh.grid();
235 :
236 : // first make sure, that all elems are assigned to a subset, since
237 : // those won't be processed later on.
238 :
239 : // num is not part of ISubsetHandler and thus causes problems, if sh has type ISubsetHandler
240 : //if(sh.template num<TElem>() != grid.num<TElem>()){
241 : for(ElemIter iter = grid.begin<TElem>();
242 0 : iter != grid.end<TElem>(); ++iter)
243 : {
244 0 : if(sh.get_subset_index(*iter) == -1)
245 0 : sh.assign_subset(*iter, si);
246 : }
247 : //}
248 : }
249 :
250 : ////////////////////////////////////////////////////////////////////////
251 : template <class TSubsetHandler>
252 0 : void AdjustSubsetsForSimulation(TSubsetHandler& sh,
253 : bool preserveExistingSubsets)
254 : {
255 : // access the grid on which sh operates.
256 0 : if(!sh.grid())
257 0 : return;
258 :
259 0 : Grid& grid = *sh.grid();
260 0 : Selector sel(grid);
261 :
262 0 : if(grid.num_volumes() > 0){
263 : // deselect all elements of lower dimension, if existing subsets are
264 : // not to be preserved.
265 0 : if(!preserveExistingSubsets){
266 : sh.assign_subset(grid.begin<Face>(), grid.end<Face>(), -1);
267 : sh.assign_subset(grid.begin<Edge>(), grid.end<Edge>(), -1);
268 : sh.assign_subset(grid.begin<Vertex>(), grid.end<Vertex>(), -1);
269 : }
270 :
271 0 : AssignUnassignedElemsToSubset<Volume>(sh, GetFirstFreeSubset(sh));
272 0 : AssignSidesToSubsets<Volume>(sh);
273 :
274 0 : SelectInterfaceElements(sel, sh, grid.begin<Face>(), grid.end<Face>());
275 0 : SelectBoundaryElements(sel, grid.begin<Face>(), grid.end<Face>());
276 0 : SelectAssociatedEdges(sel, sel.begin<Face>(), sel.end<Face>());
277 0 : SelectAssociatedVertices(sel, sel.begin<Face>(), sel.end<Face>());
278 :
279 0 : AssignSidesToSubsets<Face>(sh, &sel);
280 0 : AssignSidesToSubsets<Edge>(sh, &sel);
281 :
282 : // finally assign vertices on edge interfaces
283 0 : sel.clear<Edge>();
284 0 : sel.clear<Vertex>();
285 0 : SelectInterfaceElements(sel, sh, grid.begin<Edge>(),
286 : grid.end<Edge>(), true);
287 0 : sel.clear<Face>();
288 0 : SelectAssociatedVertices(sel, sel.begin<Edge>(), sel.end<Edge>());
289 :
290 0 : AssignSidesToSubsets<Edge>(sh, &sel);
291 : }
292 0 : else if(grid.num_faces() > 0){
293 : // deselect all elements of lower dimension, if existing subsets are
294 : // not to be preserved.
295 0 : if(!preserveExistingSubsets){
296 : sh.assign_subset(grid.begin<Edge>(), grid.end<Edge>(), -1);
297 : sh.assign_subset(grid.begin<Vertex>(), grid.end<Vertex>(), -1);
298 : }
299 :
300 0 : AssignUnassignedElemsToSubset<Face>(sh, GetFirstFreeSubset(sh));
301 0 : AssignSidesToSubsets<Face>(sh);
302 :
303 0 : SelectInterfaceElements(sel, sh, grid.begin<Edge>(), grid.end<Edge>());
304 0 : SelectBoundaryElements(sel, grid.begin<Edge>(), grid.end<Edge>());
305 0 : SelectAssociatedVertices(sel, sel.begin<Edge>(), sel.end<Edge>());
306 :
307 0 : AssignSidesToSubsets<Edge>(sh, &sel);
308 : }
309 0 : else if(grid.num_edges() > 0){
310 : // deselect all elements of lower dimension, if existing subsets are
311 : // not to be preserved.
312 0 : if(!preserveExistingSubsets){
313 : sh.assign_subset(grid.begin<Vertex>(), grid.end<Vertex>(), -1);
314 : }
315 :
316 0 : AssignUnassignedElemsToSubset<Edge>(sh, GetFirstFreeSubset(sh));
317 0 : AssignSidesToSubsets<Edge>(sh);
318 : }
319 : else{
320 0 : AssignUnassignedElemsToSubset<Vertex>(sh, GetFirstFreeSubset(sh));
321 : }
322 :
323 : // erase empty subsets again
324 0 : EraseEmptySubsets(sh);
325 0 : }
326 :
327 : ////////////////////////////////////////////////////////////////////////
328 : template <class TAAPosVRT>
329 0 : number FaceArea(ISubsetHandler& sh, int si, size_t lvl, TAAPosVRT& aaPos)
330 : {
331 : number sum = 0.;
332 0 : GridObjectCollection goc = sh.get_grid_objects_in_subset(si);
333 :
334 0 : if (goc.num<Face>(lvl) == 0) {
335 : UG_WARNING("WARNING: Given subset doesn't contain any faces on given level.");
336 : } else {
337 : typedef geometry_traits<Face>::const_iterator CIT;
338 0 : for (CIT cit = goc.faces_begin(lvl); cit != goc.faces_end(lvl); cit++)
339 0 : sum += FaceArea(*cit, aaPos);
340 : }
341 :
342 0 : return sum;
343 : }
344 :
345 :
346 : ////////////////////////////////////////////////////////////////////////
347 : // AssignAssociatedVerticesToSubset
348 : template <class TIterator>
349 : void AssignAssociatedVerticesToSubset(ISubsetHandler& sh, TIterator elemsBegin,
350 : TIterator elemsEnd, int subsetIndex)
351 : {
352 : // iterate through the elements
353 : for(;elemsBegin != elemsEnd; elemsBegin++)
354 : {
355 : typename TIterator::value_type elem = *elemsBegin;
356 : uint numVrts = elem->num_vertices();
357 : // iterate through the vertices of elem and assign them
358 : for(uint i = 0; i < numVrts; ++i)
359 : sh.assign_subset(elem->vertex(i), subsetIndex);
360 : }
361 : }
362 :
363 : ////////////////////////////////////////////////////////////////////////
364 : template <class TElem, class TSubsetHandler>
365 : void AssignAssociatedVerticesToSubsets(TSubsetHandler& sh,
366 : const ISubsetHandler& srcIndHandler)
367 : {
368 : typedef typename geometry_traits<TElem>::const_iterator iterator;
369 : for(size_t l = 0; l < sh.num_levels(); ++l){
370 : for(int si = 0; si < sh.num_subsets(); ++si){
371 : for(iterator iter = sh.template begin<TElem>(si, l);
372 : iter != sh.template end<TElem>(si, l); ++iter)
373 : {
374 : TElem* e = *iter;
375 : for(size_t i = 0; i < e->num_vertices(); ++i)
376 : {
377 : Vertex* vrt = e->vertex(i);
378 : sh.assign_subset(vrt, srcIndHandler.get_subset_index(vrt));
379 : }
380 : }
381 : }
382 : }
383 : }
384 :
385 : ////////////////////////////////////////////////////////////////////////
386 : template <class TElem, class TSubsetHandler>
387 : void AssignAssociatedEdgesToSubsets(TSubsetHandler& sh,
388 : const ISubsetHandler& srcIndHandler)
389 : {
390 : typedef typename geometry_traits<TElem>::const_iterator iterator;
391 : std::vector<Edge*> vEdges;
392 :
393 : for(size_t l = 0; l < sh.num_levels(); ++l){
394 : for(int si = 0; si < sh.num_subsets(); ++si){
395 : for(iterator iter = sh.template begin<TElem>(si, l);
396 : iter != sh.template end<TElem>(si, l); ++iter)
397 : {
398 : TElem* e = *iter;
399 : CollectEdges(vEdges, *sh.grid(), e);
400 :
401 : for(size_t i = 0; i < vEdges.size(); ++i)
402 : {
403 : Edge* edge = vEdges[i];
404 : sh.assign_subset(edge, srcIndHandler.get_subset_index(edge));
405 : }
406 : }
407 : }
408 : }
409 : }
410 :
411 : ////////////////////////////////////////////////////////////////////////
412 : template <class TElem, class TSubsetHandler>
413 : void AssignAssociatedFacesToSubsets(TSubsetHandler& sh,
414 : const ISubsetHandler& srcIndHandler)
415 : {
416 : typedef typename geometry_traits<TElem>::const_iterator iterator;
417 : std::vector<Face*> vFaces;
418 :
419 : for(size_t l = 0; l < sh.num_levels(); ++l){
420 : for(int si = 0; si < sh.num_subsets(); ++si){
421 : for(iterator iter = sh.template begin<TElem>(si, l);
422 : iter != sh.template end<TElem>(si, l); ++iter)
423 : {
424 : TElem* e = *iter;
425 : CollectFaces(vFaces, *sh.grid(), e);
426 :
427 : for(size_t i = 0; i < vFaces.size(); ++i)
428 : {
429 : Face* f = vFaces[i];
430 : sh.assign_subset(f, srcIndHandler.get_subset_index(f));
431 : }
432 : }
433 : }
434 : }
435 : }
436 :
437 : template <class TElem, class TSubsetHandlerDest, class TSubsetHandlerSrc>
438 : void AssignAssociatedSidesToSubsets(TSubsetHandlerDest& sh,
439 : const TSubsetHandlerSrc& srcIndHandler)
440 : {
441 : typedef typename geometry_traits<TElem>::const_iterator iterator;
442 : typedef typename TElem::lower_dim_base_object Side;
443 : std::vector<Side*> vSides;
444 : Grid& grid = *sh.grid();
445 :
446 : for(size_t l = 0; l < sh.num_levels(); ++l){
447 : for(int si = 0; si < sh.num_subsets(); ++si){
448 : for(iterator iter = sh.template begin<TElem>(si, l);
449 : iter != sh.template end<TElem>(si, l); ++iter)
450 : {
451 : TElem* e = *iter;
452 : CollectAssociated(vSides, grid, e);
453 :
454 : for(size_t i = 0; i < vSides.size(); ++i)
455 : {
456 : Side* s = vSides[i];
457 : sh.assign_subset(s, srcIndHandler.get_subset_index(s));
458 : }
459 : }
460 : }
461 : }
462 : }
463 :
464 : /// helper with with dummy-param for compile-time function selection.
465 : template <class TElem, class TSubsetHandlerDest, class TSubsetHandlerSrc>
466 : void AssignAssociatedLowerDimElemsToSubsets(TSubsetHandlerDest& sh,
467 : const TSubsetHandlerSrc& srcIndHandler,
468 : const Volume&)
469 : {
470 : // we have to find all associated elements of lower dimension.
471 : if(srcIndHandler.template num<Face>() > 0)
472 : AssignAssociatedFacesToSubsets<TElem>(sh, srcIndHandler);
473 : if(srcIndHandler.template num<Edge>() > 0)
474 : AssignAssociatedEdgesToSubsets<TElem>(sh, srcIndHandler);
475 : if(srcIndHandler.template num<Vertex>() > 0)
476 : AssignAssociatedVerticesToSubsets<TElem>(sh, srcIndHandler);
477 : }
478 :
479 : /// helper with with dummy-param for compile-time function selection.
480 : template <class TElem, class TSubsetHandlerDest, class TSubsetHandlerSrc>
481 : void AssignAssociatedLowerDimElemsToSubsets(TSubsetHandlerDest& sh,
482 : const TSubsetHandlerSrc& srcIndHandler,
483 : const Face&)
484 : {
485 : // we have to find all associated elements of lower dimension.
486 : if(srcIndHandler.template num<Edge>() > 0)
487 : AssignAssociatedEdgesToSubsets<TElem>(sh, srcIndHandler);
488 : if(srcIndHandler.template num<Vertex>() > 0)
489 : AssignAssociatedVerticesToSubsets<TElem>(sh, srcIndHandler);
490 : }
491 :
492 : /// helper with with dummy-param for compile-time function selection.
493 : template <class TElem, class TSubsetHandlerDest, class TSubsetHandlerSrc>
494 : void AssignAssociatedLowerDimElemsToSubsets(TSubsetHandlerDest& sh,
495 : const TSubsetHandlerSrc& srcIndHandler,
496 : const Edge&)
497 : {
498 : // we have to find all associated elements of lower dimension.
499 : if(srcIndHandler.template num<Vertex>() > 0)
500 : AssignAssociatedVerticesToSubsets<TElem>(sh, srcIndHandler);
501 : }
502 :
503 : ////////////////////////////////////////////////////////////////////////
504 : template <class TElem, class TSubsetHandlerDest, class TSubsetHandlerSrc>
505 : void AssignAssociatedLowerDimElemsToSubsets(TSubsetHandlerDest& sh,
506 : const TSubsetHandlerSrc& srcIndHandler)
507 : {
508 : AssignAssociatedLowerDimElemsToSubsets<TElem>(sh,
509 : srcIndHandler, TElem());
510 : }
511 :
512 : ////////////////////////////////////////////////////////////////////////
513 : template <typename TBaseObj>
514 : void FindSubsetGroups
515 : (
516 : std::vector<int> & minCondInd,
517 : const std::vector<bool> & isMarked,
518 : const ISubsetHandler & sh,
519 : const NeighborhoodType nbhType
520 : )
521 : {
522 : typedef typename geometry_traits<TBaseObj>::const_iterator elem_iterator;
523 :
524 : UG_ASSERT (((int) isMarked.size ()) == sh.num_subsets (), "FindSubsetGroups: array size mismatch");
525 :
526 : std::vector<TBaseObj*> neighbours;
527 :
528 : // Prepare minCondInd
529 : minCondInd.resize (sh.num_subsets ());
530 : for (size_t si = 0; si < minCondInd.size (); si++)
531 : minCondInd [si] = (isMarked [si])? si : -1;
532 :
533 : // Loop over the subsets:
534 : for (size_t si = 0; si < minCondInd.size (); si++)
535 : {
536 : int min_si;
537 :
538 : // Marked subset?
539 : if ((min_si = minCondInd [si]) < 0)
540 : continue; // no, we do not treat this subset
541 :
542 : // Yes, loop over the elements in the subdomain (in the grid level 0):
543 : GridObjectCollection goc = sh.get_grid_objects_in_subset (si);
544 : elem_iterator e_end = goc.end<TBaseObj> (0);
545 : bool is_empty = true;
546 : for (elem_iterator e_iter = goc.begin<TBaseObj> (0); e_iter != e_end; ++e_iter)
547 : {
548 : is_empty = false;
549 : // Loop over the neighbours:
550 : CollectNeighbors (neighbours, *e_iter, *sh.grid(), nbhType);
551 : for (size_t k = 0; k < neighbours.size (); k++)
552 : {
553 : int min_nbr_si;
554 : int nbr_si = sh.get_subset_index (neighbours [k]);
555 :
556 : if (nbr_si < 0 || nbr_si >= (int) minCondInd.size ())
557 : UG_THROW ("FindSubsetGroups: Illegal neighbour subset index.");
558 : if ((min_nbr_si = minCondInd [nbr_si]) < 0)
559 : continue; // we do not treat this subset
560 :
561 : // Set the same smallest index to both groups of the subsets:
562 : if (min_nbr_si < min_si)
563 : {
564 : for (size_t l = 0; l < minCondInd.size (); l++)
565 : if (minCondInd [l] == min_si)
566 : minCondInd [l] = min_nbr_si;
567 : }
568 : else if (min_nbr_si > min_si)
569 : {
570 : for (size_t l = 0; l < minCondInd.size (); l++)
571 : if (minCondInd [l] == min_nbr_si)
572 : minCondInd [l] = min_si;
573 : }
574 : }
575 : }
576 : if (is_empty) minCondInd [si] = -2;
577 : }
578 : }
579 :
580 : }// end of namespace
581 :
582 : #endif
|
