1 #ifndef DUNE_FEM_SPACE_DISCONTINUOUSGALERKIN_LOCALRESTRICTPROLONG_HH_
2 #define DUNE_FEM_SPACE_DISCONTINUOUSGALERKIN_LOCALRESTRICTPROLONG_HH_
3 
4 // dune-fem includes
5 #include <dune/fem/function/localfunction/temporary.hh>
6 #include <dune/fem/quadrature/cachingquadrature.hh>
7 #include <dune/fem/space/common/adaptationmanager.hh>
8 #include <dune/fem/space/common/localrestrictprolong.hh>
9 
10 // local includes
11 #include "declaration.hh"
12 #include "localdgmassmatrix.hh"
13 
14 
15 namespace Dune
16 {
17 
18   namespace Fem
19   {
20 
21     /** @ingroup RestrictProlongImpl
22         @{
23     **/
24 
25 
26     // DiscontinuousGalerkinLocalRestrictProlong
27     // -----------------------------------------
28 
29     template< class DiscreteFunctionSpace, bool applyInverse >
30     class DiscontinuousGalerkinLocalRestrictProlong
31     {
32       typedef DiscontinuousGalerkinLocalRestrictProlong< DiscreteFunctionSpace, applyInverse > ThisType;
33 
34     public:
35       typedef DiscreteFunctionSpace DiscreteFunctionSpaceType;
36 
37       typedef typename DiscreteFunctionSpaceType::DomainFieldType DomainFieldType;
38       typedef typename DiscreteFunctionSpaceType::RangeFieldType RangeFieldType;
39       typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
40 
41       typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
42 
43       typedef CachingQuadrature< GridPartType, 0 > QuadratureType;
44 
45       typedef LocalMassMatrix< DiscreteFunctionSpaceType, QuadratureType > LocalMassMatrixType;
46 
DiscontinuousGalerkinLocalRestrictProlong(const DiscreteFunctionSpaceType & space)47       DiscontinuousGalerkinLocalRestrictProlong ( const DiscreteFunctionSpaceType& space )
48       : localMassMatrix_( space, space.order() * 2 ),
49         weight_( -1 ),
50         temp_( space )
51       {}
52 
setFatherChildWeight(const DomainFieldType & weight)53       void setFatherChildWeight ( const DomainFieldType &weight )
54       {
55         weight_ = weight;
56       }
57 
58       //! restrict data to father
59       template< class LFFather, class LFSon, class LocalGeometry >
restrictLocal(LFFather & lfFather,const LFSon & lfSon,const LocalGeometry & geometryInFather,bool initialize) const60       void restrictLocal ( LFFather &lfFather, const LFSon &lfSon,
61                            const LocalGeometry &geometryInFather, bool initialize ) const
62       {
63         typedef ConstantLocalRestrictProlong< DiscreteFunctionSpaceType > ConstantLocalRestrictProlongType;
64         const DomainFieldType weight = (weight_ < DomainFieldType( 0 ) ? ConstantLocalRestrictProlongType::calcWeight( lfFather.entity(), lfSon.entity() ) : weight_);
65 
66         assert( weight > 0.0 );
67 
68         if( initialize )
69           lfFather.clear();
70 
71         if( applyInverse )
72         {
73           temp_.init( lfFather.entity() );
74           temp_.clear();
75         }
76 
77         typedef typename LFSon :: EntityType  EntityType ;
78         typedef typename EntityType :: Geometry   Geometry;
79         const EntityType& sonEntity = lfSon.entity();
80         const Geometry sonGeo = sonEntity.geometry();
81 
82         QuadratureType quad( sonEntity, 2*lfFather.order()+1 );
83         const int nop = quad.nop();
84         for( int qp = 0; qp < nop; ++qp )
85         {
86           RangeFieldType quadWeight = quad.weight( qp );
87 
88           // in case of non-orthonormal basis we have to
89           // apply the integration element and the
90           // inverse mass matrix later
91           if( applyInverse )
92           {
93             quadWeight *= sonGeo.integrationElement( quad.point(qp) );
94           }
95           else
96             quadWeight *= weight ;
97 
98           RangeType value;
99           lfSon.evaluate( quad[ qp ], value );
100           value *= quadWeight;
101 
102           if( applyInverse )
103             temp_.axpy( geometryInFather.global( quad.point( qp ) ), value );
104           else
105             lfFather.axpy( geometryInFather.global( quad.point( qp ) ), value );
106         }
107 
108         if( applyInverse )
109         {
110           localMassMatrix_.applyInverse( temp_ );
111           lfFather += temp_;
112         }
113       }
114       template< class LFFather >
restrictFinalize(LFFather & lfFather) const115       void restrictFinalize ( LFFather &lfFather ) const
116       {}
117 
118       template< class LFFather, class LFSon, class LocalGeometry >
prolongLocal(const LFFather & lfFather,LFSon & lfSon,const LocalGeometry & geometryInFather,bool initialize) const119       void prolongLocal ( const LFFather &lfFather, LFSon &lfSon,
120                           const LocalGeometry &geometryInFather, bool initialize ) const
121       {
122         lfSon.clear();
123 
124         typedef typename LFSon :: EntityType  EntityType ;
125         typedef typename EntityType :: Geometry   Geometry;
126         const EntityType& sonEntity = lfSon.entity();
127         const Geometry sonGeo = sonEntity.geometry();
128 
129         QuadratureType quad( sonEntity, 2*lfSon.order()+1 );
130         const int nop = quad.nop();
131         for( int qp = 0; qp < nop; ++qp )
132         {
133           RangeFieldType quadWeight = quad.weight( qp );
134 
135           // in case of non-orthonormal basis we have to
136           // apply the integration element and the
137           // inverse mass matrix later
138           if( applyInverse )
139           {
140             quadWeight *= sonGeo.integrationElement( quad.point(qp) );
141           }
142 
143           RangeType value;
144           lfFather.evaluate( geometryInFather.global( quad.point( qp ) ), value );
145           value *= quadWeight;
146           lfSon.axpy( quad[ qp ], value );
147         }
148 
149         if( applyInverse )
150         {
151           localMassMatrix_.applyInverse( sonEntity, lfSon );
152         }
153       }
154 
needCommunication() const155       bool needCommunication () const { return true; }
156 
157     protected:
158       LocalMassMatrixType localMassMatrix_;
159       DomainFieldType weight_;
160       mutable TemporaryLocalFunction< DiscreteFunctionSpace > temp_;
161     };
162 
163 
164 
165     // DefaultLocalRestrictProlong for DiscontinuousGalerkinSpace
166     // ----------------------------------------------------------
167 
168     template< class FunctionSpaceImp, class GridPartImp, int polOrd, class StorageImp >
169     class DefaultLocalRestrictProlong< DiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrd, StorageImp > >
170     : public DiscontinuousGalerkinLocalRestrictProlong< DiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrd, StorageImp >, false >
171     {
172     public:
173       typedef DiscontinuousGalerkinLocalRestrictProlong< DiscontinuousGalerkinSpace<
174         FunctionSpaceImp, GridPartImp, polOrd, StorageImp >, false  >  BaseType;
DefaultLocalRestrictProlong(const DiscontinuousGalerkinSpace<FunctionSpaceImp,GridPartImp,polOrd,StorageImp> & space)175       DefaultLocalRestrictProlong ( const DiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrd, StorageImp > & space )
176         : BaseType( space )
177       {}
178     };
179 
180     template< class FunctionSpaceImp, class GridPartImp, class StorageImp >
181     class DefaultLocalRestrictProlong< DiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, 0, StorageImp > >
182     : public ConstantLocalRestrictProlong< DiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, 0, StorageImp > >
183     {
184     public:
DefaultLocalRestrictProlong(const DiscontinuousGalerkinSpace<FunctionSpaceImp,GridPartImp,0,StorageImp> &)185       DefaultLocalRestrictProlong ( const DiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, 0, StorageImp > & )
186       {}
187     };
188 
189 
190 
191     // DefaultLocalRestrictProlong for LegendreDiscontinuousGalerkinSpace
192     // ------------------------------------------------------------------
193 
194     template< class FunctionSpaceImp, class GridPartImp, int polOrd, class StorageImp >
195     class DefaultLocalRestrictProlong< LegendreDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrd, StorageImp > >
196     : public DiscontinuousGalerkinLocalRestrictProlong< LegendreDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrd, StorageImp >, false >
197     {
198     public:
199       typedef DiscontinuousGalerkinLocalRestrictProlong<
200         LegendreDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrd,  StorageImp >, false > BaseType;
DefaultLocalRestrictProlong(const LegendreDiscontinuousGalerkinSpace<FunctionSpaceImp,GridPartImp,polOrd,StorageImp> & space)201       DefaultLocalRestrictProlong ( const LegendreDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrd, StorageImp > & space )
202         : BaseType( space )
203       {}
204     };
205 
206     template< class FunctionSpaceImp, class GridPartImp, class StorageImp >
207     class DefaultLocalRestrictProlong< LegendreDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, 0, StorageImp > >
208     : public ConstantLocalRestrictProlong< LegendreDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, 0, StorageImp > >
209     {
210     public:
DefaultLocalRestrictProlong(const LegendreDiscontinuousGalerkinSpace<FunctionSpaceImp,GridPartImp,0,StorageImp> &)211       DefaultLocalRestrictProlong ( const LegendreDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, 0, StorageImp > & )
212       {}
213     };
214 
215 
216 
217     // DefaultLocalRestrictProlong for HierarchicLegendreDiscontinuousGalerkinSpace
218     // ----------------------------------------------------------------------------
219 
220     template< class FunctionSpaceImp, class GridPartImp, int polOrd, class StorageImp >
221     class DefaultLocalRestrictProlong< HierarchicLegendreDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrd, StorageImp > >
222     : public DiscontinuousGalerkinLocalRestrictProlong< HierarchicLegendreDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrd, StorageImp >, false >
223     {
224     public:
225       typedef DiscontinuousGalerkinLocalRestrictProlong<
226         HierarchicLegendreDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrd,  StorageImp >, false > BaseType;
DefaultLocalRestrictProlong(const HierarchicLegendreDiscontinuousGalerkinSpace<FunctionSpaceImp,GridPartImp,polOrd,StorageImp> & space)227       DefaultLocalRestrictProlong ( const HierarchicLegendreDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrd, StorageImp > & space )
228         : BaseType( space )
229       {}
230     };
231 
232     template< class FunctionSpaceImp, class GridPartImp, class StorageImp >
233     class DefaultLocalRestrictProlong< HierarchicLegendreDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, 0, StorageImp > >
234     : public ConstantLocalRestrictProlong< HierarchicLegendreDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, 0, StorageImp > >
235     {
236     public:
DefaultLocalRestrictProlong(const HierarchicLegendreDiscontinuousGalerkinSpace<FunctionSpaceImp,GridPartImp,0,StorageImp> &)237       DefaultLocalRestrictProlong ( const HierarchicLegendreDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, 0, StorageImp > & )
238       {}
239     };
240 
241 
242 
243     // DefaultLocalRestrictProlong for LagrangeDiscontinuousGalerkinSpace
244     // ------------------------------------------------------------------
245 
246     template< class FunctionSpaceImp, class GridPartImp, int polOrd, class StorageImp >
247     class DefaultLocalRestrictProlong< LagrangeDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrd, StorageImp > >
248     : public DiscontinuousGalerkinLocalRestrictProlong< LagrangeDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrd, StorageImp >, true >
249     {
250     public:
251       typedef DiscontinuousGalerkinLocalRestrictProlong< LagrangeDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrd, StorageImp >, true >  BaseType ;
DefaultLocalRestrictProlong(const LagrangeDiscontinuousGalerkinSpace<FunctionSpaceImp,GridPartImp,polOrd,StorageImp> & space)252       DefaultLocalRestrictProlong ( const LagrangeDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrd, StorageImp > & space )
253         : BaseType( space )
254       {}
255     };
256 
257     template< class FunctionSpaceImp, class GridPartImp, class StorageImp >
258     class DefaultLocalRestrictProlong< LagrangeDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, 0, StorageImp > >
259     : public ConstantLocalRestrictProlong< LagrangeDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, 0, StorageImp > >
260     {
261     public:
DefaultLocalRestrictProlong(const LagrangeDiscontinuousGalerkinSpace<FunctionSpaceImp,GridPartImp,0,StorageImp> &)262       DefaultLocalRestrictProlong ( const LagrangeDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, 0, StorageImp > & )
263       {}
264     };
265 
266     ///@}
267 
268   } // namespace Fem
269 
270 } // namespace Dune
271 
272 #endif // #ifndef DUNE_FEM_SPACE_DISCONTINUOUSGALERKIN_LOCALRESTRICTPROLONG_HH_
273