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, dimworld> enBary;
132 FieldVector<ctype, dimension-1> mid(0.5);
134 const int map[3] = {1, 2, 4};
136 const Geometry geo = it->geometry();
137 if ( ! geo.type().isCube() )
return false;
140 for(
int i=0; i<dimension; ++i)
142 h[i] = (geo.corner( 0 ) - geo.corner( map[i] )).two_norm();
147 for(MacroIteratorType it = macroView.template begin<0> ();
150 const EntityType& en = *it;
151 const Geometry geo = en.geometry();
153 const FieldVector<ctype, dimworld> enBary =
154 geo.global( geoInfo.localCenter( geo.type() ));
156 typedef typename MacroGridView :: IntersectionIterator IntersectionIteratorType;
159 if ( ! geo.type().isCube() )
return false;
161 for(
int i=0; i<dimension; ++i)
163 const ctype w = (geo.corner(0) - geo.corner( map[i] )).two_norm();
164 if(
std::abs( h[i] - w ) > 1e-15 )
return false;
167 IntersectionIteratorType endnit = macroView.iend( en );
168 for(IntersectionIteratorType nit = macroView.ibegin( en );
169 nit != endnit; ++nit)
171 const typename IntersectionIteratorType::Intersection& inter=*nit;
174 if( inter.neighbor() )
176 EntityType nb = inter.outside();
177 Geometry nbGeo = nb.geometry();
179 FieldVector<ctype, dimworld> diff = nbGeo.global( geoInfo.localCenter( nbGeo.type() ));
181 assert( diff.two_norm() > 1e-15 );
182 diff /= diff.two_norm();
185 if(
std::abs(
std::abs((diff * inter.unitOuterNormal(mid))) - 1.0) > 1e-12 )
return false;
192 template <
class ctype,
int dim>
195 const FieldVector<ctype,dim-1> mid_;
196 enum { numberOfNormals = 2 * dim };
197 FieldVector<ctype,dim> refNormal_[numberOfNormals];
201 for(
int i=0; i<numberOfNormals; ++i)
204 FieldVector<ctype,dim>& refNormal = refNormal_[i];
208 int comp = ((int) i/2);
209 refNormal[comp] = ((i % 2) == 0) ? -1 : 1;
215 assert( i >= 0 && i< numberOfNormals );
216 return refNormal_[i];
227 template <
class IntersectionType>
230 if ( ! nit.type().isCube() )
return false;
232 typedef typename IntersectionType :: Entity EntityType;
233 typedef typename EntityType :: Geometry :: ctype ctype;
234 enum { dimworld = EntityType :: Geometry :: coorddimension };
240 FieldVector<ctype,dimworld> unitNormal = nit.unitOuterNormal(normals.
faceMidPoint());
244 if( unitNormal.infinity_norm() > 1e-10 )
return false;
250 static inline bool check(
const GridPartType& gridPart)
252 bool cartesian =
doCheck( gridPart.grid() , gridPart.indexSet() );
253 int val = (cartesian) ? 1 : 0;
255 val = gridPart.comm().min( val );
256 return (val == 1) ? true :
false;
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