1#ifndef DUNE_FEM_CHECKGEOMETRYAFFINITY_HH
2#define DUNE_FEM_CHECKGEOMETRYAFFINITY_HH
4#include <dune/common/fvector.hh>
5#include <dune/grid/common/gridenums.hh>
20 template <
class QuadratureType>
24 template<
class EntityType,
class ElementGeometryType >
26 const ElementGeometryType& geo,
34 QuadratureType volQuad( entity, quadOrd );
35 const int nop = volQuad.nop();
38 const double oldIntel = geo.integrationElement( volQuad.point(0) );
39 for(
int l=1; l<nop; ++l)
41 const double intel = geo.integrationElement( volQuad.point(l) );
42 if(
std::abs( oldIntel - intel ) > 1e-12 )
50 template<
class EntityType >
51 static bool checkGeometry(
const EntityType& entity,
const int quadOrd )
57 template <
class IteratorType>
59 const IteratorType& endit,
62 for(IteratorType it = begin; it != endit; ++it)
70 template <
class Gr
idPartType,
class Vector >
73 Vector& affineGeomtryVec )
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)
81 const EntityType& entity = *it ;
82 const int index = gridPart.indexSet().index( entity );
92 template <
class Gr
idPartType>
95 typedef typename GridPartType :: GridType
GridType ;
98 template <
class IndexSetType>
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;
105 typedef typename GridType :: LevelGridView MacroGridView ;
108 MacroGridView macroView = grid.levelGridView( 0 );
110 const MacroIteratorType endit = macroView.template end<0> ();
111 MacroIteratorType it = macroView.template begin<0> ();
114 if( it == endit )
return true;
117 GeometryInformationType geoInfo( indexSet );
120 if( geoInfo.geomTypes(0).size() > 1 )
return false;
123 if( ! geoInfo.geomTypes(0)[0].isCube() )
return false;
125 typedef typename GridType :: ctype ctype;
126 enum { dimension = GridType :: dimension };
127 enum { dimworld = GridType :: dimensionworld };
130 FieldVector<ctype, dimension> h(0);
131 FieldVector<ctype, dimension-1> mid(0.5);
133 const int map[3] = {1, 2, 4};
135 const Geometry geo = it->geometry();
136 if ( ! geo.type().isCube() )
return false;
139 for(
int i=0; i<dimension; ++i)
141 h[i] = (geo.corner( 0 ) - geo.corner( map[i] )).two_norm();
146 for(MacroIteratorType it = macroView.template begin<0> ();
149 const EntityType& en = *it;
150 const Geometry geo = en.geometry();
152 const FieldVector<ctype, dimworld> enBary =
153 geo.global( geoInfo.localCenter( geo.type() ));
155 typedef typename MacroGridView :: IntersectionIterator IntersectionIteratorType;
158 if ( ! geo.type().isCube() )
return false;
160 for(
int i=0; i<dimension; ++i)
162 const ctype w = (geo.corner(0) - geo.corner( map[i] )).two_norm();
163 if(
std::abs( h[i] - w ) > 1e-15 )
return false;
166 IntersectionIteratorType endnit = macroView.iend( en );
167 for(IntersectionIteratorType nit = macroView.ibegin( en );
168 nit != endnit; ++nit)
170 const typename IntersectionIteratorType::Intersection& inter=*nit;
173 if( inter.neighbor() )
175 EntityType nb = inter.outside();
176 Geometry nbGeo = nb.geometry();
178 FieldVector<ctype, dimworld> diff = nbGeo.global( geoInfo.localCenter( nbGeo.type() ));
180 assert( diff.two_norm() > 1e-15 );
181 diff /= diff.two_norm();
184 if(
std::abs(
std::abs((diff * inter.unitOuterNormal(mid))) - 1.0) > 1e-12 )
return false;
191 template <
class ctype,
int dim>
194 const FieldVector<ctype,dim-1> mid_;
195 enum { numberOfNormals = 2 * dim };
196 FieldVector<ctype,dim> refNormal_[numberOfNormals];
200 for(
int i=0; i<numberOfNormals; ++i)
203 FieldVector<ctype,dim>& refNormal = refNormal_[i];
207 int comp = ((int) i/2);
208 refNormal[comp] = ((i % 2) == 0) ? -1 : 1;
214 assert( i >= 0 && i< numberOfNormals );
215 return refNormal_[i];
226 template <
class IntersectionType>
229 if ( ! nit.type().isCube() )
return false;
231 typedef typename IntersectionType :: Entity EntityType;
232 typedef typename EntityType :: Geometry :: ctype ctype;
233 enum { dimworld = EntityType :: Geometry :: coorddimension };
239 FieldVector<ctype,dimworld> unitNormal = nit.unitOuterNormal(normals.
faceMidPoint());
243 if( unitNormal.infinity_norm() > 1e-10 )
return false;
249 static inline bool check(
const GridPartType& gridPart)
251 bool cartesian =
doCheck( gridPart.grid() , gridPart.indexSet() );
252 int val = (cartesian) ? 1 : 0;
254 val = gridPart.comm().min( val );
255 return (val == 1) ? true :
false;
Dune::Fem::Double abs(const Dune::Fem::Double &a)
Definition: double.hh:942
Definition: bindguard.hh:11
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:227
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:249
Definition: checkgeomaffinity.hh:193
ReferenceNormals()
Definition: checkgeomaffinity.hh:198
const FieldVector< ctype, dim > & referenceNormal(const int i) const
Definition: checkgeomaffinity.hh:212
const FieldVector< ctype, dim-1 > & faceMidPoint() const
Definition: checkgeomaffinity.hh:218
default implementation uses method geomTypes of given index set. Used in DiscreteFunctionSpaces.
Definition: allgeomtypes.hh:99