dune-grid  2.2.0
utility/gridinfo.hh
Go to the documentation of this file.
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