1 #ifndef DUNE_FEM_LPNORM_HH 2 #define DUNE_FEM_LPNORM_HH 3 4 #include <dune/fem/quadrature/integrator.hh> 5 6 #include "domainintegral.hh" 7 8 namespace Dune 9 { 10 11 namespace Fem 12 { 13 14 // TODO weighte LP norm might be adapted later 15 // LPNorm 16 // 17 // !!!!! It is not cleared which quadrature order have to be applied for p > 2!!!! 18 // !!!!! For p = 1 this norm does not work !!!! 19 // ------ 20 21 22 //! Quadrature Order Interface 23 struct OrderCalculatorInterface 24 { 25 virtual int operator() (const double p)=0; 26 }; 27 28 //! default Quadrature Order Calculator 29 // can be re-implemented in order to use 30 // a different type of calculation 31 // which can be sepcified in the second template argument of LPNorm 32 struct DefaultOrderCalculator : public OrderCalculatorInterface 33 { operator ()Dune::Fem::DefaultOrderCalculator34 int operator() (const double p) 35 { 36 int ret=0; 37 const double q = p / (p-1); 38 double max = std::max(p,q); 39 assert(max < std::numeric_limits<int>::max()/2. ); 40 ret = max +1; 41 return ret; 42 } 43 }; 44 45 template< class GridPart, class OrderCalculator = DefaultOrderCalculator > 46 class LPNorm : public IntegralBase < GridPart, LPNorm< GridPart, OrderCalculator > > 47 { 48 typedef LPNorm< GridPart, OrderCalculator > ThisType; 49 typedef IntegralBase< GridPart, ThisType > BaseType ; 50 51 public: 52 typedef GridPart GridPartType; 53 54 using BaseType::gridPart; 55 using BaseType::comm; 56 57 protected: 58 template< class Function > 59 struct FunctionMultiplicator; 60 61 template< class UFunction, class VFunction > 62 struct FunctionDistance; 63 64 typedef typename BaseType::EntityType EntityType; 65 typedef CachingQuadrature< GridPartType, 0 > QuadratureType; 66 typedef Integrator< QuadratureType > IntegratorType; 67 68 public: 69 /** \brief constructor 70 * \param gridPart specific gridPart for selection of entities 71 * \param p p in Lp-norm 72 * \param communicate if true global (over all ranks) norm is computed (default = true) 73 */ 74 explicit LPNorm ( const GridPartType &gridPart, const double p, const bool communicate = true ); 75 76 //! || u ||_Lp on given set of entities (partition set) 77 template< class DiscreteFunctionType, class PartitionSet > 78 typename Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type 79 norm ( const DiscreteFunctionType &u, const PartitionSet& partitionSet ) const; 80 81 //! || u ||_Lp on interior partition entities 82 template< class DiscreteFunctionType > 83 typename Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type norm(const DiscreteFunctionType & u) const84 norm ( const DiscreteFunctionType &u ) const 85 { 86 return norm( u, Partitions::interior ); 87 } 88 89 //! || u - v ||_Lp on given set of entities (partition set) 90 template< class UDiscreteFunctionType, class VDiscreteFunctionType, class PartitionSet > 91 typename Dune::FieldTraits< typename UDiscreteFunctionType::RangeFieldType >::real_type 92 distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const PartitionSet& partitionSet ) const; 93 94 //! || u - v ||_Lp on interior partition entities 95 template< class UDiscreteFunctionType, class VDiscreteFunctionType > 96 typename Dune::FieldTraits< typename UDiscreteFunctionType::RangeFieldType >::real_type distance(const UDiscreteFunctionType & u,const VDiscreteFunctionType & v) const97 distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v ) const 98 { 99 return distance( u, v, Partitions::interior ); 100 } 101 102 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType > 103 void distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const; 104 105 template< class LocalFunctionType, class ReturnType > 106 void normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const; 107 108 int order ( const int spaceOrder ) const ; 109 110 protected: 111 double p_ ; 112 OrderCalculator *orderCalc_; 113 const bool communicate_; 114 }; 115 116 117 118 119 // WeightedLPNorm 120 // -------------- 121 122 template< class WeightFunction, class OrderCalculator = DefaultOrderCalculator > 123 class WeightedLPNorm 124 : public LPNorm< typename WeightFunction::DiscreteFunctionSpaceType::GridPartType, 125 OrderCalculator > 126 { 127 typedef WeightedLPNorm< WeightFunction, OrderCalculator > ThisType; 128 typedef LPNorm< typename WeightFunction::DiscreteFunctionSpaceType::GridPartType, OrderCalculator> BaseType; 129 130 public: 131 typedef WeightFunction WeightFunctionType; 132 133 typedef typename WeightFunctionType::DiscreteFunctionSpaceType WeightFunctionSpaceType; 134 typedef typename WeightFunctionSpaceType::GridPartType GridPartType; 135 136 protected: 137 template< class Function > 138 struct WeightedFunctionMultiplicator; 139 140 typedef ConstLocalFunction< WeightFunctionType > LocalWeightFunctionType; 141 typedef typename WeightFunctionType::RangeType WeightType; 142 143 typedef typename BaseType::EntityType EntityType; 144 typedef typename BaseType::IntegratorType IntegratorType; 145 146 using BaseType::gridPart; 147 using BaseType::comm; 148 149 public: 150 using BaseType::norm; 151 using BaseType::distance; 152 153 explicit WeightedLPNorm ( const WeightFunctionType &weightFunction, const double p, const bool communicate = true ); 154 155 template< class LocalFunctionType, class ReturnType > 156 void normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const; 157 158 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType > 159 void distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const; 160 161 private: 162 LocalWeightFunctionType wfLocal_; 163 const double p_; 164 }; 165 166 167 // Implementation of LPNorm 168 // ------------------------ 169 170 template< class GridPart, class OrderCalculator > LPNorm(const GridPartType & gridPart,const double p,const bool communicate)171 inline LPNorm< GridPart, OrderCalculator >::LPNorm ( const GridPartType &gridPart, const double p, const bool communicate ) 172 : BaseType( gridPart ), 173 p_(p), 174 orderCalc_( new OrderCalculator() ), 175 communicate_( BaseType::checkCommunicateFlag( communicate ) ) 176 { 177 } 178 179 template< class GridPart, class OrderCalculator> order(const int spaceOrder) const180 inline int LPNorm< GridPart, OrderCalculator>::order(const int spaceOrder) const 181 { 182 return spaceOrder * orderCalc_->operator() (p_); 183 } 184 185 186 template< class GridPart, class OrderCalculator> 187 template< class DiscreteFunctionType, class PartitionSet > 188 inline typename Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type norm(const DiscreteFunctionType & u,const PartitionSet & partitionSet) const189 LPNorm< GridPart, OrderCalculator >::norm ( const DiscreteFunctionType &u, const PartitionSet& partitionSet ) const 190 { 191 typedef typename DiscreteFunctionType::RangeFieldType RangeFieldType; 192 typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType; 193 typedef FieldVector< RealType, 1 > ReturnType ; 194 195 // calculate integral over each element 196 ReturnType sum = BaseType :: forEach( u, ReturnType(0), partitionSet ); 197 198 // communicate_ indicates global norm 199 if( communicate_ ) 200 { 201 sum[ 0 ] = comm().sum( sum[ 0 ] ); 202 } 203 204 // return result 205 return std::pow ( sum[ 0 ], (1.0 / p_) ); 206 } 207 208 template< class GridPart, class OrderCalculator > 209 template< class UDiscreteFunctionType, class VDiscreteFunctionType, class PartitionSet > 210 inline typename Dune::FieldTraits< typename UDiscreteFunctionType::RangeFieldType >::real_type 211 LPNorm< GridPart, OrderCalculator > distance(const UDiscreteFunctionType & u,const VDiscreteFunctionType & v,const PartitionSet & partitionSet) const212 ::distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v, const PartitionSet& partitionSet ) const 213 { 214 typedef typename UDiscreteFunctionType::RangeFieldType RangeFieldType; 215 typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType; 216 typedef FieldVector< RealType, 1 > ReturnType ; 217 218 // calculate integral over each element 219 ReturnType sum = BaseType :: forEach( u, v, ReturnType(0), partitionSet ); 220 221 // communicate_ indicates global norm 222 if( communicate_ ) 223 { 224 sum[ 0 ] = comm().sum( sum[ 0 ] ); 225 } 226 227 // return result 228 return std::pow( sum[ 0 ], (1.0/p_) ); 229 } 230 231 template< class GridPart, class OrderCalculator > 232 template< class LocalFunctionType, class ReturnType > 233 inline void normLocal(const EntityType & entity,unsigned int order,const LocalFunctionType & uLocal,ReturnType & sum) const234 LPNorm< GridPart, OrderCalculator >::normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const 235 { 236 IntegratorType integrator( order ); 237 238 FunctionMultiplicator< LocalFunctionType > uLocalp( uLocal, p_ ); 239 240 integrator.integrateAdd( entity, uLocalp, sum ); 241 } 242 243 template< class GridPart, class OrderCalculator > 244 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType > 245 inline void distanceLocal(const EntityType & entity,unsigned int order,const ULocalFunctionType & uLocal,const VLocalFunctionType & vLocal,ReturnType & sum) const246 LPNorm< GridPart, OrderCalculator >::distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const 247 { 248 typedef FunctionDistance< ULocalFunctionType, VLocalFunctionType > LocalDistanceType; 249 250 IntegratorType integrator( order ); 251 252 LocalDistanceType dist( uLocal, vLocal ); 253 FunctionMultiplicator< LocalDistanceType > distp( dist, p_ ); 254 255 integrator.integrateAdd( entity, distp, sum ); 256 } 257 258 259 template< class GridPart, class OrderCalculator > 260 template< class Function > 261 struct LPNorm< GridPart, OrderCalculator >::FunctionMultiplicator 262 { 263 typedef Function FunctionType; 264 265 typedef typename FunctionType::RangeFieldType RangeFieldType; 266 typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType; 267 typedef FieldVector< RealType, 1 > RangeType ; 268 FunctionMultiplicatorDune::Fem::LPNorm::FunctionMultiplicator269 explicit FunctionMultiplicator ( const FunctionType &function, double p ) 270 : function_( function ), 271 p_(p) 272 {} 273 274 template< class Point > evaluateDune::Fem::LPNorm::FunctionMultiplicator275 void evaluate ( const Point &x, RangeType &ret ) const 276 { 277 typename FunctionType::RangeType phi; 278 function_.evaluate( x, phi ); 279 ret = std :: pow ( phi.two_norm(), p_); 280 } 281 282 private: 283 const FunctionType &function_; 284 double p_; 285 }; 286 287 288 template< class GridPart, class OrderCalculator > 289 template< class UFunction, class VFunction > 290 struct LPNorm< GridPart, OrderCalculator >::FunctionDistance 291 { 292 typedef UFunction UFunctionType; 293 typedef VFunction VFunctionType; 294 295 typedef typename UFunctionType::RangeFieldType RangeFieldType; 296 typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType; 297 typedef typename UFunctionType::RangeType RangeType; 298 typedef typename UFunctionType::JacobianRangeType JacobianRangeType; 299 FunctionDistanceDune::Fem::LPNorm::FunctionDistance300 FunctionDistance ( const UFunctionType &u, const VFunctionType &v ) 301 : u_( u ), v_( v ) 302 {} 303 304 template< class Point > evaluateDune::Fem::LPNorm::FunctionDistance305 void evaluate ( const Point &x, RangeType &ret ) const 306 { 307 RangeType phi; 308 u_.evaluate( x, ret ); 309 v_.evaluate( x, phi ); 310 ret -= phi; 311 } 312 313 template< class Point > jacobianDune::Fem::LPNorm::FunctionDistance314 void jacobian ( const Point &x, JacobianRangeType &ret ) const 315 { 316 JacobianRangeType phi; 317 u_.jacobian( x, ret ); 318 v_.jacobian( x, phi ); 319 ret -= phi; 320 } 321 322 private: 323 const UFunctionType &u_; 324 const VFunctionType &v_; 325 }; 326 327 328 // Implementation of WeightedL2Norm 329 // -------------------------------- 330 331 template< class WeightFunction, class OrderCalculator > 332 inline WeightedLPNorm< WeightFunction, OrderCalculator > WeightedLPNorm(const WeightFunctionType & weightFunction,double p,const bool communicate)333 ::WeightedLPNorm ( const WeightFunctionType &weightFunction, double p, const bool communicate ) 334 : BaseType( weightFunction.space().gridPart(), p, communicate ), 335 wfLocal_( weightFunction ), 336 p_(p) 337 { 338 static_assert( (WeightFunctionSpaceType::dimRange == 1), 339 "Weight function must be scalar." ); 340 } 341 342 343 template< class WeightFunction, class OrderCalculator > 344 template< class LocalFunctionType, class ReturnType > 345 inline void normLocal(const EntityType & entity,unsigned int order,const LocalFunctionType & uLocal,ReturnType & sum) const346 WeightedLPNorm< WeightFunction, OrderCalculator >::normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const 347 { 348 // !!!! order !!!! 349 IntegratorType integrator( order ); 350 351 wfLocal_.bind( entity ); 352 WeightedFunctionMultiplicator< LocalFunctionType > uLocal2( wfLocal_, uLocal, p_ ); 353 integrator.integrateAdd( entity, uLocal2, sum ); 354 wfLocal_.unbind(); 355 } 356 357 358 template< class WeightFunction,class OrderCalculator > 359 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType > 360 inline void distanceLocal(const EntityType & entity,unsigned int order,const ULocalFunctionType & uLocal,const VLocalFunctionType & vLocal,ReturnType & sum) const361 WeightedLPNorm< WeightFunction, OrderCalculator >::distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const 362 { 363 typedef typename BaseType::template FunctionDistance< ULocalFunctionType, VLocalFunctionType > LocalDistanceType; 364 365 // !!!! order !!!! 366 IntegratorType integrator( order ); 367 368 wfLocal_.bind( entity ); 369 LocalDistanceType dist( uLocal, vLocal ); 370 WeightedFunctionMultiplicator< LocalDistanceType > dist2( wfLocal_, dist ); 371 372 integrator.integrateAdd( entity, dist2, sum ); 373 wfLocal_.unbind(); 374 } 375 376 377 template< class WeightFunction, class OrderCalculator> 378 template< class Function > 379 struct WeightedLPNorm< WeightFunction, OrderCalculator>::WeightedFunctionMultiplicator 380 { 381 typedef Function FunctionType; 382 383 typedef typename FunctionType::RangeFieldType RangeFieldType; 384 typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType; 385 typedef FieldVector< RealType, 1 > RangeType; 386 WeightedFunctionMultiplicatorDune::Fem::WeightedLPNorm::WeightedFunctionMultiplicator387 WeightedFunctionMultiplicator ( const LocalWeightFunctionType &weightFunction, 388 const FunctionType &function, 389 double p ) 390 : weightFunction_( weightFunction ), 391 function_( function ), 392 p_(p) 393 {} 394 395 template< class Point > evaluateDune::Fem::WeightedLPNorm::WeightedFunctionMultiplicator396 void evaluate ( const Point &x, RangeType &ret ) const 397 { 398 WeightType weight; 399 weightFunction_.evaluate( x, weight ); 400 401 typename FunctionType::RangeType phi; 402 function_.evaluate( x, phi ); 403 ret = weight[ 0 ] * std::pow ( phi.two_norm(), p_); 404 } 405 406 private: 407 const LocalWeightFunctionType &weightFunction_; 408 const FunctionType &function_; 409 double p_; 410 }; 411 412 } // namespace Fem 413 414 } // namespace Dune 415 416 #endif // #ifndef DUNE_FEM_LPNORM_HH 417