ug4
misc_util_impl.hpp
Go to the documentation of this file.
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"
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 
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;
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 
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 
116 // CalculateCenter
117 template <class TIterator, class TAAPosVRT>
118 typename TAAPosVRT::ValueType
119 CalculateCenter(TIterator begin, TIterator end, TAAPosVRT& aaPos)
120 {
121  int counter = 0;
122  typename TAAPosVRT::ValueType center;
123  VecSet(center, 0);
124  for(TIterator iter = begin; iter != end; ++iter, ++counter)
125  VecAdd(center, center, CalculateCenter(*iter, aaPos));
126 
127  if(counter > 0)
128  VecScale(center, center, 1./(number)counter);
129 
130  return center;
131 }
132 
134 // FindClosestByCoordinate
135 template<class TElem, class TVertexPositionAttachmentAccessor>
136 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  if(iterBegin == iterEnd)
142  return NULL;
143 
144  typename geometry_traits<TElem>::iterator iter = iterBegin;
145  TElem* bestElem = *iter;
146  number bestDistSq = VecDistanceSq(coord, CalculateCenter(bestElem, aaPosVRT));
147  iter++;
148 
149  while(iter != iterEnd)
150  {
151  number distSq = VecDistanceSq(coord, CalculateCenter(*iter, aaPosVRT));
152  if(distSq < bestDistSq)
153  {
154  bestDistSq = distSq;
155  bestElem = *iter;
156  }
157 
158  ++iter;
159  }
160 
161  return bestElem;
162 }
163 
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 
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 
220 // EraseElements
221 template <class TElem>
222 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  while(iterBegin != iterEnd)
227  {
228  TElem* e = *iterBegin;
229  iterBegin++;
230  grid.erase(e);
231  }
232 }
233 
235 // AssignIndices
236 template <class TElem>
238  typename geometry_traits<TElem>::iterator iterEnd,
240 {
241  int index = 0;
242  while(iterBegin != iterEnd)
243  {
244  aaInt[*iterBegin] = index++;
245  iterBegin++;
246  }
247 }
248 
250 // ElementDiameter
251 template <class TElem, class TAAPos>
253  TAAPos& aaPos,
254  TElem* elem)
255 {
257  grid.associated_elements(vVert, elem);
258 
259  number max = 0.0;
260  for(size_t i = 0; i < vVert.size(); ++i)
261  for(size_t j = i+1; j < vVert.size(); ++j)
262  max = std::max(max, VecDistanceSq(aaPos[vVert[i]], aaPos[vVert[j]]));
263 
264  return max;
265 }
266 
267 template <class TAAPos>
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>
283  TAAPos& aaPos,
284  TElem* elem)
285 {
286  return std::sqrt(ElementDiameterSq(grid, aaPos, elem));
287 }
288 
289 template <class TAAPos, class TIterator>
290 number MaxElementDiameter(Grid& grid, TAAPos& aaPos,
291  TIterator iter, TIterator iterEnd)
292 {
293  number max = 0.0;
294  for(; iter != iterEnd; ++iter)
295  max = std::max(max, ElementDiameterSq(grid, aaPos, *iter));
296 
297 #ifdef UG_PARALLEL
298  // share value between all procs
300  max = com.allreduce(max, PCL_RO_MAX);
301 #endif
302 
303  return std::sqrt(max);
304 }
305 
306 template <class TAAPos, class TIterator>
307 number MinElementDiameter(Grid& grid, TAAPos& aaPos,
308  TIterator iter, TIterator iterEnd)
309 {
310  number min = std::numeric_limits<number>::max();
311  for(; iter != iterEnd; ++iter)
312  min = std::min(min, ElementDiameterSq(grid, aaPos, *iter));
313 
314 #ifdef UG_PARALLEL
315  // share value between all procs
317  min = com.allreduce(min, PCL_RO_MIN);
318 #endif
319 
320  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
Definition: pcl_process_communicator.h:70
void allreduce(const void *sendBuf, void *recBuf, int count, DataType type, ReduceOperation op) const
performs MPI_Allreduce on the processes of the communicator.
Definition: pcl_process_communicator.cpp:318
Base-class for edges.
Definition: grid_base_objects.h:397
Faces are 2-dimensional objects.
Definition: grid_base_objects.h:510
the generic attachment-accessor for access to grids attachment pipes.
Definition: grid.h:182
Manages the elements of a grid and their interconnection.
Definition: grid.h:132
void end_marking()
ends a marking sequence. Call this method when you're done with marking.
Definition: grid.cpp:1285
bool is_marked(GridObject *obj) const
returns true if the object is marked, false if not.
Definition: grid_impl.hpp:843
void mark(GridObject *obj)
marks the object. Calls are only valid between calls to Grid::begin_marking and Grid::end_marking.
Definition: grid_impl.hpp:773
void erase(GridObject *geomObj)
Definition: grid.cpp:459
void begin_marking()
begin marking.
Definition: grid.cpp:1262
void associated_elements(traits< Vertex >::secure_container &elemsOut, TElem *e)
Puts all elements of type TAss which are contained in 'e' or which contain 'e' into elemsOut.
Definition: grid_impl.hpp:466
The base class for all geometric objects, such as vertices, edges, faces, volumes,...
Definition: grid_base_objects.h:157
virtual int base_object_id() const =0
Container which holds an array of pointers.
Definition: pointer_const_array.h:84
size_t size() const
returns the size of the associated array.
Definition: pointer_const_array_impl.hpp:106
Base-class for all vertex-types.
Definition: grid_base_objects.h:231
Volumes are 3-dimensional objects.
Definition: grid_base_objects.h:754
Definition: grid_base_object_traits.h:68
number ElementDiameter(const TElem &elem, TDomain &domain)
returns the maximal distance between to element vertices
number ElementDiameterSq(const TElem &elem, TDomain &domain)
returns the maximal squared distance between to element vertices
number MinElementDiameter(TDomain &domain, int level)
returns the minimal diameter of all elements between iterBegin and iterEnd.
Definition: domain_util.cpp:113
number MaxElementDiameter(TDomain &domain, int level)
returns the maximal diameter of all elements between iterBegin and iterEnd.
Definition: domain_util.cpp:104
void AssignIndices(TIterator begin, TIterator end, TAAInt &aaInt, int baseIndex=0)
assigns indices to the elements between begin and end.
Definition: attachment_util_impl.hpp:178
UG_API void EraseElements(Grid &grid, typename geometry_traits< TElem >::iterator iterBegin, typename geometry_traits< TElem >::iterator iterEnd)
erases all elements between iterBegin and iterEnd.
Definition: misc_util_impl.hpp:222
UG_API TElem * FindClosestByCoordinate(const typename TVertexPositionAttachmentAccessor::ValueType &coord, typename geometry_traits< TElem >::iterator iterBegin, typename geometry_traits< TElem >::iterator iterEnd, TVertexPositionAttachmentAccessor &aaPosVRT)
returns the element of the given type, whose center is closest to coord.
Definition: misc_util_impl.hpp:136
TAAPos::ValueType GetDirection(TElem1 *e1, TElem2 *e2, const TAAPos &aaPos)
Returns the direction from the center of e1 to the center of e2.
Definition: misc_util_impl.hpp:326
UG_API TAAPosVRT::ValueType CalculateGridObjectCenter(const GridObject *o, TAAPosVRT &aaPosVRT)
calculates the center for arbitrary geometric object
Definition: misc_util_impl.hpp:81
UG_API size_t NumSharedVertices(Grid &grid, TElemPtr1 elem1, TElemPtr2 elem2)
returns the number of vertices that are shared by two objects
Definition: misc_util_impl.hpp:200
UG_API TElem::side * GetSharedSide(Grid &grid, TElem *e1, TElem *e2)
returns the shared side between the two elements or NULL if no such side exists.
Definition: misc_util_impl.hpp:52
bool CheckDirection(TElem1 *e1, TElem2 *e2, const TAAPos &aaPos, const typename TAAPos::ValueType &dir, number minAngle, number maxAngle)
Checks whether the center of e2 can be reached from the center of e1 in the given direction.
Definition: misc_util_impl.hpp:338
Vertex * GetVertex(Vertex *v, size_t i)
returns the i'th vertex of a vertex
Definition: grid_util_impl.hpp:99
size_t NumVertices(Vertex *elem)
Returns the number of vertices of the given geometric object.
Definition: grid_util_impl.hpp:123
#define PCL_RO_MAX
Definition: pcl_methods.h:61
#define PCL_RO_MIN
Definition: pcl_methods.h:62
static const int dim
#define UG_THROW(msg)
Definition: error.h:57
#define UG_API
Definition: ug_config.h:65
double number
Definition: types.h:124
void CalculateCenter(vector_t &centerOut, const vector_t *pointSet, size_t numPoints)
calculates the center of a point-set
Definition: math_util_impl.hpp:98
TNumber rad_to_deg(TNumber rad)
Definition: math_util_impl.hpp:58
void VecSet(vector_t &vInOut, typename vector_t::value_type s)
Set each vector component to scalar (componentwise)
Definition: math_vector_functions_common_impl.hpp:539
vector_t::value_type VecAngle(const vector_t &v1, const vector_t &v2)
returns the angle between two vectors in radiants
Definition: math_vector_functions_common_impl.hpp:401
vector_t::value_type VecDistanceSq(const vector_t &v1, const vector_t &v2)
returns the squared distance of two vector_ts.
Definition: math_vector_functions_common_impl.hpp:351
void VecAdd(vector_t &vOut, const vector_t &v1, const vector_t &v2)
adds two MathVector<N>s and stores the result in a third one
Definition: math_vector_functions_common_impl.hpp:185
void VecScale(vector_t &vOut, const vector_t &v, typename vector_t::value_type s)
scales a MathVector<N>
Definition: math_vector_functions_common_impl.hpp:252
the ug namespace
AABox< typename TAAPos::ValueType > CalculateBoundingBox(Vertex *e, TAAPos aaPos)
calculates the smallest axis aligned box that contains the given vertex
Definition: bounding_box_util.h:43
@ VOLUME
Definition: grid_base_objects.h:63
@ VERTEX
Definition: grid_base_objects.h:60
@ EDGE
Definition: grid_base_objects.h:61
@ FACE
Definition: grid_base_objects.h:62
T value_type
Definition: sparsematrix_interface.h:2