1 #ifndef DUNE_FEM_SPACE_COMMON_INTERPOLATE_HH
2 #define DUNE_FEM_SPACE_COMMON_INTERPOLATE_HH
3 
4 #include <algorithm>
5 #include <type_traits>
6 #include <utility>
7 
8 #include <dune/common/typetraits.hh>
9 
10 #include <dune/grid/common/partitionset.hh>
11 #include <dune/grid/common/rangegenerators.hh>
12 
13 #include <dune/fem/common/bindguard.hh>
14 
15 #include <dune/fem/function/common/discretefunction.hh>
16 #include <dune/fem/function/common/gridfunctionadapter.hh>
17 #include <dune/fem/function/common/localcontribution.hh>
18 #include <dune/fem/function/localfunction/const.hh>
19 #include <dune/fem/space/common/capabilities.hh>
20 #include <dune/fem/space/common/localinterpolation.hh>
21 
22 namespace Dune
23 {
24 
25   namespace Fem
26   {
27 
28     // interpolate
29     // -----------
30 
31     /**
32      * \function interpolate
33      * \ingroup  DiscreteFunctionSpace
34      * \brief    perform native interpolation of a discrete function space
35      *
36      * By definition of its degrees of freedom, each discrete function space
37      * has a native interpolation, which can be computed very quickly.
38      *
39      * For example, the native interpolation of a Lagrange discrete function
40      * space is the evaluation in its Lagrange points.
41      * An orthonormal DG space would instead perform an \f$L^2\f$-Projection.
42      *
43      * The actual implementation must locally be provided by the discrete
44      * function space through the method
45      * \code
46      * template< class LocalFunction, class LocalDofVector >
47      * void interpolate ( const LocalFunction &f, LocalDofVector &dofs ) const;
48      * \endcode
49      *
50      * \param[in]   u  grid function to interpolate
51      * \param[out]  v  discrete function to represent the interpolation
52      */
53     template< class GridFunction, class DiscreteFunction >
interpolate(const GridFunction & u,DiscreteFunction & v)54     static inline void interpolate ( const GridFunction &u, DiscreteFunction &v )
55     {
56       // just call interpolate for the all partition
57       interpolate( u, v, Partitions::all );
58     }
59 
60     template< class Function, class DiscreteFunction, unsigned int partitions >
61     static inline std::enable_if_t< !std::is_convertible< Function, HasLocalFunction >::value >
interpolate(const Function & u,DiscreteFunction & v,PartitionSet<partitions> ps)62     interpolate ( const Function &u, DiscreteFunction &v, PartitionSet< partitions > ps )
63     {
64       typedef typename DiscreteFunction :: DiscreteFunctionSpaceType :: GridPartType  GridPartType;
65       typedef GridFunctionAdapter< Function, GridPartType > GridFunctionType;
66 
67       GridFunctionType uGrid( "uGrid", u, v.space().gridPart() );
68 
69       interpolate( uGrid, v, ps );
70     }
71 
72     template< class GridFunction, class DiscreteFunction, unsigned int partitions >
73     static inline std::enable_if_t< std::is_convertible< GridFunction, HasLocalFunction >::value && Capabilities::hasInterpolation< typename DiscreteFunction::DiscreteFunctionSpaceType >::v >
interpolate(const GridFunction & u,DiscreteFunction & v,PartitionSet<partitions> ps)74     interpolate ( const GridFunction &u, DiscreteFunction &v, PartitionSet< partitions > ps )
75     {
76       ConstLocalFunction< GridFunction > uLocal( u );
77       LocalContribution< DiscreteFunction, Assembly::Set > vLocal( v );
78       LocalInterpolation< typename DiscreteFunction::DiscreteFunctionSpaceType >
79         interpolation( v.space() );
80 
81       // iterate over selected partition
82       for( const auto entity : elements( v.gridPart(), ps ) )
83       {
84         // initialize u to entity
85         auto uGuard = bindGuard( uLocal, entity );
86 
87         // bind v to entity
88         auto vGuard = bindGuard( vLocal, entity );
89 
90         // bind interpolation to entity
91         auto iGuard = bindGuard( interpolation, entity );
92 
93         // perform local interpolation
94         interpolation( uLocal, vLocal );
95       }
96     }
97 
98 
99 
100     namespace Impl
101     {
102 
103       template< class Entity, class FunctionSpace, class Weight >
104       struct WeightLocalFunction
105       {
106         typedef Entity EntityType;
107         typedef FunctionSpace FunctionSpaceType;
108 
109         typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
110         typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
111 
112         typedef typename FunctionSpaceType::DomainType DomainType;
113         typedef typename FunctionSpaceType::RangeType RangeType;
114         typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
115         typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
116 
117         typedef typename EntityType::Geometry::LocalCoordinate LocalCoordinateType;
118 
119         static constexpr int dimDomain = FunctionSpaceType::dimDomain;
120         static constexpr int dimRange = FunctionSpaceType::dimRange;
121 
WeightLocalFunctionDune::Fem::Impl::WeightLocalFunction122         explicit WeightLocalFunction ( Weight &weight, int order ) : weight_( weight ), order_(order) {}
123 
bindDune::Fem::Impl::WeightLocalFunction124         void bind ( const EntityType &entity ) { entity_ = entity; weight_.setEntity( entity ); }
unbindDune::Fem::Impl::WeightLocalFunction125         void unbind () {}
126 
entityDune::Fem::Impl::WeightLocalFunction127         const EntityType entity() const { return entity_; }
128 
129         template< class Point >
evaluateDune::Fem::Impl::WeightLocalFunction130         void evaluate ( const Point &x, RangeType &value ) const
131         {
132           const RangeFieldType weight = weight_( coordinate( x ) );
133           for( int i = 0; i < dimRange; ++i )
134             value[ i ] = weight;
135         }
136 
137         template< class Quadrature, class Values >
evaluateQuadratureDune::Fem::Impl::WeightLocalFunction138         auto evaluateQuadrature ( const Quadrature &quadrature, Values &values ) const
139           -> std::enable_if_t< std::is_same< decltype( values[ 0 ] ), RangeType & >::value >
140         {
141           for( const auto &qp : quadrature )
142             evaluate( qp, values[ qp.index() ] );
143         }
144 
orderDune::Fem::Impl::WeightLocalFunction145         int order() const
146         { return order_; }
147       private:
148         Weight &weight_;
149         int order_;
150         Entity entity_;
151       };
152 
153     } // namespace Impl
154 
155 
156 
157     // interpolate with weights
158     // ------------------------
159 
160     template< class GridFunction, class DiscreteFunction, class Weight >
interpolate(const GridFunction & u,DiscreteFunction & v,Weight && weight)161     inline static auto interpolate ( const GridFunction &u, DiscreteFunction &v, Weight &&weight )
162       -> std::enable_if_t< !std::is_const< std::remove_reference_t< Weight > >::value  >
163     {
164       DiscreteFunction w( v );
165       interpolate( u, v, std::forward< Weight >( weight ), w );
166     }
167 
168     template< class GridFunction, class DiscreteFunction, class Weight >
interpolate(const GridFunction & u,DiscreteFunction & v,Weight && weight)169     inline static auto interpolate ( const GridFunction &u, DiscreteFunction &v, Weight &&weight )
170       -> std::enable_if_t< !std::is_const< std::remove_reference_t< Weight > >::value && !std::is_base_of< HasLocalFunction, GridFunction >::value >
171     {
172       interpolate( gridFunctionAdapter( u, v.gridPart() ), v, std::forward< Weight >( weight ) );
173     }
174 
175     template< class GridFunction, class DiscreteFunction, class Weight >
interpolate(const GridFunction & u,DiscreteFunction & v,Weight && weight,DiscreteFunction & w)176     inline static auto interpolate ( const GridFunction &u, DiscreteFunction &v, Weight &&weight, DiscreteFunction &w )
177       -> std::enable_if_t< std::is_base_of< HasLocalFunction, GridFunction >::value && Capabilities::hasInterpolation< typename DiscreteFunction::DiscreteFunctionSpaceType >::v,
178                            void_t< decltype( std::declval< Weight >().setEntity( std::declval< const typename DiscreteFunction::DiscreteFunctionSpaceType::EntityType & >() ) ) > >
179     {
180       typedef typename DiscreteFunction::DiscreteFunctionSpaceType::EntityType EntityType;
181 
182       const auto &space = w.space();
183       Impl::WeightLocalFunction< EntityType, std::remove_reference_t< typename DiscreteFunction::FunctionSpaceType >, Weight > localWeight( weight, w.order() );
184       LocalInterpolation< typename DiscreteFunction::DiscreteFunctionSpaceType > interpolation( space );
185       interpolate( u, v, [ &interpolation, &localWeight ] ( const EntityType &entity, AddLocalContribution< DiscreteFunction > &w ) {
186           auto weightGuard = bindGuard( localWeight, entity );
187           auto iGuard = bindGuard( interpolation, entity );
188           interpolation( localWeight, w );
189         }, w );
190     }
191 
192     template< class GridFunction, class DiscreteFunction, class Weight >
interpolate(const GridFunction & u,DiscreteFunction & v,Weight && weight,DiscreteFunction & w)193     inline static auto interpolate ( const GridFunction &u, DiscreteFunction &v, Weight &&weight, DiscreteFunction &w )
194       -> std::enable_if_t< std::is_base_of< HasLocalFunction, GridFunction >::value && Capabilities::hasInterpolation< typename DiscreteFunction::DiscreteFunctionSpaceType >::v,
195                            void_t< decltype( std::declval< Weight >()( std::declval< typename DiscreteFunction::DiscreteFunctionSpaceType::EntityType & >(), std::declval< AddLocalContribution< DiscreteFunction > & >() ) ) > >
196     {
197       typedef typename DiscreteFunction::DofType DofType;
198 
199       v.clear();
200       w.clear();
201 
202       {
203         ConstLocalFunction< GridFunction > uLocal( u );
204         AddLocalContribution< DiscreteFunction > vLocal( v ), wLocal( w );
205         LocalInterpolation< typename DiscreteFunction::DiscreteFunctionSpaceType >
206           interpolation( v.space() );
207 
208         for( const auto &entity : v.space() )
209         {
210           auto uGuard = bindGuard( uLocal, entity );
211           auto vGuard = bindGuard( vLocal, entity );
212           auto wGuard = bindGuard( wLocal, entity );
213           auto iGuard = bindGuard( interpolation, entity );
214 
215           // interpolate u and store in v
216           interpolation( uLocal, vLocal );
217 
218           // evaluate DoF-wise weight
219           weight( entity, wLocal );
220 
221           // multiply interpolated values by weight
222           std::transform( vLocal.begin(), vLocal.end(), wLocal.begin(), vLocal.begin(), std::multiplies< DofType >() );
223         }
224       } // ensure the local contributions go out of scope, here (communication)
225 
226       std::transform( v.dbegin(), v.dend(), w.dbegin(), v.dbegin(), [] ( DofType v, DofType w ) {
227           // weights are non-negative, so cancellation cannot occur
228           return (w > DofType( 0 ) ? v / w : v);
229         } );
230     }
231 
232   } // namespace Fem
233 
234 } // namespace Dune
235 
236 #endif // #ifndef DUNE_FEM_SPACE_COMMON_INTERPOLATE_HH
237