1 #ifndef DUNE_FEMPY_FUNCTION_SIMPLEGRIDFUNCTION_HH 2 #define DUNE_FEMPY_FUNCTION_SIMPLEGRIDFUNCTION_HH 3 4 #include <cassert> 5 6 #include <limits> 7 #include <type_traits> 8 #include <utility> 9 10 #include <dune/common/classname.hh> 11 #include <dune/common/exceptions.hh> 12 13 #include <dune/python/grid/simplegridfunction.hh> 14 #include <dune/python/grid/function.hh> 15 #include <dune/fem/function/common/discretefunction.hh> 16 #include <dune/fem/space/common/functionspace.hh> 17 #include <dune/fempy/function/virtualizedgridfunction.hh> 18 #include <dune/fem/common/intersectionside.hh> 19 20 namespace Dune 21 { 22 23 namespace FemPy 24 { 25 26 // SimpleLocalFunction 27 // ------------------- 28 29 template< class GridPart, class LocalEvaluator > 30 class SimpleLocalFunction 31 { 32 typedef SimpleLocalFunction< GridPart, LocalEvaluator > This; 33 34 public: 35 typedef typename GridPart::template Codim< 0 >::EntityType EntityType; 36 typedef typename EntityType::Geometry::LocalCoordinate LocalCoordinateType; 37 typedef std::decay_t< std::result_of_t< LocalEvaluator( EntityType, LocalCoordinateType ) > > Value_; 38 typedef std::conditional_t<std::is_same_v<Value_,double>, Dune::FieldVector<double,1>, Value_ > Value; 39 40 typedef Fem::FunctionSpace< typename GridPart::ctype, typename FieldTraits< Value >::field_type, GridPart::dimensionworld, Value::dimension > FunctionSpaceType; 41 42 static const int dimDomain = FunctionSpaceType::dimDomain; 43 static const int dimRange = FunctionSpaceType::dimRange; 44 45 typedef typename FunctionSpaceType::DomainType DomainType; 46 typedef typename FunctionSpaceType::RangeType RangeType; 47 typedef typename FunctionSpaceType::RangeFieldType RangeFieldType; 48 typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType; 49 typedef typename FunctionSpaceType::HessianRangeType HessianRangeType; 50 51 template< class GridFunction, std::enable_if_t< std::is_same< This, typename GridFunction::LocalFunctionType >::value, int > = 0 > SimpleLocalFunction(const GridFunction & gridFunction)52 SimpleLocalFunction ( const GridFunction &gridFunction ) 53 : localEvaluator_( gridFunction.localEvaluator() ), order_( gridFunction.order() ) 54 {} 55 SimpleLocalFunction(const EntityType & entity,LocalEvaluator localEvaluator,int order)56 SimpleLocalFunction ( const EntityType &entity, LocalEvaluator localEvaluator, int order ) 57 : entity_( &entity ), localEvaluator_( std::move( localEvaluator ) ), order_( order ) 58 {} 59 init(const EntityType & entity)60 void init ( const EntityType &entity ) { entity_ = &entity; } 61 62 template< class Point > evaluate(const Point & x,RangeType & value) const63 void evaluate ( const Point &x, RangeType &value ) const 64 { 65 using Fem::coordinate; 66 value = localEvaluator_( entity(), coordinate( x ) ); 67 } 68 69 template< class Quadrature, class Values > evaluateQuadrature(const Quadrature & quadrature,Values & values) const70 void evaluateQuadrature ( const Quadrature &quadrature, Values &values ) const 71 { 72 for( const auto qp : quadrature ) 73 evaluate( qp, values[ qp.index() ] ); 74 } 75 76 template< class Point > jacobian(const Point & x,JacobianRangeType & jacobian) const77 void jacobian ( const Point &x, JacobianRangeType &jacobian ) const 78 { 79 DUNE_THROW( NotImplemented, "SimpleLocalFunction::jacobian not implemented" ); 80 } 81 82 template< class Quadrature, class Jacobians > jacobianQuadrature(const Quadrature & quadrature,Jacobians & jacobians) const83 void jacobianQuadrature ( const Quadrature &quadrature, Jacobians &jacobians ) const 84 { 85 for( const auto qp : quadrature ) 86 jacobian( qp, jacobians[ qp.index() ] ); 87 } 88 89 template< class Point > hessian(const Point & x,HessianRangeType & hessian) const90 void hessian ( const Point &x, HessianRangeType &hessian ) const 91 { 92 DUNE_THROW( NotImplemented, "SimpleLocalFunction::hessian not implemented" ); 93 } 94 95 template< class Quadrature, class Hessians > hessianQuadrature(const Quadrature & quadrature,Hessians & hessians) const96 void hessianQuadrature ( const Quadrature &quadrature, Hessians &hessians ) const 97 { 98 for( const auto qp : quadrature ) 99 hessian( qp, hessians[ qp.index() ] ); 100 } 101 order() const102 int order () const { return order_; } 103 entity() const104 const EntityType &entity () const { assert( entity_ ); return *entity_; } 105 106 private: 107 const EntityType *entity_ = nullptr; 108 LocalEvaluator localEvaluator_; 109 int order_; 110 }; 111 112 113 114 // IsLocalEvaluator 115 // ---------------- 116 117 template< class LE, class E, class X > 118 std::true_type __isLocalEvaluator ( const LE &, const E &, const X &, decltype( std::declval< LE >()( std::declval< E >(), std::declval< X >() ) ) * = nullptr ); 119 120 std::false_type __isLocalEvaluator ( ... ); 121 122 template< class GP, class LE > 123 struct IsLocalEvaluator 124 : public decltype( __isLocalEvaluator( std::declval< LE >(), std::declval< typename GP::template Codim< 0 >::EntityType >(), std::declval< typename GP::template Codim< 0 >::GeometryType::LocalCoordinate >() ) ) 125 {}; 126 127 128 129 // SimpleGridFunction 130 // ------------------ 131 132 template< class GridPart, class LocalEvaluator > 133 class SimpleGridFunction 134 : public Fem::Function< typename SimpleLocalFunction< GridPart, LocalEvaluator >::FunctionSpaceType, SimpleGridFunction< GridPart, LocalEvaluator > >, 135 public Fem::HasLocalFunction 136 { 137 typedef SimpleGridFunction< GridPart, LocalEvaluator > This; 138 typedef Fem::Function< typename SimpleLocalFunction< GridPart, LocalEvaluator >::FunctionSpaceType, SimpleGridFunction< GridPart, LocalEvaluator > > Base; 139 // typedef typename SimpleLocalFunction< GridPart, LocalEvaluator >::FunctionSpaceType FunctionSpaceType; 140 141 struct Space 142 { 143 typedef typename SimpleLocalFunction< GridPart, LocalEvaluator >::FunctionSpaceType FunctionSpaceType; 144 typedef GridPart GridPartType; 145 static const int dimRange = FunctionSpaceType::RangeType::dimension; SpaceDune::FemPy::SimpleGridFunction::Space146 Space(const GridPart &gridPart, int o) 147 : gp_(gridPart), o_(o) {} orderDune::FemPy::SimpleGridFunction::Space148 int order() const 149 { 150 return o_; 151 } gridPartDune::FemPy::SimpleGridFunction::Space152 const GridPart& gridPart() const 153 { 154 return gp_; 155 } 156 const GridPart &gp_; 157 int o_; 158 }; 159 public: 160 typedef Space DiscreteFunctionSpaceType; 161 typedef GridPart GridPartType; 162 typedef typename GridPart::GridViewType GridView; 163 164 typedef SimpleLocalFunction< GridPart, LocalEvaluator > LocalFunctionType; 165 166 typedef typename LocalFunctionType::EntityType EntityType; 167 168 typedef typename Base::DomainType DomainType; 169 typedef typename Base::RangeType RangeType; 170 typedef typename RangeType::field_type RangeFieldType; 171 typedef typename Base::JacobianRangeType JacobianRangeType; 172 173 static const int dimRange = RangeType::dimension; 174 SimpleGridFunction(std::string name,const GridPartType & gridPart,LocalEvaluator localEvaluator,int order)175 SimpleGridFunction ( std::string name, const GridPartType &gridPart, LocalEvaluator localEvaluator, int order) // = std::numeric_limits< int >::max() ) 176 : name_( std::move( name ) ), 177 gridPart_( gridPart ), 178 localEvaluator_( std::move( localEvaluator ) ), 179 order_( order ) 180 {} 181 SimpleGridFunction(const GridPartType & gridPart,LocalEvaluator localEvaluator,int order)182 SimpleGridFunction ( const GridPartType &gridPart, LocalEvaluator localEvaluator, int order) // = std::numeric_limits< int >::max() ) 183 : name_( "unnamed-" + className< LocalEvaluator >() ), 184 gridPart_( gridPart ), 185 localEvaluator_( std::move( localEvaluator ) ), 186 order_( order ) 187 {} 188 localFunction(const EntityType & entity) const189 LocalFunctionType localFunction ( const EntityType &entity ) const { return LocalFunctionType( entity, localEvaluator_, order_ ); } 190 name() const191 const std::string &name () const { return name_; } 192 gridPart() const193 const GridPartType &gridPart () const { return gridPart_; } 194 evaluate(const DomainType & x,RangeType & value) const195 void evaluate ( const DomainType &x, RangeType &value ) const 196 { 197 Fem::EntitySearch< GridPartType, 0 > entitySearch( gridPart() ); 198 const EntityType entity = entitySearch( x ); 199 const auto geometry = entity.geometry(); 200 localFunction( entity ).evaluate( geometry.local( x ), value ); 201 value = RangeType(0); 202 } 203 jacobian(const DomainType & x,JacobianRangeType & jacobian) const204 void jacobian ( const DomainType &x, JacobianRangeType &jacobian ) const 205 { 206 Fem::EntitySearch< GridPartType, 0 > entitySearch( gridPart() ); 207 const EntityType entity = entitySearch( x ); 208 const auto geometry = entity.geometry(); 209 localFunction( entity ).jacobian( geometry.local( x ), jacobian ); 210 jacobian = JacobianRangeType(0); 211 } 212 space() const213 const Space space() const 214 { 215 return Space(gridPart(), order_); 216 } 217 order() const218 int order () const { return order_; } 219 localEvaluator() const220 const LocalEvaluator &localEvaluator () const { return localEvaluator_; } 221 222 protected: 223 std::string name_; 224 const GridPartType &gridPart_; 225 LocalEvaluator localEvaluator_; 226 int order_; 227 }; 228 229 230 231 // IsEvaluator 232 // ----------- 233 234 template< class E, class X > 235 std::true_type __isEvaluator ( const E &, const X &, decltype( std::declval< E >()( std::declval< X >() ) ) * = nullptr ); 236 237 std::false_type __isEvaluator ( ... ); 238 239 template< class GP, class E > 240 struct IsEvaluator 241 : public decltype( __isEvaluator( std::declval< E >(), std::declval< typename GP::template Codim< 0 >::GeometryType::GlobalCoordinate >() ) ) 242 {}; 243 244 245 246 // LocalEvaluatorAdapter 247 // --------------------- 248 249 template< class Entity, class Evaluator > 250 struct LocalEvaluatorAdapter 251 { 252 typedef typename Entity::Geometry::GlobalCoordinate GlobalCoordinate; 253 typedef typename Entity::Geometry::LocalCoordinate LocalCoordinate; 254 255 typedef decltype( std::declval< Evaluator >()( std::declval< GlobalCoordinate >() ) ) Value; 256 LocalEvaluatorAdapterDune::FemPy::LocalEvaluatorAdapter257 LocalEvaluatorAdapter ( Evaluator evaluator ) : evaluator_( std::move( evaluator ) ) {} 258 operator ()Dune::FemPy::LocalEvaluatorAdapter259 Value operator () ( const GlobalCoordinate &x ) const { return evaluator_( x ); } operator ()Dune::FemPy::LocalEvaluatorAdapter260 Value operator () ( const Entity &entity, const LocalCoordinate &x ) const { return evaluator_( entity.geometry().global( x ) ); } 261 262 private: 263 Evaluator evaluator_; 264 }; 265 266 267 268 // SimpleGlobalGridFunction 269 // ------------------------ 270 271 template< class GridPart, class Evaluator > 272 class SimpleGlobalGridFunction 273 : public SimpleGridFunction< GridPart, LocalEvaluatorAdapter< typename GridPart::template Codim< 0 >::EntityType, Evaluator > > 274 { 275 typedef SimpleGlobalGridFunction< GridPart, Evaluator > This; 276 typedef SimpleGridFunction< GridPart, LocalEvaluatorAdapter< typename GridPart::template Codim< 0 >::EntityType, Evaluator > > Base; 277 278 public: 279 typedef typename Base::GridPartType GridPartType; 280 typedef typename GridPart::GridViewType GridView; 281 282 typedef typename Base::DomainType DomainType; 283 typedef typename Base::RangeType RangeType; 284 SimpleGlobalGridFunction(std::string name,const GridPartType & gridPart,Evaluator evaluator,int order=std::numeric_limits<int>::max ())285 SimpleGlobalGridFunction ( std::string name, const GridPartType &gridPart, Evaluator evaluator, int order = std::numeric_limits< int >::max() ) 286 : Base( std::move( name ), gridPart, std::move( evaluator ), order ) 287 {} 288 SimpleGlobalGridFunction(const GridPartType & gridPart,Evaluator evaluator,int order=std::numeric_limits<int>::max ())289 SimpleGlobalGridFunction ( const GridPartType &gridPart, Evaluator evaluator, int order = std::numeric_limits< int >::max() ) 290 : Base( "unnamed-" + className< Evaluator >(), gridPart, std::move( evaluator ), order ) 291 {} 292 evaluate(const DomainType & x,RangeType & value) const293 void evaluate ( const DomainType &x, RangeType &value ) const { value = localEvaluator_( x ); } 294 295 protected: 296 using Base::localEvaluator_; 297 }; 298 299 300 301 // simpleGridFunction 302 // ------------------ 303 304 template< class GridPart, class LocalEvaluator > 305 inline static std::enable_if_t< IsLocalEvaluator< GridPart, LocalEvaluator >::value, SimpleGridFunction< GridPart, LocalEvaluator > > simpleGridFunction(std::string name,const GridPart & gridPart,LocalEvaluator localEvaluator,int order)306 simpleGridFunction ( std::string name, const GridPart &gridPart, LocalEvaluator localEvaluator, int order) 307 { 308 return SimpleGridFunction< GridPart, LocalEvaluator >( std::move( name ), gridPart, std::move( localEvaluator ), order ); 309 } 310 311 template< class GridPart, class LocalEvaluator > 312 inline static std::enable_if_t< IsLocalEvaluator< GridPart, LocalEvaluator >::value, SimpleGridFunction< GridPart, LocalEvaluator > > simpleGridFunction(const GridPart & gridPart,LocalEvaluator localEvaluator,int order)313 simpleGridFunction ( const GridPart &gridPart, LocalEvaluator localEvaluator, int order) 314 { 315 return SimpleGridFunction< GridPart, LocalEvaluator >( gridPart, std::move( localEvaluator ), order ); 316 } 317 318 319 // simpleGridFunction 320 // ------------------ 321 322 template< class GridPart, class Evaluator > 323 inline static std::enable_if_t< IsEvaluator< GridPart, Evaluator >::value, SimpleGlobalGridFunction< GridPart, Evaluator > > simpleGridFunction(std::string name,const GridPart & gridPart,Evaluator evaluator,int order)324 simpleGridFunction ( std::string name, const GridPart &gridPart, Evaluator evaluator, int order) 325 { 326 return SimpleGlobalGridFunction< GridPart, Evaluator >( std::move( name ), gridPart, std::move( evaluator ), order ); 327 } 328 329 template< class GridPart, class Evaluator > 330 inline static std::enable_if_t< IsEvaluator< GridPart, Evaluator >::value, SimpleGlobalGridFunction< GridPart, Evaluator > > simpleGridFunction(const GridPart & gridPart,Evaluator evaluator,int order)331 simpleGridFunction ( const GridPart &gridPart, Evaluator evaluator, int order) 332 { 333 return SimpleGlobalGridFunction< GridPart, Evaluator >( gridPart, std::move( evaluator ), order ); 334 } 335 336 template <class GridPart, class GridFunction, 337 std::enable_if_t<std::is_class_v<typename GridFunction::LocalEvaluator>,int> v=0> getGridFunction(const GridPart & gridPart,const GridFunction & gridFunction,int order,PriorityTag<1>)338 auto getGridFunction( const GridPart &gridPart, const GridFunction &gridFunction, int order, PriorityTag<1>) 339 { 340 return simpleGridFunction(gridPart,gridFunction.localEvaluator(),order); 341 } 342 template <class GridPart, class GridFunction> getGridFunction(const GridPart & gridPart,const GridFunction & gridFunction,int order,PriorityTag<0>)343 auto getGridFunction( const GridPart &gridPart, const GridFunction &gridFunction, int order, PriorityTag<0>) 344 { 345 return gridFunction; 346 } 347 348 template <class GridPart, class GridFunction> 349 struct GetGridFunction 350 { 351 using value = decltype(getGridFunction( 352 std::declval<const GridPart&>(), 353 std::declval<const GridFunction&>(), 354 1,PriorityTag<42>())); 355 }; 356 357 template <class GridPart, class RangeType, class Functor> addVirtualizedFunctions(Functor functor)358 void addVirtualizedFunctions(Functor functor) 359 { 360 typedef typename GridPart::template Codim<0>::EntityType Entity; 361 typedef typename Entity::Geometry::LocalCoordinate Coordinate; 362 static const int dimRange = RangeType::dimension; 363 typedef std::tuple< 364 VirtualizedGridFunction< GridPart, RangeType >*, 365 Dune::Python::SimpleGridFunction< typename GridPart::GridViewType, 366 Dune::Python::detail::PyGridFunctionEvaluator< typename GridPart::GridViewType, dimRange, pybind11::function >>*, 367 Dune::Python::SimpleGridFunction< typename GridPart::GridViewType, 368 Dune::Python::detail::PyGridFunctionEvaluator< typename GridPart::GridViewType, dimRange, 369 std::function<Dune::FieldVector<double,dimRange>(const Entity&,const Coordinate&)> >>* 370 > VirtualFcts; 371 Hybrid::forEach(VirtualFcts{0,0,0},functor); 372 typedef std::tuple< 373 Dune::Python::SimpleGridFunction< typename GridPart::GridViewType, 374 Dune::Python::detail::PyGridFunctionEvaluator< typename GridPart::GridViewType, 0, pybind11::function >>*, 375 Dune::Python::SimpleGridFunction< typename GridPart::GridViewType, 376 Dune::Python::detail::PyGridFunctionEvaluator< typename GridPart::GridViewType, 0, 377 std::function<double(const Entity&,const Coordinate&)> >>* 378 > VirtualFctsDbl; 379 if constexpr (dimRange == 1) 380 Hybrid::forEach(VirtualFctsDbl{0,0},functor); 381 } 382 } // namespace FemPy 383 384 namespace Fem 385 { 386 namespace Impl 387 { 388 template< class GV, class Evaluate > 389 struct ConstLocalFunction< Dune::Python::SimpleGridFunction<GV,Evaluate> > 390 { 391 struct Type 392 { 393 typedef Dune::Python::SimpleGridFunction<GV,Evaluate> GridFunctionType; 394 typedef decltype(localFunction(std::declval<const GridFunctionType&>())) LocalFunctionType; 395 typedef typename LocalFunctionType::Element EntityType; 396 397 typedef typename LocalFunctionType::LocalCoordinate DomainType; 398 typedef typename LocalFunctionType::Value RangeType; 399 // typedef typename LocalFunctionType::JacobianRangeType JacobianRangeType; 400 // typedef typename LocalFunctionType::HessianRangeType HessianRangeType; 401 TypeDune::Fem::Impl::ConstLocalFunction::Type402 explicit Type ( const GridFunctionType &gridFunction ) 403 : gridFunction_( gridFunction ), 404 localFunction_( localFunction(gridFunction) ) 405 {} TypeDune::Fem::Impl::ConstLocalFunction::Type406 explicit Type ( const EntityType &entity, const GridFunctionType &gridFunction ) 407 : gridFunction_( gridFunction ), 408 localFunction_( localFunction(gridFunction) ) 409 { bind(entity); } 410 411 //! evaluate local function 412 template< class Point > operator ()Dune::Fem::Impl::ConstLocalFunction::Type413 RangeType operator() ( const Point &p ) const 414 { 415 return localFunction_( Fem::coordinate(p) ); 416 } 417 template< class Point > evaluateDune::Fem::Impl::ConstLocalFunction::Type418 RangeType evaluate ( const Point &p ) const 419 { 420 return localFunction_( Fem::coordinate(p) ); 421 } 422 //! evaluate local function 423 template< class Point > evaluateDune::Fem::Impl::ConstLocalFunction::Type424 void evaluate ( const Point &p, RangeType &val ) const 425 { 426 val = localFunction_( Fem::coordinate(p) ); 427 } 428 template< class Point, class JacobianRangeType > jacobianDune::Fem::Impl::ConstLocalFunction::Type429 void jacobian ( const Point &p, JacobianRangeType &val ) const 430 { 431 DUNE_THROW(NotImplemented,"SimpleLocalFunction does not provide a jacobian!"); 432 } 433 template< class Point, class HessianRangeType > hessianDune::Fem::Impl::ConstLocalFunction::Type434 void hessian ( const Point &p, HessianRangeType &val ) const 435 { 436 DUNE_THROW(NotImplemented,"SimpleLocalFunction does not provide a hessian!"); 437 } 438 bindDune::Fem::Impl::ConstLocalFunction::Type439 void bind ( const EntityType &entity ) { localFunction_.bind( entity ); } unbindDune::Fem::Impl::ConstLocalFunction::Type440 void unbind () {} 441 template <class IntersectionType> bindDune::Fem::Impl::ConstLocalFunction::Type442 void bind(const IntersectionType &intersection, IntersectionSide side) 443 { defaultIntersectionBind(localFunction_, intersection, side); } 444 gridFunctionDune::Fem::Impl::ConstLocalFunction::Type445 const GridFunctionType &gridFunction () const { return gridFunction_; } 446 447 private: 448 const GridFunctionType &gridFunction_; 449 LocalFunctionType localFunction_; 450 }; 451 }; 452 } // namespace Impl 453 } // namespace Fem 454 455 } // namespace Dune 456 457 #endif // #ifndef DUNE_FEMPY_FUNCTION_SIMPLEGRIDFUNCTION_HH 458