Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
meshutil.hpp
1 /*
2  This file is part of the Feel library
3  Copyright (C) 2001,2002,2003,2004 EPFL, INRIA and Politechnico di Milano
4 
5  This library is free software; you can redistribute it and/or
6  modify it under the terms of the GNU Lesser General Public
7  License as published by the Free Software Foundation; either
8  version 3.0 of the License, or (at your option) any later version.
9 
10  This library is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13  Lesser General Public License for more details.
14 
15  You should have received a copy of the GNU Lesser General Public
16  License along with this library; if not, write to the Free Software
17  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18 */
19 #ifndef __MESH_UTIL_BASE__
20 #define __MESH_UTIL_BASE__
21 
22 #include <vector>
23 #include <algorithm>
24 #include <set>
25 
26 #include <feel/feelcore/feel.hpp>
27 
28 #include <feel/feelmesh/entities.hpp>
29 #include <feel/feelmesh/bareitems.hpp>
30 #include <feel/feelmesh/marker.hpp>
31 
32 
34 namespace Feel
35 {
44 typedef std::map<BareFace, std::pair<size_type, size_type >, cmpBareItem<BareFace> > TempFaceContainer;
46 
48 typedef std::map<BareEdge, std::pair<size_type, size_type>, cmpBareItem<BareEdge> > TempEdgeContainer;
49 
50 template<typename MeshType>
51 struct TempEntityContainer
52 {
53  typedef typename mpl::if_<mpl::equal_to<mpl::int_<MeshType::nDim>, mpl::int_<3> >,
54  mpl::identity<BareFace>,
55  mpl::identity<BareEdge> >::type::type entity_type;
56  typedef std::map<entity_type, std::pair<size_type, size_type >, cmpBareItem<entity_type> > type;
57 };
58 template<typename Ele, int Dim>
59 struct MakeBareEntity
60 {};
61 template<typename Ele>
62 struct MakeBareEntity<Ele, 3>
63 {
64  typedef BareFace entity_type;
65  static const int numVertices = Ele::GeoBShape::numVertices;
66 
67  MakeBareEntity( Ele const& ele )
68  :
69  M_element( ele )
70  {}
71  entity_type
72  operator()() const
73  {
74  entity_type bface;
75 
76  size_type i1 = ( M_element.point( 0 ) ).id();
77  size_type i2 = ( M_element.point( 1 ) ).id();
78  size_type i3 = ( M_element.point( 2 ) ).id();
79 
80  if ( Ele::face_type::numVertices == 4 )
81  {
82  size_type i4 = ( M_element.point( 3 ) ).id();
83  bface = makeBareFace( i1, i2, i3, i4 ).first;
84  }
85 
86  else
87  {
88  bface = makeBareFace( i1, i2, i3 ).first;
89  }
90 
91  return bface;
92  }
93  entity_type
94  operator()( uint16_type j ) const
95  {
96  entity_type bface;
97  size_type i1 = M_element.fToP( j, 0 );
98  size_type i2 = M_element.fToP( j, 1 );
99  size_type i3 = M_element.fToP( j, 2 );
100  // go to global
101  i1 = ( M_element.point( i1 ) ).id();
102  i2 = ( M_element.point( i2 ) ).id();
103  i3 = ( M_element.point( i3 ) ).id();
104 
105  if ( numVertices == 4 )
106  {
107  size_type i4 = M_element.fToP( j, 3 );
108  i4 = ( M_element.point( i4 ) ).id();
109  bface = ( makeBareItem( i1, i2, i3, i4 ) ).first;
110  }
111 
112  else
113  {
114  bface = ( makeBareItem( i1, i2, i3 ) ).first;
115  }
116 
117  return bface;
118  }
119  Ele const& M_element;
120 };
121 
122 template<typename Ele>
123 struct MakeBareEntity<Ele, 2>
124 {
125  typedef BareEdge entity_type;
126 
127  MakeBareEntity( Ele const& ele )
128  :
129  M_element( ele )
130  {}
131 
132  entity_type
133  operator()() const
134  {
135  size_type i1 = ( M_element.point( 0 ) ).id();
136  size_type i2 = ( M_element.point( 1 ) ).id();
137  entity_type bface;
138  bface = makeBareEdge( i1, i2 ).first;
139  return bface;
140  }
141 
142  entity_type
143  operator()( uint16_type j ) const
144  {
145  entity_type bface;
146  size_type i1 = M_element.fToP( j, 0 );
147  size_type i2 = M_element.fToP( j, 1 );
148  // go to global
149  i1 = ( M_element.point( i1 ) ).id();
150  i2 = ( M_element.point( i2 ) ).id();
151  bface = ( makeBareItem( i1, i2 ) ).first;
152  return bface;
153  }
154  Ele const& M_element;
155 };
156 
157 template<typename Ele>
158 struct MakeBareEntity<Ele, 1>
159 {
160  typedef BarePoint entity_type;
161 
162  MakeBareEntity( Ele const& ele )
163  :
164  M_element( ele )
165  {}
166 
167  entity_type
168  operator()() const
169  {
170  size_type i1 = ( M_element.point( 0 ) ).id();
171  DVLOG(2) << "[mesh1d::updateFaces] point index in face " << i1 << "\n";
172  entity_type bface;
173  bface = makeBarePoint( i1 ).first;
174  return bface;
175  }
176  entity_type
177  operator()( uint16_type j ) const
178  {
179  entity_type bface;
180  size_type i1 = M_element.point( j ).id();
181  bface = makeBarePoint( i1 ).first;
182  return bface;
183  }
184  Ele const& M_element;
185 };
186 
187 
188 template<typename Ele, int Dim>
189 struct MakeBareEntityFromFace
190 {};
191 
192 template<typename Ele>
193 struct MakeBareEntityFromFace<Ele, 3>
194 {
195  typedef BareFace entity_type;
196  static const int numVertices = Ele::numVertices;
197 
198  MakeBareEntityFromFace( Ele const& ele )
199  :
200  M_element( ele )
201  {}
202 
203  entity_type
204  operator()() const
205  {
206  entity_type bface;
207  // go to global
208  size_type i1 = ( M_element.point( 0 ) ).id();
209  size_type i2 = ( M_element.point( 1 ) ).id();
210  size_type i3 = ( M_element.point( 2 ) ).id();
211 
212  if ( numVertices == 4 )
213  {
214  size_type i4 = ( M_element.point( 3 ) ).id();
215  bface = ( makeBareItem( i1, i2, i3, i4 ) ).first;
216  }
217 
218  else
219  {
220  bface = ( makeBareItem( i1, i2, i3 ) ).first;
221  }
222 
223  return bface;
224  }
225  Ele const& M_element;
226 };
227 
228 template<typename Ele>
229 struct MakeBareEntityFromFace<Ele, 2>
230 {
231  typedef BareEdge entity_type;
232  static const int numVertices = Ele::numVertices;
233 
234  MakeBareEntityFromFace( Ele const& ele )
235  :
236  M_element( ele )
237  {}
238 
239  entity_type
240  operator()() const
241  {
242  entity_type bface;
243  // go to global
244  size_type i1 = ( M_element.point( 0 ) ).id();
245  size_type i2 = ( M_element.point( 1 ) ).id();
246  bface = ( makeBareItem( i1, i2 ) ).first;
247  return bface;
248  }
249  Ele const& M_element;
250 };
251 
252 } // Feel
253 
254 #include <feel/feelmesh/sphere.hpp>
255 
256 namespace Feel
257 {
262 typedef std::pair<Point, Point> MeshBoundingBox;
263 
269 template<typename MeshType>
270 inline
271 MeshBoundingBox
272 boundingBox ( const MeshType& mesh )
273 {
274  // processor bounding box with no arguments
275  // computes the global bounding box
276  return processorBoundingBox( mesh );
277 }
278 
279 
283 template<typename MeshType>
284 inline
285 Sphere
286 boundingSphere ( const MeshType& mesh )
287 {
288  MeshBoundingBox bbox = boundingBox( mesh );
289 
290  const double diag = Feel::distance( bbox.second, bbox.first );
291  const Point cent = Feel::middle( bbox.second, bbox.first );
292 
293  return Sphere ( cent, .5*diag );
294 
295 }
296 
302 template<typename MeshType>
303 inline
304 MeshBoundingBox
305 processorBoundingBox ( const MeshType& mesh,
306  const size_type pid = invalid_size_type_value )
307 {
308  FEELPP_ASSERT ( mesh.numPoints() != 0 ).error( "mesh has no points" );
309 
310  Point min( 1.e30, 1.e30, 1.e30 );
311  Point max( -1.e30, -1.e30, -1.e30 );
312 
313  // By default no processor is specified and we compute
314  // the bounding box for the whole domain.
315  if ( pid == invalid_size_type_value )
316  {
317  DVLOG(2) << "[processorBoundingBox] np pid given\n";
318 
319  for ( unsigned int n=0; n<mesh.numPoints(); n++ )
320  for ( unsigned int i=0; i<mesh.dimension(); i++ )
321  {
322  min( i ) = std::min( min( i ), mesh.point( n )( i ) );
323  max( i ) = std::max( max( i ), mesh.point( n )( i ) );
324  }
325  }
326 
327  // if a specific processor id is specified then we need
328  // to only consider those elements living on that processor
329  else
330  {
331  DVLOG(2) << "[processorBoundingBox] process bounding box on pid " << pid << "\n";
332  typename MeshType::element_iterator it = mesh.beginElementWithProcessId( pid );
333  typename MeshType::element_iterator en = mesh.endElementWithProcessId( pid );
334 
335  for ( ; it != en; ++it )
336  for ( unsigned int n=0; n< MeshType::element_type::numPoints; n++ )
337  for ( unsigned int i=0; i<mesh.dimension(); i++ )
338  {
339  min( i ) = std::min( min( i ), mesh.point( n )( i ) );
340  max( i ) = std::max( max( i ), mesh.point( n )( i ) );
341  }
342  }
343 
344  for ( unsigned int i=mesh.dimension(); i< min.node().size(); i++ )
345  {
346  min( i ) = 0;
347  max( i ) = 0;
348  }
349 
350  DVLOG(2) << "[processorBoundingBox] min= " << min << "\n";
351  DVLOG(2) << "[processorBoundingBox] max= " << max << "\n";
352  const MeshBoundingBox ret_val( min, max );
353 
354  return ret_val;
355 }
356 
360 template<typename MeshType>
361 inline
362 Sphere
363 processorBoundingSphere ( const MeshType& mesh,
364  const size_type pid = invalid_size_type_value )
365 {
366  MeshBoundingBox bbox = processorBoundingBox( mesh,pid );
367 
368  const Real diag = Feel::distance( bbox.second, bbox.first );
369  const Point cent = Feel::middle( bbox.second, bbox.first );
370 
371  DVLOG(2) << "[processorBoundingSphere] processor " << mesh.comm().rank() << "\n";
372  DVLOG(2) << "[processorBoundingSphere] center " << cent << "\n";
373  DVLOG(2) << "[processorBoundingSphere] radius " << 0.5*diag << "\n";
374  return Sphere ( cent, .5*diag );
375 }
376 
377 
378 } // Feel
379 
381 #endif

Generated on Sun Oct 20 2013 08:25:02 for Feel++ by doxygen 1.8.4