Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel
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 "dof_distribution.h"
34 : #include "lib_disc/function_spaces/grid_function.h"
35 :
36 : #include "common/log.h"
37 : #include "lib_disc/domain.h"
38 : #include "lib_disc/local_finite_element/local_dof_set.h"
39 : #include "lib_disc/reference_element/reference_element.h"
40 : #include "lib_disc/reference_element/reference_element_traits.h"
41 : #include "lib_disc/reference_element/reference_mapping.h"
42 : #include "lib_disc/reference_element/reference_mapping_provider.h"
43 : #include "lib_disc/common/groups_util.h"
44 : #include "common/util/string_util.h"
45 : #include "orientation.h"
46 : #include "lib_grid/tools/periodic_boundary_manager.h"
47 :
48 : #include "lib_grid/file_io/file_io.h"
49 : #include "lib_grid/algorithms/debug_util.h"
50 :
51 : using namespace std;
52 :
53 : namespace ug{
54 :
55 : ////////////////////////////////////////////////////////////////////////////////
56 : // DoFDistribution
57 : ////////////////////////////////////////////////////////////////////////////////
58 :
59 0 : DoFDistribution::
60 : DoFDistribution(SmartPtr<MultiGrid> spMG,
61 : SmartPtr<MGSubsetHandler> spMGSH,
62 : ConstSmartPtr<DoFDistributionInfo> spDDInfo,
63 : SmartPtr<SurfaceView> spSurfView,
64 : const GridLevel& level, bool bGrouped,
65 0 : SmartPtr<DoFIndexStorage> spDoFIndexStorage)
66 : : DoFDistributionInfoProvider(spDDInfo),
67 0 : m_bGrouped(bGrouped),
68 : m_spMG(spMG),
69 0 : m_pMG(m_spMG.get()),
70 : m_spMGSH(spMGSH),
71 : m_spSurfView(spSurfView),
72 0 : m_gridLevel(level),
73 : m_spDoFIndexStorage(spDoFIndexStorage),
74 0 : m_numIndex(0)
75 : {
76 0 : if(m_spDoFIndexStorage.invalid())
77 0 : m_spDoFIndexStorage = SmartPtr<DoFIndexStorage>(new DoFIndexStorage(spMG, spDDInfo));
78 :
79 0 : check_subsets();
80 0 : m_vNumIndexOnSubset.resize(num_subsets(), 0);
81 :
82 : #ifdef UG_PARALLEL
83 : m_spAlgebraLayouts = SmartPtr<AlgebraLayouts>(new AlgebraLayouts);
84 : #endif
85 :
86 0 : reinit();
87 0 : }
88 :
89 :
90 0 : DoFDistribution::
91 0 : ~DoFDistribution() {}
92 :
93 :
94 0 : void DoFDistribution::check_subsets()
95 : {
96 : // check, that all geom objects are assigned to a subset
97 0 : if( m_spMGSH->num<Vertex>() != multi_grid()->num<Vertex>())
98 : {
99 0 : geometry_traits<Vertex>::const_iterator iter = multi_grid()->begin<Vertex>();
100 0 : geometry_traits<Vertex>::const_iterator iterEnd = multi_grid()->end<Vertex>();
101 0 : for (; iter != iterEnd; ++iter)
102 0 : if (m_spMGSH->get_subset_index(*iter) == -1)
103 0 : UG_LOGN("Missing subset info for " << ElementDebugInfo(*multi_grid(), *iter));
104 0 : UG_THROW("All Vertices "
105 : " must be assigned to a subset. The passed subset handler "
106 : " contains non-assigned elements, thus the dof distribution"
107 : " is not possible, aborting.");
108 : }
109 0 : if( m_spMGSH->num<Edge>() != multi_grid()->num<Edge>())
110 : {
111 0 : geometry_traits<Edge>::const_iterator iter = multi_grid()->begin<Edge>();
112 0 : geometry_traits<Edge>::const_iterator iterEnd = multi_grid()->end<Edge>();
113 0 : for (; iter != iterEnd; ++iter)
114 0 : if (m_spMGSH->get_subset_index(*iter) == -1)
115 0 : UG_LOGN("Missing subset info for " << ElementDebugInfo(*multi_grid(), *iter));
116 :
117 0 : UG_THROW("All Edges "
118 : " must be assigned to a subset. The passed subset handler "
119 : " contains non-assigned elements, thus the dof distribution"
120 : " is not possible, aborting.");
121 : }
122 0 : if( m_spMGSH->num<Face>() != multi_grid()->num<Face>())
123 0 : UG_THROW("All Faces "
124 : " must be assigned to a subset. The passed subset handler "
125 : " contains non-assigned elements, thus the dof distribution"
126 : " is not possible, aborting.");
127 :
128 0 : if( m_spMGSH->num<Volume>() != multi_grid()->num<Volume>())
129 0 : UG_THROW("All Volumes "
130 : " must be assigned to a subset. The passed subset handler "
131 : " contains non-assigned elements, thus the dof distribution"
132 : " is not possible, aborting.");
133 0 : }
134 :
135 0 : SurfaceView::SurfaceConstants DoFDistribution::defaultValidSurfState() const{
136 0 : if(m_gridLevel.is_level()) return SurfaceView::ALL;
137 0 : else if(m_gridLevel.is_surface()) return SurfaceView::ALL_BUT_SHADOW_COPY;
138 0 : else UG_THROW("DoFDistribution: GridLevel::type not valid.")
139 : }
140 :
141 : ////////////////////////////////////////////////////////////////////////////////
142 : // DoFDistribution: Index Access
143 : ////////////////////////////////////////////////////////////////////////////////
144 :
145 : template <typename TBaseObject>
146 0 : void DoFDistribution::
147 : add(TBaseObject* obj, const ReferenceObjectID roid, const int si)
148 : {
149 : UG_ASSERT(si >= 0, "Invalid subset index passed");
150 :
151 : // if no dofs on this subset for the roid, do nothing
152 0 : if(num_dofs(roid,si) == 0) return;
153 :
154 : // check for periodicity
155 : bool master = false;
156 0 : if(m_spMG->has_periodic_boundaries()){
157 0 : PeriodicBoundaryManager& pbm = *m_spMG->periodic_boundary_manager();
158 0 : if(pbm.is_slave(obj)) return; // ignore slaves
159 0 : if(pbm.is_master(obj)) master = true;
160 : }
161 :
162 : // compute the number of indices needed on the Geometric object
163 : size_t numNewIndex = 1;
164 0 : if(!m_bGrouped) numNewIndex = num_dofs(roid,si);
165 :
166 : // set first available index to the object. The first available index is the
167 : // first managed index plus the size of the index set. (If holes are in the
168 : // index set, this is not treated here, holes remain)
169 0 : obj_index(obj) = m_numIndex;
170 :
171 : // number of managed indices and the number of managed indices on the subset has
172 : // changed. Thus, increase the counters.
173 0 : m_numIndex += numNewIndex;
174 0 : m_vNumIndexOnSubset[si] += numNewIndex;
175 :
176 : // if obj is a master, assign all its slaves
177 0 : if(master) {
178 : typedef typename PeriodicBoundaryManager::Group<TBaseObject>::SlaveContainer SlaveContainer;
179 : typedef typename PeriodicBoundaryManager::Group<TBaseObject>::SlaveIterator SlaveIterator;
180 0 : SlaveContainer& slaves = *m_spMG->periodic_boundary_manager()->slaves(obj);
181 0 : size_t master_index = obj_index(obj);
182 0 : for(SlaveIterator iter = slaves.begin(); iter != slaves.end(); ++iter){
183 0 : obj_index(*iter) = master_index;
184 : }
185 : }
186 : }
187 :
188 : template <typename TBaseElem>
189 0 : size_t DoFDistribution::
190 : extract_inner_algebra_indices(TBaseElem* elem,
191 : std::vector<size_t>& ind) const
192 : {
193 : // get roid type and subset index
194 0 : const int si = m_spMGSH->get_subset_index(elem);
195 0 : const ReferenceObjectID roid = elem->reference_object_id();
196 :
197 : // check if dofs present
198 0 : if(num_dofs(roid,si) > 0)
199 : {
200 : // get index
201 0 : const size_t firstIndex = obj_index(elem);
202 :
203 0 : if(!m_bGrouped)
204 : {
205 0 : for(size_t fct = 0; fct < num_fct(); ++fct)
206 : {
207 : // check that function is def on subset
208 0 : if(!is_def_in_subset(fct, si)) continue;
209 :
210 : // get number of DoFs in this sub-geometric object
211 : const size_t numDoFsOnSub = num_fct_dofs(fct,roid,si);
212 :
213 : // compute index
214 0 : const size_t index = firstIndex + offset(roid,si,fct);
215 :
216 : // add dof to local indices
217 0 : for(size_t j = 0; j < numDoFsOnSub; ++j)
218 0 : ind.push_back(index+j);
219 : }
220 : }
221 : else
222 : {
223 : // add dof to local indices
224 0 : ind.push_back(firstIndex);
225 : }
226 : }
227 :
228 : // return number of indices
229 0 : return ind.size();
230 : }
231 :
232 : template<typename TBaseElem>
233 0 : void DoFDistribution::
234 : extract_inner_algebra_indices(const typename Grid::traits<TBaseElem>::secure_container& vElem,
235 : std::vector<size_t>& ind) const
236 : {
237 : // loop passed elements
238 0 : for(size_t i = 0; i < vElem.size(); ++i)
239 0 : inner_algebra_indices(vElem[i], ind, false);
240 0 : }
241 :
242 : template<typename TBaseElem>
243 : size_t DoFDistribution::_inner_algebra_indices(TBaseElem* elem,
244 : std::vector<size_t>& ind,
245 : bool bClear) const
246 : {
247 : // clear indices
248 0 : if(bClear) ind.clear();
249 :
250 : // return
251 0 : return extract_inner_algebra_indices(elem, ind);
252 : }
253 :
254 : // this function could be merged with inner_algebra_indices with additional
255 : // default paramter e.g. selectedFct==-1 -> all functions.
256 0 : size_t DoFDistribution::inner_algebra_indices_for_fct(GridObject* elem,
257 : std::vector<size_t>& ind,
258 : bool bClear, int fct) const
259 : {
260 : // clear indices
261 0 : if(bClear) ind.clear();
262 :
263 : // get roid type and subset index
264 0 : const int si = m_spMGSH->get_subset_index(elem);
265 0 : const ReferenceObjectID roid = elem->reference_object_id();
266 :
267 : // check if dofs present
268 0 : if(num_dofs(roid,si) > 0)
269 : {
270 : // get index
271 0 : const size_t firstIndex = obj_index(elem);
272 :
273 0 : if(!m_bGrouped)
274 : {
275 : // check that function is def on subset
276 0 : if(!is_def_in_subset(fct, si)) return ind.size();
277 :
278 : // get number of DoFs in this sub-geometric object
279 : const size_t numDoFsOnSub = num_fct_dofs(fct,roid,si);
280 :
281 : // compute index
282 0 : const size_t index = firstIndex + offset(roid,si,fct);
283 :
284 : // add dof to local indices
285 0 : for(size_t j = 0; j < numDoFsOnSub; ++j)
286 0 : ind.push_back(index+j);
287 : }
288 : else
289 : {
290 : // add dof to local indices
291 0 : ind.push_back(firstIndex);
292 : }
293 : }
294 :
295 : // return number of indices
296 0 : return ind.size();
297 : }
298 :
299 : template<typename TBaseElem>
300 0 : size_t DoFDistribution::_algebra_indices(TBaseElem* elem,
301 : std::vector<size_t>& ind,
302 : bool bClear) const
303 : {
304 : // clear indices
305 0 : if(bClear) ind.clear();
306 :
307 : // reference dimension
308 : static const int dim = TBaseElem::dim;
309 :
310 : // get all sub-elements and add indices on those
311 0 : if(max_dofs(VERTEX) > 0)
312 : {
313 : Grid::SecureVertexContainer vVrt;
314 0 : m_pMG->associated_elements(vVrt, elem);
315 0 : extract_inner_algebra_indices<Vertex>(vVrt, ind);
316 : }
317 0 : if(dim >= EDGE && max_dofs(EDGE) > 0)
318 : {
319 : Grid::SecureEdgeContainer vEdge;
320 0 : m_pMG->associated_elements(vEdge, elem);
321 0 : extract_inner_algebra_indices<Edge>(vEdge, ind);
322 : }
323 0 : if(dim >= FACE && max_dofs(FACE) > 0)
324 : {
325 : Grid::SecureFaceContainer vFace;
326 0 : m_pMG->associated_elements(vFace, elem);
327 0 : extract_inner_algebra_indices<Face>(vFace, ind);
328 : }
329 0 : if(dim >= VOLUME && max_dofs(VOLUME) > 0)
330 : {
331 : Grid::SecureVolumeContainer vVol;
332 0 : m_pMG->associated_elements(vVol, elem);
333 0 : extract_inner_algebra_indices<Volume>(vVol, ind);
334 : }
335 :
336 : // return number of indices
337 0 : return ind.size();
338 : }
339 :
340 : template<typename TBaseElem, typename TSubBaseElem>
341 0 : void DoFDistribution::
342 : dof_indices(TBaseElem* elem, const ReferenceObjectID roid,
343 : size_t fct, std::vector<DoFIndex>& ind,
344 : const typename Grid::traits<TSubBaseElem>::secure_container& vElem) const
345 : {
346 : // storage for offsets
347 : std::vector<size_t> vOrientOffset;
348 :
349 : // loop passed elements
350 0 : for(size_t i = 0; i < vElem.size(); ++i)
351 : {
352 : // get subelement
353 : TSubBaseElem* subElem = vElem[i];
354 :
355 : // get subset index
356 0 : const int si = m_spMGSH->get_subset_index(subElem);
357 :
358 : // check if function is defined on the subset
359 0 : if(!is_def_in_subset(fct, si)) continue;
360 :
361 : // get reference object id for subselement
362 0 : const ReferenceObjectID subRoid = subElem->reference_object_id();
363 :
364 : // check if dof given
365 0 : if(num_dofs(subRoid,si) == 0) continue;
366 :
367 : // get number of DoFs in this sub-geometric object
368 : const size_t numDoFsOnSub = num_fct_dofs(fct, subRoid, si);
369 :
370 : // get the orientation for this subelement
371 0 : ComputeOrientationOffset(vOrientOffset, elem, subElem, i, lfeid(fct));
372 :
373 : UG_ASSERT(vOrientOffset.size() == numDoFsOnSub ||
374 : vOrientOffset.empty(), "Something wrong with orientation");
375 :
376 0 : if(!m_bGrouped)
377 : {
378 : // compute index
379 0 : const size_t index = obj_index(subElem) + offset(subRoid,si,fct);
380 :
381 0 : if(vOrientOffset.empty()){
382 0 : for(size_t j = 0; j < numDoFsOnSub; ++j)
383 0 : ind.push_back(DoFIndex(index + j,0));
384 : }
385 : else{
386 0 : for(size_t j = 0; j < numDoFsOnSub; ++j)
387 0 : ind.push_back(DoFIndex(index + vOrientOffset[j],0));
388 : }
389 : }
390 : else
391 : {
392 : // compute index
393 : const size_t comp = offset(subRoid,si,fct);
394 0 : const size_t firstIndex = obj_index(subElem);
395 :
396 0 : if(vOrientOffset.empty()){
397 0 : for(size_t j = 0; j < numDoFsOnSub; ++j)
398 0 : ind.push_back(DoFIndex(firstIndex, comp + j));
399 : }
400 : else{
401 0 : for(size_t j = 0; j < numDoFsOnSub; ++j)
402 0 : ind.push_back(DoFIndex(firstIndex, comp + vOrientOffset[j]));
403 : }
404 : }
405 :
406 : } // end loop sub elements
407 0 : }
408 :
409 : template<typename TBaseElem>
410 0 : size_t DoFDistribution::_inner_dof_indices(TBaseElem* elem, size_t fct,
411 : std::vector<DoFIndex>& ind,
412 : bool bClear) const
413 : {
414 : // clear if requested
415 0 : if(bClear) ind.clear();
416 :
417 : // get subset index
418 0 : const int si = m_spMGSH->get_subset_index(elem);
419 :
420 : // check if function is defined on the subset
421 0 : if(!is_def_in_subset(fct, si)) return ind.size();
422 :
423 : // get roid type
424 0 : const ReferenceObjectID roid = elem->reference_object_id();
425 :
426 : // get number of DoFs in this sub-geometric object
427 : const size_t numDoFsOnSub = num_fct_dofs(fct,roid,si);
428 :
429 : // check if dof given
430 0 : if(numDoFsOnSub == 0) return ind.size();
431 :
432 : // Note: No orientation needed
433 0 : if(!m_bGrouped)
434 : {
435 : // compute index
436 0 : const size_t index = obj_index(elem) + offset(roid,si,fct);
437 :
438 0 : for(size_t j = 0; j < numDoFsOnSub; ++j)
439 0 : ind.push_back(DoFIndex(index+j,0));
440 : }
441 : else
442 : {
443 : // compute index
444 : const size_t comp = offset(roid,si,fct);
445 0 : const size_t firstIndex = obj_index(elem);
446 :
447 0 : for(size_t j = 0; j < numDoFsOnSub; ++j)
448 0 : ind.push_back(DoFIndex(firstIndex, comp+j));
449 : }
450 :
451 : // done
452 0 : return ind.size();
453 : }
454 :
455 : template <typename TConstraining, typename TConstrained, typename TBaseElem>
456 0 : void DoFDistribution::
457 : constrained_vertex_dof_indices(size_t fct,std::vector<DoFIndex>& ind,
458 : const typename Grid::traits<TBaseElem>::secure_container& vSubElem) const
459 : {
460 : // loop all edges
461 0 : for(size_t i = 0; i < vSubElem.size(); ++i)
462 : {
463 : // only constraining objects are of interest
464 0 : TConstraining* constrainingObj = dynamic_cast<TConstraining*>(vSubElem[i]);
465 0 : if(constrainingObj == NULL) continue;
466 :
467 : // get subset index
468 0 : const int si = m_spMGSH->get_subset_index(vSubElem[i]);
469 :
470 : // loop constraining vertices
471 0 : for(size_t j = 0; j != constrainingObj->num_constrained_vertices(); ++j)
472 : {
473 : // get vertex
474 : TConstrained* vrt = constrainingObj->constrained_vertex(j);
475 :
476 : // get roid
477 0 : const ReferenceObjectID subRoid = vrt->reference_object_id();
478 :
479 : // check if dof given
480 0 : if(num_dofs(subRoid,si) == 0) continue;
481 :
482 : // get subset index
483 0 : int si = m_spMGSH->get_subset_index(vrt);
484 :
485 : // check that function is defined on subset
486 0 : if(!is_def_in_subset(fct, si)) continue;
487 :
488 : // get number of DoFs in this sub-geometric object
489 : const size_t numDoFsOnSub = num_fct_dofs(fct, subRoid, si);
490 :
491 0 : if(!m_bGrouped)
492 : {
493 : // compute index
494 0 : const size_t index = obj_index(vrt) + offset(subRoid,si,fct);
495 :
496 0 : for (size_t k=0;k<numDoFsOnSub;k++)
497 0 : ind.push_back(DoFIndex(index + k,0));
498 : }
499 : else
500 : {
501 : // compute index
502 0 : const size_t index = obj_index(vrt);
503 : const size_t comp = offset(subRoid,si,fct);
504 :
505 : // add dof to local indices
506 0 : for(size_t k = 0; k < numDoFsOnSub; ++k)
507 0 : ind.push_back(DoFIndex(index, comp + k));
508 : }
509 : }
510 : }
511 0 : }
512 :
513 : template <typename TBaseElem,typename TConstraining, typename TConstrained, typename TSubElem>
514 0 : void DoFDistribution::
515 : constrained_edge_dof_indices(TBaseElem* elem,size_t fct,std::vector<DoFIndex>& ind,
516 : const typename Grid::traits<TSubElem>::secure_container& vSubElem) const
517 : {
518 : // storage for offsets
519 : std::vector<size_t> vOrientOffset;
520 :
521 : // loop all edges
522 0 : for(size_t i = 0; i < vSubElem.size(); ++i)
523 : {
524 : // only constraining objects are of interest
525 0 : TConstraining* constrainingObj = dynamic_cast<TConstraining*>(vSubElem[i]);
526 0 : if(constrainingObj == NULL) continue;
527 :
528 : std::vector<size_t> sortedInd;
529 0 : sort_constrained_edges<TBaseElem,TConstraining,TConstrained>(sortedInd,elem,constrainingObj,i);
530 :
531 : // get the orientation for this subelement
532 0 : ComputeOrientationOffset(vOrientOffset, elem, constrainingObj, i, lfeid(fct));
533 :
534 : // loop constraining edges
535 0 : for(size_t j = 0; j != constrainingObj->num_constrained_edges(); ++j)
536 : {
537 : // get edge
538 0 : TConstrained* edg = constrainingObj->constrained_edge(sortedInd[j]);
539 :
540 : // get roid
541 0 : const ReferenceObjectID subRoid = edg->reference_object_id();
542 :
543 : // get subset index
544 0 : int si = m_spMGSH->get_subset_index(edg);
545 :
546 :
547 : // check that function is defined on subset
548 0 : if(!is_def_in_subset(fct, si)) continue;
549 :
550 : // check if dof given
551 0 : if(num_dofs(subRoid,si) == 0) continue;
552 :
553 : const size_t numDoFsOnSub=num_fct_dofs(fct, subRoid, si);
554 :
555 0 : if(!m_bGrouped)
556 : {
557 : // compute index
558 0 : const size_t index = obj_index(edg) + offset(subRoid,si,fct);
559 :
560 0 : if(vOrientOffset.empty()){
561 0 : for(size_t k = 0; k < numDoFsOnSub; ++k)
562 0 : ind.push_back(DoFIndex(index + k,0));
563 : }
564 : else{
565 0 : for(size_t k = 0; k < numDoFsOnSub; ++k)
566 0 : ind.push_back(DoFIndex(index + vOrientOffset[k],0));
567 : }
568 : }
569 : else
570 : {
571 : // compute index
572 0 : const size_t index = obj_index(edg);
573 : const size_t comp = offset(subRoid,si,fct);
574 :
575 0 : if(vOrientOffset.empty()){
576 0 : for(size_t k = 0; k < numDoFsOnSub; ++k)
577 0 : ind.push_back(DoFIndex(index, comp + k));
578 : }
579 : else{
580 0 : for(size_t k = 0; k < numDoFsOnSub; ++k)
581 0 : ind.push_back(DoFIndex(index, comp + vOrientOffset[k]));
582 : }
583 : }
584 : }
585 : }
586 0 : }
587 :
588 : template <typename TBaseElem,typename TConstraining, typename TConstrained, typename TSubElem>
589 0 : void DoFDistribution::
590 : constrained_face_dof_indices(TBaseElem* elem,size_t fct,std::vector<DoFIndex>& ind,
591 : const typename Grid::traits<TSubElem>::secure_container& vSubElem) const
592 : {
593 : // storage for offsets
594 : std::vector<size_t> vOrientOffset;
595 :
596 : // loop all edges
597 0 : for(size_t i = 0; i < vSubElem.size(); ++i)
598 : {
599 : // only constraining objects are of interest
600 0 : TConstraining* constrainingObj = dynamic_cast<TConstraining*>(vSubElem[i]);
601 0 : if(constrainingObj == NULL) continue;
602 :
603 : std::vector<size_t> sortedInd;
604 0 : sort_constrained_faces<TBaseElem,TConstraining,TConstrained>(sortedInd,elem,constrainingObj,i);
605 :
606 : // get the orientation for this subelement
607 0 : ComputeOrientationOffset(vOrientOffset, elem, constrainingObj, i, lfeid(fct));
608 :
609 : // loop constraining edges
610 0 : for(size_t j = 0; j != constrainingObj->num_constrained_faces(); ++j)
611 : {
612 : // get face
613 0 : TConstrained* face = constrainingObj->constrained_face(sortedInd[j]);
614 :
615 : // get roid
616 0 : const ReferenceObjectID subRoid = face->reference_object_id();
617 :
618 : // get subset index
619 0 : int si = m_spMGSH->get_subset_index(face);
620 :
621 :
622 : // check that function is defined on subset
623 0 : if(!is_def_in_subset(fct, si)) continue;
624 :
625 : // check if dof given
626 0 : if(num_dofs(subRoid,si) == 0) continue;
627 :
628 : const size_t numDoFsOnSub=num_fct_dofs(fct, subRoid, si);
629 :
630 0 : if(!m_bGrouped)
631 : {
632 : // compute index
633 0 : const size_t index = obj_index(face) + offset(subRoid,si,fct);
634 :
635 0 : if(vOrientOffset.empty()){
636 0 : for(size_t k = 0; k < numDoFsOnSub; ++k)
637 0 : ind.push_back(DoFIndex(index + k,0));
638 : }
639 : else{
640 0 : for(size_t k = 0; k < numDoFsOnSub; ++k)
641 0 : ind.push_back(DoFIndex(index + vOrientOffset[k],0));
642 : }
643 : }
644 : else
645 : {
646 : // compute index
647 0 : const size_t index = obj_index(face);
648 : const size_t comp = offset(subRoid,si,fct);
649 :
650 0 : if(vOrientOffset.empty()){
651 0 : for(size_t k = 0; k < numDoFsOnSub; ++k)
652 0 : ind.push_back(DoFIndex(index, comp + k));
653 : }
654 : else{
655 0 : for(size_t k = 0; k < numDoFsOnSub; ++k)
656 0 : ind.push_back(DoFIndex(index, comp + vOrientOffset[k]));
657 : }
658 : }
659 : }
660 : }
661 0 : }
662 :
663 :
664 : template<typename TBaseElem>
665 0 : size_t DoFDistribution::_dof_indices(TBaseElem* elem, size_t fct,
666 : std::vector<DoFIndex>& ind,
667 : bool bHang, bool bClear) const
668 : {
669 : // clear indices
670 0 : if(bClear) ind.clear();
671 :
672 : // reference dimension
673 : static const int dim = TBaseElem::dim;
674 :
675 : // reference object id
676 0 : const ReferenceObjectID roid = elem->reference_object_id();
677 :
678 : // storage for (maybe needed) subelements
679 : Grid::SecureVertexContainer vCorner;
680 : Grid::SecureEdgeContainer vEdge;
681 : Grid::SecureFaceContainer vFace;
682 : Grid::SecureVolumeContainer vVol;
683 :
684 : // collect elements, if needed
685 : if(dim >= VERTEX)
686 0 : if(max_dofs(VERTEX) > 0) m_pMG->associated_elements_sorted(vCorner, elem);
687 : if(dim >= EDGE)
688 0 : if(max_dofs(EDGE) > 0 || bHang) m_pMG->associated_elements_sorted(vEdge, elem);
689 : if(dim >= FACE)
690 0 : if(max_dofs(FACE) > 0 || bHang) m_pMG->associated_elements_sorted(vFace, elem);
691 : if(dim >= VOLUME)
692 0 : if(max_dofs(VOLUME) > 0) m_pMG->associated_elements_sorted(vVol, elem);
693 :
694 : // get regular dofs on all subelements and the element itself
695 : // use specialized function for vertices (since only one position and one reference object)
696 0 : if(dim >= VERTEX && max_dofs(VERTEX) > 0) dof_indices<TBaseElem, Vertex>(elem, roid, fct, ind, vCorner);
697 0 : if(dim >= EDGE && max_dofs(EDGE) > 0) dof_indices<TBaseElem, Edge>(elem, roid, fct, ind, vEdge);
698 0 : if(dim >= FACE && max_dofs(FACE) > 0) dof_indices<TBaseElem, Face>(elem, roid, fct, ind, vFace);
699 0 : if(dim >= VOLUME && max_dofs(VOLUME) > 0) dof_indices<TBaseElem, Volume>(elem, roid, fct, ind, vVol);
700 :
701 : // If no hanging dofs are required, we're done
702 0 : if(!bHang) return ind.size();
703 :
704 : // get dofs on hanging vertices
705 0 : if(max_dofs(VERTEX > 0))
706 : {
707 0 : if(dim >= EDGE) constrained_vertex_dof_indices<ConstrainingEdge, Vertex, Edge>(fct,ind,vEdge);
708 0 : if(dim >= FACE) constrained_vertex_dof_indices<ConstrainingQuadrilateral, Vertex, Face>(fct,ind,vFace);
709 : }
710 :
711 : // get dofs on hanging edges
712 0 : if (max_dofs(EDGE) > 0){
713 0 : if(dim >= EDGE) constrained_edge_dof_indices<TBaseElem,ConstrainingEdge, Edge, Edge>(elem,fct,ind, vEdge);
714 0 : if(dim >= FACE) constrained_edge_dof_indices<TBaseElem,ConstrainingTriangle, Edge, Face>(elem,fct,ind, vFace);
715 0 : if(dim >= FACE) constrained_edge_dof_indices<TBaseElem,ConstrainingQuadrilateral, Edge, Face>(elem,fct,ind, vFace);
716 : }
717 :
718 : // get dofs on hanging faces
719 0 : if (max_dofs(FACE) > 0){
720 0 : if(dim >= FACE) constrained_face_dof_indices<TBaseElem,ConstrainingTriangle, Face, Face>(elem,fct,ind,vFace);
721 0 : if(dim >= FACE) constrained_face_dof_indices<TBaseElem,ConstrainingQuadrilateral, Face, Face>(elem,fct,ind,vFace);
722 : }
723 :
724 : // return number of indices
725 0 : return ind.size();
726 : }
727 :
728 : template<typename TBaseElem>
729 0 : void DoFDistribution::indices_on_vertex(TBaseElem* elem, const ReferenceObjectID roid,
730 : LocalIndices& ind,
731 : const Grid::SecureVertexContainer& vElem) const
732 : {
733 : // get reference object id for subelement
734 : static const ReferenceObjectID subRoid = ROID_VERTEX;
735 :
736 : // add normal dofs
737 0 : for(size_t i = 0; i < vElem.size(); ++i)
738 : {
739 : // get subset index
740 0 : const int si = m_spMGSH->get_subset_index(vElem[i]);
741 :
742 : // loop all functions
743 0 : for(size_t fct = 0; fct < num_fct(); ++fct)
744 : {
745 : // check if function is defined on the subset
746 0 : if(!is_def_in_subset(fct, si)) continue;
747 :
748 : // get number of DoFs in this sub-geometric object
749 : const size_t numDoFsOnSub = num_fct_dofs(fct,subRoid,si);
750 :
751 : // Always no orientation needed
752 0 : if(!m_bGrouped)
753 : {
754 : // compute index
755 0 : const size_t index = obj_index(vElem[i]) + offset(subRoid,si,fct);
756 :
757 : // add dof to local indices
758 0 : for(size_t j = 0; j < numDoFsOnSub; ++j)
759 0 : ind.push_back_index(fct, index+j);
760 : }
761 : else
762 : {
763 : // compute index
764 0 : const size_t index = obj_index(vElem[i]);
765 : const size_t comp = offset(subRoid,si,fct);
766 :
767 : // add dof to local indices
768 0 : for(size_t j = 0; j < numDoFsOnSub; ++j)
769 0 : ind.push_back_multi_index(fct, index, comp+j);
770 : }
771 : } // end loop functions
772 : } // end loop subelement
773 :
774 0 : }
775 :
776 : template<typename TBaseElem, typename TSubBaseElem>
777 0 : void DoFDistribution::indices(TBaseElem* elem, const ReferenceObjectID roid,
778 : LocalIndices& ind,
779 : const typename Grid::traits<TSubBaseElem>::secure_container& vElem) const
780 : {
781 : // storage for offsets
782 : std::vector<size_t> vOrientOffset;
783 :
784 : // add normal dofs
785 0 : for(size_t i = 0; i < vElem.size(); ++i)
786 : {
787 : // get subelement
788 : TSubBaseElem* subElem = vElem[i];
789 :
790 : // get subset index
791 0 : const int si = m_spMGSH->get_subset_index(subElem);
792 :
793 : // get reference object id for subselement
794 0 : const ReferenceObjectID subRoid = subElem->reference_object_id();
795 :
796 : // loop all functions
797 0 : for(size_t fct = 0; fct < num_fct(); ++fct)
798 : {
799 : // check if function is defined on the subset
800 0 : if(!is_def_in_subset(fct, si)) continue;
801 :
802 : // get number of DoFs in this sub-geometric object
803 : const size_t numDoFsOnSub = num_fct_dofs(fct,subRoid,si);
804 :
805 : // Orientation is required: Thus, we compute the offsets, that are
806 : // no longer in the usual order [0, 1, 2, ...]. Orientation is
807 : // required if there are more than 1 dof on a subelement of a
808 : // finite element and thus, when gluing two elements together,
809 : // also the dofs on the subelements have to fit in order to
810 : // guarantee continuity. This is not needed for Vertices, since there
811 : // no distinction can be made when all dofs are at the same position.
812 : // This is also not needed for the highest dimension of a finite
813 : // element, since the dofs on this geometric object must not be
814 : // identified with other dofs.
815 0 : ComputeOrientationOffset(vOrientOffset, elem, subElem, i, lfeid(fct));
816 :
817 : UG_ASSERT(vOrientOffset.size() == numDoFsOnSub ||
818 : vOrientOffset.empty(), "Something wrong with orientation");
819 :
820 0 : if(!m_bGrouped)
821 : {
822 0 : const size_t index = obj_index(subElem) + offset(subRoid,si,fct);
823 :
824 0 : if(vOrientOffset.empty()){
825 0 : for(size_t j = 0; j < numDoFsOnSub; ++j)
826 0 : ind.push_back_index(fct, index + j);
827 : }else {
828 0 : for(size_t j = 0; j < numDoFsOnSub; ++j)
829 0 : ind.push_back_index(fct, index + vOrientOffset[j]);
830 : }
831 : }
832 : else
833 : {
834 : // compute index
835 0 : const size_t index = obj_index(subElem);
836 : const size_t comp = offset(subRoid,si,fct);
837 :
838 0 : if(vOrientOffset.empty()){
839 0 : for(size_t j = 0; j < numDoFsOnSub; ++j)
840 0 : ind.push_back_multi_index(fct, index, comp + j);
841 : }else{
842 0 : for(size_t j = 0; j < numDoFsOnSub; ++j)
843 0 : ind.push_back_multi_index(fct, index, comp + vOrientOffset[j]);
844 : }
845 : }
846 : } // end loop functions
847 : } // end loop subelement
848 :
849 0 : }
850 :
851 : template <typename TConstraining, typename TConstrained, typename TBaseElem>
852 0 : void DoFDistribution::
853 : constrained_vertex_indices(LocalIndices& ind,
854 : const typename Grid::traits<TBaseElem>::secure_container& vSubElem) const
855 : {
856 : // loop all edges
857 0 : for(size_t i = 0; i < vSubElem.size(); ++i)
858 : {
859 : // only constraining objects are of interest
860 0 : TConstraining* constrainingObj = dynamic_cast<TConstraining*>(vSubElem[i]);
861 0 : if(constrainingObj == NULL) continue;
862 :
863 : // loop constraining vertices
864 0 : for(size_t j = 0; j != constrainingObj->num_constrained_vertices(); ++j)
865 : {
866 : // get vertex
867 : TConstrained* vrt = constrainingObj->constrained_vertex(j);
868 :
869 : // get roid
870 0 : const ReferenceObjectID subRoid = vrt->reference_object_id();
871 :
872 : // get subset index
873 0 : int si = m_spMGSH->get_subset_index(vrt);
874 :
875 : // loop functions
876 0 : for(size_t fct = 0; fct < num_fct(); ++fct)
877 : {
878 : // check that function is defined on subset
879 0 : if(!is_def_in_subset(fct, si)) continue;
880 :
881 0 : if(!m_bGrouped)
882 : {
883 : // compute index
884 0 : const size_t index = obj_index(vrt) + offset(subRoid,si,fct);
885 :
886 : // add dof to local indices
887 : ind.push_back_index(fct, index);
888 : }
889 : else
890 : {
891 : // compute index
892 0 : const size_t index = obj_index(vrt);
893 : const size_t comp = offset(subRoid,si,fct);
894 :
895 : // add dof to local indices
896 : ind.push_back_multi_index(fct, index, comp);
897 : }
898 : }
899 : }
900 : }
901 0 : }
902 :
903 : // sort constrained edges by association to reference object vertices of coarse edge
904 : template <typename TBaseElem,typename TConstraining, typename TConstrained>
905 0 : void DoFDistribution::
906 : sort_constrained_edges(std::vector<size_t>& sortedInd,TBaseElem* elem,TConstraining* constrainingObj,size_t objIndex) const
907 : {
908 : static const int dim = TBaseElem::dim;
909 0 : ReferenceObjectID roid = elem->reference_object_id();
910 : const DimReferenceElement<dim>& refElem
911 : = ReferenceElementProvider::get<dim>(roid);
912 : // get edge belonging to reference id vertex 0 on edge
913 0 : const size_t vertexIndex = refElem.id(1,objIndex,0,0);
914 0 : sortedInd.resize(2);
915 : Vertex* vertex0 = NULL;
916 : // get child of vertex
917 : if (dim==2){
918 : Face* baseElem = dynamic_cast<Face*>(elem);
919 0 : vertex0 = multi_grid()->template get_child<Vertex,Vertex>(baseElem->vertex(vertexIndex),0);
920 : }
921 : if (dim==3){
922 : Volume* baseElem = dynamic_cast<Volume*>(elem);
923 0 : vertex0 = multi_grid()->template get_child<Vertex,Vertex>(baseElem->vertex(vertexIndex),0);
924 : }
925 : TConstrained* edg = constrainingObj->constrained_edge(0);
926 : bool found = false;
927 0 : for (size_t k=0;k<2;k++){
928 0 : Vertex* vrt = edg->vertex(k);
929 0 : if (vrt==vertex0){
930 : found = true;
931 : break;
932 : }
933 : }
934 0 : if (found==true){
935 0 : sortedInd[0]=0;
936 0 : sortedInd[1]=1;
937 : } else {
938 0 : sortedInd[0]=1;
939 0 : sortedInd[1]=0;
940 : // check
941 : bool found2 = false;
942 : TConstrained* otherEdge = constrainingObj->constrained_edge(1);
943 0 : for (size_t k=0;k<2;k++){
944 0 : Vertex* vrt = otherEdge->vertex(k);
945 0 : if (vrt==vertex0){
946 : found2 = true;
947 : break;
948 : }
949 : }
950 0 : if (found2==false) UG_THROW("no edge found belonging to vertex 0\n");
951 : }
952 0 : }
953 :
954 : // sort constrained faces as given by association to reference object vertices of coarse faces
955 : template <typename TBaseElem,typename TConstraining, typename TConstrained>
956 0 : void DoFDistribution::
957 : sort_constrained_faces(std::vector<size_t>& sortedInd,TBaseElem* elem,TConstraining* constrainingObj,size_t objIndex) const
958 : {
959 : static const int dim = TBaseElem::dim;
960 0 : ReferenceObjectID roid = elem->reference_object_id();
961 : const DimReferenceElement<dim>& refElem
962 : = ReferenceElementProvider::get<dim>(roid);
963 0 : const size_t numVrt = constrainingObj->num_vertices();
964 0 : sortedInd.resize(4);
965 : Vertex* vrt = NULL;
966 0 : Volume* baseElem = dynamic_cast<Volume*>(elem);
967 0 : for (size_t i=0;i<numVrt;i++){
968 0 : const size_t vertexIndex = refElem.id(2,objIndex,0,i);
969 0 : vrt = multi_grid()->template get_child<Vertex,Vertex>(baseElem->vertex(vertexIndex),0);
970 : // loop constrained faces to find face corresponding to vertex
971 : bool found = false;
972 0 : for (size_t j=0;j<numVrt;j++){
973 : TConstrained* face = constrainingObj->constrained_face(j);
974 0 : for (size_t k=0;k<face->num_vertices();k++){
975 0 : if (face->vertex(k)==vrt){
976 : found = true;
977 0 : sortedInd[i] = j;
978 0 : break;
979 : }
980 : }
981 : }
982 0 : if (found==false) UG_THROW("corresponding constrained object vertex not found");
983 : }
984 : // for triangle refinement inner face is still missing, it should be constrained_face(3)
985 0 : if (numVrt==3){
986 : // check if it is not face 3
987 0 : for (size_t i=0;i<3;i++) if (sortedInd[i]==3) {
988 : bool found = false;
989 0 : for (size_t j=0;j<3;j++){
990 0 : for (size_t k=0;k<3;k++){
991 0 : if (sortedInd[k]==j){
992 : found = true;
993 : break;
994 : }
995 : }
996 0 : if (found==false){
997 0 : sortedInd[3]=j;
998 0 : return;
999 : }
1000 : }
1001 : }
1002 0 : sortedInd[3]=3;
1003 : }
1004 : }
1005 :
1006 : template <typename TBaseElem,typename TConstraining, typename TConstrained, typename TSubElem>
1007 0 : void DoFDistribution::
1008 : constrained_edge_indices(TBaseElem* elem,LocalIndices& ind,
1009 : const typename Grid::traits<TSubElem>::secure_container& vSubElem) const
1010 : {
1011 : // loop all edges
1012 0 : for(size_t i = 0; i < vSubElem.size(); ++i)
1013 : {
1014 : // only constraining objects are of interest
1015 0 : TConstraining* constrainingObj = dynamic_cast<TConstraining*>(vSubElem[i]);
1016 0 : if(constrainingObj == NULL) continue;
1017 :
1018 : std::vector<size_t> sortedInd;
1019 0 : sort_constrained_edges<TBaseElem,TConstraining,TConstrained>(sortedInd,elem,constrainingObj,i);
1020 :
1021 : // loop constraining edges
1022 0 : for(size_t j = 0; j != constrainingObj->num_constrained_edges(); ++j)
1023 : {
1024 : // get edge
1025 0 : TConstrained* edg = constrainingObj->constrained_edge(sortedInd[j]);
1026 :
1027 : // get roid
1028 0 : const ReferenceObjectID subRoid = edg->reference_object_id();
1029 :
1030 : // get subset index
1031 0 : int si = m_spMGSH->get_subset_index(edg);
1032 :
1033 : // loop functions
1034 0 : for(size_t fct = 0; fct < num_fct(); ++fct)
1035 : {
1036 : // check that function is defined on subset
1037 0 : if(!is_def_in_subset(fct, si)) continue;
1038 :
1039 0 : if(!m_bGrouped)
1040 : {
1041 : // compute index
1042 0 : const size_t index = obj_index(edg) + offset(subRoid,si,fct);
1043 :
1044 : // add dof to local indices
1045 : ind.push_back_index(fct, index);
1046 : }
1047 : else
1048 : {
1049 : // compute index
1050 0 : const size_t index = obj_index(edg);
1051 : const size_t comp = offset(subRoid,si,fct);
1052 :
1053 : // add dof to local indices
1054 : ind.push_back_multi_index(fct, index, comp);
1055 : }
1056 : }
1057 : }
1058 : }
1059 0 : }
1060 :
1061 : template <typename TBaseElem,typename TConstraining, typename TConstrained, typename TSubElem>
1062 0 : void DoFDistribution::
1063 : constrained_face_indices(TBaseElem* elem,LocalIndices& ind,
1064 : const typename Grid::traits<TSubElem>::secure_container& vSubElem) const
1065 : {
1066 : // loop all faces
1067 0 : for(size_t i = 0; i < vSubElem.size(); ++i)
1068 : {
1069 : // only constraining objects are of interest
1070 0 : TConstraining* constrainingObj = dynamic_cast<TConstraining*>(vSubElem[i]);
1071 :
1072 0 : if(constrainingObj == NULL) continue;
1073 :
1074 : std::vector<size_t> sortedInd;
1075 0 : sort_constrained_faces<TBaseElem,TConstraining,TConstrained>(sortedInd,elem,constrainingObj,i);
1076 :
1077 : // loop constraining faces
1078 0 : for(size_t j = 0; j != constrainingObj->num_constrained_faces(); ++j)
1079 : {
1080 : // get face
1081 0 : TConstrained* face = constrainingObj->constrained_face(sortedInd[j]);
1082 :
1083 : // get roid
1084 0 : const ReferenceObjectID subRoid = face->reference_object_id();
1085 :
1086 : // get subset index
1087 0 : int si = m_spMGSH->get_subset_index(face);
1088 :
1089 : // loop functions
1090 0 : for(size_t fct = 0; fct < num_fct(); ++fct)
1091 : {
1092 : // check that function is defined on subset
1093 0 : if(!is_def_in_subset(fct, si)) continue;
1094 :
1095 0 : if(!m_bGrouped)
1096 : {
1097 : // compute index
1098 0 : const size_t index = obj_index(face) + offset(subRoid,si,fct);
1099 :
1100 : // add dof to local indices
1101 : ind.push_back_index(fct, index);
1102 : }
1103 : else
1104 : {
1105 : // compute index
1106 0 : const size_t index = obj_index(face);
1107 : const size_t comp = offset(subRoid,si,fct);
1108 :
1109 : // add dof to local indices
1110 : ind.push_back_multi_index(fct, index, comp);
1111 : }
1112 : }
1113 : }
1114 : }
1115 0 : }
1116 :
1117 : template<typename TBaseElem>
1118 0 : void DoFDistribution::_indices(TBaseElem* elem, LocalIndices& ind, bool bHang) const
1119 : {
1120 : // reference dimension
1121 : static const int dim = TBaseElem::dim;
1122 :
1123 : // resize the number of functions
1124 : ind.resize_fct(num_fct());
1125 0 : for(size_t fct = 0; fct < num_fct(); ++fct) ind.clear_dof(fct);
1126 :
1127 : // storage for (maybe needed) subelements
1128 : Grid::SecureVertexContainer vCorner;
1129 : Grid::SecureEdgeContainer vEdge;
1130 : Grid::SecureFaceContainer vFace;
1131 : Grid::SecureVolumeContainer vVol;
1132 :
1133 : // collect elements, if needed
1134 : if(dim >= VERTEX)
1135 0 : if(max_dofs(VERTEX) > 0) m_pMG->associated_elements_sorted(vCorner, elem);
1136 : if(dim >= EDGE)
1137 0 : if(max_dofs(EDGE) > 0 || bHang) m_pMG->associated_elements_sorted(vEdge, elem);
1138 : if(dim >= FACE)
1139 0 : if(max_dofs(FACE) > 0 || bHang) m_pMG->associated_elements_sorted(vFace, elem);
1140 : if(dim >= VOLUME)
1141 0 : if(max_dofs(VOLUME) > 0) m_pMG->associated_elements_sorted(vVol, elem);
1142 :
1143 : // get reference object id
1144 0 : const ReferenceObjectID roid = elem->reference_object_id();
1145 :
1146 : // get regular dofs on all subelements and the element itself
1147 : // use specialized function for vertices (since only one position and one reference object)
1148 0 : if(dim >= VERTEX && max_dofs(VERTEX) > 0) indices_on_vertex<TBaseElem>(elem, roid, ind, vCorner);
1149 0 : if(dim >= EDGE && max_dofs(EDGE) > 0) indices<TBaseElem, Edge>(elem, roid, ind, vEdge);
1150 0 : if(dim >= FACE && max_dofs(FACE) > 0) indices<TBaseElem, Face>(elem, roid, ind, vFace);
1151 0 : if(dim >= VOLUME && max_dofs(VOLUME) > 0) indices<TBaseElem, Volume>(elem, roid, ind, vVol);
1152 :
1153 : // If no hanging dofs are required, we're done
1154 0 : if(!bHang) return;
1155 :
1156 : // get dofs on hanging vertices
1157 0 : if (max_dofs(VERTEX) > 0)
1158 : {
1159 0 : if(dim >= EDGE) constrained_vertex_indices<ConstrainingEdge, Vertex, Edge>(ind, vEdge);
1160 0 : if(dim >= FACE) constrained_vertex_indices<ConstrainingQuadrilateral, Vertex, Face>(ind, vFace);
1161 : }
1162 :
1163 : // get dofs on hanging edges
1164 0 : if (max_dofs(EDGE) > 0){
1165 0 : if(dim >= EDGE) constrained_edge_indices<TBaseElem,ConstrainingEdge, Edge, Edge>(elem,ind, vEdge);
1166 0 : if(dim >= FACE) constrained_edge_indices<TBaseElem,ConstrainingTriangle, Edge, Face>(elem,ind, vFace);
1167 0 : if(dim >= FACE) constrained_edge_indices<TBaseElem,ConstrainingQuadrilateral, Edge, Face>(elem,ind, vFace);
1168 : }
1169 :
1170 : // get dofs on hanging faces
1171 0 : if (max_dofs(FACE) > 0){
1172 0 : if(dim >= FACE) constrained_face_indices<TBaseElem,ConstrainingTriangle, Face, Face>(elem,ind, vFace);
1173 0 : if(dim >= FACE) constrained_face_indices<TBaseElem,ConstrainingQuadrilateral, Face, Face>(elem,ind, vFace);
1174 : }
1175 :
1176 : // we're done
1177 : return;
1178 : }
1179 :
1180 :
1181 : template <typename TBaseElem>
1182 0 : void DoFDistribution::
1183 : changable_indices(std::vector<size_t>& vIndex,
1184 : const std::vector<TBaseElem*>& vElem) const
1185 : {
1186 : // Get connected indices
1187 0 : for(size_t i = 0; i < vElem.size(); ++i)
1188 : {
1189 : // Get Vertices of adjacent edges
1190 0 : TBaseElem* elem = vElem[i];
1191 : UG_ASSERT(m_spMGSH->get_subset_index(elem) >= 0, "Must have subset");
1192 :
1193 : // get adjacent index
1194 0 : const size_t adjInd = obj_index(elem);
1195 :
1196 : UG_ASSERT(adjInd < (size_t)1e10, "adjInd = " << adjInd); // <-- I get an adjInd of (size_t) (-1) here!
1197 :
1198 : // add to index list
1199 0 : vIndex.push_back(adjInd);
1200 : }
1201 0 : }
1202 :
1203 : ///////////////////////////////////////////////////////////////////////////////
1204 : // forwarding fcts
1205 : ///////////////////////////////////////////////////////////////////////////////
1206 :
1207 0 : void DoFDistribution::indices(Vertex* elem, LocalIndices& ind, bool bHang) const {_indices<Vertex>(elem, ind, bHang);}
1208 0 : void DoFDistribution::indices(Edge* elem, LocalIndices& ind, bool bHang) const {_indices<Edge>(elem, ind, bHang);}
1209 0 : void DoFDistribution::indices(Face* elem, LocalIndices& ind, bool bHang) const {_indices<Face>(elem, ind, bHang);}
1210 0 : void DoFDistribution::indices(Volume* elem, LocalIndices& ind, bool bHang) const {_indices<Volume>(elem, ind, bHang);}
1211 0 : void DoFDistribution::indices(GridObject* elem, LocalIndices& ind, bool bHang) const
1212 : {
1213 0 : switch(elem->base_object_id())
1214 : {
1215 0 : case VERTEX: return indices(static_cast<Vertex*>(elem), ind, bHang);
1216 0 : case EDGE: return indices(static_cast<Edge*>(elem), ind, bHang);
1217 0 : case FACE: return indices(static_cast<Face*>(elem), ind, bHang);
1218 0 : case VOLUME: return indices(static_cast<Volume*>(elem), ind, bHang);
1219 0 : default: UG_THROW("Geometric Base element not found.");
1220 : }
1221 : }
1222 :
1223 0 : size_t DoFDistribution::dof_indices(Vertex* elem, size_t fct, std::vector<DoFIndex>& ind, bool bHang, bool bClear) const {return _dof_indices<Vertex>(elem, fct, ind, bHang, bClear);}
1224 0 : size_t DoFDistribution::dof_indices(Edge* elem, size_t fct, std::vector<DoFIndex>& ind, bool bHang, bool bClear) const {return _dof_indices<Edge>(elem, fct, ind, bHang, bClear);}
1225 0 : size_t DoFDistribution::dof_indices(Face* elem, size_t fct, std::vector<DoFIndex>& ind, bool bHang, bool bClear) const {return _dof_indices<Face>(elem, fct, ind, bHang, bClear);}
1226 0 : size_t DoFDistribution::dof_indices(Volume* elem, size_t fct, std::vector<DoFIndex>& ind, bool bHang, bool bClear) const {return _dof_indices<Volume>(elem, fct, ind, bHang, bClear);}
1227 0 : size_t DoFDistribution::dof_indices(GridObject* elem, size_t fct, std::vector<DoFIndex>& ind, bool bHang, bool bClear) const
1228 : {
1229 0 : switch(elem->base_object_id())
1230 : {
1231 0 : case VERTEX: return dof_indices(static_cast<Vertex*>(elem), fct, ind, bHang, bClear);
1232 0 : case EDGE: return dof_indices(static_cast<Edge*>(elem), fct, ind, bHang, bClear);
1233 0 : case FACE: return dof_indices(static_cast<Face*>(elem), fct, ind, bHang, bClear);
1234 0 : case VOLUME: return dof_indices(static_cast<Volume*>(elem), fct, ind, bHang, bClear);
1235 0 : default: UG_THROW("Geometric Base element not found.");
1236 : }
1237 : }
1238 :
1239 0 : size_t DoFDistribution::inner_dof_indices(Vertex* elem, size_t fct, std::vector<DoFIndex>& ind, bool bHang) const {return _inner_dof_indices<Vertex>(elem, fct, ind, bHang);}
1240 0 : size_t DoFDistribution::inner_dof_indices(Edge* elem, size_t fct, std::vector<DoFIndex>& ind, bool bHang) const {return _inner_dof_indices<Edge>(elem, fct, ind, bHang);}
1241 0 : size_t DoFDistribution::inner_dof_indices(Face* elem, size_t fct, std::vector<DoFIndex>& ind, bool bHang) const {return _inner_dof_indices<Face>(elem, fct, ind, bHang);}
1242 0 : size_t DoFDistribution::inner_dof_indices(Volume* elem, size_t fct, std::vector<DoFIndex>& ind, bool bHang) const {return _inner_dof_indices<Volume>(elem, fct, ind, bHang);}
1243 0 : size_t DoFDistribution::inner_dof_indices(GridObject* elem, size_t fct, std::vector<DoFIndex>& ind, bool bClear) const
1244 : {
1245 0 : switch(elem->base_object_id())
1246 : {
1247 0 : case VERTEX: return inner_dof_indices(static_cast<Vertex*>(elem), fct, ind, bClear);
1248 0 : case EDGE: return inner_dof_indices(static_cast<Edge*>(elem), fct, ind, bClear);
1249 0 : case FACE: return inner_dof_indices(static_cast<Face*>(elem), fct, ind, bClear);
1250 0 : case VOLUME: return inner_dof_indices(static_cast<Volume*>(elem), fct, ind, bClear);
1251 0 : default: UG_THROW("Geometric Base element not found.");
1252 : }
1253 : }
1254 :
1255 :
1256 0 : size_t DoFDistribution::algebra_indices(Vertex* elem, std::vector<size_t>& ind, bool bClear) const {return _algebra_indices<Vertex>(elem, ind, bClear);}
1257 0 : size_t DoFDistribution::algebra_indices(Edge* elem, std::vector<size_t>& ind, bool bClear) const {return _algebra_indices<Edge>(elem, ind, bClear);}
1258 0 : size_t DoFDistribution::algebra_indices(Face* elem, std::vector<size_t>& ind, bool bClear) const {return _algebra_indices<Face>(elem, ind, bClear);}
1259 0 : size_t DoFDistribution::algebra_indices(Volume* elem, std::vector<size_t>& ind, bool bClear) const {return _algebra_indices<Volume>(elem, ind, bClear);}
1260 0 : size_t DoFDistribution::algebra_indices(GridObject* elem, std::vector<size_t>& ind, bool bClear) const
1261 : {
1262 0 : switch(elem->base_object_id())
1263 : {
1264 0 : case VERTEX: return algebra_indices(static_cast<Vertex*>(elem), ind, bClear);
1265 0 : case EDGE: return algebra_indices(static_cast<Edge*>(elem), ind, bClear);
1266 0 : case FACE: return algebra_indices(static_cast<Face*>(elem), ind, bClear);
1267 0 : case VOLUME: return algebra_indices(static_cast<Volume*>(elem), ind, bClear);
1268 0 : default: UG_THROW("Geometric Base element not found.");
1269 : }
1270 : }
1271 :
1272 0 : size_t DoFDistribution::inner_algebra_indices(Vertex* elem, std::vector<size_t>& ind, bool bClear) const {return _inner_algebra_indices<Vertex>(elem, ind, bClear);}
1273 0 : size_t DoFDistribution::inner_algebra_indices(Edge* elem, std::vector<size_t>& ind, bool bClear) const {return _inner_algebra_indices<Edge>(elem, ind, bClear);}
1274 0 : size_t DoFDistribution::inner_algebra_indices(Face* elem, std::vector<size_t>& ind, bool bClear) const {return _inner_algebra_indices<Face>(elem, ind, bClear);}
1275 0 : size_t DoFDistribution::inner_algebra_indices(Volume* elem, std::vector<size_t>& ind, bool bClear) const {return _inner_algebra_indices<Volume>(elem, ind, bClear);}
1276 0 : size_t DoFDistribution::inner_algebra_indices(GridObject* elem, std::vector<size_t>& ind, bool bClear) const
1277 : {
1278 0 : switch(elem->base_object_id())
1279 : {
1280 0 : case VERTEX: return inner_algebra_indices(static_cast<Vertex*>(elem), ind, bClear);
1281 0 : case EDGE: return inner_algebra_indices(static_cast<Edge*>(elem), ind, bClear);
1282 0 : case FACE: return inner_algebra_indices(static_cast<Face*>(elem), ind, bClear);
1283 0 : case VOLUME: return inner_algebra_indices(static_cast<Volume*>(elem), ind, bClear);
1284 0 : default: UG_THROW("Geometric Base element not found.");
1285 : }
1286 : }
1287 :
1288 :
1289 : ////////////////////////////////////////////////////////////////////////////////
1290 : // Resize management
1291 : ////////////////////////////////////////////////////////////////////////////////
1292 :
1293 0 : void DoFDistribution::manage_grid_function(IGridFunction& gridFct)
1294 : {
1295 : // if already registered, we're done
1296 0 : if(std::find(m_vpGridFunction.begin(), m_vpGridFunction.end(), &gridFct)
1297 : != m_vpGridFunction.end())
1298 : return;
1299 :
1300 : // add to managed functions
1301 0 : m_vpGridFunction.push_back(&gridFct);
1302 : }
1303 :
1304 0 : void DoFDistribution::unmanage_grid_function(IGridFunction& gridFct)
1305 : {
1306 0 : m_vpGridFunction.erase(
1307 0 : std::remove(m_vpGridFunction.begin(), m_vpGridFunction.end(), &gridFct)
1308 : , m_vpGridFunction.end());
1309 0 : }
1310 :
1311 0 : void DoFDistribution::permute_values(const std::vector<size_t>& vIndNew)
1312 : {
1313 : // swap values of handled grid functions
1314 0 : for(size_t i = 0; i < m_vpGridFunction.size(); ++i)
1315 0 : m_vpGridFunction[i]->permute_values(vIndNew);
1316 0 : }
1317 :
1318 0 : void DoFDistribution::copy_values(const std::vector<std::pair<size_t, size_t> >& vIndexMap,
1319 : bool bDisjunct)
1320 : {
1321 : // swap values of handled grid functions
1322 0 : for(size_t i = 0; i < m_vpGridFunction.size(); ++i)
1323 0 : m_vpGridFunction[i]->copy_values(vIndexMap, bDisjunct);
1324 0 : }
1325 :
1326 0 : void DoFDistribution::resize_values(size_t newSize)
1327 : {
1328 : // swap values of handled grid functions
1329 0 : for(size_t i = 0; i < m_vpGridFunction.size(); ++i)
1330 0 : m_vpGridFunction[i]->resize_values(newSize);
1331 0 : }
1332 :
1333 : ////////////////////////////////////////////////////////////////////////////////
1334 : // Init DoFs
1335 : ////////////////////////////////////////////////////////////////////////////////
1336 :
1337 :
1338 : template <typename TBaseElem>
1339 0 : void DoFDistribution::reinit()
1340 : {
1341 : typedef typename traits<TBaseElem>::iterator iterator;
1342 : static const int dim = TBaseElem::dim;
1343 :
1344 : // check if indices in the dimension
1345 0 : if(max_dofs(dim) == 0) return;
1346 :
1347 : // SURFACE
1348 0 : if(grid_level().type() == GridLevel::SURFACE){
1349 :
1350 : // in order to also cater for some seldom occurring parallel cases, we have
1351 : // to perform a slightly cumbersome iteration here. The basic idea is the following:
1352 : // Each PURE_SURFACE element requires a dof. Furthermore all SHADOWING elements
1353 : // which are not SHADOWED require a dof. SHADOWED_NONCOPY elements also require
1354 : // a dof and SHADOWED_COPY elements have to copy their dofs from their SHADOWING
1355 : // element.
1356 : // In parallel, however, a problem occurs. The associated SHADOWING element of
1357 : // a SHADOW_COPY element does not necessarily exist on the same process
1358 : // (should occur on vertices only).
1359 : // We thus can't simply iterate over PURE_SURFACE | SHADOWING and pass values on
1360 : // to parents. Instead we have to iterate over all candidates and also give a
1361 : // dof to SHADOW_COPY elements which don't have children.
1362 :
1363 : const SurfaceView& sv = *m_spSurfView;
1364 : MultiGrid& mg = *m_spMG;
1365 :
1366 0 : for(int si = 0; si < num_subsets(); ++si)
1367 : {
1368 : // skip if no dofs shall exist on the given subset
1369 0 : if(max_dofs(dim, si) == 0) continue;
1370 :
1371 : // iterate over all surface elements (including shadows)
1372 : iterator iter = begin<TBaseElem>(si, SurfaceView::ALL);
1373 : iterator iterEnd = end<TBaseElem>(si, SurfaceView::ALL);
1374 :
1375 0 : for(; iter != iterEnd; ++iter){
1376 : TBaseElem* elem = *iter;
1377 0 : if(sv.is_contained(elem, grid_level(), SurfaceView::SHADOW_RIM_COPY)){
1378 0 : if(mg.num_children<TBaseElem>(elem) > 0){
1379 : TBaseElem* child = mg.get_child<TBaseElem>(elem, 0);
1380 0 : if(sv.is_contained(child, grid_level(), SurfaceView::SURFACE_RIM))
1381 0 : continue;
1382 : }
1383 : }
1384 :
1385 : // create a dof and copy it down to SHADOW_COPY parents
1386 0 : const ReferenceObjectID roid = elem->reference_object_id();
1387 0 : add(elem, roid, si);
1388 :
1389 0 : TBaseElem* p = dynamic_cast<TBaseElem*>(mg.get_parent(elem));
1390 0 : while(p && sv.is_contained(p, grid_level(), SurfaceView::SHADOW_RIM_COPY)){
1391 0 : obj_index(p) = obj_index(elem);
1392 0 : p = dynamic_cast<TBaseElem*>(mg.get_parent(p));
1393 : }
1394 : }
1395 : } // end subset
1396 : }
1397 :
1398 : // LEVEL
1399 0 : else if(grid_level().type() == GridLevel::LEVEL){
1400 :
1401 0 : for(int si = 0; si < num_subsets(); ++si)
1402 : {
1403 : // skip if no dofs to be distributed
1404 0 : if(max_dofs(dim, si) == 0) continue;
1405 :
1406 : // get iterators of elems
1407 0 : iterator iter = begin<TBaseElem>(si);
1408 0 : iterator iterEnd = end<TBaseElem>(si);
1409 :
1410 : // loop elems
1411 0 : for(; iter != iterEnd; ++iter)
1412 : {
1413 : // get vertex
1414 : TBaseElem* elem = *iter;
1415 :
1416 : // get roid
1417 0 : const ReferenceObjectID roid = elem->reference_object_id();
1418 :
1419 : // add element
1420 0 : add(elem, roid, si);
1421 : }
1422 : }
1423 :
1424 : }
1425 : else{
1426 0 : UG_THROW("DoFDistribution: GridLevel-Type"<<grid_level().type()<<" not supported");
1427 : }
1428 : }
1429 :
1430 0 : void DoFDistribution::reinit()
1431 : {
1432 0 : m_numIndex = 0;
1433 0 : m_vNumIndexOnSubset.resize(0);
1434 0 : m_vNumIndexOnSubset.resize(num_subsets(), 0);
1435 :
1436 0 : if(max_dofs(VERTEX)) reinit<Vertex>();
1437 0 : if(max_dofs(EDGE)) reinit<Edge>();
1438 0 : if(max_dofs(FACE)) reinit<Face>();
1439 0 : if(max_dofs(VOLUME)) reinit<Volume>();
1440 :
1441 : #ifdef UG_PARALLEL
1442 : reinit_layouts_and_communicator();
1443 : #endif
1444 0 : }
1445 :
1446 :
1447 : #ifdef UG_PARALLEL
1448 : void DoFDistribution::reinit_layouts_and_communicator()
1449 : {
1450 : pcl::ProcessCommunicator commWorld;
1451 :
1452 : // -----------------------------------
1453 : // CREATE PROCESS COMMUNICATOR
1454 : // -----------------------------------
1455 : // The idea of local processes is to exclude processes from
1456 : // e.g. norm computation that does not have a grid on a given
1457 : // level. If no DoFs exist on the level, that level is excluded from
1458 : // norm computations. In those cases the process votes false for
1459 : // the subcommunicator.
1460 :
1461 : // choose if this process participates
1462 : bool participate = !commWorld.empty() && (num_indices() > 0);
1463 :
1464 : // create process communicator for interprocess layouts
1465 : layouts()->proc_comm() = commWorld.create_sub_communicator(participate);
1466 :
1467 : // -----------------------------------
1468 : // CREATE INDEX LAYOUTS ON LEVEL
1469 : // -----------------------------------
1470 :
1471 : reinit_index_layout(layouts()->master(), INT_H_MASTER);
1472 : reinit_index_layout(layouts()->slave(), INT_H_SLAVE);
1473 : reinit_index_layout(layouts()->vertical_slave(), INT_V_SLAVE);
1474 :
1475 : // vertical layouts for ghosts
1476 : if(grid_level().ghosts()){
1477 : reinit_index_layout(layouts()->vertical_master(), INT_V_MASTER);
1478 : }else{
1479 : layouts()->vertical_master().clear();
1480 : }
1481 : }
1482 :
1483 : void DoFDistribution::reinit_index_layout(IndexLayout& layout, int keyType)
1484 : {
1485 : // clear layout
1486 : layout.clear();
1487 :
1488 : // add the index from grid layouts
1489 : if(max_dofs(VERTEX)) add_indices_from_layouts<Vertex>(layout, keyType);
1490 : if(max_dofs(EDGE)) add_indices_from_layouts<Edge>(layout, keyType);
1491 : if(max_dofs(FACE)) add_indices_from_layouts<Face>(layout, keyType);
1492 : if(max_dofs(VOLUME)) add_indices_from_layouts<Volume>(layout, keyType);
1493 :
1494 : // touching an interface means creation. Thus we remove the empty interfaces
1495 : // to avoid storage, communication (should not happen any longer) etc...
1496 : pcl::RemoveEmptyInterfaces(layout);
1497 : }
1498 :
1499 : template <typename TBaseElem>
1500 : void DoFDistribution::
1501 : add_indices_from_layouts(IndexLayout& indexLayout,int keyType)
1502 : {
1503 : // get the grid layout map
1504 : GridLayoutMap& layoutMap = m_spMG->distributed_grid_manager()->grid_layout_map();
1505 :
1506 : // check if layout present
1507 : if(!layoutMap.has_layout<TBaseElem>(keyType)) return;
1508 :
1509 : // types
1510 : typedef typename GridLayoutMap::Types<TBaseElem>::Layout::LevelLayout TLayout;
1511 : typedef typename TLayout::iterator InterfaceIterator;
1512 : typedef typename TLayout::Interface ElemInterface;
1513 : typedef typename ElemInterface::iterator ElemIterator;
1514 : typedef IndexLayout::Interface IndexInterface;
1515 :
1516 : // vector for algebra indices
1517 : std::vector<size_t> vIndex;
1518 :
1519 : // find levels to loop
1520 : int lvlTo = (grid_level().top() ? (multi_grid()->num_levels()-1) : grid_level().level());
1521 : int lvlFrom = (grid_level().is_surface() ? 0 : lvlTo);
1522 :
1523 : // loop all level
1524 : for(int lvl = lvlFrom; lvl <= lvlTo; ++lvl)
1525 : {
1526 : // get element layout
1527 : TLayout& elemLayout = layoutMap.get_layout<TBaseElem>(keyType).layout_on_level(lvl);
1528 :
1529 : // iterate over all grid element interfaces
1530 : for(InterfaceIterator iIter = elemLayout.begin();
1531 : iIter != elemLayout.end(); ++iIter)
1532 : {
1533 : // get a grid element interface
1534 : ElemInterface& elemInterface = elemLayout.interface(iIter);
1535 :
1536 : // get a corresponding index interface
1537 : IndexInterface& indexInterface = indexLayout.interface(elemLayout.proc_id(iIter));
1538 :
1539 : // iterate over entries in the grid element interface
1540 : for(ElemIterator eIter = elemInterface.begin();
1541 : eIter != elemInterface.end(); ++eIter)
1542 : {
1543 : // get the grid element
1544 : typename ElemInterface::Element elem = elemInterface.get_element(eIter);
1545 :
1546 : // check if element is a surface element
1547 : // Using SURFACE will also admit SHADOW elements from a level distribution.
1548 : // SHADOW_RIM elements must not be included; otherwise, any communication
1549 : // between hMaster and hSlaves would be done twice (also by the corresponding
1550 : // shadowing element, which references the same DoF), which is why we use
1551 : // SURFACE here instead of, e.g., ALL.
1552 : if(!m_spSurfView->is_contained(elem, grid_level(), SurfaceView::SURFACE))
1553 : continue;
1554 :
1555 : // add the indices to the interface
1556 : inner_algebra_indices(elem, vIndex);
1557 : for(size_t i = 0; i < vIndex.size(); ++i)
1558 : indexInterface.push_back(vIndex[i]);
1559 : }
1560 : }
1561 : }
1562 : }
1563 : #endif
1564 :
1565 : ////////////////////////////////////////////////////////////////////////////////
1566 : // DoF Statistic
1567 : ////////////////////////////////////////////////////////////////////////////////
1568 :
1569 :
1570 : template <typename TBaseElem>
1571 0 : void DoFDistribution::sum_dof_count(DoFCount& cnt) const
1572 : {
1573 : typedef typename traits<TBaseElem>::const_iterator iterator;
1574 :
1575 : const SurfaceView& sv = *m_spSurfView;
1576 0 : DistributedGridManager* spDstGrdMgr = sv.subset_handler()->grid()->distributed_grid_manager();
1577 :
1578 : // get iterators of elems
1579 0 : iterator iter = begin<TBaseElem>();
1580 0 : iterator iterEnd = end<TBaseElem>();
1581 :
1582 : // loop elems
1583 0 : for(; iter != iterEnd; ++iter){
1584 : TBaseElem* elem = *iter;
1585 0 : const ReferenceObjectID roid = elem->reference_object_id();
1586 0 : const int si = sv.subset_handler()->get_subset_index(elem);
1587 :
1588 : // Surface State
1589 0 : SurfaceView::SurfaceState SurfaceState = sv.surface_state(elem, grid_level());
1590 :
1591 : // skip shadow-rim-copy: those are already numbered by their SURFACE_RIM
1592 0 : if(SurfaceState.contains(SurfaceView::MG_SHADOW_RIM_COPY))
1593 0 : continue;
1594 :
1595 : // Interface State
1596 : byte InterfaceState = ES_NONE;
1597 0 : if(spDstGrdMgr) InterfaceState = spDstGrdMgr->get_status(elem);
1598 :
1599 : // loop all functions
1600 0 : for(size_t fct = 0; fct < num_fct(); ++fct){
1601 :
1602 : const size_t numDoF = num_fct_dofs(fct,roid,si);
1603 :
1604 0 : cnt.add(fct, si, SurfaceState, InterfaceState, numDoF);
1605 : }
1606 : }
1607 0 : }
1608 :
1609 0 : DoFCount DoFDistribution::dof_count() const
1610 : {
1611 : PROFILE_FUNC();
1612 0 : DoFCount cnt(grid_level(), dof_distribution_info());
1613 :
1614 0 : if(max_dofs(VERTEX)) sum_dof_count<Vertex>(cnt);
1615 0 : if(max_dofs(EDGE)) sum_dof_count<Edge>(cnt);
1616 0 : if(max_dofs(FACE)) sum_dof_count<Face>(cnt);
1617 0 : if(max_dofs(VOLUME)) sum_dof_count<Volume>(cnt);
1618 :
1619 0 : return cnt;
1620 : }
1621 :
1622 :
1623 : ////////////////////////////////////////////////////////////////////////////////
1624 : // Permutation of indices
1625 : ////////////////////////////////////////////////////////////////////////////////
1626 :
1627 : template <typename TBaseElem>
1628 0 : void DoFDistribution::
1629 : get_connections(std::vector<std::vector<size_t> >& vvConnection) const
1630 : {
1631 : // dimension of Base Elem
1632 : static const int dim = TBaseElem::dim;
1633 :
1634 : // Adjacent geometric objects
1635 : std::vector<Vertex*> vVrts;
1636 : std::vector<Edge*> vEdges;
1637 : std::vector<Face*> vFaces;
1638 : std::vector<Volume*> vVols;
1639 :
1640 : // Iterators
1641 : typedef typename traits<TBaseElem>::const_iterator const_iterator;
1642 0 : const_iterator iterEnd = end<TBaseElem>();
1643 0 : const_iterator iter = begin<TBaseElem>();
1644 :
1645 : // loop elem
1646 0 : for(; iter != iterEnd; ++iter)
1647 : {
1648 : // Get elem
1649 : TBaseElem* elem = *iter;
1650 : std::vector<size_t> vIndex;
1651 :
1652 : // Get connected elements
1653 0 : if(dim >= VERTEX && max_dofs(VERTEX) > 0) {
1654 : collect_associated(vVrts, elem);
1655 0 : changable_indices<Vertex>(vIndex, vVrts);
1656 : }
1657 0 : if(dim >= EDGE && max_dofs(EDGE) > 0) {
1658 : collect_associated(vEdges, elem);
1659 0 : changable_indices<Edge>(vIndex, vEdges);
1660 : }
1661 0 : if(dim >= FACE && max_dofs(FACE) > 0) {
1662 : collect_associated(vFaces, elem);
1663 0 : changable_indices<Face>(vIndex, vFaces);
1664 : }
1665 0 : if(dim >= VOLUME && max_dofs(VOLUME) > 0) {
1666 : collect_associated(vVols, elem);
1667 0 : changable_indices<Volume>(vIndex, vVols);
1668 : }
1669 :
1670 : // remove doubles
1671 0 : std::sort(vIndex.begin(), vIndex.end());
1672 0 : vIndex.erase(std::unique(vIndex.begin(), vIndex.end()), vIndex.end());
1673 :
1674 : // add coupling to adjacency graph
1675 : std::vector<size_t>::iterator it;
1676 :
1677 0 : for(size_t i = 0; i < vIndex.size(); ++i){
1678 0 : for(size_t j = i; j < vIndex.size(); ++j)
1679 : {
1680 0 : const size_t iIndex = vIndex[i];
1681 0 : const size_t jIndex = vIndex[j];
1682 :
1683 : // add connection (iIndex->jIndex)
1684 0 : it = find(vvConnection[iIndex].begin(), vvConnection[iIndex].end(), jIndex);
1685 0 : if(it == vvConnection[iIndex].end()) vvConnection[iIndex].push_back(jIndex);
1686 :
1687 :
1688 : // add the opposite direction (adjInd -> index)
1689 0 : it = find(vvConnection[jIndex].begin(), vvConnection[jIndex].end(), iIndex);
1690 0 : if(it == vvConnection[jIndex].end()) vvConnection[jIndex].push_back(iIndex);
1691 : }
1692 : }
1693 : }
1694 0 : }
1695 :
1696 0 : void DoFDistribution::
1697 : get_connections(std::vector<std::vector<size_t> >& vvConnection) const
1698 : {
1699 : vvConnection.clear();
1700 :
1701 : // if no subset given, we're done
1702 0 : if(num_subsets() == 0) return;
1703 :
1704 : // check that in all subsets same number of functions and at least one
1705 : // if this is not the case for non-grouped DoFs, we cannot allow reordering
1706 0 : if(!grouped())
1707 : {
1708 : size_t numDoFs = 0;
1709 0 : for(int si = 0; si < num_subsets(); ++si)
1710 0 : for(int i = 0; i < NUM_REFERENCE_OBJECTS; ++i)
1711 : {
1712 : const ReferenceObjectID roid = (ReferenceObjectID) i;
1713 0 : if(num_dofs(roid,si) == 0) continue;
1714 0 : if(numDoFs == 0) {numDoFs = num_dofs(roid,si); continue;}
1715 0 : if(num_dofs(roid,si) != numDoFs)
1716 0 : UG_THROW("DoFDistribution::get_connections: "
1717 : "Currently only implemented iff same number of DoFs on"
1718 : " all geometric objects in all subsets: \n"
1719 : "num_dofs("<<roid<<","<<si<<")="<<num_dofs(roid,si)<<
1720 : ", but previously found "<<numDoFs);
1721 : }
1722 : }
1723 :
1724 : // clear neighbors
1725 0 : vvConnection.resize(num_indices());
1726 :
1727 0 : get_connections<Vertex>(vvConnection);
1728 0 : get_connections<Edge>(vvConnection);
1729 0 : get_connections<Face>(vvConnection);
1730 0 : get_connections<Volume>(vvConnection);
1731 : }
1732 :
1733 : template <typename TBaseElem>
1734 0 : void DoFDistribution::permute_indices(const std::vector<size_t>& vNewInd)
1735 : {
1736 : // loop Vertices
1737 : typedef typename traits<TBaseElem>::const_iterator const_iterator;
1738 :
1739 0 : const_iterator iterEnd = end<TBaseElem>(SurfaceView::ALL);
1740 0 : const_iterator iter = begin<TBaseElem>(SurfaceView::ALL);
1741 0 : for(; iter != iterEnd; ++iter)
1742 : {
1743 : // get vertex
1744 : TBaseElem* elem = *iter;
1745 :
1746 : // get current (old) index
1747 0 : const size_t oldIndex = obj_index(elem);
1748 :
1749 : // replace old index by new one
1750 0 : obj_index(elem) = vNewInd[oldIndex];
1751 : }
1752 0 : }
1753 :
1754 0 : void DoFDistribution::permute_indices(const std::vector<size_t>& vNewInd)
1755 : {
1756 0 : if(max_dofs(VERTEX)) permute_indices<Vertex>(vNewInd);
1757 0 : if(max_dofs(EDGE)) permute_indices<Edge>(vNewInd);
1758 0 : if(max_dofs(FACE)) permute_indices<Face>(vNewInd);
1759 0 : if(max_dofs(VOLUME)) permute_indices<Volume>(vNewInd);
1760 :
1761 : #ifdef UG_PARALLEL
1762 : reinit_layouts_and_communicator();
1763 : #endif
1764 :
1765 : // permute indices in associated vectors
1766 0 : permute_values(vNewInd);
1767 0 : }
1768 :
1769 : } // end namespace ug
|
