dune-fem 2.8.0
Loading...
Searching...
No Matches
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
19namespace Dune
20{
21 namespace Fem
22 {
23
24 template <class GridPartType>
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
const EntityType & entity() const
Definition: vtxprojection.hh:81
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
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::DiscreteFunctionSpaceType::GridPartType GridPartType
Definition: vtxprojection.hh:246
RType RangeType
Definition: vtxprojection.hh:242
VtxProjection()
Constructor.
Definition: vtxprojection.hh:248
DType DomainType
Definition: vtxprojection.hh:241
RType::RangeFieldType RangeFieldType
Definition: vtxprojection.hh:244
Dune::FieldTraits< RangeFieldType >::real_type RealType
Definition: vtxprojection.hh:245
void operator()(const DomainType &f, RangeType &discFunc, WeightType &weight) const
apply projection
Definition: vtxprojection.hh:252
DomainType::RangeFieldType DomainFieldType
Definition: vtxprojection.hh:243
actual interface class for quadratures
Definition: quadrature.hh:405
Definition: common/localinterpolation.hh:22