Line data Source code
1 : /*
2 : * Copyright (c) 2009-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__LIB_GRID__MISC_UTIL__IMPL__
34 : #define __H__LIB_GRID__MISC_UTIL__IMPL__
35 :
36 : #include "misc_util.h"
37 : #include "lib_grid/grid/grid_util.h"
38 : #include "vertex_util.h"
39 : #include "edge_util.h"
40 : #include "face_util.h"
41 : #include "volume_util.h"
42 :
43 : #include <limits>
44 :
45 : namespace ug
46 : {
47 :
48 : ////////////////////////////////////////////////////////////////////////////////////////////
49 : template <class TElem>
50 : UG_API
51 : typename TElem::side*
52 : GetSharedSide(Grid& grid, TElem* e1, TElem* e2)
53 : {
54 : if(!TElem::HAS_SIDES)
55 : return NULL;
56 :
57 : typedef typename TElem::side side_t;
58 : typename Grid::traits<side_t>::secure_container sides1;
59 : typename Grid::traits<side_t>::secure_container sides2;
60 :
61 : grid.associated_elements(sides1, e1);
62 : grid.associated_elements(sides2, e2);
63 :
64 : for(size_t i1 = 0; i1 < sides1.size(); ++i1){
65 : for(size_t i2 = 0; i2 < sides2.size(); ++i2){
66 : if(sides1[i1] == sides2[i2])
67 : return sides1[i1];
68 : }
69 : }
70 :
71 : return NULL;
72 : }
73 :
74 :
75 : ////////////////////////////////////////////////////////////////////////
76 : // CalculateGridObjectCenter
77 : template<class TAAPosVRT>
78 : UG_API
79 : inline
80 : typename TAAPosVRT::ValueType
81 : CalculateGridObjectCenter(const GridObject* o, TAAPosVRT& aaPosVRT)
82 : {
83 : switch(o->base_object_id()){
84 : case VERTEX: return CalculateCenter(static_cast<const Vertex*>(o), aaPosVRT);
85 : case EDGE: return CalculateCenter(static_cast<const Edge*>(o), aaPosVRT);
86 : case FACE: return CalculateCenter(static_cast<const Face*>(o), aaPosVRT);
87 : case VOLUME: return CalculateCenter(static_cast<const Volume*>(o), aaPosVRT);
88 : default:
89 : UG_THROW("Unknown geometric-object type.");
90 : }
91 : }
92 :
93 : template<class TAAPosVRT, class TAAWeightVRT>
94 : UG_API
95 : inline
96 : typename TAAPosVRT::ValueType
97 : CalculateGridObjectCenter(const GridObject* o, TAAPosVRT& aaPosVRT,
98 : TAAWeightVRT& aaWeight)
99 : {
100 : switch(o->base_object_id()){
101 : case VERTEX:
102 : return CalculateCenter(static_cast<const Vertex*>(o), aaPosVRT, aaWeight);
103 : case EDGE:
104 : return CalculateCenter(static_cast<const Edge*>(o), aaPosVRT, aaWeight);
105 : case FACE:
106 : return CalculateCenter(static_cast<const Face*>(o), aaPosVRT, aaWeight);
107 : case VOLUME:
108 : return CalculateCenter(static_cast<const Volume*>(o), aaPosVRT, aaWeight);
109 : default:
110 : UG_THROW("Unknown geometric-object type.");
111 : }
112 : }
113 :
114 :
115 : ////////////////////////////////////////////////////////////////////////
116 : // CalculateCenter
117 : template <class TIterator, class TAAPosVRT>
118 : typename TAAPosVRT::ValueType
119 0 : CalculateCenter(TIterator begin, TIterator end, TAAPosVRT& aaPos)
120 : {
121 : int counter = 0;
122 : typename TAAPosVRT::ValueType center;
123 : VecSet(center, 0);
124 0 : for(TIterator iter = begin; iter != end; ++iter, ++counter)
125 0 : VecAdd(center, center, CalculateCenter(*iter, aaPos));
126 :
127 0 : if(counter > 0)
128 0 : VecScale(center, center, 1./(number)counter);
129 :
130 0 : return center;
131 : }
132 :
133 : ////////////////////////////////////////////////////////////////////////
134 : // FindClosestByCoordinate
135 : template<class TElem, class TVertexPositionAttachmentAccessor>
136 0 : TElem* FindClosestByCoordinate(const typename TVertexPositionAttachmentAccessor::ValueType& coord,
137 : typename geometry_traits<TElem>::iterator iterBegin,
138 : typename geometry_traits<TElem>::iterator iterEnd,
139 : TVertexPositionAttachmentAccessor& aaPosVRT)
140 : {
141 0 : if(iterBegin == iterEnd)
142 : return NULL;
143 :
144 : typename geometry_traits<TElem>::iterator iter = iterBegin;
145 : TElem* bestElem = *iter;
146 0 : number bestDistSq = VecDistanceSq(coord, CalculateCenter(bestElem, aaPosVRT));
147 : iter++;
148 :
149 0 : while(iter != iterEnd)
150 : {
151 0 : number distSq = VecDistanceSq(coord, CalculateCenter(*iter, aaPosVRT));
152 0 : if(distSq < bestDistSq)
153 : {
154 : bestDistSq = distSq;
155 : bestElem = *iter;
156 : }
157 :
158 : ++iter;
159 : }
160 :
161 : return bestElem;
162 : }
163 :
164 : ////////////////////////////////////////////////////////////////////////
165 : template<class vector_t, class TIterator, class TAAPos>
166 : void CalculateBoundingBox(vector_t& vMinOut, vector_t& vMaxOut,
167 : TIterator begin, TIterator end,
168 : TAAPos& aaPos)
169 : {
170 : if(begin == end){
171 : VecSet(vMinOut, 0);
172 : VecSet(vMaxOut, 0);
173 : return;
174 : }
175 :
176 : vMinOut = vMaxOut = aaPos[GetVertex(*begin, 0)];
177 :
178 : const int dim = vector_t::Size;
179 :
180 : // iterate over all elements and find min and max values
181 : vector_t tmin, tmax;
182 : for(TIterator iter = begin; iter != end; ++iter){
183 : typename TIterator::value_type elem = *iter;
184 :
185 : for(size_t i_vrt = 0; i_vrt < NumVertices(elem); ++i_vrt){
186 : for(int i = 0; i < dim; ++i){
187 : vector_t& v = aaPos[GetVertex(elem, i_vrt)];
188 : if(v[i] < vMinOut[i])
189 : vMinOut[i] = v[i];
190 : else if(v[i] > vMaxOut[i])
191 : vMaxOut[i] = v[i];
192 : }
193 : }
194 : }
195 : }
196 :
197 : ////////////////////////////////////////////////////////////////////////
198 : // NumSharedVertices
199 : template <class TElemPtr1, class TElemPtr2>
200 : size_t NumSharedVertices(Grid& grid, TElemPtr1 elem1, TElemPtr2 elem2)
201 : {
202 : grid.begin_marking();
203 : // first mark all vertices of elem1
204 : for(size_t i = 0; i < elem1->num_vertices(); ++i)
205 : grid.mark(elem1->vertex(i));
206 :
207 : // now count how many of vertex 2 are marked.
208 : size_t counter = 0;
209 : for(size_t i = 0; i < elem2->num_vertices(); ++i){
210 : if(grid.is_marked(elem2->vertex(i)))
211 : ++counter;
212 : }
213 :
214 : grid.end_marking();
215 :
216 : return counter;
217 : }
218 :
219 : ////////////////////////////////////////////////////////////////////////
220 : // EraseElements
221 : template <class TElem>
222 0 : void EraseElements(Grid& grid, typename geometry_traits<TElem>::iterator iterBegin,
223 : typename geometry_traits<TElem>::iterator iterEnd)
224 : {
225 : // be careful to not invalidate the iterators.
226 0 : while(iterBegin != iterEnd)
227 : {
228 : TElem* e = *iterBegin;
229 : iterBegin++;
230 0 : grid.erase(e);
231 : }
232 0 : }
233 :
234 : ////////////////////////////////////////////////////////////////////////
235 : // AssignIndices
236 : template <class TElem>
237 : void AssignIndices(typename geometry_traits<TElem>::iterator iterBegin,
238 : typename geometry_traits<TElem>::iterator iterEnd,
239 : Grid::AttachmentAccessor<TElem, AInt>& aaInt)
240 : {
241 : int index = 0;
242 0 : while(iterBegin != iterEnd)
243 : {
244 0 : aaInt[*iterBegin] = index++;
245 : iterBegin++;
246 : }
247 : }
248 :
249 : ////////////////////////////////////////////////////////////////////////
250 : // ElementDiameter
251 : template <class TElem, class TAAPos>
252 0 : number ElementDiameterSq(Grid& grid,
253 : TAAPos& aaPos,
254 : TElem* elem)
255 : {
256 : PointerConstArray<Vertex*> vVert;
257 : grid.associated_elements(vVert, elem);
258 :
259 0 : number max = 0.0;
260 0 : for(size_t i = 0; i < vVert.size(); ++i)
261 0 : for(size_t j = i+1; j < vVert.size(); ++j)
262 0 : max = std::max(max, VecDistanceSq(aaPos[vVert[i]], aaPos[vVert[j]]));
263 :
264 0 : return max;
265 : }
266 :
267 : template <class TAAPos>
268 : number ElementDiameterSq(Grid& grid,
269 : TAAPos& aaPos,
270 : GridObject* elem)
271 : {
272 : switch(elem->base_object_id()){
273 : case VERTEX: return ElementDiameterSq(grid, aaPos, static_cast<Vertex*>(elem));
274 : case EDGE: return ElementDiameterSq(grid, aaPos, static_cast<Edge*>(elem));
275 : case FACE: return ElementDiameterSq(grid, aaPos, static_cast<Face*>(elem));
276 : case VOLUME: return ElementDiameterSq(grid, aaPos, static_cast<Volume*>(elem));
277 : default: UG_THROW("ElementDiameterSq: Element type not found.")
278 : }
279 : }
280 :
281 : template <class TElem, class TAAPos>
282 : number ElementDiameter(Grid& grid,
283 : TAAPos& aaPos,
284 : TElem* elem)
285 : {
286 0 : return std::sqrt(ElementDiameterSq(grid, aaPos, elem));
287 : }
288 :
289 : template <class TAAPos, class TIterator>
290 0 : number MaxElementDiameter(Grid& grid, TAAPos& aaPos,
291 : TIterator iter, TIterator iterEnd)
292 : {
293 0 : number max = 0.0;
294 0 : for(; iter != iterEnd; ++iter)
295 0 : max = std::max(max, ElementDiameterSq(grid, aaPos, *iter));
296 :
297 : #ifdef UG_PARALLEL
298 : // share value between all procs
299 : pcl::ProcessCommunicator com;
300 : max = com.allreduce(max, PCL_RO_MAX);
301 : #endif
302 :
303 0 : return std::sqrt(max);
304 : }
305 :
306 : template <class TAAPos, class TIterator>
307 0 : number MinElementDiameter(Grid& grid, TAAPos& aaPos,
308 : TIterator iter, TIterator iterEnd)
309 : {
310 0 : number min = std::numeric_limits<number>::max();
311 0 : for(; iter != iterEnd; ++iter)
312 0 : min = std::min(min, ElementDiameterSq(grid, aaPos, *iter));
313 :
314 : #ifdef UG_PARALLEL
315 : // share value between all procs
316 : pcl::ProcessCommunicator com;
317 : min = com.allreduce(min, PCL_RO_MIN);
318 : #endif
319 :
320 0 : return std::sqrt(min);
321 : }
322 :
323 :
324 : template <class TElem1, class TElem2, class TAAPos>
325 : typename TAAPos::ValueType
326 : GetDirection (TElem1* e1, TElem2* e2, const TAAPos& aaPos)
327 : {
328 : typedef typename TAAPos::ValueType vector_t;
329 :
330 : vector_t c1 = CalculateCenter (e1, aaPos);
331 : vector_t c2 = CalculateCenter (e2, aaPos);
332 :
333 : c2 -= c1;
334 : return c2;
335 : }
336 :
337 : template <class TElem1, class TElem2, class TAAPos>
338 : bool CheckDirection (TElem1* e1,
339 : TElem2* e2,
340 : const TAAPos& aaPos,
341 : const typename TAAPos::ValueType& dir,
342 : number minAngle,
343 : number maxAngle)
344 : {
345 : typedef typename TAAPos::ValueType vector_t;
346 :
347 : const vector_t v = GetDirection (e1, e2, aaPos);
348 : const number angle = rad_to_deg(VecAngle(v, dir));
349 : return (minAngle <= angle) && (maxAngle >= angle);
350 : }
351 :
352 : }// end of namespace
353 :
354 : #endif
|
