Line data Source code
1 : /*
2 : * Copyright (c) 2012-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Martin Scherer
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 PERIODIC_IDENTIFIER_IMPL_HPP_
34 : #define PERIODIC_IDENTIFIER_IMPL_HPP_
35 :
36 : // include declarations
37 : #include "periodic_boundary_manager.h"
38 : #include "lib_grid/algorithms/debug_util.h"
39 : #include "lib_grid/grid_objects/grid_dim_traits.h"
40 : #include "common/assert.h"
41 : #include "common/error.h"
42 : #include "pcl/pcl_base.h"
43 :
44 : #include <boost/mpl/map.hpp>
45 : #include <boost/mpl/at.hpp>
46 :
47 : #include <algorithm>
48 :
49 : namespace ug {
50 :
51 : template <class TAAPos>
52 : template <class TElem>
53 0 : bool ParallelShiftIdentifier<TAAPos>::match_impl(TElem* e1, TElem* e2) const {
54 0 : if (e1 == e2)
55 : return false;
56 :
57 0 : AttachmentType c1 = CalculateCenter(e1, m_aaPos),
58 0 : c2 = CalculateCenter(e2, m_aaPos), diff, error;
59 : bool result = false;
60 :
61 : VecSubtract(diff, c1, c2);
62 : VecSubtract(error, diff, m_shift);
63 : number len = VecLengthSq(error);
64 0 : if (std::abs(len) < 10E-8)
65 : result = true;
66 : else // check for opposite shift
67 : {
68 : VecSubtract(error, diff, m_shift_opposite);
69 : len = VecLengthSq(error);
70 0 : if (std::abs(len) < 10E-8)
71 : result = true;
72 : }
73 :
74 : return result;
75 : }
76 :
77 : template <class TElem>
78 0 : void PeriodicBoundaryManager::identify(TElem* e1, TElem* e2,
79 : IIdentifier& ident) {
80 : typedef typename
81 : Grid::traits<typename TElem::side>::secure_container container;
82 :
83 : // determine masters
84 0 : TElem* m1 = master(e1);
85 0 : TElem* m2 = master(e2);
86 :
87 0 : if (m1 || m2) {
88 : // if m1 == m2, there's nothing to do
89 0 : if (m1 != m2) {
90 0 : if (!m1) { // m2 is master
91 0 : make_slave(group(m2), e1);
92 0 : } else if (!m2) { // m1 is master
93 0 : make_slave(group(m1), e2);
94 : } else if (m1 && m2) { // both are distinct masters
95 0 : merge_groups(group(m1), group(m2));
96 : } else {
97 : UG_THROW("should never get here")
98 : }
99 : }
100 : } else {
101 : // create new group with e1 as master
102 0 : Group<TElem>* g = new Group<TElem>(e1);
103 : // set group of master
104 0 : set_group(g, e1);
105 : // make e2 slave
106 0 : make_slave(g, e2);
107 : }
108 :
109 : // while elements have sides, recursively identify sub type elements
110 : if (TElem::HAS_SIDES) {
111 : container sides1, sides2;
112 : // collect sides and identify them
113 0 : m_pGrid->associated_elements<TElem>(sides1, e1);
114 0 : m_pGrid->associated_elements<TElem>(sides2, e2);
115 0 : for (size_t i = 0; i < sides1.size(); ++i) {
116 0 : for (size_t j = 0; j < sides2.size(); ++j) {
117 0 : if(ident.match(sides1[i], sides2[j])) {
118 0 : identify(sides1[i], sides2[j], ident);
119 : }
120 : }
121 : }
122 : }
123 0 : }
124 :
125 : // is considered periodic if it is a master or a slave
126 : template <class TElem>
127 0 : bool PeriodicBoundaryManager::is_periodic(TElem* e) const {
128 0 : if(e)
129 0 : return get_periodic_status_accessor<TElem>()[e] != P_NOT_PERIODIC;
130 : else
131 0 : UG_THROW("null pointer is never periodic.");
132 : }
133 :
134 : // gets master of e, may be null
135 : template <class TElem>
136 0 : TElem* PeriodicBoundaryManager::master(TElem* e) const {
137 0 : if (group(e)) {
138 0 : return group(e)->m_master;
139 : }
140 : return NULL;
141 : }
142 :
143 : // gets slaves of e
144 : template <class TElem>
145 : typename PeriodicBoundaryManager::Group<TElem>::SlaveContainer*
146 0 : PeriodicBoundaryManager::slaves(TElem* e) const {
147 0 : if (group(e)) {
148 0 : return &group(e)->get_slaves();
149 : }
150 : return NULL;
151 : }
152 :
153 : template <class TElem>
154 0 : void PeriodicBoundaryManager::print_identification() const {
155 : typedef typename ElementStorage<TElem>::SectionContainer::iterator Iterator;
156 : typedef typename Group<TElem>::SlaveIterator SlaveIter;
157 :
158 0 : for (Iterator elem = m_pGrid->begin<TElem>(); elem != m_pGrid->end<TElem>();
159 : ++elem) {
160 : // log masters and their slaves
161 0 : if (!(master(*elem) == *elem))
162 0 : continue;
163 :
164 : Group<TElem>* g = group(*elem);
165 : UG_ASSERT(g, "group not valid")
166 0 : UG_LOG("group of " << (*elem)->reference_object_id() << "\tlevel: " <<
167 : m_pGrid->get_level(*elem) << "\tmaster: " <<
168 : GetGridObjectCenter(*m_pGrid, g->m_master) << "\tslaves: ");
169 0 : for (SlaveIter slave = g->get_slaves().begin();
170 0 : slave != g->get_slaves().end(); ++slave) {
171 0 : TElem* e = *slave;
172 0 : UG_LOG(GetGridObjectCenter(*m_pGrid, e) << ", ")
173 : UG_ASSERT(m_pGrid->get_level(*elem) == m_pGrid->get_level(e),
174 : "wrong level in group")
175 : }
176 : UG_LOG(std::endl)
177 : }
178 0 : }
179 :
180 : /**
181 : * TParent should be only of type Vertex, Edge, Face.
182 : * Volumes are not meant to be periodic.
183 : *
184 : * If replacesParent is true, e is meant to replace pParent
185 : */
186 : template <class TElem>
187 0 : void PeriodicBoundaryManager::replace_parent(TElem* e, TElem* pParent) {
188 :
189 0 : if (is_master(pParent)) {
190 : // if parent is master, set newly created item as new master
191 : Group<TElem>* g = group(pParent);
192 0 : g->m_master = e;
193 0 : set_group(g, e);
194 0 : } else if(is_slave(pParent)) { // slave
195 : // iterate over parent-groups slave list and replace pointer to parent
196 0 : typename Group<TElem>::SlaveContainer* slaveCon = slaves(pParent);
197 0 : for(typename Group<TElem>::SlaveContainer::iterator i = slaveCon->begin();
198 0 : i != slaveCon->end(); ++i)
199 : {
200 0 : if(*i == pParent)
201 0 : *i = e;
202 : }
203 0 : set_group(group(pParent), e);
204 : }
205 0 : }
206 :
207 : template <class TElem, class TParent>
208 0 : void PeriodicBoundaryManager::handle_creation(TElem* e, TParent* pParent) {
209 :
210 : typedef typename Group<TParent>::SlaveContainer ParentSlaveContainer;
211 :
212 0 : if (!pParent)
213 : return;
214 :
215 0 : if (!is_periodic(pParent))
216 : return;
217 :
218 : UG_ASSERT(m_pGrid,
219 : "PeriodicBoundaryManager has to operate on a grid. No grid assigned.");
220 0 : MultiGrid& mg = *m_pGrid;
221 :
222 0 : mg.begin_marking();
223 : if(TElem::BASE_OBJECT_ID != VERTEX){
224 : Grid::vertex_traits::secure_container vrts;
225 : mg.associated_elements(vrts, e);
226 0 : for(size_t i = 0; i < vrts.size(); ++i){
227 : Group<Vertex>* grp = group(vrts[i]);
228 0 : if(!grp)
229 0 : continue;
230 0 : mg.mark(grp->m_master);
231 : mg.mark(grp->get_slaves().begin(), grp->get_slaves().end());
232 : }
233 : }
234 :
235 0 : if (is_master<TParent>(pParent)) {
236 : // create new group for e, with e as master
237 0 : Group<TElem>* newGroup = new Group<TElem>(e);
238 : UG_ASSERT(group(e) == NULL, "element already has a group.")
239 0 : set_group(newGroup, e);
240 :
241 : // iterate over slaves of parent
242 0 : ParentSlaveContainer* parentSlaves = slaves(pParent);
243 0 : for (typename ParentSlaveContainer::iterator i_slave =
244 0 : parentSlaves->begin(); i_slave != parentSlaves->end();
245 : ++i_slave) {
246 0 : TParent* parentSlave = *i_slave;
247 : UG_ASSERT(parentSlave, "parent slave not valid")
248 :
249 : // iterate over all children of parentSlave and make them slaves of e
250 : // only matching children of type TElem are considered.
251 : // note that not all children already have to exist at this point.
252 : // children which are created later on are handled by the 'else' section.
253 0 : for (size_t i_child = 0;
254 0 : i_child < mg.num_children<TElem>(parentSlave); ++i_child) {
255 : TElem* c = mg.get_child<TElem>(parentSlave, i_child);
256 : UG_ASSERT(!(e->base_object_id() == VERTEX)
257 : || (mg.num_children<TElem>(parentSlave) == 1),
258 : "At most one 1 vertex-child is currently allowed");
259 : // We use a special case for vertices here, since the position
260 : // attachment has not yet been set (this is currently done after
261 : // all observers have been informed about the creation of the new
262 : // vertex). Note that this is no problem for edges or faces, since
263 : // associated vertices have been properly initialized at that time.
264 0 : if((e->base_object_id() == VERTEX)) {
265 0 : make_slave(newGroup, c);
266 : } else {
267 : Grid::vertex_traits::secure_container vrts;
268 : mg.associated_elements(vrts, c);
269 : bool allMarked = true;
270 0 : for(size_t i = 0; i < vrts.size(); ++i){
271 0 : if(!mg.is_marked(vrts[i])){
272 : allMarked = false;
273 : break;
274 : }
275 : }
276 :
277 0 : if(allMarked)
278 0 : make_slave(newGroup, c);
279 : }
280 : }
281 : }
282 : } else {
283 : // find the associated master of e by checking for a matching child
284 : // in the set of children of the parents master.
285 : // If no appropriate child already exists at this point, no master can
286 : // be found. We simply leave the new element as it is, since it will
287 : // be found as a slave when the associated master will be created later
288 : // on (the 'if(is_master...)' case applies then.
289 :
290 : // attention: parentGroup may be null
291 : Group<TParent>* parentGroup = group<TParent>(pParent);
292 :
293 0 : if (parentGroup == NULL) {
294 0 : get_periodic_status_accessor<TElem>()[e] = P_SLAVE_MASTER_UNKNOWN;
295 0 : return;
296 : }
297 :
298 0 : TParent* parentMaster = parentGroup->m_master;
299 : bool master_found = false;
300 0 : for (size_t i_child = 0; i_child < mg.num_children<TElem>(parentMaster);
301 : ++i_child) {
302 : TElem* c = mg.get_child<TElem>(parentMaster, i_child);
303 : UG_ASSERT(!(e->base_object_id() == VERTEX)
304 : || (mg.num_children<TElem>(parentMaster) == 1),
305 : "At most one 1 vertex-child is currently allowed");
306 : // We use a special case for vertices here, since the position
307 : // attachment has not yet been set (this is currently done after
308 : // all observers have been informed about the creation of the new
309 : // vertex). Note that this is no problem for edges or faces, since
310 : // associated vertices have been properly initialized at that time.
311 0 : if((e->base_object_id() == VERTEX)) {
312 0 : make_slave(group(c), e);
313 : master_found = true;
314 : break;
315 : } else {
316 : Grid::vertex_traits::secure_container vrts;
317 : mg.associated_elements(vrts, c);
318 : bool allMarked = true;
319 0 : for(size_t i = 0; i < vrts.size(); ++i){
320 0 : if(!mg.is_marked(vrts[i])){
321 : allMarked = false;
322 : break;
323 : }
324 : }
325 :
326 0 : if(allMarked){
327 0 : make_slave(group(c), e);
328 : master_found = true;
329 : break;
330 : }
331 : }
332 : }
333 :
334 : // wait until a master for e will be known
335 : if (!master_found) {
336 0 : get_periodic_status_accessor<TElem>()[e] = P_SLAVE_MASTER_UNKNOWN;
337 : }
338 : }
339 :
340 0 : mg.end_marking();
341 : }
342 :
343 : template <class TElem>
344 0 : void PeriodicBoundaryManager::handle_creation_cast_wrapper(TElem* e,
345 : GridObject* pParent, bool replacesParent) {
346 : // we can only identify periodic elements, which have a periodic parent
347 0 : if(!pParent)
348 : return;
349 :
350 0 : if(replacesParent){
351 0 : replace_parent(e, static_cast<TElem*>(pParent));
352 : }
353 : else{
354 0 : switch (pParent->base_object_id()) {
355 0 : case VERTEX:
356 0 : handle_creation(e, static_cast<Vertex*>(pParent));
357 0 : break;
358 0 : case EDGE:
359 0 : handle_creation(e, static_cast<Edge*>(pParent));
360 0 : break;
361 0 : case FACE:
362 0 : handle_creation(e, static_cast<Face*>(pParent));
363 0 : break;
364 : // ignore volumes, as these are not meant to be periodic
365 : case VOLUME:
366 : break;
367 0 : default:
368 0 : UG_THROW("no handling for parent type: " << pParent->base_object_id())
369 : }
370 : }
371 : }
372 :
373 : /// handles deletion of element type
374 : template <class TElem>
375 0 : void PeriodicBoundaryManager::handle_deletion(TElem* e, TElem* replacedBy) {
376 0 : if((!is_periodic(e)) || replacedBy)
377 : return;
378 :
379 0 : if (is_master(e)) {
380 : // delete a master completely...
381 : if (replacedBy) {
382 : UG_ASSERT(!group(replacedBy),
383 : "replacing element is already in group")
384 : make_master(group(e), replacedBy);
385 : } else {
386 0 : remove_group(group(e));
387 : }
388 : } else { // slave
389 : UG_ASSERT(is_slave(e), "e should be a slave")
390 0 : if(!remove_slave(e))
391 0 : UG_THROW("old slave not removed.")
392 : }
393 : }
394 :
395 : template <class TElem>
396 0 : void PeriodicBoundaryManager::make_slave(Group<TElem>* g, TElem* slave) {
397 : UG_ASSERT(g, "invalid group")
398 : UG_ASSERT(slave, "invalid slave")
399 : UG_ASSERT(group(slave) != g, "trying to add a duplicate slave!")
400 : UG_ASSERT(m_pGrid->get_level(g->m_master) == m_pGrid->get_level(slave),
401 : "level of slave and group mismatch")
402 :
403 : // if slave is already engrouped, remove it first
404 0 : if(group(slave)) {
405 : UG_LOG("slave already engrouped. removing it from former group\n")
406 0 : if(!remove_slave(slave)) {
407 0 : UG_THROW("slave could not be removed")
408 : }
409 : }
410 :
411 : // add slave to group and set group attachment/status
412 0 : g->add_slave(slave);
413 0 : set_group(g, slave);
414 0 : }
415 :
416 : template <class TElem>
417 : void PeriodicBoundaryManager::make_master(Group<TElem>* g, TElem* new_master) {
418 : // group has a master, reset its group
419 : if(g->m_master) {
420 : set_group<TElem>(NULL, g->m_master);
421 : }
422 : g->m_master = new_master;
423 : }
424 :
425 : template <class TElem>
426 0 : bool PeriodicBoundaryManager::remove_slave(TElem* e) {
427 0 : if (slaves(e)) {
428 0 : typename Group<TElem>::SlaveContainer& s = *slaves(e);
429 0 : typename Group<TElem>::SlaveIterator pos = std::find(s.begin(), s.end(), e);
430 0 : if (pos != s.end()) {
431 : // s.erase(pos);
432 : // set_group<TElem>(NULL, e);
433 :
434 : // todo this is the only slave left, remove whole group?!
435 0 : if(s.size() == 1) {
436 : UG_LOGN("delete last slave, remove_group")
437 0 : remove_group(group(e));
438 : }
439 : // todo what if e was the last slave of of group?
440 : // if(s.empty()) {
441 : // UG_THROW("empty group leftover....")
442 : // }
443 0 : return true;
444 : }
445 : }
446 : return false;
447 : }
448 :
449 : template <class TElem>
450 0 : void PeriodicBoundaryManager::remove_group(Group<TElem>* g) {
451 : UG_ASSERT(g, "should remove invalid group")
452 :
453 : // reset group pointers of all group members to NULL
454 0 : set_group<TElem>(NULL, g->m_master);
455 : typename Group<TElem>::SlaveContainer& s = g->get_slaves();
456 : typename Group<TElem>::SlaveIterator iter;
457 0 : for (iter = s.begin(); iter != s.end(); ++iter) {
458 0 : set_group<TElem>(NULL, *iter);
459 : }
460 :
461 0 : delete g;
462 0 : }
463 :
464 : /**
465 : * merges g1 in g0 and deletes g1 afterwards
466 : */
467 : template <class TElem>
468 0 : void PeriodicBoundaryManager::merge_groups(Group<TElem>* g0, Group<TElem>* g1) {
469 : UG_ASSERT(g0 && g1, "groups not valid")
470 : UG_ASSERT(g0 != g1, "groups are equal")
471 :
472 : typedef typename Group<TElem>::SlaveContainer SlaveContainer;
473 : typedef typename Group<TElem>::SlaveIterator SlaveIterator;
474 :
475 : SlaveContainer& slaves_g1 = g1->get_slaves();
476 :
477 : // insert slaves of g1 at the end of slaves of g0
478 0 : for (SlaveIterator iter = slaves_g1.begin(); iter != slaves_g1.end();
479 : ++iter) {
480 0 : TElem* e = *iter;
481 : UG_ASSERT(e, "slave not valid")
482 : UG_ASSERT(e != g0->m_master, "slave of g1 is master of g0!")
483 : // we do not use make_slave() here, because g1 will be deleted at the end
484 0 : g0->add_slave(e);
485 0 : set_group(g0, e);
486 : }
487 :
488 : // make old master a slave of group g1
489 : // we do not use make_slave() here, because g1 will be deleted at the end
490 0 : g0->add_slave(g1->m_master);
491 0 : set_group(g0, g1->m_master);
492 :
493 : // remove old group
494 0 : delete g1;
495 0 : }
496 :
497 : template <class TElem>
498 : bool PeriodicBoundaryManager::is_slave(TElem* e) const {
499 0 : PeriodicStatus p = get_periodic_status_accessor<TElem>()[e];
500 : return (p == P_SLAVE || p == P_SLAVE_MASTER_UNKNOWN);
501 : }
502 :
503 : template <class TElem>
504 : bool PeriodicBoundaryManager::is_master(TElem* e) const {
505 0 : return get_periodic_status_accessor<TElem>()[e] == P_MASTER;
506 : }
507 :
508 : template <class TElem>
509 : PeriodicBoundaryManager::Group<TElem>* PeriodicBoundaryManager::group(
510 : TElem* e) const {
511 : UG_ASSERT(e, "element not valid.")
512 0 : return get_group_accessor<TElem>()[e];
513 : }
514 :
515 : template <class TElem>
516 0 : void PeriodicBoundaryManager::set_group(Group<TElem>* g, TElem* e) {
517 : UG_ASSERT(e, "element not valid for attachment access.")
518 : // set group pointer of element e
519 0 : get_group_accessor<TElem>()[e] = g;
520 :
521 : // set periodic status
522 0 : if (g == NULL) {
523 0 : get_periodic_status_accessor<TElem>()[e] = P_NOT_PERIODIC;
524 : } else {
525 0 : if (g->m_master == e)
526 0 : get_periodic_status_accessor<TElem>()[e] = P_MASTER;
527 : else
528 0 : get_periodic_status_accessor<TElem>()[e] = P_SLAVE;
529 : }
530 0 : }
531 :
532 : template <class TElem, class TIterator>
533 0 : void PeriodicBoundaryManager::check_elements_periodicity(
534 : TIterator begin,
535 : TIterator end,
536 : typename Group<TElem>::unique_pairs& s,
537 : ISubsetHandler* sh) {
538 :
539 : typedef typename Group<TElem>::SlaveContainer Container;
540 : typedef typename Group<TElem>::SlaveIterator SlaveIter;
541 : typedef typename ElementStorage<TElem>::SectionContainer::const_iterator SecContainerIter;
542 :
543 0 : for (SecContainerIter iter = begin; iter != end; ++iter) {
544 : TElem* e = *iter;
545 0 : if(! e)
546 : continue;
547 0 : if(!is_periodic(e)) {
548 : // lookup subset name of element
549 : const char* sh_name = "";
550 0 : if(sh) {
551 : int element_si = sh->get_subset_index(e);
552 0 : sh_name = sh->get_subset_name(element_si);
553 : }
554 0 : UG_THROW("Element in subset '" << sh_name
555 : << "' is not periodic after identification: "
556 : << GetGridObjectCenter(*get_grid(), e)
557 : << "\nCheck your geometry for symmetry!\n")
558 : }
559 :
560 0 : if (master(e) == e) {
561 0 : Container* _slaves = slaves(e);
562 0 : if(! _slaves)
563 0 : UG_THROW("masters slave storage is not valid.")
564 :
565 0 : if(_slaves->empty())
566 0 : UG_THROW("master has no slaves")
567 :
568 0 : for (SlaveIter i = _slaves->begin(); i != _slaves->end(); ++i) {
569 0 : TElem* slave = *i;
570 0 : typename Group<TElem>::master_slave_pair p = std::make_pair(e, slave);
571 : bool inserted = (s.insert(p)).second;
572 0 : if(! inserted)
573 0 : UG_THROW("master/slave pair already exists.");
574 : }
575 : }
576 : }
577 0 : }
578 :
579 : template <class TDomain>
580 0 : void IdentifySubsets(TDomain& dom, const char* sName1, const char* sName2) {
581 : // get subset handler from domain
582 : typedef typename TDomain::subset_handler_type subset_handler_type;
583 :
584 0 : subset_handler_type& sh = *dom.subset_handler();
585 :
586 0 : int si1 = sh.get_subset_index(sName1);
587 0 : int si2 = sh.get_subset_index(sName2);
588 :
589 : #ifdef UG_PARALLEL
590 : //if(pcl::NumProcs() > 1)
591 : // UG_THROW("sorry, in real parallel environment periodic bnds are not impled yet.");
592 : #endif
593 :
594 0 : if (si1 == -1)
595 0 : UG_THROW("IdentifySubsets: given subset name " << sName1 << " does not exist");
596 0 : if (si2 == -1)
597 0 : UG_THROW("IdentifySubsets: given subset name " << sName2 << " does not exist");
598 :
599 0 : IdentifySubsets(dom, si1, si2);
600 0 : }
601 :
602 : /// performs geometric ident of periodic elements and master slave
603 : template <class TDomain>
604 0 : void IdentifySubsets(TDomain& dom, int sInd1, int sInd2) {
605 :
606 : #ifdef UG_PARALLEL
607 : //if(pcl::NumProcs() > 1)
608 : // UG_THROW("sorry, in real parallel environment periodic bnds are not impled yet.");
609 : #endif
610 :
611 0 : if (sInd1 == -1 || sInd2 == -1) {
612 0 : UG_THROW("IdentifySubsets: at least one invalid subset given!")
613 : }
614 :
615 0 : if(sInd1 == sInd2) {
616 0 : UG_THROW("IdentifySubsets: can not identify two identical subsets!")
617 : }
618 :
619 : // ensure grid has support for periodic boundaries
620 0 : if (!dom.grid()->has_periodic_boundaries()) {
621 0 : dom.grid()->set_periodic_boundaries(true);
622 : }
623 :
624 0 : PeriodicBoundaryManager& pbm = *dom.grid()->periodic_boundary_manager();
625 :
626 : typedef typename TDomain::position_type position_type;
627 : typedef typename TDomain::position_accessor_type position_accessor_type;
628 :
629 : // get subset handler from domain
630 : typedef typename TDomain::subset_handler_type subset_handler_type;
631 :
632 0 : subset_handler_type& sh = *dom.subset_handler();
633 :
634 : // get aaPos from domain
635 : position_accessor_type& aaPos = dom.position_accessor();
636 :
637 : // create parallel shift identifier to match subset elements
638 : ParallelShiftIdentifier<position_accessor_type> ident(aaPos);
639 :
640 : // shift vector between subsets
641 : position_type shift;
642 :
643 : // collect all geometric objects (even for all levels in case of multi grid)
644 0 : GridObjectCollection goc1 = sh.get_grid_objects_in_subset(sInd1);
645 0 : GridObjectCollection goc2 = sh.get_grid_objects_in_subset(sInd2);
646 :
647 0 : if(goc1.num<Vertex>() != goc2.num<Vertex>()) {
648 0 : UG_THROW("IdentifySubsets: Given subsets have different number of vertices."
649 : "\nnum# in " << sh.get_subset_name(sInd1) << ": " << goc1.num<Vertex>() <<
650 : "\nnum# in " << sh.get_subset_name(sInd2) << ": " << goc2.num<Vertex>())
651 : }
652 :
653 0 : if(goc1.num<Edge>() != goc2.num<Edge>()) {
654 0 : UG_THROW("IdentifySubsets: Given subsets have different number of edges."
655 : "\nnum# in " << sh.get_subset_name(sInd1) << ": " << goc1.num<Edge>() <<
656 : "\nnum# in " << sh.get_subset_name(sInd2) << ": " << goc2.num<Edge>())
657 : }
658 :
659 0 : if(goc1.num<Face>() != goc2.num<Face>()) {
660 0 : UG_THROW("IdentifySubsets: Given subsets have different number of faces."
661 : "\nnum# in " << sh.get_subset_name(sInd1) << ": " << goc1.num<Face>() <<
662 : "\nnum# in " << sh.get_subset_name(sInd2) << ": " << goc2.num<Face>())
663 : }
664 :
665 : // map start type of recursion dependent to TDomain
666 : // in 3d start with faces, in 2d with edges, in 1d with vertices
667 : // namespace mpl = boost::mpl;
668 : // typedef mpl::map<mpl::pair<Domain1d, Vertex>,
669 : // mpl::pair<Domain2d, Edge>,
670 : // mpl::pair<Domain3d, Face> > m;
671 :
672 : // typedef typename mpl::at<m, TDomain>::type TElem;
673 : typedef typename grid_dim_traits<TDomain::dim>::side_type TElem;
674 : typedef typename ElementStorage<TElem>::SectionContainer::iterator gocIter;
675 :
676 : // calculate shift vector for top level
677 0 : position_type c1 = CalculateCenter(goc1.begin<TElem>(0), goc1.end<TElem>(0),
678 : aaPos);
679 0 : position_type c2 = CalculateCenter(goc2.begin<TElem>(0), goc2.end<TElem>(0),
680 : aaPos);
681 :
682 : VecSubtract(shift, c1, c2);
683 : ident.set_shift(shift);
684 :
685 : // for each level of multi grid. In case of simple grid only one iteration
686 0 : for (size_t lvl = 0; lvl < goc1.num_levels(); lvl++) {
687 : // identify corresponding elements for second subset. A element is considered
688 : // to have symmetric element in second subset if there exists a shift vector between them.
689 : // todo use kd-tree for fast lookup of shifted elements
690 : for (gocIter iter1 = goc1.begin<TElem>(lvl);
691 0 : iter1 != goc1.end<TElem>(lvl); ++iter1)
692 : for (gocIter iter2 = goc2.begin<TElem>(lvl);
693 0 : iter2 != goc2.end<TElem>(lvl); ++iter2) {
694 0 : if(ident.match(*iter1, *iter2)) {
695 0 : pbm.identify(*iter1, *iter2, ident);
696 : }
697 : }
698 : }
699 :
700 : // ensure periodic identification has been performed correctly
701 0 : pbm.check_periodicity(goc1, goc2, &sh);
702 0 : }
703 :
704 : } // end namespace ug
705 :
706 : #endif /* PERIODIC_IDENTIFIER_IMPL_HPP_ */
|
