/usr/include/dune/grid/utility/gridinfo.hh is in libdune-grid-dev 2.3.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 | // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_GRID_UTILITY_GRIDINFO_HH
#define DUNE_GRID_UTILITY_GRIDINFO_HH
#include <algorithm>
#include <cstddef>
#include <functional>
#include <limits>
#include <map>
#include <ostream>
#include <string>
#include <vector>
#include <dune/common/classname.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/forloop.hh>
#include <dune/common/fvector.hh>
#include <dune/geometry/type.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/multilineargeometry.hh>
#include <dune/grid/common/mcmgmapper.hh>
namespace Dune {
//! Structure to hold statistical information about one type of entity
template<class ctype>
struct EntityInfo {
//! number of entities in the set
std::size_t count;
//! minimum volume of all entities in the set.
/**
* If the set is empty, this is \f$+\infty\f$. If the volume of the
* entities cannot be determined (some instance of the none GeometryType)
* this is NaN.
*/
ctype volumeMin;
//! maximum volume of all entities in the set.
/**
* If the set is empty, this is \f$-\infty\f$. If the volume of the
* entities cannot be determined (some instance of the none GeometryType)
* this is NaN.
*/
ctype volumeMax;
//! sum of volumes of all entities in the set.
/**
* If the set is empty, this is 0. If the volume of the entities cannot
* be determined (some instance of the none GeometryType) this is NaN.
*/
ctype volumeSum;
//! initialize the structure
/**
* This assumes an empty set of entities so that information can be added
* later: \c count is set to 0, \c volumeMin to +infinity, and \c
* volumeMax to -infinity.
*/
EntityInfo() :
count(0), volumeMin(std::numeric_limits<ctype>::infinity()),
volumeMax(-std::numeric_limits<ctype>::infinity()), volumeSum(0)
{ }
};
//! Comparison object to sort GeometryType by majorly dimension
/**
* This differs from the standard GeometryType::operator<() in that it puts
* more emphasis on the dimension than the isNone() property. If it is used
* as the comparison object in a map where the key is of type GeometryType,
* all entries of one dimension will be lumped together.
*/
struct GridViewInfoGTCompare :
public std::binary_function<GeometryType, GeometryType, bool>
{
//! compare two GeometryTypes
inline bool operator()(const GeometryType &a, const GeometryType &b) const
{
return a.dim() < b.dim() ||
(a.dim() == b.dim() && (a.isNone() < b.isNone() ||
(a.isNone() == b.isNone() && (a.id() >> 1) < (b.id() >> 1))));
// topologyId is set to 0 for None, so no harm im comparing them even if
// isNone()==true
}
};
//! structure to hold information about a certain GridView.
/**
* This is a map from GeometryType to EntityInfo structures. The entries in
* the map are sorted in dimension-first, isNone()-second, and
* topologyId-last order.
*/
template<class ctype>
struct GridViewInfo :
public std::map<GeometryType, EntityInfo<ctype>, GridViewInfoGTCompare>
{
//! name of the grid class this information was extracted from
std::string gridName;
//! name of the class of the GridView this information was extracted from
std::string gridViewName;
//! name of the partition this information was extracted from
/**
* May be empty if not applicable (serial grids, for instance) or may be a
* combination of partitions such as "interior+border".
*/
std::string partitionName;
//! print the information contained in this object
/**
* \param stream Stream object to print to.
* \param prefix Prefix to print in front of each line.
*
* Sample output:
* \verbatim
* prefix>
* \endverbatim
*
* If \c gridName, \c gridViewName, or \c partitionName is emtpy, the
* corresponding line is not printed and no extra indentation is added for
* the subsequent lines.
*/
void print(std::ostream &stream, std::string prefix) const {
if(!gridName.empty()) {
stream << prefix << gridName << ":\n";
prefix += " ";
}
if(!gridViewName.empty()) {
stream << prefix << gridViewName << ":\n";
prefix += " ";
}
if(!partitionName.empty()) {
stream << prefix << partitionName << ":\n";
prefix += " ";
}
typedef typename GridViewInfo::const_iterator Iterator;
std::size_t dim = ~0;
const Iterator &end = this->end();
for(Iterator it = this->begin(); it != end; ++it) {
if(it->first.dim() != dim) {
dim = it->first.dim();
stream << prefix << "Dim = " << dim << ":\n";
}
stream << prefix << " " << it->first << ": Count = "
<< it->second.count << ", Volume range = "
<< "(" << it->second.volumeMin << ".."
<< it->second.volumeMax << "), Total volume = "
<< it->second.volumeSum << "\n";
}
}
};
//! write a GridViewInfo object
/**
* \relates GridViewInfo
*
* This is equivalent to callinf info.print(stream, "").
*/
template<class ctype>
std::ostream &operator<<(std::ostream &stream,
const GridViewInfo<ctype> &info)
{
info.print(stream, "");
return stream;
}
#ifndef DOXYGEN
//! operation for ForLoop Internally used by fillGridViewInfoSerial
template<int codim>
struct FillGridInfoOperation {
template<class Entity, class Mapper, class Visited, class RefElem>
static void apply(const Entity &e, const Mapper &mapper, Visited &visited,
const typename Entity::Geometry &geo,
const RefElem &refelem,
GridViewInfo<typename Entity::ctype> &gridViewInfo)
{
typedef typename Entity::ctype ctype;
static const std::size_t dimw = Entity::Geometry::dimensionworld;
static const std::size_t dim = Entity::dimension;
std::vector<FieldVector<ctype, dimw> > coords;
for(int i = 0; i < refelem.size(codim); ++i) {
int index = mapper.map(e, i, codim);
if(visited[index])
continue;
visited[index] = true;
GeometryType gt = refelem.type(i, codim);
coords.clear();
coords.resize( refelem.size(i, codim, dim) );
for(std::size_t corner = 0; corner < coords.size(); ++corner)
coords[ corner ] = geo.corner( refelem.subEntity( i, codim, corner, dim ) );
MultiLinearGeometry<ctype, dim-codim, dimw> mygeo(gt, coords);
ctype volume = mygeo.volume();
EntityInfo<ctype> &ei = gridViewInfo[mygeo.type()];
ei.volumeMin = std::min(ei.volumeMin, volume);
ei.volumeMax = std::max(ei.volumeMax, volume);
ei.volumeSum += volume;
}
}
};
template<int dimgrid>
struct MCMGNonElementLayout {
bool contains(GeometryType gt) const { return gt.dim() < dimgrid; }
};
#endif // !DOXYGEN
//! fill a GridViewInfo structure from a serial grid
/**
* If used on a parallel grid, it will gather information for entities of
* all partitions on each rank locally.
*/
template<class GV>
void fillGridViewInfoSerial(const GV &gv,
GridViewInfo<typename GV::ctype> &gridViewInfo)
{
typedef typename GV::ctype ctype;
static const std::size_t dim = GV::dimension;
typedef typename GV::template Codim<0>::Iterator EIterator;
typedef typename GV::template Codim<0>::Geometry EGeometry;
typedef typename GV::IndexSet IndexSet;
typedef typename GridViewInfo<ctype>::iterator InfoIterator;
typedef ReferenceElements<ctype, dim> RefElems;
MultipleCodimMultipleGeomTypeMapper<GV, MCMGNonElementLayout> mapper(gv);
std::vector<bool> visited(mapper.size(), false);
gridViewInfo.gridName = className<typename GV::Grid>();
gridViewInfo.gridViewName = className<GV>();
gridViewInfo.partitionName = "";
gridViewInfo.clear();
const EIterator &eend = gv.template end<0>();
for(EIterator eit = gv.template begin<0>(); eit != eend; ++eit) {
ctype volume = eit->geometry().volume();
EntityInfo<ctype> &ei = gridViewInfo[eit->type()];
ei.volumeMin = std::min(ei.volumeMin, volume);
ei.volumeMax = std::max(ei.volumeMax, volume);
ei.volumeSum += volume;
if(!eit->type().isNone()) {
const EGeometry &geo = eit->geometry();
ForLoop<FillGridInfoOperation, 1, dim>::
apply(*eit, mapper, visited, geo, RefElems::general(eit->type()),
gridViewInfo);
}
}
GeometryType gt;
gt.makeNone(dim);
if(gridViewInfo.count(gt) > 0) {
for(std::size_t codim = 0; codim < dim; ++codim) {
gt.makeNone(dim-codim);
EntityInfo<ctype> & ei = gridViewInfo[gt];
ei.volumeMin = ei.volumeMax = ei.volumeSum =
std::numeric_limits<ctype>::quiet_NaN();
}
gt.makeNone(0);
EntityInfo<ctype> & ei = gridViewInfo[gt];
ei.volumeMin = ei.volumeMax = ei.volumeSum = 0;
}
const InfoIterator &end = gridViewInfo.end();
const IndexSet &is = gv.indexSet();
for(InfoIterator it = gridViewInfo.begin(); it != end; ++it) {
it->second.count = is.size(it->first);
if(it->second.count == 0)
DUNE_THROW(Exception, "Found Entities of geomentry type " <<
it->first << " while iterating through the grid, but "
"indexSet.size() == 0 for that geometry type");
}
}
} // namespace Dune
#endif // DUNE_GRID_UTILITY_GRIDINFO_HH
|