dune-fem  2.8-git
vtxprojection.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_VTXPROJECTION_HH
2 #define DUNE_FEM_VTXPROJECTION_HH
3 
4 #include <cstddef>
5 #include <cmath>
6 
7 #include <algorithm>
8 #include <functional>
9 #include <vector>
10 
11 #include <dune/grid/common/rangegenerators.hh>
12 
18 
19 namespace Dune
20 {
21  namespace Fem
22  {
23 
24  template <class GridPartType>
25  struct WeightDefault {
26  typedef typename GridPartType::template Codim<0>::EntityType EntityType;
27  typedef FieldVector<double,GridPartType::dimensionworld> WorldDomainType;
28  typedef FieldVector<double,GridPartType::dimension> DomainType;
29  void setEntity(const EntityType& en) {
30  en_ = en;
31  volume_ = en.geometry().volume();
32  baryCenter_ = en.geometry().center();
33  }
34  double operator()(const DomainType& point) {
35  // return volume_;
36  auto tau = en_.geometry().global(point);
37  tau -= baryCenter_;
38  return tau.two_norm() / volume_;
39  }
40  double volume_;
43  };
44 
46  {
47  template< class DiscreteFunction, class Intersection >
49  {
50  typedef typename DiscreteFunction::LocalFunctionType LocalFunctionType;
51 
52  typedef typename LocalFunctionType::EntityType EntityType;
53  typedef typename LocalFunctionType::FunctionSpaceType FunctionSpaceType;
54 
55  typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
56  typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
57 
58  typedef typename FunctionSpaceType::DomainType DomainType;
59  typedef typename FunctionSpaceType::RangeType RangeType;
60 
61  typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
62  typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
63 
64  typedef typename EntityType::Geometry::LocalCoordinate LocalCoordinateType;
65 
66  static const int dimDomain = FunctionSpaceType::dimDomain;
67  static const int dimRange = FunctionSpaceType::dimRange;
68 
69  typedef typename Intersection::LocalGeometry GeometryType;
70 
71  explicit OutsideLocalFunction ( const DiscreteFunction &df ) : order_(df.order()), localFunction_( df ) {}
72 
73  void init ( const EntityType &outside, const GeometryType &geoIn, const GeometryType &geoOut )
74  {
75  localFunction_.init( outside );
76  geoIn_ = &geoIn;
77  geoOut_ = &geoOut;
78  entity_ = outside;
79  }
80 
81  const EntityType &entity() const { return entity_; }
82 
83  int order() const { return order_; }
84 
85  template< class Point >
86  void evaluate ( const Point &x, RangeType &value ) const
87  {
88  localFunction_.evaluate( geoOut_->global( geoIn_->local( coordinate( x ) ) ), value );
89  };
90 
91  template< class Point >
92  void jacobian ( const Point &x, JacobianRangeType &jacobian ) const
93  {
94  DUNE_THROW( NotImplemented, "Vertex projection weights do not provide Jacobians." );
95  }
96 
97  template< class Point >
98  void hessian ( const Point &x, HessianRangeType &hessian ) const
99  {
100  DUNE_THROW( NotImplemented, "Vertex projection weights do not provide Hessians." );
101  }
102 
103  template< class Quadrature, class Values >
104  void evaluateQuadrature ( const Quadrature &quadrature, Values &values ) const
105  {
106  for( const auto &qp : quadrature )
107  doEvaluate( qp, values[ qp.index() ] );
108  }
109 
110  private:
111  template< class Point >
112  void doEvaluate ( const Point &x, RangeType &value ) const
113  {
114  evaluate( x, value );
115  };
116 
117  template< class Point >
118  void doEvaluate ( const Point &x, JacobianRangeType &value ) const
119  {
120  jacobian( x, value );
121  };
122 
123  template< class Point >
124  void doEvaluate ( const Point &x, HessianRangeType &value ) const
125  {
126  hessian( x, value );
127  };
128 
129  int order_;
130  LocalFunctionType localFunction_;
131  const GeometryType *geoIn_ = nullptr;
132  const GeometryType *geoOut_ = nullptr;
133  EntityType entity_;
134  };
135 
136 
137  template< class DiscreteFunction >
138  static void makeConforming ( DiscreteFunction &u )
139  {
140  typedef typename DiscreteFunction::GridPartType GridPartType;
141 
142  typedef typename GridPartType::IntersectionType IntersectionType;
143  typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
144  typedef typename DiscreteFunction::DiscreteFunctionSpaceType
145  DiscreteFunctionSpaceType;
146 
147  if( GridPartType::Traits::conforming || !Fem::GridPartCapabilities::hasGrid< GridPartType >::v )
148  return;
149 
150  const auto &blockMapper = u.space().blockMapper();
151 
152  std::vector< typename DiscreteFunction::DofType > ldu, ldw;
153  const std::size_t localBlockSize = DiscreteFunctionSpaceType::localBlockSize;
154  const std::size_t maxNumBlocks = blockMapper.maxNumDofs();
155  ldu.reserve( maxNumBlocks * localBlockSize );
156  ldw.reserve( maxNumBlocks * localBlockSize );
157 
158  std::vector< bool > onSubEntity;
159  onSubEntity.reserve( maxNumBlocks );
160 
162 
163  LocalInterpolation< DiscreteFunctionSpaceType > interpolation( u.space() );
164 
165  for( const EntityType &inside : u.space() )
166  {
167  for( const IntersectionType &intersection : intersections( u.gridPart(), inside ) )
168  {
169  if( intersection.conforming() || !intersection.neighbor() )
170  continue;
171 
172  // skip this intersection, if inside is finer than outside
173  const EntityType outside = intersection.outside();
174  if( inside.level() <= outside.level() )
175  continue;
176 
177  // initialize "outside local function"
178  // note: intersection geometries must live until end of intersection loop!
179  const auto geoIn = intersection.geometryInInside();
180  const auto geoOut = intersection.geometryInOutside();
181  uOutside.init( outside, geoIn, geoOut );
182 
183  // resize local DoF vectors
184  const std::size_t numBlocks = blockMapper.numDofs( inside );
185  ldu.resize( numBlocks * localBlockSize );
186  ldw.resize( numBlocks * localBlockSize );
187 
188  // interpolate "outside" values
189  auto guard = bindGuard( interpolation, inside );
190  interpolation( uOutside, ldw );
191 
192  // fetch inside local DoFs
193  u.getLocalDofs( inside, ldu );
194 
195  // patch DoFs assigned to the face
196  blockMapper.onSubEntity( inside, intersection.indexInInside(), 1, onSubEntity );
197  for( std::size_t i = 0; i < numBlocks; ++i )
198  {
199  if( !onSubEntity[ i ] )
200  continue;
201  for( std::size_t j = 0; j < localBlockSize; ++j )
202  ldu[ i*localBlockSize + j ] = ldw[ i*localBlockSize + j ];
203  }
204 
205  // write back inside local DoFs
206  u.setLocalDofs( inside, ldu );
207  }
208  }
209  }
210 
211  template< class Function, class DiscreteFunction, class Weight >
212  static auto project ( const Function &f, DiscreteFunction &u, Weight &weight )
213  -> std::enable_if_t<std::is_void< decltype( interpolate(f,u,weight,u )) >::value>
214  // -> std::enable_if_t<std::is_void< decltype( interpolate(f,u,weight,std::declval<std::remove_cv<std::remove_reference<DiscreteFunction>>&>()) ) >::value>
215  {
216  DiscreteFunction w ( "weight", u.space() );
217  interpolate( f, u, weight, w );
218  makeConforming( u );
219  }
220  template< class Function, class DiscreteFunction >
221  static auto project ( const Function &f, DiscreteFunction &u )
222  -> std::enable_if_t< std::is_void< decltype( project( f,u,std::declval<WeightDefault<typename DiscreteFunction::GridPartType>&>() ) ) >::value>
223  {
225  project(f, u, weight);
226  }
227  };
228 
229  /*======================================================================*/
236  /*======================================================================*/
237  template < typename DType, typename RType >
238  class VtxProjection : public Operator< DType, RType >
239  {
240  public:
241  typedef DType DomainType;
242  typedef RType RangeType;
243  typedef typename DomainType :: RangeFieldType DomainFieldType;
244  typedef typename RType :: RangeFieldType RangeFieldType;
245  typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
246  typedef typename RType :: DiscreteFunctionSpaceType :: GridPartType GridPartType;
249 
251  template <class WeightType>
252  void operator() (const DomainType& f, RangeType& discFunc,
253  WeightType& weight) const
254  {
255  VtxProjectionImpl::project(f,discFunc,weight);
256  }
258  void operator() (const DomainType& f, RangeType& discFunc) const
259  {
260  VtxProjectionImpl::project(f,discFunc);
261  }
262 
263  private:
264  };
265 
266  } // namespace Fem
267 
268 } // name space Dune
269 
270 #endif // #ifndef DUNE_FEM_VTXPROJECTION_HH
static void interpolate(const GridFunction &u, DiscreteFunction &v)
perform native interpolation of a discrete function space
Definition: common/interpolate.hh:54
Definition: bindguard.hh:11
static const Point & coordinate(const Point &x)
Definition: coordinate.hh:14
static auto bindGuard(Object &object, Args &&... args) -> std::enable_if_t< isBindable< Object, Args... >::value, BindGuard< Object > >
Definition: bindguard.hh:67
Abstract class representing a function.
Definition: common/function.hh:50
specialize with 'false' if grid part has no underlying dune grid (default=true)
Definition: gridpart/common/capabilities.hh:18
abstract operator
Definition: operator.hh:34
Definition: vtxprojection.hh:25
double operator()(const DomainType &point)
Definition: vtxprojection.hh:34
FieldVector< double, GridPartType::dimensionworld > WorldDomainType
Definition: vtxprojection.hh:27
double volume_
Definition: vtxprojection.hh:40
GridPartType::template Codim< 0 >::EntityType EntityType
Definition: vtxprojection.hh:26
WorldDomainType baryCenter_
Definition: vtxprojection.hh:41
FieldVector< double, GridPartType::dimension > DomainType
Definition: vtxprojection.hh:28
void setEntity(const EntityType &en)
Definition: vtxprojection.hh:29
EntityType en_
Definition: vtxprojection.hh:42
Definition: vtxprojection.hh:46
static auto project(const Function &f, DiscreteFunction &u) -> std::enable_if_t< std::is_void< decltype(project(f, u, std::declval< WeightDefault< typename DiscreteFunction::GridPartType > & >())) >::value >
Definition: vtxprojection.hh:221
static auto project(const Function &f, DiscreteFunction &u, Weight &weight) -> std::enable_if_t< std::is_void< decltype(interpolate(f, u, weight, u)) >::value >
Definition: vtxprojection.hh:212
static void makeConforming(DiscreteFunction &u)
Definition: vtxprojection.hh:138
void init(const EntityType &outside, const GeometryType &geoIn, const GeometryType &geoOut)
Definition: vtxprojection.hh:73
static const int dimRange
Definition: vtxprojection.hh:67
FunctionSpaceType::JacobianRangeType JacobianRangeType
Definition: vtxprojection.hh:61
void evaluate(const Point &x, RangeType &value) const
Definition: vtxprojection.hh:86
int order() const
Definition: vtxprojection.hh:83
const EntityType & entity() const
Definition: vtxprojection.hh:81
void evaluateQuadrature(const Quadrature &quadrature, Values &values) const
Definition: vtxprojection.hh:104
FunctionSpaceType::RangeType RangeType
Definition: vtxprojection.hh:59
EntityType::Geometry::LocalCoordinate LocalCoordinateType
Definition: vtxprojection.hh:64
LocalFunctionType::FunctionSpaceType FunctionSpaceType
Definition: vtxprojection.hh:53
LocalFunctionType::EntityType EntityType
Definition: vtxprojection.hh:52
DiscreteFunction::LocalFunctionType LocalFunctionType
Definition: vtxprojection.hh:50
FunctionSpaceType::DomainFieldType DomainFieldType
Definition: vtxprojection.hh:55
Intersection::LocalGeometry GeometryType
Definition: vtxprojection.hh:69
static const int dimDomain
Definition: vtxprojection.hh:66
void jacobian(const Point &x, JacobianRangeType &jacobian) const
Definition: vtxprojection.hh:92
void hessian(const Point &x, HessianRangeType &hessian) const
Definition: vtxprojection.hh:98
FunctionSpaceType::HessianRangeType HessianRangeType
Definition: vtxprojection.hh:62
FunctionSpaceType::RangeFieldType RangeFieldType
Definition: vtxprojection.hh:56
OutsideLocalFunction(const DiscreteFunction &df)
Definition: vtxprojection.hh:71
FunctionSpaceType::DomainType DomainType
Definition: vtxprojection.hh:58
The Projection class which average discontinuous function using the interpolation of the space (e....
Definition: vtxprojection.hh:239
RType RangeType
Definition: vtxprojection.hh:242
VtxProjection()
Constructor.
Definition: vtxprojection.hh:248
DType DomainType
Definition: vtxprojection.hh:241
void operator()(const DomainType &f, RangeType &discFunc, WeightType &weight) const
apply projection
Definition: vtxprojection.hh:252
RType ::RangeFieldType RangeFieldType
Definition: vtxprojection.hh:244
DomainType ::RangeFieldType DomainFieldType
Definition: vtxprojection.hh:243
RType ::DiscreteFunctionSpaceType ::GridPartType GridPartType
Definition: vtxprojection.hh:246
Dune::FieldTraits< RangeFieldType >::real_type RealType
Definition: vtxprojection.hh:245
actual interface class for quadratures
Definition: quadrature.hh:405
Definition: common/localinterpolation.hh:22