dune-grid
2.2.0
|
00001 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*- 00002 // vi: set et ts=8 sw=2 sts=2: 00003 00004 #ifndef DUNE_GRID_UTILITY_GRIDINFO_HH 00005 #define DUNE_GRID_UTILITY_GRIDINFO_HH 00006 00007 #include <algorithm> 00008 #include <cstddef> 00009 #include <functional> 00010 #include <limits> 00011 #include <map> 00012 #include <ostream> 00013 #include <string> 00014 #include <vector> 00015 00016 #include <dune/common/classname.hh> 00017 #include <dune/common/exceptions.hh> 00018 #include <dune/common/forloop.hh> 00019 #include <dune/common/fvector.hh> 00020 #include <dune/geometry/type.hh> 00021 #include <dune/geometry/referenceelements.hh> 00022 00023 #include <dune/geometry/mockgeometry.hh> 00024 00025 #include <dune/grid/common/mcmgmapper.hh> 00026 00027 namespace Dune { 00028 00030 template<class ctype> 00031 struct EntityInfo { 00033 std::size_t count; 00035 00040 ctype volumeMin; 00042 00047 ctype volumeMax; 00048 00050 00054 ctype volumeSum; 00055 00057 00062 EntityInfo() : 00063 count(0), volumeMin(std::numeric_limits<ctype>::infinity()), 00064 volumeMax(-std::numeric_limits<ctype>::infinity()), volumeSum(0) 00065 { } 00066 }; 00067 00069 00075 struct GridViewInfoGTCompare : 00076 public std::binary_function<GeometryType, GeometryType, bool> 00077 { 00079 inline bool operator()(const GeometryType &a, const GeometryType &b) const 00080 { 00081 return a.dim() < b.dim() || 00082 (a.dim() == b.dim() && (a.isNone() < b.isNone() || 00083 (a.isNone() == b.isNone() && (a.id() >> 1) < (b.id() >> 1)))); 00084 // topologyId is set to 0 for None, so no harm im comparing them even if 00085 // isNone()==true 00086 } 00087 }; 00088 00090 00095 template<class ctype> 00096 struct GridViewInfo : 00097 public std::map<GeometryType, EntityInfo<ctype>, GridViewInfoGTCompare> 00098 { 00100 std::string gridName; 00102 std::string gridViewName; 00104 00108 std::string partitionName; 00109 00111 00124 void print(std::ostream &stream, std::string prefix) const { 00125 if(!gridName.empty()) { 00126 stream << prefix << gridName << ":\n"; 00127 prefix += " "; 00128 } 00129 if(!gridViewName.empty()) { 00130 stream << prefix << gridViewName << ":\n"; 00131 prefix += " "; 00132 } 00133 if(!partitionName.empty()) { 00134 stream << prefix << partitionName << ":\n"; 00135 prefix += " "; 00136 } 00137 00138 typedef typename GridViewInfo::const_iterator Iterator; 00139 std::size_t dim = ~0; 00140 const Iterator &end = this->end(); 00141 for(Iterator it = this->begin(); it != end; ++it) { 00142 if(it->first.dim() != dim) { 00143 dim = it->first.dim(); 00144 stream << prefix << "Dim = " << dim << ":\n"; 00145 } 00146 stream << prefix << " " << it->first << ": Count = " 00147 << it->second.count << ", Volume range = " 00148 << "(" << it->second.volumeMin << ".." 00149 << it->second.volumeMax << "), Total volume = " 00150 << it->second.volumeSum << "\n"; 00151 } 00152 } 00153 }; 00154 00156 00161 template<class ctype> 00162 std::ostream &operator<<(std::ostream &stream, 00163 const GridViewInfo<ctype> &info) 00164 { 00165 info.print(stream, ""); 00166 return stream; 00167 } 00168 00169 #ifndef DOXYGEN 00170 00171 template<int codim> 00172 struct FillGridInfoOperation { 00173 template<class Entity, class Mapper, class Visited, class RefElem> 00174 static void apply(const Entity &e, const Mapper &mapper, Visited &visited, 00175 const typename Entity::Geometry &geo, 00176 const RefElem &refelem, 00177 GridViewInfo<typename Entity::ctype> &gridViewInfo) 00178 { 00179 typedef typename Entity::ctype ctype; 00180 static const std::size_t dimw = Entity::dimensionworld; 00181 static const std::size_t dim = Entity::dimension; 00182 std::vector<FieldVector<ctype, dimw> > coords; 00183 for(int i = 0; i < refelem.size(codim); ++i) { 00184 int index = mapper.map(e, i, codim); 00185 if(visited[index]) 00186 continue; 00187 visited[index] = true; 00188 00189 GeometryType gt = refelem.type(i, codim); 00190 coords.clear(); 00191 coords.resize( refelem.size(i, codim, dim) ); 00192 for(std::size_t corner = 0; corner < coords.size(); ++corner) 00193 coords[ corner ] = geo.corner( refelem.subEntity( i, codim, corner, dim ) ); 00194 MockGeometry<ctype, dim-codim, dimw> mygeo(gt, coords); 00195 00196 ctype volume = mygeo.volume(); 00197 EntityInfo<ctype> &ei = gridViewInfo[mygeo.type()]; 00198 ei.volumeMin = std::min(ei.volumeMin, volume); 00199 ei.volumeMax = std::max(ei.volumeMax, volume); 00200 ei.volumeSum += volume; 00201 } 00202 } 00203 }; 00204 00205 template<int dimgrid> 00206 struct MCMGNonElementLayout { 00207 bool contains(GeometryType gt) const { return gt.dim() < dimgrid; } 00208 }; 00209 #endif // !DOXYGEN 00210 00212 00216 template<class GV> 00217 void fillGridViewInfoSerial(const GV &gv, 00218 GridViewInfo<typename GV::ctype> &gridViewInfo) 00219 { 00220 typedef typename GV::ctype ctype; 00221 static const std::size_t dim = GV::dimension; 00222 typedef typename GV::template Codim<0>::Iterator EIterator; 00223 typedef typename GV::template Codim<0>::Geometry EGeometry; 00224 typedef typename GV::IntersectionIterator IIterator; 00225 typedef typename GV::IndexSet IndexSet; 00226 00227 typedef typename GridViewInfo<ctype>::iterator InfoIterator; 00228 00229 typedef GenericReferenceElements<ctype, dim> RefElems; 00230 00231 MultipleCodimMultipleGeomTypeMapper<GV, MCMGNonElementLayout> mapper(gv); 00232 std::vector<bool> visited(mapper.size(), false); 00233 00234 gridViewInfo.gridName = className<typename GV::Grid>(); 00235 gridViewInfo.gridViewName = className<GV>(); 00236 gridViewInfo.partitionName = ""; 00237 gridViewInfo.clear(); 00238 00239 const EIterator &eend = gv.template end<0>(); 00240 for(EIterator eit = gv.template begin<0>(); eit != eend; ++eit) { 00241 ctype volume = eit->geometry().volume(); 00242 EntityInfo<ctype> &ei = gridViewInfo[eit->type()]; 00243 ei.volumeMin = std::min(ei.volumeMin, volume); 00244 ei.volumeMax = std::max(ei.volumeMax, volume); 00245 ei.volumeSum += volume; 00246 00247 if(!eit->type().isNone()) { 00248 const EGeometry &geo = eit->geometry(); 00249 ForLoop<FillGridInfoOperation, 1, dim>:: 00250 apply(*eit, mapper, visited, geo, RefElems::general(eit->type()), 00251 gridViewInfo); 00252 } 00253 } 00254 00255 GeometryType gt; 00256 gt.makeNone(dim); 00257 if(gridViewInfo.count(gt) > 0) { 00258 for(std::size_t codim = 0; codim < dim; ++codim) { 00259 gt.makeNone(dim-codim); 00260 EntityInfo<ctype> & ei = gridViewInfo[gt]; 00261 ei.volumeMin = ei.volumeMax = ei.volumeSum = 00262 std::numeric_limits<ctype>::quiet_NaN(); 00263 } 00264 gt.makeNone(0); 00265 EntityInfo<ctype> & ei = gridViewInfo[gt]; 00266 ei.volumeMin = ei.volumeMax = ei.volumeSum = 0; 00267 } 00268 00269 const InfoIterator &end = gridViewInfo.end(); 00270 const IndexSet &is = gv.indexSet(); 00271 for(InfoIterator it = gridViewInfo.begin(); it != end; ++it) { 00272 it->second.count = is.size(it->first); 00273 if(it->second.count == 0) 00274 DUNE_THROW(Exception, "Found Entities of geomentry type " << 00275 it->first << " while iterating through the grid, but " 00276 "indexSet.size() == 0 for that geometry type"); 00277 } 00278 00279 } 00280 00281 } // namespace Dune 00282 00283 00284 #endif // DUNE_GRID_UTILITY_GRIDINFO_HH