1 #ifndef DUNE_FEM_L2NORM_HH 2 #define DUNE_FEM_L2NORM_HH 3 4 #include <dune/fem/quadrature/cachingquadrature.hh> 5 #include <dune/fem/quadrature/integrator.hh> 6 7 #include <dune/fem/misc/domainintegral.hh> 8 9 namespace Dune 10 { 11 12 namespace Fem 13 { 14 15 // L2Norm 16 // ------ 17 18 template< class GridPart > 19 class L2Norm : public IntegralBase< GridPart, L2Norm< GridPart > > 20 { 21 typedef IntegralBase< GridPart, L2Norm< GridPart > > BaseType ; 22 typedef L2Norm< GridPart > ThisType; 23 24 public: 25 typedef GridPart GridPartType; 26 27 using BaseType :: gridPart ; 28 using BaseType :: comm ; 29 30 template< class Function > 31 struct FunctionSquare; 32 33 template< class UFunction, class VFunction > 34 struct FunctionDistance; 35 36 protected: 37 typedef typename BaseType::EntityType EntityType; 38 typedef CachingQuadrature< GridPartType, 0 > QuadratureType; 39 typedef Integrator< QuadratureType > IntegratorType; 40 41 const unsigned int order_; 42 const bool communicate_; 43 public: 44 /** \brief constructor 45 * \param gridPart specific gridPart for selection of entities 46 * \param order order of integration quadrature (default = 2*space.order()) 47 * \param communicate if true global (over all ranks) norm is computed (default = true) 48 */ 49 explicit L2Norm ( const GridPartType &gridPart, 50 const unsigned int order = 0, 51 const bool communicate = true ); 52 53 //! || u ||_L2 on given set of entities (partition set) 54 template< class DiscreteFunctionType, class PartitionSet > 55 typename Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type 56 norm ( const DiscreteFunctionType &u, const PartitionSet& partitionSet ) const; 57 58 //! || u ||_L2 on interior partition entities 59 template< class DiscreteFunctionType > 60 typename Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type norm(const DiscreteFunctionType & u) const61 norm ( const DiscreteFunctionType &u ) const 62 { 63 return norm( u, Partitions::interior ); 64 } 65 66 //! || u - v ||_L2 on given set of entities (partition set) 67 template< class UDiscreteFunctionType, class VDiscreteFunctionType, class PartitionSet > 68 typename Dune::FieldTraits< typename UDiscreteFunctionType::RangeFieldType >::real_type 69 distance ( const UDiscreteFunctionType &u, 70 const VDiscreteFunctionType &v, 71 const PartitionSet& partitionSet ) const; 72 73 //! || u - v ||_L2 on interior partition entities 74 template< class UDiscreteFunctionType, class VDiscreteFunctionType > 75 typename Dune::FieldTraits< typename UDiscreteFunctionType::RangeFieldType >::real_type distance(const UDiscreteFunctionType & u,const VDiscreteFunctionType & v) const76 distance ( const UDiscreteFunctionType &u, const VDiscreteFunctionType &v ) const 77 { 78 return distance( u, v, Partitions::interior ); 79 } 80 81 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType > 82 void distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const; 83 84 template< class LocalFunctionType, class ReturnType > 85 void normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const; 86 }; 87 88 89 // Implementation of L2Norm 90 // ------------------------ 91 92 template< class GridPart > L2Norm(const GridPartType & gridPart,const unsigned int order,const bool communicate)93 inline L2Norm< GridPart >::L2Norm ( const GridPartType &gridPart, const unsigned int order, const bool communicate ) 94 : BaseType( gridPart ), 95 order_( order ), 96 communicate_( BaseType::checkCommunicateFlag( communicate ) ) 97 { 98 } 99 100 101 template< class GridPart > 102 template< class DiscreteFunctionType, class PartitionSet > 103 typename Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type norm(const DiscreteFunctionType & u,const PartitionSet & partitionSet) const104 L2Norm< GridPart >::norm ( const DiscreteFunctionType &u, const PartitionSet& partitionSet ) const 105 { 106 typedef typename Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type RealType; 107 typedef FieldVector< RealType, 1 > ReturnType ; 108 109 // calculate integral over each element 110 ReturnType sum = BaseType :: forEach( u, ReturnType(0), partitionSet, order_ ); 111 112 // communicate_ indicates global norm 113 if( communicate_ ) 114 { 115 sum[ 0 ] = comm().sum( sum[ 0 ] ); 116 } 117 118 // return result, e.g. sqrt of calculated sum 119 return std::sqrt( sum[ 0 ] ); 120 } 121 122 123 template< class GridPart > 124 template< class UDiscreteFunctionType, class VDiscreteFunctionType, class PartitionSet > 125 typename Dune::FieldTraits< typename UDiscreteFunctionType::RangeFieldType >::real_type 126 L2Norm< GridPart > distance(const UDiscreteFunctionType & u,const VDiscreteFunctionType & v,const PartitionSet & partitionSet) const127 ::distance ( const UDiscreteFunctionType &u, 128 const VDiscreteFunctionType &v, 129 const PartitionSet& partitionSet ) const 130 { 131 typedef typename Dune::FieldTraits< typename UDiscreteFunctionType::RangeFieldType >::real_type RealType; 132 typedef FieldVector< RealType, 1 > ReturnType ; 133 134 // calculate integral over each element 135 ReturnType sum = BaseType :: forEach( u, v, ReturnType(0), partitionSet, order_ ); 136 137 // communicate_ indicates global norm 138 if( communicate_ ) 139 { 140 sum[ 0 ] = comm().sum( sum[ 0 ] ); 141 } 142 143 // return result, e.g. sqrt of calculated sum 144 return std::sqrt( sum[ 0 ] ); 145 } 146 147 template< class GridPart > 148 template< class LocalFunctionType, class ReturnType > normLocal(const EntityType & entity,unsigned int order,const LocalFunctionType & uLocal,ReturnType & sum) const149 inline void L2Norm< GridPart >::normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const 150 { 151 // evaluate norm locally 152 IntegratorType integrator( order ); 153 154 FunctionSquare< LocalFunctionType > uLocal2( uLocal ); 155 156 integrator.integrateAdd( entity, uLocal2, sum ); 157 } 158 159 template< class GridPart > 160 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType > 161 inline void distanceLocal(const EntityType & entity,unsigned int order,const ULocalFunctionType & uLocal,const VLocalFunctionType & vLocal,ReturnType & sum) const162 L2Norm< GridPart >::distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const 163 { 164 // evaluate norm locally 165 IntegratorType integrator( order ); 166 167 typedef FunctionDistance< ULocalFunctionType, VLocalFunctionType > LocalDistanceType; 168 169 LocalDistanceType dist( uLocal, vLocal ); 170 FunctionSquare< LocalDistanceType > dist2( dist ); 171 172 integrator.integrateAdd( entity, dist2, sum ); 173 } 174 175 176 template< class GridPart > 177 template< class Function > 178 struct L2Norm< GridPart >::FunctionSquare 179 { 180 typedef Function FunctionType; 181 182 typedef typename FunctionType::RangeFieldType RangeFieldType; 183 typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType; 184 typedef FieldVector< RealType, 1 > RangeType; 185 FunctionSquareDune::Fem::L2Norm::FunctionSquare186 explicit FunctionSquare ( const FunctionType &function ) 187 : function_( function ) 188 {} 189 190 template< class Point > evaluateDune::Fem::L2Norm::FunctionSquare191 void evaluate ( const Point &x, RangeType &ret ) const 192 { 193 typename FunctionType::RangeType phi; 194 function_.evaluate( x, phi ); 195 ret = phi.two_norm2(); 196 } 197 198 private: 199 const FunctionType &function_; 200 }; 201 202 203 template< class GridPart > 204 template< class UFunction, class VFunction > 205 struct L2Norm< GridPart >::FunctionDistance 206 { 207 typedef UFunction UFunctionType; 208 typedef VFunction VFunctionType; 209 210 typedef typename UFunctionType::RangeFieldType RangeFieldType; 211 typedef typename UFunctionType::RangeType RangeType; 212 typedef typename UFunctionType::JacobianRangeType JacobianRangeType; 213 FunctionDistanceDune::Fem::L2Norm::FunctionDistance214 FunctionDistance ( const UFunctionType &u, const VFunctionType &v ) 215 : u_( u ), v_( v ) 216 {} 217 218 template< class Point > evaluateDune::Fem::L2Norm::FunctionDistance219 void evaluate ( const Point &x, RangeType &ret ) const 220 { 221 RangeType phi; 222 u_.evaluate( x, ret ); 223 v_.evaluate( x, phi ); 224 ret -= phi; 225 } 226 227 template< class Point > jacobianDune::Fem::L2Norm::FunctionDistance228 void jacobian ( const Point &x, JacobianRangeType &ret ) const 229 { 230 JacobianRangeType phi; 231 u_.jacobian( x, ret ); 232 v_.jacobian( x, phi ); 233 ret -= phi; 234 } 235 236 private: 237 const UFunctionType &u_; 238 const VFunctionType &v_; 239 }; 240 241 242 243 // WeightedL2Norm 244 // -------------- 245 246 template< class WeightFunction > 247 class WeightedL2Norm 248 : public L2Norm< typename WeightFunction::DiscreteFunctionSpaceType::GridPartType > 249 { 250 typedef WeightedL2Norm< WeightFunction > ThisType; 251 typedef L2Norm< typename WeightFunction::DiscreteFunctionSpaceType::GridPartType > BaseType; 252 253 public: 254 typedef WeightFunction WeightFunctionType; 255 256 typedef typename WeightFunctionType::DiscreteFunctionSpaceType WeightFunctionSpaceType; 257 typedef typename WeightFunctionSpaceType::GridPartType GridPartType; 258 259 protected: 260 template< class Function > 261 struct WeightedFunctionSquare; 262 263 typedef ConstLocalFunction< WeightFunctionType > LocalWeightFunctionType; 264 typedef typename WeightFunctionType::RangeType WeightType; 265 266 typedef typename BaseType::EntityType EntityType; 267 typedef typename BaseType::IntegratorType IntegratorType; 268 269 using BaseType::gridPart; 270 using BaseType::comm; 271 272 public: 273 274 using BaseType::norm; 275 using BaseType::distance; 276 277 explicit WeightedL2Norm ( const WeightFunctionType &weightFunction, const unsigned int order = 0, const bool communicate = true ); 278 279 template< class LocalFunctionType, class ReturnType > 280 void normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const; 281 282 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType > 283 void distanceLocal ( const EntityType &entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType &sum ) const; 284 285 private: 286 LocalWeightFunctionType wfLocal_; 287 }; 288 289 290 291 292 // Implementation of WeightedL2Norm 293 // -------------------------------- 294 295 template< class WeightFunction > 296 inline WeightedL2Norm< WeightFunction > WeightedL2Norm(const WeightFunctionType & weightFunction,const unsigned int order,const bool communicate)297 ::WeightedL2Norm ( const WeightFunctionType &weightFunction, const unsigned int order, const bool communicate ) 298 : BaseType( weightFunction.space().gridPart(), order, communicate ), 299 wfLocal_( weightFunction ) 300 { 301 static_assert( (WeightFunctionSpaceType::dimRange == 1), 302 "Wight function must be scalar." ); 303 } 304 305 306 template< class WeightFunction > 307 template< class LocalFunctionType, class ReturnType > 308 inline void normLocal(const EntityType & entity,unsigned int order,const LocalFunctionType & uLocal,ReturnType & sum) const309 WeightedL2Norm< WeightFunction >::normLocal ( const EntityType &entity, unsigned int order, const LocalFunctionType &uLocal, ReturnType &sum ) const 310 { 311 // !!!! order !!!! 312 IntegratorType integrator( order ); 313 314 wfLocal_.bind( entity ); 315 WeightedFunctionSquare< LocalFunctionType > uLocal2( wfLocal_, uLocal ); 316 integrator.integrateAdd( entity, uLocal2, sum ); 317 wfLocal_.unbind(); 318 } 319 320 321 template< class WeightFunction > 322 template< class ULocalFunctionType, class VLocalFunctionType, class ReturnType > 323 inline void distanceLocal(const EntityType & entity,unsigned int order,const ULocalFunctionType & uLocal,const VLocalFunctionType & vLocal,ReturnType & sum) const324 WeightedL2Norm< WeightFunction >::distanceLocal ( const EntityType& entity, unsigned int order, const ULocalFunctionType &uLocal, const VLocalFunctionType &vLocal, ReturnType& sum ) const 325 { 326 typedef typename BaseType::template FunctionDistance< ULocalFunctionType, VLocalFunctionType > LocalDistanceType; 327 328 // !!!! order !!!! 329 IntegratorType integrator( order ); 330 331 wfLocal_.bind( entity ); 332 LocalDistanceType dist( uLocal, vLocal ); 333 WeightedFunctionSquare< LocalDistanceType > dist2( wfLocal_, dist ); 334 335 integrator.integrateAdd( entity, dist2, sum ); 336 wfLocal_.unbind(); 337 } 338 339 340 template< class WeightFunction > 341 template< class Function > 342 struct WeightedL2Norm< WeightFunction >::WeightedFunctionSquare 343 { 344 typedef Function FunctionType; 345 346 typedef typename FunctionType::RangeFieldType RangeFieldType; 347 typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType; 348 typedef FieldVector< RealType, 1 > RangeType; 349 WeightedFunctionSquareDune::Fem::WeightedL2Norm::WeightedFunctionSquare350 WeightedFunctionSquare ( const LocalWeightFunctionType &weightFunction, 351 const FunctionType &function ) 352 : weightFunction_( weightFunction ), 353 function_( function ) 354 {} 355 356 template< class Point > evaluateDune::Fem::WeightedL2Norm::WeightedFunctionSquare357 void evaluate ( const Point &x, RangeType &ret ) const 358 { 359 WeightType weight; 360 weightFunction_.evaluate( x, weight ); 361 362 typename FunctionType::RangeType phi; 363 function_.evaluate( x, phi ); 364 ret = weight[ 0 ] * (phi * phi); 365 } 366 367 private: 368 const LocalWeightFunctionType &weightFunction_; 369 const FunctionType &function_; 370 }; 371 372 } // end namespace Fem 373 374 } // end namespace Dune 375 376 #endif // #ifndef DUNE_FEM_L2NORM_HH 377