/usr/include/dune/pdelab/common/functionutilities.hh is in libdune-pdelab-dev 2.0.0-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 | // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=8 sw=2 sts=2:
#ifndef DUNE_PDELAB_COMMON_FUNCTIONUTILITIES_HH
#define DUNE_PDELAB_COMMON_FUNCTIONUTILITIES_HH
#include <limits>
#include <ostream>
#include <dune/common/debugstream.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/type.hh>
#include <dune/grid/common/gridenums.hh>
#include <dune/grid/utility/hierarchicsearch.hh>
namespace Dune {
namespace PDELab {
//! \addtogroup PDELab_Function Function
//! \ingroup PDELab
//! \{
//! Integrate a GridFunction
/**
* \code
#include <dune/pdelab/common/functionutilities.hh>
* \endcode
*
* Integrate a GridFunction over the domain given by the GridFunction's
* GridView. In the parallel case, this function integrates over the
* Interior_Partition only. If the accumulated result over all processors
* result is required, use something like
* \code
integrateGridFunction(gf, sum);
sum = gf.getGridView().comm().sum(sum);
* \endcode
*
* \tparam GF Type of the GridFunction.
* \param gf The GridFunction object.
* \param sum Resulting integral. There is no need to clear this
* variable before calling this function.
* \param qorder Quadrature order to use. If the GridFunction is
* element-wise polynomial, then this is the order of the
* highest-order monom needed to represent the function.
*/
template<typename GF>
void integrateGridFunction(const GF& gf,
typename GF::Traits::RangeType& sum,
unsigned qorder = 1) {
typedef typename GF::Traits::GridViewType GV;
typedef typename GV::template Codim<0>::
template Partition<Interior_Partition>::Iterator EIterator;
typedef typename GV::template Codim<0>::Geometry Geometry;
typedef typename GF::Traits::RangeType Range;
typedef typename GF::Traits::DomainFieldType DF;
static const int dimD = GF::Traits::dimDomain;
typedef Dune::QuadratureRule<DF,dimD> QR;
typedef Dune::QuadratureRules<DF,dimD> QRs;
typedef typename QR::const_iterator QIterator;
sum = 0;
Range val;
const EIterator eend = gf.getGridView().template end<0,
Interior_Partition>();
for(EIterator eit = gf.getGridView().template begin<0,
Interior_Partition>(); eit != eend; ++eit) {
const Geometry& geo = eit->geometry();
Dune::GeometryType gt = geo.type();
const QR& rule = QRs::rule(gt,qorder);
const QIterator qend = rule.end();
for (QIterator qit=rule.begin(); qit != qend; ++qit)
{
// evaluate the given grid functions at integration point
gf.evaluate(*eit,qit->position(),val);
// accumulate error
val *= qit->weight() * geo.integrationElement(qit->position());
sum += val;
}
}
}
//! Evaluate a GridFunction at a certain global coordinate
/**
* This class should work correctly even in a parallel setup. We look for
* the entity on containing the global coordinate on all ranks, then find
* the rank that actually has the entity. The GridFunction is then
* evaluated only for that rank, and the result is communicated to all the
* other ranks.
*
* \tparam GF Type of the GridFunction to evaluate.
*/
template<typename GF>
class GridFunctionProbe {
typedef typename GF::Traits::GridViewType GV;
typedef typename GV::template Codim<0>::EntityPointer EPtr;
typedef typename GF::Traits::DomainType Domain;
typedef typename GF::Traits::RangeType Range;
public:
//! Constructor
/**
* Construct the object and find the rank and the entity containing the
* global coordinate \c xg. The is currently not facility to recompute
* this information later, e.g. after a loadBalance() on the grid.
*
* \param gf The GridFunction to probe, either as a reference, a
* pointer, or a shared_ptr.
* \param xg The global coordinate the evaluate at.
*/
template<class GFHandle>
GridFunctionProbe(const GFHandle& gf, const Domain& xg)
{
setGridFunction(gf);
xl = 0;
evalRank = gfp->getGridView().comm().size();
int myRank = gfp->getGridView().comm().rank();
try {
e.reset(new EPtr
(HierarchicSearch<typename GV::Grid, typename GV::IndexSet>
(gfp->getGridView().grid(), gfp->getGridView().indexSet()).
template findEntity<Interior_Partition>(xg)));
// make sure only interior entities are accepted
if((*e)->partitionType() == InteriorEntity)
evalRank = myRank;
}
catch(const Dune::GridError&) { /* do nothing */ }
evalRank = gfp->getGridView().comm().min(evalRank);
if(myRank == evalRank)
xl = (*e)->geometry().local(xg);
else
e.reset();
if(myRank == 0 && evalRank == gfp->getGridView().comm().size())
dwarn << "Warning: GridFunctionProbe at (" << xg << ") is outside "
<< "the grid" << std::endl;
}
//! Set a new GridFunction
/**
* This takes the GridFunction as a refence. The referenced object must
* be valid for as long a the GridFunctionProbe is evaluated, or until
* setGridFunction() is called again.
*/
void setGridFunction(const GF &gf) {
gfsp.reset();
gfp = &gf;
}
//! Set a new GridFunction
/**
* This takes the GridFunction as a pointer. Ownership of the
* GridFunction object is transferred to the GridFunctionProbe. The
* GridFunction object will be deleted by the GridFunctionProbe on
* destruction, or the next time setGridFunction is called.
*/
void setGridFunction(const GF *gf) {
gfsp.reset(gf);
gfp = gf;
}
//! Set a new GridFunction
/**
* This takes the GridFunction as a shared_ptr. Ownership of the
* GridFunction object is shared with other places in the program, that
* hold a shared_ptr to the same GridFunction object. The GridFunction
* object will be deleted by the GridFunctionProbe on destruction, or
* the next time setGridFunction is called, if this GridFunctionProbe
* holds the last reference at that time.
*/
void setGridFunction(const Dune::shared_ptr<const GF> &gf) {
gfsp = gf;
gfp = &*gf;
}
//! evaluate the GridFunction and broadcast result to all ranks
/**
* \param val Store the result here.
*
* \note If the GridFunctionProbe is outside the grid NaN will be stored
* in val.
*/
void eval_all(Range& val) const {
typedef typename GF::Traits::RangeFieldType RF;
if(evalRank == gfp->getGridView().comm().size())
val = std::numeric_limits<RF>::quiet_NaN();
else {
if(gfp->getGridView().comm().rank() == evalRank)
gfp->evaluate(**e, xl, val);
gfp->getGridView().comm().broadcast(&val,1,evalRank);
}
}
//! evaluate the GridFunction and communicate result to the given rank
/**
* \param val Store the result here. This variable should be considered
* undefined on any rank other than the one given in \c rank
* after the function returns.
* \param rank Rank of the process to communicate the result to.
*
* \note CollectiveCommunication does not provide any direct
* process-to-process communication, so currently this function is
* identical with eval_all(), with the \c rank parameter ignored.
* \note If the GridFunctionProbe is outside the grid NaN will be stored
* in val.
*/
void eval(Range& val, int rank = 0) const {
eval_all(val);
}
private:
Dune::shared_ptr<const GF> gfsp;
const GF *gfp;
Dune::shared_ptr<EPtr> e;
Domain xl;
int evalRank;
};
//! \} Function
} // namespace PDELab
} // namespace Dune
#endif // DUNE_PDELAB_COMMON_FUNCTIONUTILITIES_HH
|