1 #ifndef DUNE_FEM_LFESPACE_LOCALRESTRICTPROLONG_HH
2 #define DUNE_FEM_LFESPACE_LOCALRESTRICTPROLONG_HH
3 
4 #include <dune/common/dynmatrix.hh>
5 #include <dune/geometry/referenceelements.hh>
6 #include <dune/fem/function/localfunction/localfunction.hh>
7 #include <dune/fem/space/localfiniteelement/space.hh>
8 #include <dune/fem/space/common/localinterpolation.hh>
9 #include <dune/fem/space/common/localrestrictprolong.hh>
10 
11 namespace Dune
12 {
13 
14   namespace Fem
15   {
16 
17     namespace Impl
18     {
19       template< class LocalGeometry, class LF>
20       struct FatherWrapper
21       {
22         typedef std::remove_reference_t< LF > LocalFunctionType;
23         typedef std::remove_reference_t< LocalGeometry > LocalGeometryType;
24         struct Traits
25         {
26           typedef typename LocalFunctionType::DomainType DomainType;
27           typedef typename LocalFunctionType::RangeType RangeType;
28         };
29         typedef typename LocalFunctionType::EntityType EntityType;
30         typedef typename LocalFunctionType::FunctionSpaceType FunctionSpaceType;
31         typedef typename LocalFunctionType::DomainType DomainType;
32         typedef typename LocalFunctionType::RangeType RangeType;
33         typedef typename LocalFunctionType::JacobianRangeType JacobianRangeType;
34         typedef typename LocalFunctionType::HessianRangeType HessianRangeType;
35 
FatherWrapperDune::Fem::Impl::FatherWrapper36         FatherWrapper ( const LocalGeometryType &localGeo, const LocalFunctionType &lfFather)
37           : localGeo_(localGeo), lfFather_(lfFather) {}
38         template <class Point>
evaluateDune::Fem::Impl::FatherWrapper39         void evaluate ( const Point &x, RangeType &y ) const
40         {
41           lfFather_.evaluate( localGeo_.global(coordinate(x)), y);
42         }
43         template <class Quadrature, class RangeArray>
evaluateQuadratureDune::Fem::Impl::FatherWrapper44         void evaluateQuadrature( const Quadrature& quadrature, RangeArray& values ) const
45         {
46           const unsigned int nop = quadrature.nop();
47           values.resize( nop );
48           for( unsigned int qp=0; qp<nop; ++qp)
49             evaluate( quadrature[ qp ], values[ qp ]);
50         }
51 
entityDune::Fem::Impl::FatherWrapper52         const EntityType& entity () const { return lfFather_.entity(); }
orderDune::Fem::Impl::FatherWrapper53         const unsigned int order () const { return lfFather_.order(); }
54 
55       private:
56         const LocalGeometryType &localGeo_;
57         const LocalFunctionType &lfFather_;
58       };
59       template< class BasisFunctionSet, class LF>
60       struct SonsWrapper
61       {
62         typedef std::remove_reference_t< LF > LocalFunctionType;
63         typedef std::remove_reference_t< BasisFunctionSet > BasisFunctionSetType;
64         struct Traits
65         {
66           typedef typename LocalFunctionType::DomainType DomainType;
67           typedef typename LocalFunctionType::RangeType RangeType;
68         };
69         typedef typename LocalFunctionType::FunctionSpaceType FunctionSpaceType;
70         typedef typename LocalFunctionType::DomainType DomainType;
71         typedef typename LocalFunctionType::RangeType RangeType;
72         typedef typename LocalFunctionType::JacobianRangeType JacobianRangeType;
73         typedef typename LocalFunctionType::HessianRangeType HessianRangeType;
74         typedef typename LocalFunctionType::EntityType EntityType;
75         typedef typename EntityType::LocalGeometry LocalGeometryType;
76 
77         template <class LFFather>
SonsWrapperDune::Fem::Impl::SonsWrapper78         SonsWrapper ( const LFFather &father,
79           const std::vector< EntityType >& childEntities,
80           const std::vector< BasisFunctionSetType >& childBasisSets,
81           const std::vector< std::vector<double> >& childDofs )
82           : father_(father.entity())
83           , order_(father.order())
84           , childEntities_(childEntities), childBasisSets_(childBasisSets), childDofs_(childDofs)
85         {}
86         template <class Point>
evaluateDune::Fem::Impl::SonsWrapper87         void evaluate ( const Point &x, RangeType &val ) const
88         {
89           val = RangeType(0);
90           RangeType tmp;
91           double weight = 0;
92           for (unsigned int i=0; i<childEntities_.size();++i)
93           {
94             const auto &refSon = Dune::ReferenceElements< typename LocalGeometryType::ctype, LocalGeometryType::mydimension >
95               ::general( childEntities_[i].type() );
96             auto y = childEntities_[i].geometryInFather().local(coordinate(x));
97             if( refSon.checkInside( y ) )
98             {
99               childBasisSets_[i].evaluateAll( y, childDofs_[i], tmp );
100               val += tmp;
101               weight += 1.;
102             }
103           }
104           assert( weight > 0); // weight==0 would mean that point was found in none of the children
105           val /= weight;
106         }
107         template <class Quadrature, class RangeArray>
evaluateQuadratureDune::Fem::Impl::SonsWrapper108         void evaluateQuadrature( const Quadrature& quadrature, RangeArray& values ) const
109         {
110           const unsigned int nop = quadrature.nop();
111           values.resize( nop );
112           for( unsigned int qp=0; qp<nop; ++qp)
113             evaluate( quadrature[ qp ], values[ qp ]);
114         }
entityDune::Fem::Impl::SonsWrapper115         const EntityType& entity () const { return father_; }
orderDune::Fem::Impl::SonsWrapper116         const unsigned int order () const { return order_; }
117 
118       private:
119         const EntityType &father_;
120         unsigned int order_;
121         const std::vector< EntityType >& childEntities_;
122         const std::vector< BasisFunctionSetType >& childBasisSets_;
123         const std::vector< std::vector<double> >& childDofs_;
124       };
125 
126       // a detailed description is given in MR308
127       // https://gitlab.dune-project.org/dune-fem/dune-fem/merge_requests/308
128       template< class LFESpace >
129       struct DefaultLocalRestrictProlongLFE
130       {
131         typedef LFESpace DiscreteFunctionSpaceType;
132         typedef DefaultLocalRestrictProlongLFE< DiscreteFunctionSpaceType > ThisType;
133 
134         typedef typename DiscreteFunctionSpaceType::DomainFieldType DomainFieldType;
135         typedef typename DiscreteFunctionSpaceType::BasisFunctionSetType BasisFunctionSetType;
136         typedef typename DiscreteFunctionSpaceType::EntityType EntityType;
137         typedef typename EntityType::LocalGeometry LocalGeometryType;
138         typedef typename EntityType::EntitySeed EntitySeedType;
139 
DefaultLocalRestrictProlongLFEDune::Fem::Impl::DefaultLocalRestrictProlongLFE140         DefaultLocalRestrictProlongLFE (const DiscreteFunctionSpaceType &space)
141         : space_( space ),
142           interpolation_( space_ ),
143           childSeeds_(0), childDofs_(0)
144         {}
145 
146         /** \brief explicit set volume ratio of son and father
147          *
148          *  \param[in]  weight  volume of son / volume of father
149          *
150          *  \note If this ratio is set, it is assume to be constant.
151          */
setFatherChildWeightDune::Fem::Impl::DefaultLocalRestrictProlongLFE152         void setFatherChildWeight ( const DomainFieldType &weight )
153         {
154         }
155 
156         //! restrict data to father
157         template< class LFFather, class LFSon, class LocalGeometry >
restrictLocalDune::Fem::Impl::DefaultLocalRestrictProlongLFE158         void restrictLocal ( LFFather &lfFather, const LFSon &lfSon,
159                              const LocalGeometry &geometryInFather, bool initialize ) const
160         {
161           assert( lfFather.numDofs() == lfSon.numDofs() );
162 
163           if (initialize)
164           {
165             childSeeds_.resize(0);
166             childDofs_.resize(0);
167           }
168 
169           childSeeds_.push_back(lfSon.entity().seed());
170           childDofs_.push_back(std::vector<double>(lfSon.size()));
171           std::vector<double>& chDofs = childDofs_.back();
172           const unsigned int numDofs = lfSon.size();
173           for (unsigned int i=0;i<numDofs; ++i )
174             chDofs[i] = lfSon[i];
175         }
176 
177         template <class LFFather>
restrictFinalizeDune::Fem::Impl::DefaultLocalRestrictProlongLFE178         void restrictFinalize( LFFather &lfFather ) const
179         {
180           std::vector< EntityType > childEntities(childSeeds_.size());
181           std::vector< BasisFunctionSetType > childBasisSets(childSeeds_.size());
182           const unsigned int cSize = childSeeds_.size();
183           for (unsigned int i=0; i<cSize; ++i)
184           {
185             childEntities[i] = space_.gridPart().entity( childSeeds_[i] );
186             childBasisSets[i] = space_.basisFunctionSet( childEntities[i] );
187           }
188 
189 
190           // bind interpolation object to father
191           auto guard = bindGuard( interpolation_, lfFather.entity() );
192 
193           typedef Impl::SonsWrapper<BasisFunctionSetType, LFFather> SonsWrapperType;
194           SonsWrapperType sonsWrapper( lfFather, childEntities, childBasisSets, childDofs_ );
195 
196           // call interpolation
197           interpolation_( sonsWrapper, lfFather );
198         }
199 
200         //! prolong data to children
201         template< class LFFather, class LFSon, class LocalGeometry >
prolongLocalDune::Fem::Impl::DefaultLocalRestrictProlongLFE202         void prolongLocal ( const LFFather &lfFather, LFSon &lfSon,
203                             const LocalGeometry &geometryInFather, bool initialize ) const
204         {
205           assert( lfFather.numDofs() == lfSon.numDofs() );
206 
207           typedef Impl::FatherWrapper<LocalGeometry,LFFather> FatherWrapperType;
208           FatherWrapperType fatherWrapper(geometryInFather,lfFather);
209 
210           // bind interpolation object to father
211           auto guard = bindGuard( interpolation_, lfSon.entity() );
212 
213           // call interpolation
214           interpolation_( fatherWrapper, lfSon );
215         }
216 
217         //! do discrete functions need a communication after restriction / prolongation?
needCommunicationDune::Fem::Impl::DefaultLocalRestrictProlongLFE218         bool needCommunication () const { return true; }
219 
220       protected:
221         template< class Entity >
calcWeightDune::Fem::Impl::DefaultLocalRestrictProlongLFE222         static DomainFieldType calcWeight ( const Entity &father, const Entity &son )
223         {
224           return son.geometry().volume() / father.geometry().volume();
225         }
226         const DiscreteFunctionSpaceType &space_;
227         mutable LocalInterpolation< DiscreteFunctionSpaceType > interpolation_;
228 
229         mutable std::vector< EntitySeedType > childSeeds_;
230         mutable std::vector< std::vector<double> > childDofs_;
231 
232         mutable std::vector< double > tmpDofs_;
233       };
234     } // namespce Impl
235 
236     template< class LFEMap, class FunctionSpace, class Storage >
237     class DefaultLocalRestrictProlong< LocalFiniteElementSpace< LFEMap, FunctionSpace, Storage > >
238     : public Impl::DefaultLocalRestrictProlongLFE
239                < LocalFiniteElementSpace< LFEMap, FunctionSpace, Storage > >
240     {
241     public:
242       typedef Impl::DefaultLocalRestrictProlongLFE
243                < LocalFiniteElementSpace< LFEMap, FunctionSpace, Storage > > BaseType;
244       using BaseType::BaseType;
245     };
246     template< class LFEMap, class FunctionSpace, class Storage >
247     class DefaultLocalRestrictProlong< DiscontinuousLocalFiniteElementSpace< LFEMap, FunctionSpace, Storage > >
248     : public Impl::DefaultLocalRestrictProlongLFE
249                < DiscontinuousLocalFiniteElementSpace< LFEMap, FunctionSpace, Storage > >
250 
251     {
252     public:
253       typedef Impl::DefaultLocalRestrictProlongLFE
254                < DiscontinuousLocalFiniteElementSpace< LFEMap, FunctionSpace, Storage > > BaseType;
255       using BaseType::BaseType;
256     };
257   } // namespace Fem
258 
259 } // namespace Dune
260 
261 #endif // #ifndef DUNE_FEM_LOCALRESTRICTPROLONG_HH
262