dune-fem  2.6-git
checkgeomaffinity.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_CHECKGEOMETRYAFFINITY_HH
2 #define DUNE_FEM_CHECKGEOMETRYAFFINITY_HH
3 
4 #include <dune/common/fvector.hh>
5 #include <dune/grid/common/gridenums.hh>
6 
8 
9 namespace Dune
10 {
11 
12  namespace Fem
13  {
14 
20  template <class QuadratureType>
22  {
23  // return true if geometry is affine
24  template< class EntityType, class ElementGeometryType >
25  static bool checkGeometry( const EntityType& entity,
26  const ElementGeometryType& geo,
27  const int quadOrd )
28  {
29  // if method tells that geometry is not affine
30  // then check it carefully
31  if( ! geo.affine() )
32  {
33  // get quadrature of desired order
34  QuadratureType volQuad( entity, quadOrd );
35  const int nop = volQuad.nop();
36 
37  // check all integration elements against the first
38  const double oldIntel = geo.integrationElement( volQuad.point(0) );
39  for(int l=1; l<nop; ++l)
40  {
41  const double intel = geo.integrationElement( volQuad.point(l) );
42  if( std::abs( oldIntel - intel ) > 1e-12 )
43  return false;
44  }
45  }
46  return true ;
47  }
48 
49  // return true if geometry is affine
50  template< class EntityType >
51  static bool checkGeometry( const EntityType& entity, const int quadOrd )
52  {
53  return checkGeometry( entity, entity.geometry(), quadOrd );
54  }
55 
57  template <class IteratorType>
58  static inline bool checkAffinity(const IteratorType& begin,
59  const IteratorType& endit,
60  const int quadOrd)
61  {
62  for(IteratorType it = begin; it != endit; ++it)
63  {
64  if( ! checkGeometry( *it, quadOrd ) ) return false ;
65  }
66  return true;
67  }
68 
70  template <class GridPartType, class Vector >
71  static inline void checkElementAffinity(const GridPartType& gridPart,
72  const int quadOrd,
73  Vector& affineGeomtryVec )
74  {
75  typedef typename GridPartType :: template Codim< 0 > :: IteratorType IteratorType;
76  typedef typename GridPartType :: template Codim< 0 > :: EntityType EntityType;
77  const IteratorType endit = gridPart.template end<0> ();
78  affineGeomtryVec.resize( gridPart.indexSet().size( 0 ) );
79  for(IteratorType it = gridPart.template begin<0>(); it != endit; ++it)
80  {
81  const EntityType& entity = *it ;
82  const int index = gridPart.indexSet().index( entity );
83  affineGeomtryVec[ index ] = checkGeometry( entity, quadOrd );
84  }
85 
86  //for( size_t i=0; i<affineGeomtryVec.size(); ++ i)
87  // std::cout << "geo is " << affineGeomtryVec[ i ] << std::endl;
88  }
89  };
90 
92  template <class GridPartType>
94  {
95  typedef typename GridPartType :: GridType GridType ;
96  protected:
98  template <class IndexSetType>
99  static inline bool doCheck(const GridType& grid, const IndexSetType& indexSet)
100  {
101  typedef typename GridType :: template Codim<0> :: LevelIterator MacroIteratorType;
102  typedef typename GridType :: template Codim<0> :: Entity EntityType;
103  typedef typename GridType :: template Codim<0> :: Geometry Geometry;
104 
105  typedef typename GridType :: LevelGridView MacroGridView ;
106 
107  // get macro grid view
108  MacroGridView macroView = grid.levelGridView( 0 );
109 
110  const MacroIteratorType endit = macroView.template end<0> ();
111  MacroIteratorType it = macroView.template begin<0> ();
112 
113  // empty grids are considerd as cartesian
114  if( it == endit ) return true;
115 
116  typedef AllGeomTypes< IndexSetType, GridType> GeometryInformationType;
117  GeometryInformationType geoInfo( indexSet );
118 
119  // if we have more than one geometry Type return false
120  if( geoInfo.geomTypes(0).size() > 1 ) return false;
121 
122  // if this type is not cube return false
123  if( ! geoInfo.geomTypes(0)[0].isCube() ) return false;
124 
125  typedef typename GridType :: ctype ctype;
126  enum { dimension = GridType :: dimension };
127  enum { dimworld = GridType :: dimensionworld };
128 
129  // grid width
130  FieldVector<ctype, dimension> h(0);
131  FieldVector<ctype, dimworld> enBary;
132  FieldVector<ctype, dimension-1> mid(0.5);
133 
134  const int map[3] = {1, 2, 4};
135  {
136  const Geometry geo = it->geometry();
137  if ( ! geo.type().isCube() ) return false;
138 
139  // calculate grid with
140  for(int i=0; i<dimension; ++i)
141  {
142  h[i] = (geo.corner( 0 ) - geo.corner( map[i] )).two_norm();
143  }
144  }
145 
146  // loop over all macro elements
147  for(MacroIteratorType it = macroView.template begin<0> ();
148  it != endit; ++it)
149  {
150  const EntityType& en = *it;
151  const Geometry geo = en.geometry();
152 
153  const FieldVector<ctype, dimworld> enBary =
154  geo.global( geoInfo.localCenter( geo.type() ));
155 
156  typedef typename MacroGridView :: IntersectionIterator IntersectionIteratorType;
157 
158  // if geometry is not cube, it's not a cartesian grid
159  if ( ! geo.type().isCube() ) return false;
160 
161  for(int i=0; i<dimension; ++i)
162  {
163  const ctype w = (geo.corner(0) - geo.corner( map[i] )).two_norm();
164  if( std::abs( h[i] - w ) > 1e-15 ) return false;
165  }
166 
167  IntersectionIteratorType endnit = macroView.iend( en );
168  for(IntersectionIteratorType nit = macroView.ibegin( en );
169  nit != endnit; ++nit)
170  {
171  const typename IntersectionIteratorType::Intersection& inter=*nit;
172  if( ! checkIntersection(inter) ) return false;
173 
174  if( inter.neighbor() )
175  {
176  EntityType nb = inter.outside();
177  Geometry nbGeo = nb.geometry();
178 
179  FieldVector<ctype, dimworld> diff = nbGeo.global( geoInfo.localCenter( nbGeo.type() ));
180  diff -= enBary;
181  assert( diff.two_norm() > 1e-15 );
182  diff /= diff.two_norm();
183 
184  // scalar product should be 1 or -1
185  if( std::abs(std::abs((diff * inter.unitOuterNormal(mid))) - 1.0) > 1e-12 ) return false;
186  }
187  }
188  }
189  return true;
190  }
191 
192  template <class ctype, int dim>
194  {
195  const FieldVector<ctype,dim-1> mid_;
196  enum { numberOfNormals = 2 * dim };
197  FieldVector<ctype,dim> refNormal_[numberOfNormals];
198  public:
199  ReferenceNormals () : mid_(0.5)
200  {
201  for(int i=0; i<numberOfNormals; ++i)
202  {
203  // get i-th reference normal
204  FieldVector<ctype,dim>& refNormal = refNormal_[i];
205  // reset normal
206  refNormal = 0;
207  // set one component
208  int comp = ((int) i/2);
209  refNormal[comp] = ((i % 2) == 0) ? -1 : 1;
210  }
211  }
212 
213  const FieldVector<ctype,dim>& referenceNormal(const int i) const
214  {
215  assert( i >= 0 && i< numberOfNormals );
216  return refNormal_[i];
217  }
218 
219  const FieldVector<ctype,dim-1>& faceMidPoint() const
220  {
221  return mid_;
222  }
223  };
224 
225  public:
226  // check that element is provided following the DUNE reference cube
227  template <class IntersectionType>
228  static bool checkIntersection(const IntersectionType& nit)
229  {
230  if ( ! nit.type().isCube() ) return false;
231 
232  typedef typename IntersectionType :: Entity EntityType;
233  typedef typename EntityType :: Geometry :: ctype ctype;
234  enum { dimworld = EntityType :: Geometry :: coorddimension };
235 
236  // get reference normals
237  static const ReferenceNormals<ctype,dimworld> normals;
238 
239  // get current normal
240  FieldVector<ctype,dimworld> unitNormal = nit.unitOuterNormal(normals.faceMidPoint());
241  unitNormal -= normals.referenceNormal( nit.indexInInside() );
242 
243  // if normals are not equal grid is not cartesian
244  if( unitNormal.infinity_norm() > 1e-10 ) return false;
245 
246  return true;
247  }
248 
250  static inline bool check(const GridPartType& gridPart)
251  {
252  bool cartesian = doCheck( gridPart.grid() , gridPart.indexSet() );
253  int val = (cartesian) ? 1 : 0;
254  // take global minimum
255  val = gridPart.comm().min( val );
256  return (val == 1) ? true : false;
257  }
258  };
259 
261 
262  } // namespace Fem
263 
264 } // namespace Dune
265 
266 #endif // #ifndef DUNE_FEM_CHECKGEOMETRYAFFINITY_HH
Definition: bindguard.hh:11
Double abs(const Double &a)
Definition: double.hh:872
Helper class to check affinity of the grid's geometries.
Definition: checkgeomaffinity.hh:22
static bool checkGeometry(const EntityType &entity, const ElementGeometryType &geo, const int quadOrd)
Definition: checkgeomaffinity.hh:25
static void checkElementAffinity(const GridPartType &gridPart, const int quadOrd, Vector &affineGeomtryVec)
check whether all geometry mappings are affine
Definition: checkgeomaffinity.hh:71
static bool checkGeometry(const EntityType &entity, const int quadOrd)
Definition: checkgeomaffinity.hh:51
static bool checkAffinity(const IteratorType &begin, const IteratorType &endit, const int quadOrd)
check whether all geometry mappings are affine
Definition: checkgeomaffinity.hh:58
Helper class to check whether grid is structured or not.
Definition: checkgeomaffinity.hh:94
static bool checkIntersection(const IntersectionType &nit)
Definition: checkgeomaffinity.hh:228
GridPartType ::GridType GridType
Definition: checkgeomaffinity.hh:95
static bool doCheck(const GridType &grid, const IndexSetType &indexSet)
check whether this is a cartesian grid or not
Definition: checkgeomaffinity.hh:99
static bool check(const GridPartType &gridPart)
check whether all the is grid is a cartesian grid
Definition: checkgeomaffinity.hh:250
Definition: checkgeomaffinity.hh:194
ReferenceNormals()
Definition: checkgeomaffinity.hh:199
const FieldVector< ctype, dim > & referenceNormal(const int i) const
Definition: checkgeomaffinity.hh:213
const FieldVector< ctype, dim-1 > & faceMidPoint() const
Definition: checkgeomaffinity.hh:219
default implementation uses method geomTypes of given index set. Used in DiscreteFunctionSpaces.
Definition: allgeomtypes.hh:99