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