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 : #include <vector>
34 : #include <stack>
35 : #include <queue>
36 : #include <map>
37 :
38 : #include "subset_util.h"
39 : #include "geom_obj_util/geom_obj_util.h"
40 : #include "polychain_util.h"
41 : #include "orientation_util.h"
42 : #include "subset_color_util.h"
43 :
44 : using namespace std;
45 :
46 : namespace ug
47 : {
48 :
49 : ////////////////////////////////////////////////////////////////////////
50 0 : int GetFirstFreeSubset(const ISubsetHandler& sh)
51 : {
52 0 : for(int i = 0; i < sh.num_subsets(); ++i){
53 0 : if(!(sh.contains_vertices(i) || sh.contains_edges(i)
54 0 : || sh.contains_faces(i) || sh.contains_volumes(i)))
55 : {
56 0 : return i;
57 : }
58 : }
59 :
60 : return sh.num_subsets();
61 : }
62 :
63 :
64 : ////////////////////////////////////////////////////////////////////////
65 : // EraseEmptySubsets
66 : /// Erases all subsets which do not contain any geometric objects
67 0 : void EraseEmptySubsets(ISubsetHandler& sh)
68 : {
69 : int si = 0;
70 0 : while(si < sh.num_subsets()){
71 0 : if(!( sh.contains_vertices(si)
72 0 : || sh.contains_edges(si)
73 0 : || sh.contains_faces(si)
74 0 : || sh.contains_volumes(si)))
75 : {
76 : // the subset is empty
77 0 : sh.erase_subset(si);
78 : }
79 : else
80 0 : ++si;
81 : }
82 0 : }
83 :
84 : ////////////////////////////////////////////////////////////////////////
85 : // AssignGridToSubset
86 0 : void AssignGridToSubset(Grid& g, ISubsetHandler& sh, int subsetInd)
87 : {
88 : UG_ASSERT(&g == sh.grid(), "Specified subset-handler has to operate on the specified grid!");
89 :
90 : sh.assign_subset(g.begin<Vertex>(), g.end<Vertex>(), subsetInd);
91 : sh.assign_subset(g.begin<Edge>(), g.end<Edge>(), subsetInd);
92 : sh.assign_subset(g.begin<Face>(), g.end<Face>(), subsetInd);
93 : sh.assign_subset(g.begin<Volume>(), g.end<Volume>(), subsetInd);
94 0 : }
95 :
96 : ////////////////////////////////////////////////////////////////////////
97 : // AssignSelectionToSubset
98 0 : void AssignSelectionToSubset(ISelector& sel, ISubsetHandler& sh, int subsetInd)
99 : {
100 : UG_ASSERT(sel.grid() == sh.grid(), "Specified selector and subset-handler "
101 : "have to operate on the same grid!");
102 :
103 0 : GridObjectCollection selGoc = sel.get_grid_objects();
104 :
105 0 : for(size_t i = 0; i < selGoc.num_levels(); ++i){
106 : sh.assign_subset(selGoc.begin<Vertex>(i),
107 : selGoc.end<Vertex>(i), subsetInd);
108 : sh.assign_subset(selGoc.begin<Edge>(i),
109 : selGoc.end<Edge>(i), subsetInd);
110 : sh.assign_subset(selGoc.begin<Face>(i),
111 : selGoc.end<Face>(i), subsetInd);
112 : sh.assign_subset(selGoc.begin<Volume>(i),
113 : selGoc.end<Volume>(i), subsetInd);
114 : }
115 0 : }
116 :
117 : ////////////////////////////////////////////////////////////////////////
118 : // AssignInterfaceEdgesToSubsets
119 0 : void AssignFaceInterfaceEdgesToSubsets(Grid& grid, SubsetHandler& sh)
120 : {
121 : //TODO: This algorithm can be improved.
122 : // In the moment it runs with arbitrarily O(numSubs^2*averageNumFacesPerSub)
123 :
124 : // get the first free edge-subset
125 0 : int edgeSubsetIndex = GetMaxSubsetIndex<Edge>(sh) + 1;
126 :
127 : // use this vector to collect adjacent faces
128 : vector<Edge*> vEdges;
129 : vector<Face*> vFaces;
130 :
131 : // iterate through all subsets
132 0 : for(int j = -1; j < sh.num_subsets(); ++j)
133 : {
134 : // find interface-edges to the other subsets.
135 0 : for(int i = j + 1; i < sh.num_subsets(); ++i)
136 : {
137 : // iterate through all faces of subset i and check
138 : // whether some edges are adjacent to subset j.
139 : for(FaceIterator iter = sh.begin<Face>(i);
140 0 : iter != sh.end<Face>(i); ++iter)
141 : {
142 : Face* f = *iter;
143 0 : CollectEdges(vEdges, grid, f);
144 0 : for(uint k = 0; k < vEdges.size(); ++k)
145 : {
146 0 : Edge* e = vEdges[k];
147 :
148 0 : if(sh.get_subset_index(e) == -1)
149 : {
150 0 : CollectFaces(vFaces, grid, e);
151 :
152 0 : if(vFaces.size() == 1)
153 0 : sh.assign_subset(e, edgeSubsetIndex);
154 : else
155 : {
156 0 : for(uint l = 0; l < vFaces.size(); ++l)
157 : {
158 0 : if(sh.get_subset_index(vFaces[l]) == j)
159 : {
160 0 : sh.assign_subset(e, edgeSubsetIndex);
161 : break;
162 : }
163 : }
164 : }
165 : }
166 : }
167 : }
168 :
169 : // if we assigned edges, we will proceed with the next subset
170 0 : if(edgeSubsetIndex < (int)sh.num_subsets())
171 0 : if(sh.num_elements<Edge>(edgeSubsetIndex) > 0)
172 0 : edgeSubsetIndex++;
173 : }
174 : }
175 0 : }
176 :
177 : ////////////////////////////////////////////////////////////////////////
178 : // AssignVolumeInterfaceFacesToSubsets
179 0 : void AssignVolumeInterfaceFacesToSubsets(Grid& grid, SubsetHandler& sh)
180 : {
181 : //TODO: This algorithm can be improved.
182 : // In the moment it runs with arbitrarily O(numSubs^2*averageNumVolumesPerSub)
183 :
184 : // get the first free face-subset
185 0 : int faceSubsetIndex = GetMaxSubsetIndex<Face>(sh) + 1;
186 :
187 : // use this vector to collect adjacent faces and volumes
188 : vector<Face*> vFaces;
189 : vector<Volume*> vVolumes;
190 :
191 : // iterate through all subsets
192 0 : for(int j = -1; j < (int)sh.num_subsets(); ++j)
193 : {
194 : // find interface-edges to the other subsets.
195 0 : for(int i = j + 1; i < (int)sh.num_subsets(); ++i)
196 : {
197 : // iterate through all volumes of subset i and check
198 : // whether some faces are adjacent to subset j.
199 : for(VolumeIterator iter = sh.begin<Volume>(i);
200 0 : iter != sh.end<Volume>(i); ++iter)
201 : {
202 : Volume* v = *iter;
203 0 : CollectFaces(vFaces, grid, v);
204 0 : for(uint k = 0; k < vFaces.size(); ++k)
205 : {
206 0 : Face* f = vFaces[k];
207 :
208 0 : if(sh.get_subset_index(f) == -1)
209 : {
210 0 : CollectVolumes(vVolumes, grid, f);
211 :
212 0 : if(vVolumes.size() == 1)
213 0 : sh.assign_subset(f, faceSubsetIndex);
214 : else
215 : {
216 0 : for(uint l = 0; l < vVolumes.size(); ++l)
217 : {
218 0 : if(sh.get_subset_index(vVolumes[l]) == j)
219 : {
220 0 : sh.assign_subset(f, faceSubsetIndex);
221 : break;
222 : }
223 : }
224 : }
225 : }
226 : }
227 : }
228 :
229 : // if we assigned edges, we will proceed with the next subset
230 0 : if(faceSubsetIndex < (int)sh.num_subsets())
231 0 : if(sh.num_elements<Face>(faceSubsetIndex) > 0)
232 0 : faceSubsetIndex++;
233 : }
234 : }
235 0 : }
236 :
237 :
238 : ////////////////////////////////////////////////////////////////////////
239 : /** Helper method that copies subset-indices based on the subset-dimension
240 : * property of each element. Subsets are only assigned to previously unassigned elements.*/
241 : template <class TIterator>
242 0 : static void CopySubsetFromHigherDimNbr(
243 : ISubsetHandler& sh,
244 : TIterator elemsBegin,
245 : TIterator elemsEnd,
246 : AChar aDimension)
247 : {
248 : typedef typename PtrToValueType<typename TIterator::value_type>::base_type TElem;
249 :
250 0 : UG_COND_THROW(!sh.grid(), "The subset-handler has to operate on a grid!");
251 0 : Grid& grid = *sh.grid();
252 : MultiElementAttachmentAccessor<AChar> aaDim(grid, aDimension, true, true, true, true);
253 :
254 : Grid::edge_traits::secure_container edges;
255 : Grid::face_traits::secure_container faces;
256 : Grid::volume_traits::secure_container vols;
257 :
258 0 : for(TIterator iter = elemsBegin; iter != elemsEnd; ++iter){
259 : TElem* elem = *iter;
260 0 : if(sh.get_subset_index(elem) != -1)
261 0 : continue;
262 :
263 0 : const int dim = aaDim[elem];
264 0 : switch(dim){
265 : case 1:{
266 : // get the edge with dim 1 and with the lowest subset index
267 : int si = sh.num_subsets();
268 : grid.associated_elements(edges, elem);
269 0 : for(size_t i = 0; i < edges.size(); ++i){
270 : Edge* e = edges[i];
271 : int esi = sh.get_subset_index(e);
272 0 : if((aaDim[e] == 1) && (esi < si))
273 : si = esi;
274 : }
275 0 : sh.assign_subset(elem, si);
276 : }break;
277 :
278 : case 2:{
279 : // get the face with dim 2 and with the lowest subset index
280 : int si = sh.num_subsets();
281 : grid.associated_elements(faces, elem);
282 0 : for(size_t i = 0; i < faces.size(); ++i){
283 : Face* e = faces[i];
284 : int esi = sh.get_subset_index(e);
285 0 : if((aaDim[e] == 2) && (esi < si))
286 : si = esi;
287 : }
288 0 : sh.assign_subset(elem, si);
289 : }break;
290 :
291 : case 3:{
292 : // get the volume with dim 3 and with the lowest subset index
293 : int si = sh.num_subsets();
294 : grid.associated_elements(vols, elem);
295 0 : for(size_t i = 0; i < vols.size(); ++i){
296 : Volume* e = vols[i];
297 : int esi = sh.get_subset_index(e);
298 0 : if((aaDim[e] == 3) && (esi < si))
299 : si = esi;
300 : }
301 0 : sh.assign_subset(elem, si);
302 : }break;
303 : }
304 : }
305 0 : }
306 :
307 0 : void CopySubsetIndicesToSides(
308 : ISubsetHandler& sh,
309 : GridObjectCollection goc,
310 : bool toUnassignedOnly)
311 : {
312 0 : UG_COND_THROW(!sh.grid(), "The subset-handler has to operate on a grid!");
313 0 : Grid& grid = *sh.grid();
314 :
315 0 : if(toUnassignedOnly){
316 : AChar aDimension;
317 0 : grid.attach_to_all(aDimension);
318 0 : ComputeLocalSubsetDimensions(sh, aDimension, true);
319 0 : for(int i = 0; i < (int)goc.num_levels(); ++i)
320 0 : CopySubsetFromHigherDimNbr(sh, goc.begin<Vertex>(i), goc.end<Vertex>(i), aDimension);
321 0 : for(int i = 0; i < (int)goc.num_levels(); ++i)
322 0 : CopySubsetFromHigherDimNbr(sh, goc.begin<Edge>(i), goc.end<Edge>(i), aDimension);
323 0 : for(int i = 0; i < (int)goc.num_levels(); ++i)
324 0 : CopySubsetFromHigherDimNbr(sh, goc.begin<Face>(i), goc.end<Face>(i), aDimension);
325 0 : grid.detach_from_all(aDimension);
326 : }
327 : else{
328 0 : for(int i = 0; i < (int)goc.num_levels(); ++i){
329 0 : CopySubsetIndicesToSides(sh, goc.begin<Volume>(i), goc.end<Volume>(i), false);
330 0 : CopySubsetIndicesToSides(sh, goc.begin<Face>(i), goc.end<Face>(i), false);
331 0 : CopySubsetIndicesToSides(sh, goc.begin<Edge>(i), goc.end<Edge>(i), false);
332 : CopySubsetIndicesToSides(sh, goc.begin<Vertex>(i), goc.end<Vertex>(i), false);
333 : }
334 : }
335 0 : }
336 :
337 0 : void CopySubsetIndicesToSides(ISubsetHandler& sh, ISelector& sel,
338 : bool toUnassignedOnly)
339 : {
340 0 : CopySubsetIndicesToSides (sh, sel.get_grid_objects(), toUnassignedOnly);
341 0 : }
342 :
343 :
344 : ////////////////////////////////////////////////////////////////////////
345 0 : void CopySubsetIndicesToSides(ISubsetHandler& sh, bool toUnassignedOnly)
346 : {
347 0 : UG_COND_THROW(!sh.grid(), "The subset-handler has to operate on a grid!");
348 0 : Grid& grid = *sh.grid();
349 0 : CopySubsetIndicesToSides(sh, grid.get_grid_objects(), toUnassignedOnly);
350 0 : }
351 :
352 :
353 : ////////////////////////////////////////////////////////////////////////
354 : // AdjustSubsetsForLgmNg
355 0 : void AdjustSubsetsForLgmNg(Grid& grid, SubsetHandler& sh,
356 : bool keepExistingInterfaceSubsets)
357 : {
358 0 : if(grid.num_volumes() > 0)
359 : {
360 : // the 3d case
361 : vector<Volume*> vVolumes;
362 :
363 : // first make sure that all volumes are assigned to subsets,
364 : // and that no empty volume-subset is between filled ones.
365 0 : MakeSubsetsConsecutive<Volume>(sh);
366 :
367 : // the first index for new subsets.
368 0 : int newVolSubsetIndex = GetMaxSubsetIndex<Volume>(sh) + 1;
369 : // iterate through all volumes and assign each, that is not assigned
370 : // to a subset, to the new one.
371 : for(VolumeIterator iter = grid.volumes_begin();
372 0 : iter != grid.volumes_end(); ++iter)
373 : {
374 0 : if(sh.get_subset_index(*iter) == -1)
375 0 : sh.assign_subset(*iter, newVolSubsetIndex);
376 : }
377 :
378 : // assign all edges to subset -1
379 0 : sh.assign_subset(grid.edges_begin(), grid.edges_end(), -1);
380 :
381 0 : if(keepExistingInterfaceSubsets){
382 : // now we'll assign all faces which are no interface
383 : // elements, to subset -1.
384 : for(FaceIterator iter = grid.faces_begin();
385 0 : iter != grid.faces_end(); ++iter)
386 : {
387 : Face* f = *iter;
388 0 : if(sh.get_subset_index(f) != -1)
389 : {
390 0 : CollectVolumes(vVolumes, grid, f);
391 :
392 0 : if(vVolumes.size() > 1)
393 : {
394 : // check if there are different adjacent volumes.
395 0 : int si = sh.get_subset_index(vVolumes[0]);
396 : bool gotOne = false;
397 0 : for(uint i = 0; i < vVolumes.size(); ++i)
398 : {
399 0 : if(sh.get_subset_index(vVolumes[i]) != si)
400 : {
401 : gotOne = true;
402 : break;
403 : }
404 : }
405 :
406 0 : if(!gotOne)
407 : {
408 : // no, they are all from the same subset.
409 0 : sh.assign_subset(f, -1);
410 : }
411 : }
412 : }
413 : }
414 : }
415 : else{
416 : // since we don't have to keep existing interface subsets we
417 : // can assign all faces to subset -1
418 : sh.assign_subset(grid.faces_begin(), grid.faces_end(), -1);
419 : }
420 :
421 : // make sure that there are no empty face-subsets between filled ones.
422 : // since we may not swap subsets (otherwise the volume-sequence would be destroyed)
423 : // we have to copy the elements.
424 0 : for(int i = 0; i < sh.num_subsets(); ++i){
425 0 : if(sh.num<Face>(i) == 0){
426 : // find the next that has faces
427 : int next = sh.num_subsets();
428 0 : for(int j = i+1; j < sh.num_subsets(); ++j){
429 0 : if(sh.num<Face>(j) > 0){
430 : next = j;
431 : break;
432 : }
433 : }
434 :
435 : // if a filled one has been found, we'll copy the elements
436 0 : if(next < sh.num_subsets()){
437 : // assign all faces in next to subset i
438 : sh.assign_subset(sh.begin<Face>(next), sh.end<Face>(next), i);
439 : }
440 : else{
441 : // we're done
442 : break;
443 : }
444 : }
445 : }
446 :
447 : // now we have to assign the interface faces to subsets.
448 0 : AssignVolumeInterfaceFacesToSubsets(grid, sh);
449 :
450 : // we now have to make sure that all face-subsets are regular manifolds.
451 0 : for(int i = 0; i < sh.num_subsets(); ++i){
452 0 : int firstFree = GetMaxSubsetIndex<Face>(sh) + 1;
453 0 : SplitIrregularManifoldSubset(sh, i, firstFree, true);
454 : }
455 :
456 : // fix orientation of all face subsets
457 0 : for(int i = 0; i < sh.num_subsets(); ++i){
458 0 : if(sh.num<Face>(i) != 0){
459 0 : FixFaceOrientation(grid, sh.begin<Face>(i), sh.end<Face>(i));
460 : }
461 : }
462 0 : }
463 : else
464 : {
465 : // the 2d - case.
466 : vector<Face*> vFaces;
467 :
468 : // first make sure that all faces are assigned to subsets,
469 : // and that no empty face-subset is between filled ones.
470 0 : MakeSubsetsConsecutive<Face>(sh);
471 :
472 : // the first index for new subsets.
473 0 : int newFaceSubsetIndex = GetMaxSubsetIndex<Face>(sh) + 1;
474 :
475 : // iterate through all faces and assign each, that is not assigned
476 : // to a subset, to the new one.
477 : for(FaceIterator iter = grid.faces_begin();
478 0 : iter != grid.faces_end(); ++iter)
479 : {
480 0 : if(sh.get_subset_index(*iter) == -1)
481 0 : sh.assign_subset(*iter, newFaceSubsetIndex);
482 : }
483 :
484 : // now we'll assign all edges which are no interface
485 : // elements, to subset -1.
486 0 : if(keepExistingInterfaceSubsets){
487 : for(EdgeIterator iter = grid.edges_begin();
488 0 : iter != grid.edges_end(); ++iter)
489 : {
490 : Edge* e = *iter;
491 0 : if(sh.get_subset_index(e) != -1)
492 : {
493 0 : CollectFaces(vFaces, grid, e);
494 :
495 0 : if(vFaces.size() > 1)
496 : {
497 : // check if there are different adjacent volumes.
498 0 : int si = sh.get_subset_index(vFaces[0]);
499 : bool gotOne = false;
500 0 : for(uint i = 0; i < vFaces.size(); ++i)
501 : {
502 0 : if(sh.get_subset_index(vFaces[i]) != si)
503 : {
504 : gotOne = true;
505 : break;
506 : }
507 : }
508 :
509 0 : if(!gotOne)
510 : {
511 : // no, they are all from the same subset.
512 0 : sh.assign_subset(e, -1);
513 : }
514 : }
515 : }
516 : }
517 : }
518 : else{
519 : // since we don't have to keep existing interface subsets we can
520 : // assign all edges to subset -1.
521 0 : sh.assign_subset(grid.edges_begin(), grid.edges_end(), -1);
522 : }
523 :
524 : // make sure that there are no empty edge-subsets between filled ones.
525 : // since we may not swap subsets (otherwise the face-sequence would be destroyed)
526 : // we have to copy the elements.
527 0 : for(int i = 0; i < sh.num_subsets(); ++i){
528 0 : if(sh.num<Edge>(i) == 0){
529 : // find the next that has faces
530 : int next = sh.num_subsets();
531 0 : for(int j = i+1; j < sh.num_subsets(); ++j){
532 0 : if(sh.num<Edge>(j) > 0){
533 : next = j;
534 : break;
535 : }
536 : }
537 :
538 : // if a filled one has been found, we'll copy the elements
539 0 : if(next < sh.num_subsets()){
540 : // assign all edges in next to subset i
541 0 : sh.assign_subset(sh.begin<Edge>(next), sh.end<Edge>(next), i);
542 : }
543 : else{
544 : // we're done
545 : break;
546 : }
547 : }
548 : }
549 :
550 : // now we have to assign the interface edges to subsets.
551 0 : AssignFaceInterfaceEdgesToSubsets(grid, sh);
552 :
553 : // make sure that all edge-subsets are regular poly-chains
554 0 : int firstFree = GetMaxSubsetIndex<Edge>(sh) + 1;
555 0 : for(int i = 0; i < sh.num_subsets(); ++i){
556 0 : if(SplitIrregularPolyChain(sh, i, firstFree))
557 0 : ++firstFree;
558 : }
559 :
560 : // Since ug3 does not support loops, we have to remove those.
561 0 : for(int i = 0; i < sh.num_subsets(); ++i){
562 0 : size_t chainType = GetPolyChainType(grid, sh.begin<Edge>(i),
563 0 : sh.end<Edge>(i),
564 0 : IsInSubset(sh, i));
565 0 : if(chainType == PCT_CLOSED){
566 : // since the chain is regular (see loop above) and since it is
567 : // closed, it is enough to simply move the first edge of the
568 : // chain to a new subset.
569 0 : sh.assign_subset(*sh.begin<Edge>(i), firstFree++);
570 : }
571 : }
572 :
573 : // fix orientation
574 0 : for(int i = 0; i < sh.num_subsets(); ++i){
575 0 : if(sh.num<Edge>(i) != 0){
576 0 : FixEdgeOrientation(grid, sh.begin<Edge>(i), sh.end<Edge>(i));
577 : }
578 : }
579 :
580 0 : }
581 0 : }
582 :
583 : ////////////////////////////////////////////////////////////////////////
584 0 : bool SplitIrregularManifoldSubset(SubsetHandler& sh, int srcIndex,
585 : int targetIndex, bool strictSplitting)
586 : {
587 : // Begin with the first face in the subset-handler and find
588 : // associated faces which build a regular manifold.
589 : // if the subset is empty, there's nothing to do.
590 0 : if(sh.empty<Face>(srcIndex))
591 : return false;
592 :
593 : // get the grid behind the subset-handler
594 0 : if(!sh.grid()){
595 : UG_LOG("ERROR in SplitIrregularManifoldSubset: No Grid associated");
596 : UG_LOG(" with the given SubsetHandler.\n");
597 0 : return false;
598 : }
599 :
600 0 : Grid& grid = *sh.grid();
601 :
602 : // edges are required
603 0 : if(!grid.option_is_enabled(FACEOPT_AUTOGENERATE_EDGES)){
604 : UG_LOG("WARNING in SplitIrregularManifoldSubset: Autoenabling FACEOPT_AUTOGENERATE_EDGES.\n");
605 0 : grid.enable_options(FACEOPT_AUTOGENERATE_EDGES);
606 : }
607 :
608 : // here we'll store some temporary results
609 : vector<Edge*> edges;
610 : vector<Face*> faces;
611 :
612 : // the queue is used to keep a list of elements which have to be checked.
613 : // begin with the first face.
614 : queue<Face*> queFaces;
615 0 : queFaces.push(*sh.begin<Face>(srcIndex));
616 :
617 : // We have to mark all faces which have once been pushed to the stack
618 0 : grid.begin_marking();
619 0 : grid.mark(queFaces.front());
620 :
621 : // If we notice that a face would cause a manifold to be irregular,
622 : // we'll mark it and push it to this vector.
623 : vector<Face*> vIrregFaces;
624 :
625 : // This counter tells us whether all faces have been processed.
626 : // If so, the surface is already a regular manifold.
627 : size_t numProcessed = 0;
628 :
629 : // now while there are unprocessed manifold faces
630 0 : while(!queFaces.empty())
631 : {
632 0 : Face* curFace = queFaces.front();
633 : queFaces.pop();
634 0 : ++numProcessed;
635 :
636 : // collect associated faces of associated edges
637 : CollectAssociated(edges, grid, curFace);
638 0 : for(size_t i_edge = 0; i_edge < edges.size(); ++i_edge)
639 : {
640 : // collect associated faces
641 0 : CollectAssociated(faces, grid, edges[i_edge]);
642 :
643 : // we only have to continue if we found exactly one face
644 : // in the same subset other than curFace.
645 : // (Do not check marks here, since we otherwise may find
646 : // unwanted neighbors).
647 0 : Face* nbr = NULL;
648 : size_t numFound = 0;
649 0 : for(size_t i_face = 0; i_face < faces.size(); ++i_face){
650 0 : Face* f = faces[i_face];
651 0 : if((f != curFace) &&
652 0 : ( (strictSplitting && (sh.get_subset_index(f) != -1))
653 0 : || (!strictSplitting && (sh.get_subset_index(f) == srcIndex))))
654 : {
655 : // we found an associated face in the same subset
656 0 : nbr = f;
657 0 : numFound++;
658 : }
659 : }
660 :
661 : // If we found exactly one valid neighbor
662 0 : if(numFound == 1){
663 : // if the neighbor was not already processed we'll schedule it
664 0 : if(!grid.is_marked(nbr)){
665 : grid.mark(nbr);
666 : queFaces.push(nbr);
667 : }
668 : }
669 0 : else if(numFound > 1){
670 : // here we have to make sure that no face which builds an
671 : // irregular manifold with another face can be identified
672 : // as a valid neighbor through another element.
673 : // Thus we'll mark them. Since they wouldn't be added to the
674 : // subset at targetInd later on (since they are marked) we have
675 : // to additionally push them into a vector.
676 0 : for(size_t i_face = 0; i_face < faces.size(); ++i_face){
677 0 : Face* f = faces[i_face];
678 0 : if((f != curFace) && (sh.get_subset_index(f) == srcIndex)){
679 0 : if(!grid.is_marked(f)){
680 : grid.mark(f);
681 0 : vIrregFaces.push_back(f);
682 : }
683 : }
684 : }
685 : }
686 : }
687 : }
688 :
689 : // if all faces of the subset have been processed, then the whole subset is
690 : // a regular manifold.
691 0 : if(numProcessed == sh.num<Face>(srcIndex)){
692 0 : grid.end_marking();
693 : return false;
694 : }
695 :
696 : // if not, we'll iterate over all faces of the subset and assign unmarked ones
697 : // to the subset at targetIndex.
698 : for(FaceIterator iter = sh.begin<Face>(srcIndex);
699 0 : iter != sh.end<Face>(srcIndex);)
700 : {
701 : // Take care with the iterator, since the element may be assigned
702 : // to another subset.
703 : Face* f = *iter;
704 : ++iter;
705 :
706 0 : if(!grid.is_marked(f))
707 0 : sh.assign_subset(f, targetIndex);
708 : }
709 :
710 : // additionally we have to assign all faces in vIrregFaces to the new subset
711 0 : for(size_t i = 0; i < vIrregFaces.size(); ++i)
712 0 : sh.assign_subset(vIrregFaces[i], targetIndex);
713 :
714 : // The subset was split - return true.
715 0 : grid.end_marking();
716 : return true;
717 0 : }
718 :
719 : ////////////////////////////////////////////////////////////////////////
720 : // SeparateFaceSubsetsByNormal
721 : // separates faces by orthogonal axis-aligned normals.
722 0 : void SeparateFaceSubsetsByNormal(Grid& grid, SubsetHandler& sh,
723 : APosition aPos, ANormal* paNorm,
724 : int applyToSubset)
725 : {
726 0 : vector<vector3> vNormals(6);
727 : vNormals[0] = vector3(1, 0, 0);
728 : vNormals[1] = vector3(0, 1, 0);
729 : vNormals[2] = vector3(0, 0, 1);
730 : vNormals[3] = vector3(-1, 0, 0);
731 : vNormals[4] = vector3(0, -1, 0);
732 : vNormals[5] = vector3(0, 0, -1);
733 :
734 0 : SeparateFaceSubsetsByNormal(grid, sh, vNormals, aPos, paNorm,
735 : applyToSubset);
736 0 : }
737 :
738 : ////////////////////////////////////////////////////////////////////////
739 : // SeparateFaceSubsetsByNormal
740 : // separates subset by the given normals.
741 0 : void SeparateFaceSubsetsByNormal(Grid& grid, SubsetHandler& sh,
742 : std::vector<vector3> vNormals,
743 : APosition aPos, ANormal* paNorm,
744 : int applyToSubset)
745 : {
746 : // iterate through all existing face-subsets.
747 : // iterate through all faces of the subset
748 : // find the normal in vNormals that is the closest to the face-normal.
749 : // add the faces to the matching subset.
750 :
751 0 : if(!grid.has_vertex_attachment(aPos))
752 0 : return;
753 :
754 0 : if(vNormals.empty())
755 : return;
756 :
757 : // the position accessor
758 : Grid::VertexAttachmentAccessor<APosition> aaPos(grid, aPos);
759 :
760 : // the normal accessor
761 : vector3 tmpNorm;// this is only required if paNorm == NULL.
762 : Grid::FaceAttachmentAccessor<ANormal> aaNorm;
763 0 : if(paNorm)
764 : {
765 0 : if(grid.has_face_attachment(*paNorm))
766 : aaNorm.access(grid, *paNorm);
767 : }
768 :
769 : // only separate faces of subsets that already existed at the beginning.
770 0 : int numFaceSubsets = GetMaxSubsetIndex<Face>(sh) + 1;
771 : int nextFreeSubset = numFaceSubsets;
772 :
773 : // if a subset was specified, we'll only apply the separation to this subset
774 : int i = 0;
775 : int iMax = numFaceSubsets - 1;
776 0 : if(applyToSubset >= 0){
777 : i = applyToSubset;
778 : iMax = applyToSubset;
779 : }
780 :
781 0 : for(; i <= iMax; ++i)
782 : {
783 : bool firstFace = true;
784 : // this vector holds associated subset indices for each normal.
785 0 : vector<int> vSubsetIndices(vNormals.size(), -1);
786 :
787 : // iterate through the faces of the active subset
788 : FaceIterator iter = sh.begin<Face>(i);
789 0 : while(iter != sh.end<Face>(i))
790 : {
791 : Face* f = *iter;
792 : iter++;// increment iterator here, since it would else be invalidated later on.
793 :
794 : vector3* pN;// pointer to the normal of the face.
795 0 : if(aaNorm.valid())
796 : pN = &aaNorm[f];
797 : else
798 : {
799 : // calculate normal on the fly
800 0 : CalculateNormal(tmpNorm, f, aaPos);
801 : pN = &tmpNorm;
802 : }
803 :
804 : // find the index of the closest normal
805 : int closestInd = -1;
806 : number closestDist = 100;// 4 is the max-distance-square for normized normals.
807 0 : for(uint j = 0; j < (uint)vNormals.size(); ++j)
808 : {
809 0 : number nDist = VecDistanceSq(*pN, vNormals[j]);
810 0 : if(nDist < closestDist)
811 : {
812 : closestDist = nDist;
813 0 : closestInd = j;
814 : }
815 : }
816 :
817 : // get the index of the matching subset.
818 0 : if(vSubsetIndices[closestInd] == -1)
819 : {
820 : // all faces with the same normal as the first face
821 : // will stay in the subset
822 0 : if(firstFace)
823 : {
824 : firstFace = false;
825 0 : vSubsetIndices[closestInd] = i;
826 : }
827 : else
828 : {
829 : // choose an empty subset
830 0 : vSubsetIndices[closestInd] = nextFreeSubset++;
831 : }
832 : }
833 :
834 : // assign the face
835 0 : if(vSubsetIndices[closestInd] != i)
836 0 : sh.assign_subset(f, vSubsetIndices[closestInd]);
837 : }
838 0 : }
839 :
840 : }
841 :
842 : /*
843 : void SeparateVolumesByFaceSubsets(Grid& grid, SubsetHandler& sh)
844 : {
845 : // first we'll assign all volumes to subset -1.
846 : sh.assign_subset(grid.volumes_begin(), grid.volumes_end(), -1);
847 :
848 : // we'll keep all unassigned volumes in a selector.
849 : Selector sel(grid);
850 : sel.select(grid.volumes_begin(), grid.volumes_end());
851 :
852 : // those vectors will be used to gather element neighbours.
853 : vector<Face*> vFaces;
854 : vector<Volume*> vVolumes;
855 :
856 : // this stack contains all volumes that we still have to check for neighbours.
857 : stack<Volume*> stkVols;
858 :
859 : // will be used to store sides of volumes
860 : FaceDescriptor fd;
861 :
862 : // now - while there are unassigned volumes.
863 : int subsetIndex = 0;
864 : while(!sel.empty())
865 : {
866 : // choose the volume with which we want to start
867 : // TODO: if material-points are supplied, this should be the
868 : // the volume that contains the i-th material point.
869 : stkVols.push(*sel.begin<Volume>());
870 : while(!stkVols.empty())
871 : {
872 : Volume* v = stkVols.top();
873 : stkVols.pop();
874 : // if the volume is unselected it has already been processed.
875 : if(!sel.is_selected(v))
876 : continue;
877 : sel.deselect(v);
878 :
879 : // assign v to its new subset
880 : sh.assign_subset(v, subsetIndex);
881 :
882 : // check neighbour-volumes, whether they belong to the same subset.
883 : // iterate through the sides of the volume
884 : for(uint i = 0; i < v->num_faces(); ++i)
885 : {
886 : v->face(i, fd);
887 :
888 : // get the corresponding face
889 : Face* f = grid.get_face(fd);
890 :
891 : bool bSeparator = false;
892 : // if it belongs to a subset other that -1, it is a separator.
893 : if(f)
894 : bSeparator = (sh.get_subset_index(f) != -1);
895 :
896 : // if fd is not corresponding to a separator, we'll add all connected volumes
897 : if(!bSeparator)
898 : {
899 : CollectVolumes(vVolumes, grid, fd);
900 :
901 : // add all volumes that are still selected (v is not selected anymore).
902 : for(uint j = 0; j < vVolumes.size(); ++j)
903 : {
904 : if(sel.is_selected(vVolumes[j]))
905 : stkVols.push(vVolumes[j]);
906 : }
907 : }
908 : }
909 : }
910 : // the stack is empty. increase subset index.
911 : subsetIndex++;
912 : }
913 : }
914 : */
915 0 : void AssignRegionToSubset(Grid& grid, ISubsetHandler& shVolsOut,
916 : const ISubsetHandler& shFaces,
917 : Volume* proxyVol, int newSubsetIndex)
918 : {
919 : // we'll utilize a stack to gather all volumes that have to
920 : // be examined.
921 : stack<Volume*> stkVols;
922 : stkVols.push(proxyVol);
923 :
924 : // this vector will be used to collect neighbours
925 : vector<Volume*> vVols;
926 :
927 : // while there are volumes on the stack we'll go on
928 0 : while(!stkVols.empty()){
929 0 : Volume* v = stkVols.top();
930 : stkVols.pop();
931 :
932 0 : shVolsOut.assign_subset(v, newSubsetIndex);
933 :
934 : // check all neighbours and decide whether they have to be
935 : // pushed to the stack.
936 : // First we'll have to check for each side of v whether we may traverse it
937 0 : for(uint i = 0; i < v->num_faces(); ++i){
938 0 : Face* f = grid.get_face(v, i);
939 0 : if(f){
940 : // check whether f lies in a subset
941 : // if so we may not traverse it and we'll continue with the next face.
942 0 : if(shFaces.get_subset_index(f) != -1)
943 0 : continue;
944 :
945 : // we may traverse it. get the neighbour-volumes
946 0 : CollectVolumes(vVols, grid, f);
947 : }
948 : else{
949 : // no associated face existed. We'll use GetNeighbours to find
950 : // neighbouring volumes. Since no face can separate them, they
951 : // belong to the region, too.
952 0 : GetNeighbours(vVols, grid, v, i);
953 : }
954 :
955 : // add the volumes in vVols to the stack - but only if they don't
956 : // already belong to the region.
957 0 : for(size_t j = 0; j < vVols.size(); ++j){
958 0 : if(vVols[j] != v){
959 0 : if(shVolsOut.get_subset_index(vVols[j]) != newSubsetIndex)
960 : stkVols.push(vVols[j]);
961 : }
962 : }
963 : }
964 : }
965 0 : }
966 :
967 : ////////////////////////////////////////////////////////////////////////
968 0 : bool SeparateRegions(Grid& grid, ISubsetHandler& shVolsOut,
969 : const ISubsetHandler& shFaces,
970 : const MarkerPointManager& mpm,
971 : int firstSubsetIndex)
972 : {
973 0 : if(!grid.has_vertex_attachment(aPosition))
974 : return false;
975 :
976 : Grid::VertexAttachmentAccessor<APosition> aaPos(grid, aPosition);
977 :
978 0 : for(size_t i = 0; i < mpm.num_markers(); ++i)
979 : {
980 0 : const vector3& pos = mpm.get_marker(i).pos;
981 :
982 : //TODO: extend to volumes in general. Add a PointIsInsideVolume method.
983 : // find the tetrahedron that contains the marker point
984 : for(TetrahedronIterator iter = grid.begin<Tetrahedron>();
985 0 : iter != grid.end<Tetrahedron>(); ++iter)
986 : {
987 0 : if(PointIsInsideTetrahedron(pos, *iter, aaPos)){
988 : // assign region to subset
989 0 : int si = firstSubsetIndex + i;
990 0 : AssignRegionToSubset(grid, shVolsOut, shFaces,
991 : *iter, si);
992 : // set subset name
993 0 : shVolsOut.subset_info(si).name = mpm.get_marker(i).name;
994 : // we don't have to check other volumes.
995 : break;
996 : }
997 : }
998 : }
999 :
1000 : // done
1001 : return true;
1002 : }
1003 :
1004 0 : void AssignInnerAndBoundarySubsets(Grid& grid, ISubsetHandler& shOut,
1005 : int inSubset, int bndSubset)
1006 : {
1007 0 : if(grid.num<Volume>() > 0){
1008 : // assign volumes to inSubset
1009 : for(VolumeIterator iter = grid.begin<Volume>();
1010 0 : iter != grid.end<Volume>(); ++iter)
1011 : {
1012 0 : if(shOut.get_subset_index(*iter) == -1)
1013 0 : shOut.assign_subset(*iter, inSubset);
1014 : }
1015 :
1016 : vector<Edge*> vEdges;
1017 : // assign volume-boundary elements to bndSubset
1018 0 : for(FaceIterator iter = grid.begin<Face>(); iter != grid.end<Face>(); ++iter)
1019 : {
1020 : Face* f = *iter;
1021 0 : if(shOut.get_subset_index(f) == -1){
1022 0 : if(IsVolumeBoundaryFace(grid, f)){
1023 : // assign the face to the boundary subset.
1024 0 : shOut.assign_subset(f, bndSubset);
1025 :
1026 : // assign associated vertices and edges to the boundary subset, too.
1027 0 : for(size_t i = 0; i < f->num_vertices(); ++i){
1028 0 : if(shOut.get_subset_index(f->vertex(i)) == -1)
1029 0 : shOut.assign_subset(f->vertex(i), bndSubset);
1030 : }
1031 :
1032 0 : CollectEdges(vEdges, grid, f);
1033 0 : for(size_t i = 0; i < vEdges.size(); ++i){
1034 0 : if(shOut.get_subset_index(vEdges[i]) == -1)
1035 0 : shOut.assign_subset(vEdges[i], bndSubset);
1036 : }
1037 : }
1038 : else{
1039 : // assing the face to the inner subset
1040 0 : shOut.assign_subset(f, inSubset);
1041 : }
1042 : }
1043 : }
1044 0 : }
1045 : else{
1046 : for(FaceIterator iter = grid.begin<Face>();
1047 0 : iter != grid.end<Face>(); ++iter)
1048 : {
1049 0 : if(shOut.get_subset_index(*iter) == -1)
1050 0 : shOut.assign_subset(*iter, inSubset);
1051 : }
1052 : }
1053 :
1054 : // assign edges and vertices
1055 : for(EdgeIterator iter = grid.begin<Edge>();
1056 0 : iter != grid.end<Edge>(); ++iter)
1057 : {
1058 : Edge* e = *iter;
1059 0 : if(shOut.get_subset_index(e) == -1){
1060 0 : if(IsBoundaryEdge2D(grid, e))
1061 : {
1062 0 : shOut.assign_subset(e, bndSubset);
1063 0 : for(size_t i = 0; i < 2; ++i)
1064 : {
1065 0 : if(shOut.get_subset_index(e->vertex(i)) == -1)
1066 0 : shOut.assign_subset(e->vertex(i), bndSubset);
1067 : }
1068 : }
1069 : else
1070 0 : shOut.assign_subset(e, inSubset);
1071 : }
1072 : }
1073 :
1074 : // assign unassigned vertices
1075 : for(VertexIterator iter = grid.begin<Vertex>();
1076 0 : iter != grid.end<Vertex>(); ++iter)
1077 : {
1078 0 : if(shOut.get_subset_index(*iter) == -1)
1079 0 : shOut.assign_subset(*iter, inSubset);
1080 : }
1081 0 : }
1082 :
1083 :
1084 : ////////////////////////////////////////////////////////////////////////
1085 0 : vector3 GetColorFromStandardPalette(int index)
1086 : {
1087 0 : return GetColorFromDefaultPalette (index);
1088 : }
1089 :
1090 : ////////////////////////////////////////////////////////////////////////
1091 0 : void AssignSubsetColors(ISubsetHandler& sh)
1092 : {
1093 0 : AssignDefaultSubsetColors (sh);
1094 0 : }
1095 :
1096 :
1097 : ////////////////////////////////////////////////////////////////////////
1098 : template <class TElem>
1099 0 : void AssignSidesToSubsets(ISubsetHandler& sh, ISelector* psel)
1100 : {
1101 : typedef typename TElem::lower_dim_base_object Side;
1102 : typedef typename geometry_traits<Side>::iterator SideIter;
1103 : typedef std::vector<int> IntString;
1104 :
1105 : // access the grid on which sh operates.
1106 0 : if(!sh.grid())
1107 0 : return;
1108 :
1109 0 : Grid& grid = *sh.grid();
1110 :
1111 : // we'll use those marks to check whether a subset has already been
1112 : // processed in an iteration. A subset is considered to be marked, if
1113 : // marks[subInd] == curMarker.
1114 : // Note that curMarker is increased in each iteration.
1115 : // Note also, that marks has to be resized whenever a new subset is added.
1116 : size_t curMark = 1;
1117 0 : vector<size_t> marks(sh.num_subsets(), 0);
1118 :
1119 : // here we'll check whether a subset already contains elements with
1120 : // a given neighborhood subset-constellation.
1121 : // Note that the key will be calculated from the subset-constellation.
1122 : // we're using double to allow more values
1123 : map<IntString, int> subsetMap;
1124 :
1125 : // this vector will be used to collect all associated subsets of an element.
1126 : // we'll use the subset-mark mechanism to avoid duplicate entries.
1127 : IntString skey;
1128 :
1129 : // used to collect neighbors of type TElem
1130 : vector<TElem*> elems;
1131 :
1132 : // Now iterate over all elements of type Side
1133 : for(SideIter iterSide = grid.begin<Side>();
1134 0 : iterSide != grid.end<Side>(); ++iterSide)
1135 : {
1136 : Side* side = *iterSide;
1137 :
1138 0 : if(sh.get_subset_index(side) != -1)
1139 : continue;
1140 :
1141 : // check whether we are working on a selected side
1142 : bool isSelected = false;
1143 0 : if(psel && psel->is_selected(side))
1144 : isSelected = true;
1145 :
1146 : // increase curMark
1147 0 : ++curMark;
1148 :
1149 : // collect all associated elements of side
1150 : CollectAssociated(elems, grid, side);
1151 :
1152 : // collect the subsets in which they lie. Note that we won't push the
1153 : // same subset index twice to skey.
1154 : skey.clear();
1155 :
1156 : // this subset index shall be used if the key is new
1157 0 : int useSubsetInd = GetFirstFreeSubset(sh);
1158 :
1159 0 : for(size_t i_elems = 0; i_elems < elems.size(); ++i_elems){
1160 : bool nbrIsSelected = false;
1161 0 : if(psel && psel->is_selected(elems[i_elems]))
1162 : nbrIsSelected = true;
1163 :
1164 : // either the element is not selected or it is selected and its neighbor
1165 : // is selected, too.
1166 0 : if(!isSelected || nbrIsSelected){
1167 0 : int si = sh.get_subset_index(elems[i_elems]);
1168 : // check if the subset is marked
1169 0 : if(marks[si] != curMark){
1170 : // no - mark it and add si to subsets
1171 0 : marks[si] = curMark;
1172 0 : skey.push_back(si);
1173 : }
1174 : }
1175 : }
1176 :
1177 0 : if(skey.size() == 1){
1178 0 : if(elems.size() > 1){
1179 : // we've got an inner element. Assign the subset-index
1180 : // of those inner neighbor elements.
1181 0 : useSubsetInd = skey[0];
1182 : }
1183 : else{
1184 : // An outer interface element.
1185 : // add -2 to skey to make sure that they go into a separate subset
1186 : // (-1 might occur as a normal subset index)
1187 0 : skey.push_back(-2);
1188 : }
1189 : }
1190 :
1191 : // The subset indices in skey have to be sorted before they can be
1192 : // used as unique key
1193 0 : sort(skey.begin(), skey.end());
1194 :
1195 : // get the subset index from the subsetMap, if one exists.
1196 : map<IntString, int>::iterator tMapIter = subsetMap.find(skey);
1197 0 : if(tMapIter != subsetMap.end()){
1198 : // got one
1199 0 : sh.assign_subset(side, tMapIter->second);
1200 : }
1201 : else{
1202 : // create a new entry and resize the marks array
1203 0 : subsetMap[skey] = useSubsetInd;
1204 0 : sh.assign_subset(side, useSubsetInd);
1205 0 : if(sh.num_subsets() > (int)marks.size())
1206 0 : marks.resize(sh.num_subsets(), 0);
1207 : }
1208 : }
1209 0 : }
1210 :
1211 : // template specialization
1212 : template void AssignSidesToSubsets<Edge>(ISubsetHandler&, ISelector*);
1213 : template void AssignSidesToSubsets<Face>(ISubsetHandler&, ISelector*);
1214 : template void AssignSidesToSubsets<Volume>(ISubsetHandler&, ISelector*);
1215 :
1216 :
1217 0 : void AssignSubsetsByElementType(ISubsetHandler& sh)
1218 : {
1219 0 : if(!sh.grid())
1220 : return;
1221 :
1222 0 : AssignSubsetsByElementType(sh, sh.grid()->get_grid_objects());
1223 : }
1224 :
1225 :
1226 0 : void AssignSubsetsByElementType(ISubsetHandler& sh, GridObjectCollection g)
1227 : {
1228 : int subsetInd = 0;
1229 :
1230 0 : if(g.num<RegularVertex>() > 0){
1231 : sh.assign_subset(g.begin<RegularVertex>(), g.end<RegularVertex>(), subsetInd);
1232 0 : sh.subset_info(subsetInd++).name = "RegularVertex";
1233 : }
1234 :
1235 0 : if(g.num<ConstrainedVertex>() > 0){
1236 : sh.assign_subset(g.begin<ConstrainedVertex>(), g.end<ConstrainedVertex>(), subsetInd);
1237 0 : sh.subset_info(subsetInd++).name = "HangingVertex";
1238 : }
1239 :
1240 0 : if(g.num<RegularEdge>() > 0){
1241 : sh.assign_subset(g.begin<RegularEdge>(), g.end<RegularEdge>(), subsetInd);
1242 0 : sh.subset_info(subsetInd++).name = "RegularEdge";
1243 : }
1244 :
1245 0 : if(g.num<ConstrainingEdge>() > 0){
1246 : sh.assign_subset(g.begin<ConstrainingEdge>(), g.end<ConstrainingEdge>(), subsetInd);
1247 0 : sh.subset_info(subsetInd++).name = "ConstrainingEdge";
1248 : }
1249 :
1250 0 : if(g.num<ConstrainedEdge>() > 0){
1251 : sh.assign_subset(g.begin<ConstrainedEdge>(), g.end<ConstrainedEdge>(), subsetInd);
1252 0 : sh.subset_info(subsetInd++).name = "ConstrainedEdge";
1253 : }
1254 :
1255 0 : if(g.num<Triangle>() > 0){
1256 : sh.assign_subset(g.begin<Triangle>(), g.end<Triangle>(), subsetInd);
1257 0 : sh.subset_info(subsetInd++).name = "Triangle";
1258 : }
1259 :
1260 0 : if(g.num<ConstrainingTriangle>() > 0){
1261 : sh.assign_subset(g.begin<ConstrainingTriangle>(), g.end<ConstrainingTriangle>(), subsetInd);
1262 0 : sh.subset_info(subsetInd++).name = "ConstrainingTriangle";
1263 : }
1264 :
1265 0 : if(g.num<ConstrainedTriangle>() > 0){
1266 : sh.assign_subset(g.begin<ConstrainedTriangle>(), g.end<ConstrainedTriangle>(), subsetInd);
1267 0 : sh.subset_info(subsetInd++).name = "ConstrainedTriangle";
1268 : }
1269 :
1270 0 : if(g.num<Quadrilateral>() > 0){
1271 : sh.assign_subset(g.begin<Quadrilateral>(), g.end<Quadrilateral>(), subsetInd);
1272 0 : sh.subset_info(subsetInd++).name = "Quadrilateral";
1273 : }
1274 :
1275 0 : if(g.num<ConstrainingQuadrilateral>() > 0){
1276 : sh.assign_subset(g.begin<ConstrainingQuadrilateral>(), g.end<ConstrainingQuadrilateral>(), subsetInd);
1277 0 : sh.subset_info(subsetInd++).name = "ConstrainingQuadrilateral";
1278 : }
1279 :
1280 0 : if(g.num<ConstrainedQuadrilateral>() > 0){
1281 : sh.assign_subset(g.begin<ConstrainedQuadrilateral>(), g.end<ConstrainedQuadrilateral>(), subsetInd);
1282 0 : sh.subset_info(subsetInd++).name = "ConstrainedQuadrilateral";
1283 : }
1284 :
1285 0 : if(g.num<Tetrahedron>() > 0){
1286 : sh.assign_subset(g.begin<Tetrahedron>(), g.end<Tetrahedron>(), subsetInd);
1287 0 : sh.subset_info(subsetInd++).name = "Tetrahedron";
1288 : }
1289 :
1290 0 : if(g.num<Pyramid>() > 0){
1291 : sh.assign_subset(g.begin<Pyramid>(), g.end<Pyramid>(), subsetInd);
1292 0 : sh.subset_info(subsetInd++).name = "Pyramid";
1293 : }
1294 :
1295 0 : if(g.num<Prism>() > 0){
1296 : sh.assign_subset(g.begin<Prism>(), g.end<Prism>(), subsetInd);
1297 0 : sh.subset_info(subsetInd++).name = "Prism";
1298 : }
1299 :
1300 0 : if(g.num<Hexahedron>() > 0){
1301 : sh.assign_subset(g.begin<Hexahedron>(), g.end<Hexahedron>(), subsetInd);
1302 0 : sh.subset_info(subsetInd++).name = "Hexahedron";
1303 : }
1304 :
1305 0 : if(g.num<Octahedron>() > 0){
1306 : sh.assign_subset(g.begin<Octahedron>(), g.end<Octahedron>(), subsetInd);
1307 0 : sh.subset_info(subsetInd++).name = "Octahedron";
1308 : }
1309 0 : }
1310 :
1311 : template <class TElem>
1312 : static char GetSmallestLocalSubsetDimension(
1313 : typename Grid::traits<TElem>::secure_container& nbrs,
1314 : MultiElementAttachmentAccessor<AChar>& aaDim)
1315 : {
1316 : int d = -1;
1317 0 : for(size_t i = 0; i < nbrs.size(); ++i){
1318 : TElem* e = nbrs[i];
1319 0 : char ed = aaDim[e];
1320 0 : if((ed != -1 ) && ((d == -1) || (ed < d)))
1321 0 : d = ed;
1322 : }
1323 0 : return d;
1324 : }
1325 :
1326 : template <class TElem>
1327 0 : static bool NbrIsInSubset(
1328 : ISubsetHandler& sh,
1329 : typename Grid::traits<TElem>::secure_container& nbrs,
1330 : int si)
1331 : {
1332 0 : for(size_t i = 0; i < nbrs.size(); ++i){
1333 0 : if(sh.get_subset_index(nbrs[i]) == si)
1334 : return true;
1335 : }
1336 : return false;
1337 : }
1338 :
1339 0 : void ComputeLocalSubsetDimensions(
1340 : ISubsetHandler& sh,
1341 : AChar aDimension,
1342 : bool includeUnassigned)
1343 : {
1344 0 : UG_COND_THROW(!sh.grid(), "The subset-handler has to operate on a grid!");
1345 0 : Grid& grid = *sh.grid();
1346 0 : grid.attach_to_all(aDimension);
1347 : MultiElementAttachmentAccessor<AChar> aaDim(grid, aDimension, true, true, true, true);
1348 :
1349 :
1350 0 : for(VolumeIterator ivol = grid.begin<Volume>(); ivol != grid.end<Volume>(); ++ivol){
1351 0 : if(sh.get_subset_index(*ivol) == -1)
1352 0 : aaDim[*ivol] = -1;
1353 : else
1354 0 : aaDim[*ivol] = 3;
1355 : }
1356 :
1357 : Grid::volume_traits::secure_container vols;
1358 0 : for(FaceIterator iface = grid.begin<Face>(); iface != grid.end<Face>(); ++iface){
1359 : Face* f = *iface;
1360 : int si = sh.get_subset_index(f);
1361 0 : if(si == -1){
1362 0 : if(includeUnassigned){
1363 : grid.associated_elements(vols, f);
1364 0 : aaDim[f] = GetSmallestLocalSubsetDimension<Volume>(vols, aaDim);
1365 : }
1366 : else
1367 0 : aaDim[f] = -1;
1368 : }
1369 : else{
1370 0 : aaDim[f] = 2;
1371 : grid.associated_elements(vols, f);
1372 0 : if(NbrIsInSubset<Volume>(sh, vols, si))
1373 0 : aaDim[f] = 3;
1374 : }
1375 : }
1376 :
1377 : Grid::face_traits::secure_container faces;
1378 0 : for(EdgeIterator iedge = grid.begin<Edge>(); iedge != grid.end<Edge>(); ++iedge){
1379 : Edge* e = *iedge;
1380 : int si = sh.get_subset_index(e);
1381 0 : if(si == -1){
1382 0 : if(includeUnassigned){
1383 : grid.associated_elements(faces, e);
1384 0 : aaDim[e] = GetSmallestLocalSubsetDimension<Face>(faces, aaDim);
1385 : }
1386 : else
1387 0 : aaDim[e] = -1;
1388 : }
1389 : else{
1390 0 : aaDim[e] = 1;
1391 : grid.associated_elements(vols, e);
1392 0 : if(NbrIsInSubset<Volume>(sh, vols, si))
1393 0 : aaDim[e] = 3;
1394 : else{
1395 : grid.associated_elements(faces, e);
1396 0 : if(NbrIsInSubset<Face>(sh, faces, si))
1397 0 : aaDim[e] = 2;
1398 : }
1399 : }
1400 : }
1401 :
1402 : Grid::edge_traits::secure_container edges;
1403 0 : for(VertexIterator ivrt = grid.begin<Vertex>(); ivrt != grid.end<Vertex>(); ++ivrt){
1404 : Vertex* v = *ivrt;
1405 : int si = sh.get_subset_index(v);
1406 0 : if(si == -1){
1407 0 : if(includeUnassigned){
1408 : grid.associated_elements(edges, v);
1409 0 : aaDim[v] = GetSmallestLocalSubsetDimension<Edge>(edges, aaDim);
1410 : }
1411 : else
1412 0 : aaDim[v] = -1;
1413 : }
1414 : else{
1415 0 : aaDim[v] = 0;
1416 : grid.associated_elements(vols, v);
1417 0 : if(NbrIsInSubset<Volume>(sh, vols, si))
1418 0 : aaDim[v] = 3;
1419 : else{
1420 : grid.associated_elements(faces, v);
1421 0 : if(NbrIsInSubset<Face>(sh, faces, si))
1422 0 : aaDim[v] = 2;
1423 : else{
1424 : grid.associated_elements(edges, v);
1425 0 : if(NbrIsInSubset<Edge>(sh, edges, si))
1426 0 : aaDim[v] = 1;
1427 : }
1428 : }
1429 : }
1430 : }
1431 0 : }
1432 :
1433 : }// end of namespace
1434 :
|
