1 #ifndef DUNE_FEM_INTEGRAL_HH 2 #define DUNE_FEM_INTEGRAL_HH 3 4 #include <type_traits> 5 6 #include <dune/grid/common/rangegenerators.hh> 7 8 #include <dune/fem/quadrature/cachingquadrature.hh> 9 #include <dune/fem/quadrature/integrator.hh> 10 11 #include <dune/fem/function/common/gridfunctionadapter.hh> 12 13 #include <dune/common/bartonnackmanifcheck.hh> 14 #include <dune/fem/misc/bartonnackmaninterface.hh> 15 16 namespace Dune 17 { 18 19 namespace Fem 20 { 21 // IntegralBase 22 // ---------- 23 24 template< class GridPart, class NormImplementation > 25 class IntegralBase 26 : public BartonNackmanInterface< IntegralBase< GridPart, NormImplementation >, 27 NormImplementation > 28 { 29 typedef BartonNackmanInterface< IntegralBase< GridPart, NormImplementation >, 30 NormImplementation > BaseType ; 31 typedef IntegralBase< GridPart, NormImplementation > ThisType; 32 33 public: 34 typedef GridPart GridPartType; 35 36 protected: 37 using BaseType :: asImp ; 38 39 typedef typename GridPartType::template Codim< 0 >::EntityType EntityType; 40 41 template <bool uDiscrete, bool vDiscrete> 42 struct ForEachCaller 43 { 44 template <class UDiscreteFunctionType, class VDiscreteFunctionType, class ReturnType, class PartitionSet> forEachDune::Fem::IntegralBase::ForEachCaller45 static ReturnType forEach ( const ThisType &norm, const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, 46 const ReturnType &initialValue, 47 const PartitionSet& partitionSet, 48 unsigned int order ) 49 { 50 static_assert( uDiscrete && vDiscrete, "Distance can only be calculated between GridFunctions" ); 51 52 ReturnType sum( 0 ); 53 ConstLocalFunction< UDiscreteFunctionType > uLocal( u ); 54 ConstLocalFunction< VDiscreteFunctionType > vLocal( v ); 55 for( const EntityType &entity : elements( norm.gridPart_, partitionSet ) ) 56 { 57 uLocal.bind( entity ); 58 vLocal.bind( entity ); 59 const unsigned int uOrder = uLocal.order(); 60 const unsigned int vOrder = vLocal.order(); 61 const unsigned int orderLocal = (order == 0 ? 2*std::max( uOrder, vOrder ) : order); 62 norm.distanceLocal( entity, orderLocal, uLocal, vLocal, sum ); 63 uLocal.unbind(); 64 vLocal.unbind(); 65 } 66 return sum; 67 } 68 }; 69 70 // this specialization creates a grid function adapter 71 template <bool vDiscrete> 72 struct ForEachCaller<false, vDiscrete> 73 { 74 template <class F, 75 class VDiscreteFunctionType, 76 class ReturnType, 77 class PartitionSet> 78 static forEachDune::Fem::IntegralBase::ForEachCaller79 ReturnType forEach ( const ThisType& norm, 80 const F& f, const VDiscreteFunctionType& v, 81 const ReturnType& initialValue, 82 const PartitionSet& partitionSet, 83 const unsigned int order ) 84 { 85 typedef GridFunctionAdapter< F, GridPartType> GridFunction ; 86 GridFunction u( "Integral::adapter" , f , v.space().gridPart(), v.space().order() ); 87 88 return ForEachCaller< true, vDiscrete > :: 89 forEach( norm, u, v, initialValue, partitionSet, order ); 90 } 91 }; 92 93 // this specialization simply switches arguments 94 template <bool uDiscrete> 95 struct ForEachCaller<uDiscrete, false> 96 { 97 template <class UDiscreteFunctionType, 98 class F, 99 class ReturnType, 100 class PartitionSet> 101 static forEachDune::Fem::IntegralBase::ForEachCaller102 ReturnType forEach ( const ThisType& norm, 103 const UDiscreteFunctionType& u, 104 const F& f, 105 const ReturnType& initialValue, 106 const PartitionSet& partitionSet, 107 const unsigned int order ) 108 { 109 return ForEachCaller< false, uDiscrete > :: 110 forEach( norm, f, u, initialValue, partitionSet, order ); 111 } 112 }; 113 114 template< class DiscreteFunctionType, class ReturnType, class PartitionSet > 115 ReturnType forEach ( const DiscreteFunctionType &u, const ReturnType &initialValue, 116 const PartitionSet& partitionSet, 117 unsigned int order = 0 ) const 118 { 119 static_assert( (std::is_base_of<Fem::HasLocalFunction, DiscreteFunctionType>::value), 120 "Norm only implemented for quantities implementing a local function!" ); 121 122 ReturnType sum( 0 ); 123 ConstLocalFunction< DiscreteFunctionType > uLocal( u ); 124 for( const EntityType &entity : elements( gridPart_, partitionSet ) ) 125 { 126 uLocal.bind( entity ); 127 const unsigned int orderLocal = (order == 0 ? 2*uLocal.order() : order); 128 normLocal( entity, orderLocal, uLocal, sum ); 129 uLocal.unbind(); 130 } 131 return sum; 132 } 133 134 template< class UDiscreteFunctionType, class VDiscreteFunctionType, class ReturnType, class PartitionSet > 135 ReturnType forEach ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, 136 const ReturnType &initialValue, const PartitionSet& partitionSet, 137 unsigned int order = 0 ) const 138 { 139 enum { uDiscrete = std::is_convertible<UDiscreteFunctionType, HasLocalFunction>::value }; 140 enum { vDiscrete = std::is_convertible<VDiscreteFunctionType, HasLocalFunction>::value }; 141 142 // call forEach depending on which argument is a grid function, 143 // i.e. has a local function 144 return ForEachCaller< uDiscrete, vDiscrete > :: 145 forEach( *this, u, v, initialValue, partitionSet, order ); 146 } 147 148 public: IntegralBase(const GridPartType & gridPart)149 explicit IntegralBase ( const GridPartType &gridPart ) 150 : gridPart_( gridPart ) 151 {} 152 153 protected: 154 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType > distanceLocal(const EntityType & entity,unsigned int order,const ULocalFunctionType & uLocal,const VLocalFunctionType & vLocal,ReturnType & sum) const155 void distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const 156 { 157 CHECK_AND_CALL_INTERFACE_IMPLEMENTATION( asImp().distanceLocal( entity, order, uLocal, vLocal, sum ) ); 158 } 159 160 template< class LocalFunctionType, class ReturnType > normLocal(const EntityType & entity,unsigned int order,const LocalFunctionType & uLocal,ReturnType & sum) const161 void normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const 162 { 163 CHECK_AND_CALL_INTERFACE_IMPLEMENTATION( asImp().normLocal( entity, order, uLocal, sum ) ); 164 } 165 gridPart() const166 const GridPartType &gridPart () const { return gridPart_; } 167 comm() const168 const typename GridPartType::CollectiveCommunicationType& comm () const 169 { 170 return gridPart().comm(); 171 } 172 checkCommunicateFlag(bool communicate) const173 bool checkCommunicateFlag( bool communicate ) const 174 { 175 #ifndef NDEBUG 176 bool commMax = communicate; 177 assert( communicate == comm().max( commMax ) ); 178 #endif 179 return communicate; 180 } 181 182 private: 183 const GridPartType &gridPart_; 184 }; 185 186 // Integral 187 // ------ 188 189 template< class GridPart > 190 class Integral : public IntegralBase< GridPart, Integral< GridPart > > 191 { 192 typedef IntegralBase< GridPart, Integral< GridPart > > BaseType ; 193 typedef Integral< GridPart > ThisType; 194 195 public: 196 typedef GridPart GridPartType; 197 198 using BaseType :: gridPart ; 199 using BaseType :: comm ; 200 201 protected: 202 203 template< class UFunction, class VFunction > 204 struct FunctionDistance; 205 206 typedef typename BaseType::EntityType EntityType; 207 typedef CachingQuadrature< GridPartType, 0 > QuadratureType; 208 209 const unsigned int order_; 210 const bool communicate_; 211 public: 212 /** \brief constructor 213 * \param gridPart specific gridPart for selection of entities 214 * \param order order of integration quadrature (default = 2*space.order()) 215 * \param communicate if true global (over all ranks) norm is computed (default = true) 216 */ 217 explicit Integral ( const GridPartType &gridPart, 218 const unsigned int order = 0, 219 const bool communicate = true ); 220 221 222 //! || u ||_L1 on given set of entities (partition set) 223 template< class DiscreteFunctionType, class PartitionSet > 224 typename DiscreteFunctionType::RangeType 225 norm ( const DiscreteFunctionType &u, const PartitionSet& partitionSet ) const; 226 227 //! || u ||_L1 on interior partition entities 228 template< class DiscreteFunctionType > 229 typename DiscreteFunctionType::RangeType norm(const DiscreteFunctionType & u) const230 norm ( const DiscreteFunctionType &u ) const 231 { 232 return norm( u, Partitions::interior ); 233 } 234 235 //! || u - v ||_L2 on given set of entities (partition set) 236 template< class UDiscreteFunctionType, class VDiscreteFunctionType, class PartitionSet > 237 typename UDiscreteFunctionType::RangeType 238 distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const PartitionSet& partitionSet ) const; 239 240 //! || u - v ||_L2 on interior partition entities 241 template< class UDiscreteFunctionType, class VDiscreteFunctionType > 242 typename UDiscreteFunctionType::RangeType distance(const UDiscreteFunctionType & u,const VDiscreteFunctionType & v) const243 distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v ) const 244 { 245 return distance( u, v, Partitions::interior ); 246 } 247 248 template< class LocalFunctionType, class ReturnType > 249 void normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const; 250 251 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType > 252 void distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const; 253 }; 254 255 256 257 // Implementation of Integral 258 // ------------------------ 259 260 template< class GridPart > Integral(const GridPartType & gridPart,const unsigned int order,const bool communicate)261 inline Integral< GridPart >::Integral ( const GridPartType &gridPart, const unsigned int order, const bool communicate ) 262 : BaseType( gridPart ), 263 order_( order ), 264 communicate_( BaseType::checkCommunicateFlag( communicate ) ) 265 {} 266 267 268 template< class GridPart > 269 template< class DiscreteFunctionType, class PartitionSet > 270 inline typename DiscreteFunctionType::RangeType norm(const DiscreteFunctionType & u,const PartitionSet & partitionSet) const271 Integral< GridPart >::norm ( const DiscreteFunctionType &u, const PartitionSet& partitionSet ) const 272 { 273 typedef typename DiscreteFunctionType::RangeType RangeType; 274 typedef RangeType ReturnType ; 275 276 // calculate integral over each element 277 ReturnType sum = BaseType :: forEach( u, ReturnType(0), partitionSet, order_ ); 278 279 // communicate_ indicates global norm 280 if( communicate_ ) 281 { 282 sum = comm().sum( sum ); 283 } 284 285 return sum; 286 } 287 288 289 template< class GridPart > 290 template< class UDiscreteFunctionType, class VDiscreteFunctionType, class PartitionSet > 291 inline typename UDiscreteFunctionType::RangeType 292 Integral< GridPart > distance(const UDiscreteFunctionType & u,const VDiscreteFunctionType & v,const PartitionSet & partitionSet) const293 ::distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const PartitionSet& partitionSet ) const 294 { 295 typedef typename UDiscreteFunctionType::RangeType RangeType; 296 typedef RangeType ReturnType ; 297 298 // calculate integral over each element 299 ReturnType sum = BaseType :: forEach( u, v, ReturnType(0), partitionSet, order_ ); 300 301 // communicate_ indicates global norm 302 if( communicate_ ) 303 { 304 sum = comm().sum( sum ); 305 } 306 307 return sum; 308 } 309 310 template< class GridPart > 311 template< class LocalFunctionType, class ReturnType > 312 inline void normLocal(const EntityType & entity,unsigned int order,const LocalFunctionType & uLocal,ReturnType & sum) const313 Integral< GridPart >::normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const 314 { 315 Integrator< QuadratureType > integrator( order ); 316 317 integrator.integrateAdd( entity, uLocal, sum ); 318 } 319 320 template< class GridPart > 321 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType > 322 inline void distanceLocal(const EntityType & entity,unsigned int order,const ULocalFunctionType & uLocal,const VLocalFunctionType & vLocal,ReturnType & sum) const323 Integral< GridPart >::distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const 324 { 325 Integrator< QuadratureType > integrator( order ); 326 327 typedef FunctionDistance< ULocalFunctionType, VLocalFunctionType > LocalDistanceType; 328 329 LocalDistanceType dist( uLocal, vLocal ); 330 331 integrator.integrateAdd( entity, dist, sum ); 332 } 333 334 template< class GridPart > 335 template< class UFunction, class VFunction > 336 struct Integral< GridPart >::FunctionDistance 337 { 338 typedef UFunction UFunctionType; 339 typedef VFunction VFunctionType; 340 341 typedef typename UFunctionType::RangeFieldType RangeFieldType; 342 typedef typename UFunctionType::RangeType RangeType; 343 typedef typename UFunctionType::JacobianRangeType JacobianRangeType; 344 FunctionDistanceDune::Fem::Integral::FunctionDistance345 FunctionDistance ( const UFunctionType &u, const VFunctionType &v ) 346 : u_( u ), v_( v ) 347 {} 348 349 template< class Point > evaluateDune::Fem::Integral::FunctionDistance350 void evaluate ( const Point &x, RangeType &ret ) const 351 { 352 RangeType phi; 353 u_.evaluate( x, ret ); 354 v_.evaluate( x, phi ); 355 ret -= phi; 356 } 357 358 template< class Point > jacobianDune::Fem::Integral::FunctionDistance359 void jacobian ( const Point &x, JacobianRangeType &ret ) const 360 { 361 JacobianRangeType phi; 362 u_.jacobian( x, ret ); 363 v_.jacobian( x, phi ); 364 ret -= phi; 365 } 366 367 private: 368 const UFunctionType &u_; 369 const VFunctionType &v_; 370 }; 371 372 } // namespace Fem 373 374 } // namespace Dune 375 376 #endif // #ifndef DUNE_FEM_INTEGRAL_HH 377