2 #ifndef DUNE_PDELAB_LOCALOPERATOR_DARCYFEM_HH
3 #define DUNE_PDELAB_LOCALOPERATOR_DARCYFEM_HH
5 #include <dune/common/fvector.hh>
6 #include <dune/geometry/referenceelements.hh>
21 template<
typename P,
typename T,
typename X>
24 Dune::PDELab::GridFunctionTraits<
25 typename T::Traits::GridViewType,
26 typename T::Traits::FiniteElementType::Traits::LocalBasisType
27 ::Traits::RangeFieldType,
28 T::Traits::FiniteElementType::Traits::LocalBasisType::Traits
31 typename T::Traits::FiniteElementType::Traits
32 ::LocalBasisType::Traits::RangeFieldType,
33 T::Traits::FiniteElementType::Traits::LocalBasisType::Traits
35 DarcyVelocityFromHeadFEM<P,T,X> >
38 using LBTraits =
typename GFS::Traits::FiniteElementType::
39 Traits::LocalBasisType::Traits;
43 using LView =
typename X::template LocalView<LFSCache>;
47 typename GFS::Traits::GridViewType,
48 typename LBTraits::RangeFieldType,
51 typename LBTraits::RangeFieldType,
52 LBTraits::dimDomain> >;
66 : pgfs(stackobject_to_shared_ptr(gfs))
67 , pxg(stackobject_to_shared_ptr(x_))
68 , pp(stackobject_to_shared_ptr(
p))
84 std::vector<typename Traits::RangeFieldType> xl(lfs.size());
85 lview.bind(lfs_cache);
90 auto geo =
e.geometry();
93 const typename Traits::ElementType::Geometry::JacobianInverseTransposed
94 JgeoIT(geo.jacobianInverseTransposed(x));
97 std::vector<typename LBTraits::JacobianType> J(lfs.size());
98 lfs.finiteElement().localBasis().evaluateJacobian(x,J);
102 for(
unsigned int i = 0; i < lfs.size(); ++i) {
105 JgeoIT.umv(J[i][0], gradphi);
108 minusgrad.axpy(-xl[i], gradphi);
112 auto inside_cell_center_local = referenceElement(geo).position(0,0);
113 typename P::Traits::PermTensorType A(pp->A(
e,inside_cell_center_local));
120 return pgfs->gridView();
124 std::shared_ptr<const GFS> pgfs;
125 std::shared_ptr<X> pxg;
126 std::shared_ptr<const P> pp;
128 mutable LFSCache lfs_cache;
132 #endif // DUNE_PDELAB_LOCALOPERATOR_DARCYFEM_HH