Line data Source code
1 : /*
2 : * Copyright (c) 2012-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_position_util.h"
34 : #include "approximation_space.h"
35 : #include "lib_disc/domain.h"
36 : #include "lib_disc/domain_util.h"
37 : #include "lib_disc/domain_traits.h"
38 :
39 : #include "lib_disc/local_finite_element/local_finite_element_provider.h"
40 : #include "lib_disc/local_finite_element/local_dof_set.h"
41 : #include "lib_disc/reference_element/reference_mapping_provider.h"
42 : #include "lib_disc/reference_element/reference_element_util.h"
43 :
44 : namespace ug{
45 :
46 : ////////////////////////////////////////////////////////////////////////////////
47 : // DoFPositions
48 : ////////////////////////////////////////////////////////////////////////////////
49 :
50 : template <int refDim, int dim>
51 0 : bool InnerDoFPositionElem(std::vector<MathVector<dim> >& vPos, const ReferenceObjectID roid,
52 : const std::vector<MathVector<dim> >& vVertPos, const LFEID& lfeID)
53 : {
54 : // get local dof set
55 : const DimLocalDoFSet<refDim>& lds =
56 0 : LocalFiniteElementProvider::get_dofs<refDim>(roid, lfeID);
57 :
58 : // create a reference mapping
59 : DimReferenceMapping<refDim,dim>& map =
60 : ReferenceMappingProvider::get<refDim,dim>(roid, vVertPos);
61 :
62 : // clear pos
63 : vPos.clear();
64 :
65 : // bool flag if position is exact, or no exact position available for shapes
66 : bool bExact = true;
67 :
68 : // loop all shape functions
69 0 : for(size_t sh = 0; sh < lds.num_sh(); ++sh)
70 : {
71 : // check if dof in interior
72 0 : if(lds.local_dof(sh).dim() != refDim) continue;
73 :
74 : // get local position
75 : MathVector<refDim> locPos;
76 0 : bExact &= lds.position(sh, locPos);
77 :
78 : // map to global position
79 : MathVector<dim> globPos;
80 0 : map.local_to_global(globPos, locPos);
81 :
82 : // add
83 0 : vPos.push_back(globPos);
84 : }
85 :
86 : // return if positions are given exactly
87 0 : return bExact;
88 : };
89 :
90 : template <int dim>
91 0 : bool InnerDoFPositionVertex(std::vector<MathVector<dim> >& vPos, const ReferenceObjectID roid,
92 : const std::vector<MathVector<dim> >& vVertPos, const LFEID& lfeID)
93 : {
94 : UG_ASSERT(vVertPos.size() == 1, "Vertex should have only on inner Vertex");
95 :
96 : // get local dof set
97 0 : const CommonLocalDoFSet& lds = LocalFiniteElementProvider::get_dofs(lfeID);
98 :
99 : // clear pos
100 : vPos.clear();
101 :
102 : // loop all shape functions
103 0 : for(int sh = 0; sh < lds.num_dof(ROID_VERTEX); ++sh)
104 0 : vPos.push_back(vVertPos[0]);
105 :
106 : // return if positions are given exactly
107 0 : return true;
108 : };
109 :
110 :
111 : template <int dim>
112 0 : bool InnerDoFPosition(std::vector<MathVector<dim> >& vPos, const ReferenceObjectID roid,
113 : const std::vector<MathVector<dim> >& vCornerCoord, const LFEID& lfeID)
114 : {
115 0 : switch(ReferenceElementDimension(roid))
116 : {
117 0 : case VERTEX: return InnerDoFPositionVertex<dim>(vPos, roid, vCornerCoord, lfeID);
118 0 : case EDGE: return InnerDoFPositionElem<1,dim>(vPos, roid, vCornerCoord, lfeID);
119 0 : case FACE: return InnerDoFPositionElem<2,dim>(vPos, roid, vCornerCoord, lfeID);
120 0 : case VOLUME: return InnerDoFPositionElem<3,dim>(vPos, roid, vCornerCoord, lfeID);
121 0 : default: UG_THROW("Base Object type not found.");
122 : }
123 : }
124 :
125 : template <typename TDomain>
126 0 : bool InnerDoFPosition(std::vector<MathVector<TDomain::dim> >& vPos, GridObject* elem,
127 : const TDomain& domain, const LFEID& lfeID)
128 : {
129 : // reference object id
130 0 : const ReferenceObjectID roid = elem->reference_object_id();
131 :
132 : // get the vertices
133 : std::vector<MathVector<TDomain::dim> > vVertPos;
134 0 : switch(elem->base_object_id())
135 : {
136 : case VERTEX: CollectCornerCoordinates(vVertPos, *static_cast<Vertex*>(elem), domain, true); break;
137 : case EDGE: CollectCornerCoordinates(vVertPos, *static_cast<Edge*>(elem), domain, true); break;
138 : case FACE: CollectCornerCoordinates(vVertPos, *static_cast<Face*>(elem), domain, true); break;
139 : case VOLUME: CollectCornerCoordinates(vVertPos, *static_cast<Volume*>(elem), domain, true); break;
140 0 : default: UG_THROW("Base Object type not found.");
141 : }
142 :
143 : // forward
144 0 : return InnerDoFPosition<TDomain::dim>(vPos, roid, vVertPos, lfeID);
145 0 : }
146 :
147 :
148 :
149 : template <int refDim, int dim>
150 0 : bool DoFPositionElem(std::vector<MathVector<dim> >& vPos, const ReferenceObjectID roid,
151 : const std::vector<MathVector<dim> >& vVertPos, const LFEID& lfeID)
152 : {
153 : // get local shape function set
154 : const DimLocalDoFSet<refDim>& lds
155 0 : = LocalFiniteElementProvider::get_dofs<refDim>(roid, lfeID);
156 :
157 : // create a reference mapping
158 : const DimReferenceMapping<refDim,dim>& map =
159 : ReferenceMappingProvider::get<refDim,dim>(roid, vVertPos);
160 :
161 : // clear pos
162 0 : vPos.resize(lds.num_sh());
163 :
164 : // bool flag if position is exact, or no exact position available for shapes
165 : bool bExact = true;
166 :
167 : // loop all shape functions
168 0 : for(size_t sh = 0; sh < lds.num_sh(); ++sh)
169 : {
170 : // get local position
171 : MathVector<refDim> locPos;
172 0 : bExact &= lds.position(sh, locPos);
173 :
174 : // map to global position
175 0 : map.local_to_global(vPos[sh], locPos);
176 : }
177 :
178 : // return if positions are given exactly
179 0 : return bExact;
180 : };
181 :
182 : template <int dim>
183 0 : bool DoFPositionVertex(std::vector<MathVector<dim> >& vPos, const ReferenceObjectID roid,
184 : const std::vector<MathVector<dim> >& vVertPos, const LFEID& lfeID)
185 : {
186 : // get local dof set
187 0 : const CommonLocalDoFSet& lds = LocalFiniteElementProvider::get_dofs(lfeID);
188 :
189 : // clear pos
190 : vPos.clear();
191 :
192 : // loop all shape functions
193 0 : for(size_t co = 0; co < vVertPos.size(); ++co)
194 0 : for(int sh = 0; sh < lds.num_dof(ROID_VERTEX); ++sh)
195 0 : vPos.push_back(vVertPos[co]);
196 :
197 : // return if positions are given exactly
198 0 : return true;
199 : };
200 :
201 : template <int dim>
202 0 : bool DoFPosition(std::vector<MathVector<dim> >& vPos, const ReferenceObjectID roid,
203 : const std::vector<MathVector<dim> >& vCornerCoord, const LFEID& lfeID)
204 : {
205 0 : switch(ReferenceElementDimension(roid))
206 : {
207 0 : case VERTEX: return DoFPositionVertex<dim>(vPos, roid, vCornerCoord, lfeID);
208 0 : case EDGE: return DoFPositionElem<1,dim>(vPos, roid, vCornerCoord, lfeID);
209 0 : case FACE: return DoFPositionElem<2,dim>(vPos, roid, vCornerCoord, lfeID);
210 0 : case VOLUME: return DoFPositionElem<3,dim>(vPos, roid, vCornerCoord, lfeID);
211 0 : default: UG_THROW("Base Object type not found.");
212 : }
213 : }
214 :
215 : template <typename TDomain>
216 0 : bool DoFPosition(std::vector<MathVector<TDomain::dim> >& vPos, GridObject* elem,
217 : const TDomain& domain, const LFEID& lfeID)
218 : {
219 : // reference object id
220 0 : const ReferenceObjectID roid = elem->reference_object_id();
221 :
222 : // get the vertices
223 : std::vector<MathVector<TDomain::dim> > vVertPos;
224 0 : switch(elem->base_object_id())
225 : {
226 : case VERTEX: CollectCornerCoordinates(vVertPos, *static_cast<Vertex*>(elem), domain, true); break;
227 : case EDGE: CollectCornerCoordinates(vVertPos, *static_cast<Edge*>(elem), domain, true); break;
228 : case FACE: CollectCornerCoordinates(vVertPos, *static_cast<Face*>(elem), domain, true); break;
229 : case VOLUME: CollectCornerCoordinates(vVertPos, *static_cast<Volume*>(elem), domain, true); break;
230 0 : default: UG_THROW( "Base Object type not found.");
231 : }
232 :
233 : // forward
234 0 : return DoFPosition<TDomain::dim>(vPos, roid, vVertPos, lfeID);
235 0 : }
236 :
237 :
238 : ////////////////////////////////////////////////////////////////////////////////
239 : // ShapesAtGlobalPosition
240 : ////////////////////////////////////////////////////////////////////////////////
241 :
242 : template <int dim>
243 0 : void ShapesAtGlobalPositionVertex(std::vector<std::vector<number> >& vvShape,
244 : const std::vector<MathVector<dim> >& vGlobPos,
245 : const LFEID& lfeID)
246 : {
247 : // get local position of DoF
248 0 : std::vector<MathVector<0> > vLocPos(vGlobPos.size(), 0.0);
249 :
250 : // evaluate coarse shape fct at fine local point
251 : try{
252 : const LocalShapeFunctionSet<0>& lsfs =
253 : LocalFiniteElementProvider::get<0>(ROID_VERTEX, lfeID);
254 0 : lsfs.shapes(vvShape, vLocPos);
255 : }
256 0 : UG_CATCH_THROW("ShapesAtGlobalPosition: Cannot evalute shapes.")
257 0 : }
258 :
259 : template <int refDim, int dim>
260 0 : void ShapesAtGlobalPositionElem(std::vector<std::vector<number> >& vvShape,
261 : const std::vector<MathVector<dim> >& vGlobPos,
262 : const ReferenceObjectID roid,
263 : const std::vector<MathVector<dim> >& vCornerCoord,
264 : const LFEID& lfeID)
265 : {
266 : // get local position of DoF
267 0 : std::vector<MathVector<refDim> > vLocPos(vGlobPos.size(), 0.0);
268 : try{
269 : DimReferenceMapping<refDim, dim>& map =
270 : ReferenceMappingProvider::get<refDim, dim>(roid, vCornerCoord);
271 0 : map.global_to_local(vLocPos, vGlobPos);
272 : }
273 0 : UG_CATCH_THROW("ShapesAtGlobalPosition: Cannot find elem-local Positions.");
274 :
275 : // evaluate coarse shape fct at fine local point
276 : try{
277 : const LocalShapeFunctionSet<refDim>& lsfs =
278 : LocalFiniteElementProvider::get<refDim>(roid, lfeID);
279 0 : lsfs.shapes(vvShape, vLocPos);
280 : }
281 0 : UG_CATCH_THROW("ShapesAtGlobalPosition: Cannot evalute shapes.")
282 0 : }
283 :
284 : template <int dim>
285 0 : void ShapesAtGlobalPosition(std::vector<std::vector<number> >& vvShape,
286 : const std::vector<MathVector<dim> >& vGlobPos,
287 : const ReferenceObjectID roid,
288 : const std::vector<MathVector<dim> >& vCornerCoord,
289 : const LFEID& lfeID)
290 : {
291 0 : switch(ReferenceElementDimension(roid))
292 : {
293 0 : case VERTEX: return ShapesAtGlobalPositionVertex<dim>(vvShape, vGlobPos, lfeID);
294 0 : case EDGE: return ShapesAtGlobalPositionElem<1,dim>(vvShape, vGlobPos, roid, vCornerCoord, lfeID);
295 0 : case FACE: return ShapesAtGlobalPositionElem<2,dim>(vvShape, vGlobPos, roid, vCornerCoord, lfeID);
296 0 : case VOLUME: return ShapesAtGlobalPositionElem<3,dim>(vvShape, vGlobPos, roid, vCornerCoord, lfeID);
297 0 : default: UG_THROW("Base Object type not found.");
298 : }
299 : }
300 :
301 : template <typename TDomain>
302 0 : void ShapesAtGlobalPosition(std::vector<std::vector<number> >& vvShape,
303 : const std::vector<MathVector<TDomain::dim> >& vGlobPos,
304 : GridObject* elem, const TDomain& domain, const LFEID& lfeID)
305 : {
306 0 : const int baseDim = elem->base_object_id();
307 0 : if(baseDim == VERTEX)
308 0 : return ShapesAtGlobalPositionVertex<TDomain::dim>(vvShape, vGlobPos, lfeID);
309 :
310 : // get the vertices
311 : std::vector<MathVector<TDomain::dim> > vCornerCoord;
312 0 : switch(baseDim)
313 : {
314 : case EDGE: CollectCornerCoordinates(vCornerCoord, *static_cast<Edge*>(elem), domain, true); break;
315 : case FACE: CollectCornerCoordinates(vCornerCoord, *static_cast<Face*>(elem), domain, true); break;
316 : case VOLUME: CollectCornerCoordinates(vCornerCoord, *static_cast<Volume*>(elem), domain, true); break;
317 0 : default: UG_THROW( "Base Object type not found.");
318 : }
319 :
320 : // reference object id
321 0 : const ReferenceObjectID roid = elem->reference_object_id();
322 :
323 : // forward
324 0 : return ShapesAtGlobalPosition<TDomain::dim>(vvShape, vGlobPos, roid, vCornerCoord, lfeID);
325 0 : }
326 :
327 :
328 : ////////////////////////////////////////////////////////////////////////////////
329 : // Extract Positions
330 : ////////////////////////////////////////////////////////////////////////////////
331 :
332 : template <typename TDomain>
333 0 : void ExtractPositionsVertex(ConstSmartPtr<TDomain> domain,
334 : ConstSmartPtr<DoFDistribution> dd,
335 : std::vector<MathVector<TDomain::dim> >& vPos)
336 : {
337 : // get position accessor
338 : const typename TDomain::position_accessor_type& aaPos = domain->position_accessor();
339 :
340 : // iterator
341 : typename DoFDistribution::traits<Vertex>::const_iterator iter, iterEnd;
342 :
343 : // algebra indices vector
344 : std::vector<size_t> ind;
345 :
346 : // get iterators
347 0 : iter = dd->begin<Vertex>(SurfaceView::ALL);
348 0 : iterEnd = dd->end<Vertex>(SurfaceView::ALL);
349 :
350 : // loop all vertices
351 0 : for(;iter != iterEnd; ++iter)
352 : {
353 : // get vertex
354 : Vertex* v = *iter;
355 :
356 : // load indices associated with vertex
357 0 : dd->inner_algebra_indices(v, ind);
358 :
359 : // write position
360 0 : for(size_t i = 0; i < ind.size(); ++i)
361 : {
362 0 : const size_t index = ind[i];
363 : vPos[index] = aaPos[v];
364 : }
365 : }
366 0 : }
367 :
368 : template <typename TDomain, typename TBaseElem>
369 0 : void ExtractPositionsElem(ConstSmartPtr<TDomain> domain,
370 : ConstSmartPtr<DoFDistribution> dd,
371 : std::vector<MathVector<TDomain::dim> >& vPos)
372 : {
373 : typename DoFDistribution::traits<TBaseElem>::const_iterator iter, iterEnd;
374 :
375 : // vector for positions
376 : std::vector<MathVector<TDomain::dim> > vElemPos;
377 :
378 : // algebra indices vector
379 : std::vector<DoFIndex> ind;
380 :
381 : // loop all subsets
382 0 : for(int si = 0; si < dd->num_subsets(); ++si)
383 : {
384 : // get iterators
385 0 : iter = dd->begin<TBaseElem>(si, SurfaceView::ALL);
386 0 : iterEnd = dd->end<TBaseElem>(si, SurfaceView::ALL);
387 :
388 : // loop all elements
389 0 : for(;iter != iterEnd; ++iter)
390 : {
391 : // get vertex
392 : TBaseElem* elem = *iter;
393 :
394 : // loop all functions
395 0 : for(size_t fct = 0; fct < dd->num_fct(); ++fct)
396 : {
397 : // skip non-used function
398 0 : if(!dd->is_def_in_subset(fct,si)) continue;
399 :
400 : // load indices associated with element function
401 0 : dd->inner_dof_indices(elem, fct, ind);
402 :
403 : // load positions associated with element and function
404 0 : InnerDoFPosition(vElemPos, elem, *(const_cast<TDomain*>(domain.get())),
405 : dd->local_finite_element_id(fct));
406 :
407 : // check correct size
408 : UG_ASSERT(ind.size() == vElemPos.size(), "Num MultiIndex ("<<ind.size()
409 : <<") and Num Position ("<<vElemPos.size()<<") must match."
410 : "GeomObject dim="<<geometry_traits<TBaseElem>::BASE_OBJECT_ID);
411 :
412 : // write position
413 0 : for(size_t sh = 0; sh < ind.size(); ++sh)
414 : {
415 0 : const size_t index = ind[sh][0];
416 : vPos[index] = vElemPos[sh];
417 : }
418 : }
419 : }
420 : }
421 0 : }
422 :
423 : template <typename TDomain>
424 0 : void ExtractPositions(ConstSmartPtr<TDomain> domain,
425 : ConstSmartPtr<DoFDistribution> dd,
426 : std::vector<MathVector<TDomain::dim> >& vPos)
427 : {
428 : // number of total dofs
429 0 : int nr = dd->num_indices();
430 :
431 : // resize positions
432 0 : vPos.resize(nr);
433 :
434 : // extract for all element types
435 0 : if(dd->max_dofs(VERTEX)) ExtractPositionsVertex<TDomain>(domain, dd, vPos);
436 0 : if(dd->max_dofs(EDGE)) ExtractPositionsElem<TDomain, Edge>(domain, dd, vPos);
437 0 : if(dd->max_dofs(FACE)) ExtractPositionsElem<TDomain, Face>(domain, dd, vPos);
438 0 : if(dd->max_dofs(VOLUME)) ExtractPositionsElem<TDomain, Volume>(domain, dd, vPos);
439 0 : }
440 :
441 :
442 : template <typename TDomain, typename TBaseElem>
443 0 : void ExtractAlgebraIndices2(ConstSmartPtr<TDomain> domain,
444 : ConstSmartPtr<DoFDistribution> dd,
445 : std::vector<size_t>& fctIndex)
446 : {
447 : typename DoFDistribution::traits<TBaseElem>::const_iterator iter, iterEnd;
448 :
449 : // algebra indices vector
450 : std::vector<size_t> ind;
451 :
452 : // loop all subsets
453 0 : for(int si = 0; si < dd->num_subsets(); ++si)
454 : {
455 : // get iterators
456 0 : iter = dd->begin<TBaseElem>(si, SurfaceView::ALL);
457 0 : iterEnd = dd->end<TBaseElem>(si, SurfaceView::ALL);
458 :
459 : // loop all elements
460 0 : for(;iter != iterEnd; ++iter)
461 : {
462 : // get element
463 : TBaseElem* elem = *iter;
464 :
465 : // loop all functions
466 0 : for(size_t fct = 0; fct < dd->num_fct(); ++fct)
467 : {
468 : // skip non-used function
469 0 : if(!dd->is_def_in_subset(fct,si)) continue;
470 :
471 : // load indices associated with element function
472 0 : dd->inner_algebra_indices_for_fct(elem, ind, true, fct);
473 :
474 0 : for(size_t sh = 0; sh < ind.size(); ++sh)
475 : {
476 0 : const size_t index = ind[sh];
477 0 : if(index < fctIndex.size())
478 0 : fctIndex[index] = fct;
479 : }
480 : }
481 : }
482 : }
483 0 : }
484 :
485 :
486 : template <typename TDomain>
487 0 : void ExtractAlgebraIndices(ConstSmartPtr<TDomain> domain,
488 : ConstSmartPtr<DoFDistribution> dd,
489 : std::vector<size_t> &fctIndex)
490 : {
491 : // number of total dofs
492 0 : int nr = dd->num_indices();
493 :
494 : // resize positions
495 0 : fctIndex.resize(nr);
496 :
497 : // extract for all element types
498 0 : if(dd->max_dofs(VERTEX)) ExtractAlgebraIndices2<TDomain, Vertex>(domain, dd, fctIndex);
499 0 : if(dd->max_dofs(EDGE)) ExtractAlgebraIndices2<TDomain, Edge>(domain, dd, fctIndex);
500 0 : if(dd->max_dofs(FACE)) ExtractAlgebraIndices2<TDomain, Face>(domain, dd, fctIndex);
501 0 : if(dd->max_dofs(VOLUME)) ExtractAlgebraIndices2<TDomain, Volume>(domain, dd, fctIndex);
502 0 : }
503 :
504 :
505 : ////////////////////////////////////////////////////////////////////////////////
506 : // Extract (Positions, Index) Pairs
507 : ////////////////////////////////////////////////////////////////////////////////
508 :
509 : template <typename TDomain>
510 0 : void ExtractPositionsVertex(ConstSmartPtr<TDomain> domain,
511 : ConstSmartPtr<DoFDistribution> dd,
512 : std::vector<std::pair<MathVector<TDomain::dim>, size_t> >& vPosPair)
513 : {
514 : // get position accessor
515 : const typename TDomain::position_accessor_type& aaPos = domain->position_accessor();
516 :
517 : // resize positions
518 0 : vPosPair.resize(dd->num_indices());
519 :
520 : typedef DoFDistribution::traits<Vertex>::const_iterator const_iterator;
521 :
522 : // loop all vertices
523 : const_iterator iter = dd->begin<Vertex>(SurfaceView::ALL);
524 : const_iterator iterEnd = dd->end<Vertex>(SurfaceView::ALL);
525 :
526 : // algebra indices vector
527 : std::vector<size_t> ind;
528 :
529 0 : for(;iter != iterEnd; ++iter)
530 : {
531 : // get vertex
532 : Vertex* v = *iter;
533 :
534 : // load indices associated with vertex
535 0 : dd->inner_algebra_indices(v, ind);
536 :
537 : // write position and index
538 0 : for(size_t i = 0; i < ind.size(); ++i)
539 : {
540 0 : const size_t index = ind[i];
541 : vPosPair[index].first = aaPos[v];
542 0 : vPosPair[index].second = index;
543 : }
544 : }
545 0 : }
546 :
547 :
548 : template <typename TDomain, typename TBaseElem>
549 0 : void ExtractPositionsElem(ConstSmartPtr<TDomain> domain,
550 : ConstSmartPtr<DoFDistribution> dd,
551 : std::vector<std::pair<MathVector<TDomain::dim>, size_t> >& vPosPair)
552 : {
553 : typename DoFDistribution::traits<TBaseElem>::const_iterator iter, iterEnd;
554 :
555 : // vector for positions
556 : std::vector<MathVector<TDomain::dim> > vElemPos;
557 :
558 : // algebra indices vector
559 : std::vector<DoFIndex> ind;
560 :
561 : // loop all subsets
562 0 : for(int si = 0; si < dd->num_subsets(); ++si)
563 : {
564 : // get iterators
565 0 : iter = dd->begin<TBaseElem>(si, SurfaceView::ALL);
566 0 : iterEnd = dd->end<TBaseElem>(si, SurfaceView::ALL);
567 :
568 : // loop all elements
569 0 : for(;iter != iterEnd; ++iter)
570 : {
571 : // get vertex
572 : TBaseElem* elem = *iter;
573 :
574 : // loop all functions
575 0 : for(size_t fct = 0; fct < dd->num_fct(); ++fct)
576 : {
577 : // skip non-used function
578 0 : if(!dd->is_def_in_subset(fct,si)) continue;
579 :
580 : // load indices associated with element function
581 0 : dd->inner_dof_indices(elem, fct, ind);
582 :
583 : // load positions associated with element and function
584 0 : InnerDoFPosition(vElemPos, elem, *(const_cast<TDomain*>(domain.get())),
585 : dd->local_finite_element_id(fct));
586 :
587 : // check correct size
588 : UG_ASSERT(ind.size() == vElemPos.size(), "Num MultiIndex ("<<ind.size()
589 : <<") and Num Position ("<<vElemPos.size()<<") must match."
590 : "GeomObject dim="<<geometry_traits<TBaseElem>::BASE_OBJECT_ID);
591 :
592 : // write position
593 0 : for(size_t sh = 0; sh < ind.size(); ++sh)
594 : {
595 0 : const size_t index = ind[sh][0];
596 : vPosPair[index].first = vElemPos[sh];
597 0 : vPosPair[index].second = index;
598 :
599 : }
600 : }
601 : }
602 : }
603 0 : }
604 :
605 : template <typename TDomain>
606 0 : void ExtractPositions(ConstSmartPtr<TDomain> domain,
607 : ConstSmartPtr<DoFDistribution> dd,
608 : std::vector<std::pair<MathVector<TDomain::dim>, size_t> >& vPosPair)
609 : {
610 : // number of total dofs
611 0 : int nr = dd->num_indices();
612 :
613 : // resize positions
614 0 : vPosPair.resize(nr);
615 :
616 : // extract for all element types
617 0 : if(dd->max_dofs(VERTEX)) ExtractPositionsVertex<TDomain>(domain, dd, vPosPair);
618 0 : if(dd->max_dofs(EDGE)) ExtractPositionsElem<TDomain, Edge>(domain, dd, vPosPair);
619 0 : if(dd->max_dofs(FACE)) ExtractPositionsElem<TDomain, Face>(domain, dd, vPosPair);
620 0 : if(dd->max_dofs(VOLUME)) ExtractPositionsElem<TDomain, Volume>(domain, dd, vPosPair);
621 0 : }
622 :
623 : ////////////////////////////////////////////////////////////////////////////////
624 : // Extract (Positions, Index) Pairs for a single component
625 : ////////////////////////////////////////////////////////////////////////////////
626 :
627 :
628 : template <typename TDomain, typename TBaseElem>
629 0 : void ExtractPositionsElem(ConstSmartPtr<TDomain> domain,
630 : ConstSmartPtr<DoFDistribution> dd,
631 : const size_t fct,
632 : std::vector<std::pair<MathVector<TDomain::dim>, size_t> >& vPosPair)
633 : {
634 : typename DoFDistribution::traits<TBaseElem>::const_iterator iter, iterEnd;
635 :
636 : // vector for positions
637 : std::vector<MathVector<TDomain::dim> > vElemPos;
638 :
639 : // algebra indices vector
640 : std::vector<DoFIndex> ind;
641 :
642 : // a pair
643 : std::pair<MathVector<TDomain::dim>, size_t> pair;
644 :
645 : // loop all subsets
646 0 : for(int si = 0; si < dd->num_subsets(); ++si)
647 : {
648 : // get iterators
649 0 : iter = dd->begin<TBaseElem>(si, SurfaceView::ALL);
650 0 : iterEnd = dd->end<TBaseElem>(si, SurfaceView::ALL);
651 :
652 : // skip non-used function
653 0 : if(!dd->is_def_in_subset(fct,si)) continue;
654 :
655 : // loop all elements
656 0 : for(;iter != iterEnd; ++iter)
657 : {
658 : // get vertex
659 : TBaseElem* elem = *iter;
660 :
661 : // load indices associated with element function
662 0 : dd->inner_dof_indices(elem, fct, ind);
663 :
664 : // load positions associated with element and function
665 0 : InnerDoFPosition(vElemPos, elem, *(const_cast<TDomain*>(domain.get())),
666 : dd->local_finite_element_id(fct));
667 :
668 : // check correct size
669 : UG_ASSERT(ind.size() == vElemPos.size(), "Num MultiIndex ("<<ind.size()
670 : <<") and Num Position ("<<vElemPos.size()<<") must match."
671 : "GeomObject dim="<<geometry_traits<TBaseElem>::BASE_OBJECT_ID);
672 :
673 : // write position
674 0 : for(size_t sh = 0; sh < ind.size(); ++sh)
675 : {
676 0 : const size_t index = ind[sh][0];
677 :
678 : pair.first = vElemPos[sh];
679 0 : pair.second = index;
680 :
681 0 : vPosPair.push_back(pair);
682 : }
683 : }
684 : }
685 0 : }
686 :
687 : template <typename TDomain>
688 0 : void ExtractPositions(ConstSmartPtr<TDomain> domain,
689 : ConstSmartPtr<DoFDistribution> dd,
690 : const size_t fct,
691 : std::vector<std::pair<MathVector<TDomain::dim>, size_t> >& vPosPair)
692 : {
693 : // resize positions
694 : vPosPair.clear();
695 :
696 : // extract for all element types
697 0 : if(dd->max_dofs(VERTEX)) ExtractPositionsElem<TDomain, Vertex>(domain, dd, fct, vPosPair);
698 0 : if(dd->max_dofs(EDGE)) ExtractPositionsElem<TDomain, Edge>(domain, dd, fct, vPosPair);
699 0 : if(dd->max_dofs(FACE)) ExtractPositionsElem<TDomain, Face>(domain, dd, fct, vPosPair);
700 0 : if(dd->max_dofs(VOLUME)) ExtractPositionsElem<TDomain, Volume>(domain, dd, fct, vPosPair);
701 0 : }
702 :
703 : ////////////////////////////////////////////////////////////////////////////////
704 : // Checks correct DoF Positions
705 : ////////////////////////////////////////////////////////////////////////////////
706 :
707 : template <typename TDomain, typename TBaseElem>
708 0 : bool CheckDoFElem(ConstSmartPtr<TDomain> domain,
709 : ConstSmartPtr<DoFDistribution> dd,
710 : std::vector<MathVector<TDomain::dim> >& vPos)
711 : {
712 : typename DoFDistribution::traits<TBaseElem>::const_iterator iter, iterEnd;
713 :
714 : bool bRes = true;
715 :
716 : // vector for positions
717 : std::vector<MathVector<TDomain::dim> > vElemPos;
718 :
719 : // algebra indices vector
720 : std::vector<DoFIndex> ind;
721 :
722 : // loop all subsets
723 0 : for(int si = 0; si < dd->num_subsets(); ++si)
724 : {
725 : // get iterators
726 0 : iter = dd->begin<TBaseElem>(si, SurfaceView::ALL);
727 0 : iterEnd = dd->end<TBaseElem>(si, SurfaceView::ALL);
728 :
729 : // loop all elements
730 0 : for(;iter != iterEnd; ++iter)
731 : {
732 : // get element
733 : TBaseElem* elem = *iter;
734 :
735 : // loop all functions
736 0 : for(size_t fct = 0; fct < dd->num_fct(); ++fct)
737 : {
738 : // skip non-used function
739 0 : if(!dd->is_def_in_subset(fct,si)) continue;
740 :
741 : // load indices associated with element function
742 0 : dd->inner_dof_indices(elem, fct, ind);
743 :
744 : // load positions associated with element and function
745 0 : InnerDoFPosition(vElemPos, elem, *(const_cast<TDomain*>(domain.get())),
746 : dd->local_finite_element_id(fct));
747 :
748 : bool bWrite = false;
749 0 : for(size_t sh = 0; sh < ind.size(); ++sh)
750 : {
751 0 : size_t index = ind[sh][0];
752 :
753 0 : if(vPos[index] != MathVector<TDomain::dim>(-1)){
754 0 : if(VecDistance(vPos[index], vElemPos[sh]) < 1e-10) continue;
755 :
756 0 : if(!bWrite)
757 : UG_LOG(" **** inner_multi_index (start) ******\n")
758 : bWrite = true;
759 : bRes = false;
760 0 : UG_LOG("CheckDoFPositions "<<sh<<": inner_dof_indices: index: "
761 : <<index<<" at "<<vElemPos[sh]<<", but previously: "
762 : <<vPos[index]<<"\n");
763 : }
764 :
765 : vPos[index] = vElemPos[sh];
766 :
767 : }
768 :
769 0 : if(bWrite){
770 : std::vector<MathVector<TDomain::dim> > vVertPos;
771 : CollectCornerCoordinates(vVertPos, elem, *domain, true);
772 0 : UG_LOG("From Elem "<<elem<<" ("<<elem->reference_object_id()<<"):\n")
773 0 : for(size_t i = 0; i < vVertPos.size(); ++i)
774 0 : UG_LOG("with corner "<<i<<": "<<vVertPos[i]<<"\n");
775 : UG_LOG(" **** inner_multi_index (end) ******\n")
776 0 : }
777 :
778 : /////////////////////////
779 : // load indices associated with element function
780 0 : dd->dof_indices(elem, fct, ind);
781 :
782 : // load positions associated with element and function
783 0 : DoFPosition(vElemPos, elem, *(const_cast<TDomain*>(domain.get())),
784 : dd->local_finite_element_id(fct));
785 :
786 : // write position
787 : bWrite = false;
788 0 : for(size_t sh = 0; sh < ind.size(); ++sh)
789 : {
790 0 : size_t index = ind[sh][0];
791 :
792 0 : if(vPos[index] != MathVector<TDomain::dim>(-1)){
793 0 : if(VecDistance(vPos[index], vElemPos[sh]) < 1e-10) continue;
794 :
795 0 : if(!bWrite)
796 : UG_LOG(" **** multi_index (start) ******\n")
797 : bWrite = true;
798 : bRes = false;
799 0 : UG_LOG("CheckDoFPositions "<<sh<<": dof_indices: index: "
800 : <<index<<" at "<<vElemPos[sh]<<", but previously: "
801 : <<vPos[index]<<"\n");
802 : }
803 :
804 : vPos[index] = vElemPos[sh];
805 : }
806 :
807 0 : if(bWrite){
808 : std::vector<MathVector<TDomain::dim> > vVertPos;
809 : CollectCornerCoordinates(vVertPos, elem, *domain, true);
810 0 : UG_LOG("From Elem "<<elem<<" ("<<elem->reference_object_id()<<"):\n")
811 0 : for(size_t i = 0; i < vVertPos.size(); ++i)
812 0 : UG_LOG("with corner "<<i<<": "<<vVertPos[i]<<"\n");
813 : UG_LOG(" **** multi_index (end) ******\n")
814 0 : }
815 : }
816 : }
817 : }
818 :
819 0 : return bRes;
820 0 : }
821 :
822 : template <typename TDomain>
823 0 : bool CheckDoFPositions(ConstSmartPtr<TDomain> domain,
824 : ConstSmartPtr<DoFDistribution> dd)
825 : {
826 : // number of total dofs
827 0 : int nr = dd->num_indices();
828 0 : std::vector<MathVector<TDomain::dim> > vPos(nr, -1);
829 :
830 : // extract for all element types
831 : bool bRes = true;
832 0 : if(dd->max_dofs(VERTEX)) bRes &= CheckDoFElem<TDomain, Vertex>(domain, dd, vPos);
833 0 : if(dd->max_dofs(EDGE)) bRes &= CheckDoFElem<TDomain, Edge>(domain, dd, vPos);
834 0 : if(dd->max_dofs(FACE)) bRes &= CheckDoFElem<TDomain, Face>(domain, dd, vPos);
835 0 : if(dd->max_dofs(VOLUME)) bRes &= CheckDoFElem<TDomain, Volume>(domain, dd, vPos);
836 0 : return bRes;
837 0 : }
838 :
839 : #ifdef UG_DIM_1
840 : template bool InnerDoFPosition<Domain1d>(std::vector<MathVector<1> >& vPos, GridObject* elem, const Domain1d& domain, const LFEID& lfeID);
841 : template bool InnerDoFPosition<1>(std::vector<MathVector<1> >& vPos, const ReferenceObjectID roid,const std::vector<MathVector<1> >& vCornerCoord, const LFEID& lfeID);
842 : template bool DoFPosition<Domain1d>(std::vector<MathVector<1> >& vPos, GridObject* elem, const Domain1d& domain, const LFEID& lfeID);
843 : template bool DoFPosition<1>(std::vector<MathVector<1> >& vPos, const ReferenceObjectID roid,const std::vector<MathVector<1> >& vCornerCoord, const LFEID& lfeID);
844 : template void ExtractPositions(ConstSmartPtr<Domain1d> domain, ConstSmartPtr<DoFDistribution> dd, std::vector<MathVector<Domain1d::dim> >& vPos);
845 : template void ExtractPositions(ConstSmartPtr<Domain1d> domain, ConstSmartPtr<DoFDistribution> dd, std::vector<std::pair<MathVector<Domain1d::dim>, size_t> >& vPos);
846 : template void ExtractPositions(ConstSmartPtr<Domain1d> domain, ConstSmartPtr<DoFDistribution> dd, const size_t fct, std::vector<std::pair<MathVector<Domain1d::dim>, size_t> >& vPos);
847 : template bool CheckDoFPositions(ConstSmartPtr<Domain1d> domain, ConstSmartPtr<DoFDistribution> dd);
848 : template void ExtractAlgebraIndices<Domain1d>(ConstSmartPtr<Domain1d> domain, ConstSmartPtr<DoFDistribution> dd, std::vector<size_t> &fctIndex);
849 : template void ShapesAtGlobalPosition<1>(std::vector<std::vector<number> >& vvShape, const std::vector<MathVector<1> >& vGlobPos, const ReferenceObjectID roid,const std::vector<MathVector<1> >& vCornerCoord, const LFEID& lfeID);
850 : template void ShapesAtGlobalPosition<Domain1d>(std::vector<std::vector<number> >& vvShape, const std::vector<MathVector<1> >& vGlobPos,GridObject* elem, const Domain1d& domain, const LFEID& lfeID);
851 : #endif
852 : #ifdef UG_DIM_2
853 : template bool InnerDoFPosition<Domain2d>(std::vector<MathVector<2> >& vPos, GridObject* elem, const Domain2d& domain, const LFEID& lfeID);
854 : template bool InnerDoFPosition<2>(std::vector<MathVector<2> >& vPos, const ReferenceObjectID roid,const std::vector<MathVector<2> >& vCornerCoord, const LFEID& lfeID);
855 : template bool DoFPosition<Domain2d>(std::vector<MathVector<2> >& vPos, GridObject* elem, const Domain2d& domain, const LFEID& lfeID);
856 : template bool DoFPosition<2>(std::vector<MathVector<2> >& vPos, const ReferenceObjectID roid,const std::vector<MathVector<2> >& vCornerCoord, const LFEID& lfeID);
857 : template void ExtractPositions(ConstSmartPtr<Domain2d> domain, ConstSmartPtr<DoFDistribution> dd, std::vector<MathVector<Domain2d::dim> >& vPos);
858 : template void ExtractPositions(ConstSmartPtr<Domain2d> domain, ConstSmartPtr<DoFDistribution> dd, std::vector<std::pair<MathVector<Domain2d::dim>, size_t> >& vPos);
859 : template void ExtractPositions(ConstSmartPtr<Domain2d> domain, ConstSmartPtr<DoFDistribution> dd, const size_t fct, std::vector<std::pair<MathVector<Domain2d::dim>, size_t> >& vPos);
860 : template bool CheckDoFPositions(ConstSmartPtr<Domain2d> domain, ConstSmartPtr<DoFDistribution> dd);
861 : template void ExtractAlgebraIndices<Domain2d>(ConstSmartPtr<Domain2d> domain, ConstSmartPtr<DoFDistribution> dd, std::vector<size_t> &fctIndex);
862 : template void ShapesAtGlobalPosition<2>(std::vector<std::vector<number> >& vvShape, const std::vector<MathVector<2> >& vGlobPos, const ReferenceObjectID roid,const std::vector<MathVector<2> >& vCornerCoord, const LFEID& lfeID);
863 : template void ShapesAtGlobalPosition<Domain2d>(std::vector<std::vector<number> >& vvShape, const std::vector<MathVector<2> >& vGlobPos,GridObject* elem, const Domain2d& domain, const LFEID& lfeID);
864 : #endif
865 : #ifdef UG_DIM_3
866 : template bool InnerDoFPosition<Domain3d>(std::vector<MathVector<3> >& vPos, GridObject* elem, const Domain3d& domain, const LFEID& lfeID);
867 : template bool InnerDoFPosition<3>(std::vector<MathVector<3> >& vPos, const ReferenceObjectID roid,const std::vector<MathVector<3> >& vCornerCoord, const LFEID& lfeID);
868 : template bool DoFPosition<Domain3d>(std::vector<MathVector<3> >& vPos, GridObject* elem, const Domain3d& domain, const LFEID& lfeID);
869 : template bool DoFPosition<3>(std::vector<MathVector<3> >& vPos, const ReferenceObjectID roid,const std::vector<MathVector<3> >& vCornerCoord, const LFEID& lfeID);
870 : template void ExtractPositions(ConstSmartPtr<Domain3d> domain, ConstSmartPtr<DoFDistribution> dd, std::vector<MathVector<Domain3d::dim> >& vPos);
871 : template void ExtractPositions(ConstSmartPtr<Domain3d> domain, ConstSmartPtr<DoFDistribution> dd, std::vector<std::pair<MathVector<Domain3d::dim>, size_t> >& vPos);
872 : template void ExtractPositions(ConstSmartPtr<Domain3d> domain, ConstSmartPtr<DoFDistribution> dd, const size_t fct, std::vector<std::pair<MathVector<Domain3d::dim>, size_t> >& vPos);
873 : template bool CheckDoFPositions(ConstSmartPtr<Domain3d> domain, ConstSmartPtr<DoFDistribution> dd);
874 : template void ExtractAlgebraIndices<Domain3d>(ConstSmartPtr<Domain3d> domain, ConstSmartPtr<DoFDistribution> dd, std::vector<size_t> &fctIndex);
875 : template void ShapesAtGlobalPosition<3>(std::vector<std::vector<number> >& vvShape, const std::vector<MathVector<3> >& vGlobPos, const ReferenceObjectID roid,const std::vector<MathVector<3> >& vCornerCoord, const LFEID& lfeID);
876 : template void ShapesAtGlobalPosition<Domain3d>(std::vector<std::vector<number> >& vvShape, const std::vector<MathVector<3> >& vGlobPos,GridObject* elem, const Domain3d& domain, const LFEID& lfeID);
877 : #endif
878 :
879 : } // end namespace ug
|
