1 #ifndef DUNE_FEM_SPACE_LAGRANGE_SPACE_HH
2 #define DUNE_FEM_SPACE_LAGRANGE_SPACE_HH
3 
4 // C++ includes
5 #include <vector>
6 
7 // dune-common includes
8 #include <dune/common/exceptions.hh>
9 
10 // dune-geometry includes
11 #include <dune/geometry/type.hh>
12 
13 // dune-fem includes
14 #include <dune/fem/common/hybrid.hh>
15 
16 #include <dune/fem/space/basisfunctionset/default.hh>
17 #include <dune/fem/space/common/localinterpolation.hh>
18 #include <dune/fem/space/common/allgeomtypes.hh>
19 #include <dune/fem/space/common/basesetlocalkeystorage.hh>
20 #include <dune/fem/space/common/defaultcommhandler.hh>
21 #include <dune/fem/space/common/discretefunctionspace.hh>
22 #include <dune/fem/space/common/functionspace.hh>
23 #include <dune/fem/space/mapper/indexsetdofmapper.hh>
24 #include <dune/fem/space/mapper/nonblockmapper.hh>
25 #include <dune/fem/space/shapefunctionset/proxy.hh>
26 #include <dune/fem/space/shapefunctionset/selectcaching.hh>
27 #include <dune/fem/space/shapefunctionset/vectorial.hh>
28 #include <dune/fem/storage/singletonlist.hh>
29 
30 // local includes
31 #include "adaptmanager.hh"
32 #include "capabilities.hh"
33 #include "interpolation.hh"
34 #include "lagrangepoints.hh"
35 #include "shapefunctionset.hh"
36 #include "storage.hh"
37 
38 
39 namespace Dune
40 {
41 
42   namespace Fem
43   {
44 
45     // Forward declaration
46     // -------------------
47 
48     template< class FunctionSpace, class GridPart, int maxPolOrder, class Storage = CachingStorage >
49     class LagrangeDiscreteFunctionSpace;
50 
51 
52 
53     // LagrangeDiscreteFunctionSpaceTraits
54     // -----------------------------------
55 
56     template< class FunctionSpace, class GridPart, int maxPolOrder, class Storage >
57     struct LagrangeDiscreteFunctionSpaceTraits
58     {
59       typedef LagrangeDiscreteFunctionSpace< FunctionSpace, GridPart, maxPolOrder, Storage > DiscreteFunctionSpaceType;
60 
61       // we take 6 as the maximal polynomial order for the dynamic space
62       static const int maxPolynomialOrder = ( maxPolOrder < 0 ) ? -maxPolOrder : maxPolOrder;
63       static const int minPolynomialOrder = ( maxPolOrder < 0 ) ? 1 : maxPolOrder;
64 
65       typedef GridPart GridPartType;
66       typedef GridFunctionSpace< GridPartType, FunctionSpace > FunctionSpaceType;
67 
68       typedef IndexSetDofMapper< GridPartType,
69                                  std::conditional_t< Dune::Capabilities::isCartesian< typename GridPart::GridType >::v,
70                                                      DefaultLocalDofMapping< GridPart >,
71                                                      LagrangeLocalDofMapping< GridPart > > > BlockMapperType;
72 
73       typedef Hybrid::IndexRange< int, FunctionSpaceType::dimRange > LocalBlockIndices;
74 
75       static const int codimension = 0;
76 
77     private:
78       typedef typename GridPartType::template Codim< codimension >::EntityType EntityType;
79       static const int dimLocal = GridPartType::dimension;
80       typedef typename FunctionSpaceType::ScalarFunctionSpaceType ScalarFunctionSpaceType;
81       typedef typename ToNewDimDomainFunctionSpace< ScalarFunctionSpaceType, dimLocal >::Type ShapeFunctionSpaceType;
82 
83     public:
84       typedef LagrangeShapeFunctionSet< ShapeFunctionSpaceType, maxPolynomialOrder > LagrangeShapeFunctionSetType;
85       typedef SelectCachingShapeFunctionSet< LagrangeShapeFunctionSetType, Storage > ScalarShapeFunctionSetType;
86 
87       typedef ShapeFunctionSetProxy< ScalarShapeFunctionSetType > ScalarShapeFunctionSetProxyType;
88       typedef VectorialShapeFunctionSet< ScalarShapeFunctionSetProxyType, typename FunctionSpaceType::RangeType > ShapeFunctionSetType;
89 
90       typedef Dune::Fem::DefaultBasisFunctionSet< EntityType, ShapeFunctionSetType > BasisFunctionSetType;
91 
92       template< class DiscreteFunction, class Operation = Dune::Fem::DFCommunicationOperation::Add >
93       struct CommDataHandle
94       {
95         // type of data handle
96         typedef Dune::Fem::DefaultCommunicationHandler< DiscreteFunction, Operation > Type;
97         // type of operation to perform on scatter
98         typedef Operation OperationType;
99       };
100     };
101 
102 
103 
104     /** \addtogroup LagrangeDiscreteFunctionSpace
105      *
106      *  Provides access to bse function sets for different element types in
107      *  one grid and size of function space and maps from local to global dof
108      *  number.
109      *
110      *  \note This space can only be used with special index sets. If you want
111      *  to use the LagrangeDiscreteFunctionSpace with an index set only
112      *  supporting the index set interface you will have to use the
113      *  IndexSetWrapper class to provide the required functionality.
114      *
115      *  \note For adaptive calculations one has to use index sets that are
116      *  capable of adaption (i.e. the method adaptive returns true). See also
117      *  AdaptiveLeafIndexSet.
118      */
119 
120     // LagrangeDiscreteFunctionSpace
121     // -----------------------------
122 
123     /** \class   LagrangeDiscreteFunctionSpace
124      *  \ingroup LagrangeDiscreteFunctionSpace
125      *  \brief   Lagrange discrete function space
126      */
127 
128     template< class FunctionSpace, class GridPart, int maxPolOrder, class Storage >
129     class LagrangeDiscreteFunctionSpace
130     : public DiscreteFunctionSpaceDefault< LagrangeDiscreteFunctionSpaceTraits< FunctionSpace, GridPart, maxPolOrder, Storage > >
131     {
132 
133       typedef LagrangeDiscreteFunctionSpace< FunctionSpace, GridPart, maxPolOrder, Storage > ThisType;
134       typedef DiscreteFunctionSpaceDefault< LagrangeDiscreteFunctionSpaceTraits< FunctionSpace, GridPart, maxPolOrder, Storage > > BaseType;
135 
136     public:
137       typedef typename BaseType::Traits Traits;
138       static const int maxPolynomialOrder = Traits :: maxPolynomialOrder;
139       static const int minPolynomialOrder = Traits :: minPolynomialOrder;
140 
141       static_assert( (maxPolynomialOrder > 0), "LagrangeDiscreteFunctionSpace only defined for polOrder > 0" );
142 
143       typedef typename BaseType::FunctionSpaceType FunctionSpaceType;
144 
145       typedef typename BaseType::GridPartType GridPartType;
146       typedef typename BaseType::GridType GridType;
147       typedef typename BaseType::IndexSetType IndexSetType;
148       typedef typename BaseType::EntityType EntityType;
149       typedef typename BaseType::IntersectionType IntersectionType;
150 
151       typedef typename BaseType::Traits::ShapeFunctionSetType ShapeFunctionSetType;
152       typedef typename BaseType::BasisFunctionSetType BasisFunctionSetType;
153 
154       typedef typename BaseType::BlockMapperType BlockMapperType;
155 
156       typedef LagrangePointSet< GridPartType, maxPolynomialOrder > LagrangePointSetType;
157       typedef LagrangeLocalInterpolation< GridPartType, maxPolynomialOrder, BasisFunctionSetType >  LocalInterpolationType;
158       typedef LocalInterpolationType InterpolationImplType;
159 
160       typedef LocalInterpolationWrapper< ThisType >  InterpolationType;
161 
162     private:
163       typedef typename Traits::ScalarShapeFunctionSetType ScalarShapeFunctionSetType;
164       typedef BaseSetLocalKeyStorage< ScalarShapeFunctionSetType > ScalarShapeFunctionSetStorageType;
165 
166       typedef CompiledLocalKeyContainer< LagrangePointSetType, minPolynomialOrder, maxPolynomialOrder > LagrangePointSetContainerType;
167       typedef typename LagrangePointSetContainerType::LocalKeyStorageType LocalKeyStorageType;
168 
169       typedef LagrangeMapperSingletonKey< GridPartType, LocalKeyStorageType >
170         MapperSingletonKeyType;
171       typedef LagrangeMapperSingletonFactory< MapperSingletonKeyType, BlockMapperType >
172         BlockMapperSingletonFactoryType;
173       typedef SingletonList< MapperSingletonKeyType, BlockMapperType, BlockMapperSingletonFactoryType > BlockMapperProviderType;
174 
175       // static const InterfaceType defaultInterface = InteriorBorder_InteriorBorder_Interface;
176       static const InterfaceType defaultInterface = GridPart::indexSetInterfaceType;
177       static const CommunicationDirection defaultDirection = ForwardCommunication;
178 
179 
180       template < int pOrd >
181       struct Initialize
182       {
183         struct ScalarShapeFunctionSetFactory
184         {
185         public:
createObjectDune::Fem::LagrangeDiscreteFunctionSpace::Initialize::ScalarShapeFunctionSetFactory186           static ScalarShapeFunctionSetType *createObject ( const GeometryType &type )
187           {
188             typedef typename Traits :: LagrangeShapeFunctionSetType  LagrangeShapeFunctionSetType;
189             return new ScalarShapeFunctionSetType( type, LagrangeShapeFunctionSetType( type, pOrd ) );
190           }
191 
deleteObjectDune::Fem::LagrangeDiscreteFunctionSpace::Initialize::ScalarShapeFunctionSetFactory192           static void deleteObject ( ScalarShapeFunctionSetType *object ) { delete object; }
193         };
194 
applyDune::Fem::LagrangeDiscreteFunctionSpace::Initialize195         static void apply ( ScalarShapeFunctionSetStorageType& scalarShapeFunctionSets,
196                             const int polynomialOrder,
197                             const GeometryType &type )
198         {
199           if( polynomialOrder == pOrd )
200           {
201             typedef SingletonList< const GeometryType, ScalarShapeFunctionSetType, ScalarShapeFunctionSetFactory > SingletonProviderType;
202             scalarShapeFunctionSets.template insert< SingletonProviderType >( type );
203           }
204         }
205       };
206 
207 
208     public:
209       ///////////////////////
210       // Interface methods //
211       ///////////////////////
212 
213       using BaseType::order;
214 
LagrangeDiscreteFunctionSpace(GridPartType & gridPart,const InterfaceType commInterface,const CommunicationDirection commDirection)215       explicit LagrangeDiscreteFunctionSpace ( GridPartType &gridPart,
216                                                const InterfaceType commInterface,
217                                                const CommunicationDirection commDirection )
218       : LagrangeDiscreteFunctionSpace( gridPart, minPolynomialOrder, commInterface, commDirection )
219       {}
220 
LagrangeDiscreteFunctionSpace(GridPartType & gridPart,const int polOrder=minPolynomialOrder,const InterfaceType commInterface=defaultInterface,const CommunicationDirection commDirection=defaultDirection)221       explicit LagrangeDiscreteFunctionSpace ( GridPartType &gridPart,
222                                                const int polOrder = minPolynomialOrder,
223                                                const InterfaceType commInterface = defaultInterface,
224                                                const CommunicationDirection commDirection = defaultDirection )
225       : BaseType( gridPart, commInterface, commDirection ),
226         blockMapper_(),
227         lagrangePointSetContainer_( gridPart ),
228         // when min == max polynomial order the dynamic choice is off
229         polynomialOrder_( ( minPolynomialOrder == maxPolynomialOrder ) ? minPolynomialOrder : polOrder )
230       {
231         const IndexSetType &indexSet = gridPart.indexSet();
232 
233         AllGeomTypes< IndexSetType, GridType > allGeometryTypes( indexSet );
234         const std::vector< GeometryType > &geometryTypes = allGeometryTypes.geomTypes( 0 );
235 
236         // create shape function sets
237         for( unsigned int i = 0; i < geometryTypes.size(); ++i )
238         {
239           Fem::ForLoop< Initialize, minPolynomialOrder, maxPolynomialOrder >::
240             apply( scalarShapeFunctionSets_, polynomialOrder_, geometryTypes[ i ] );
241         }
242 
243         MapperSingletonKeyType key( gridPart, lagrangePointSetContainer_.compiledLocalKeys( polynomialOrder_ ), polynomialOrder_ );
244         blockMapper_.reset( &BlockMapperProviderType::getObject( key ) );
245         assert( blockMapper_ );
246       }
247 
248       /** \copydoc Dune::Fem::DiscreteFunctionSpaceInterface::type */
type() const249       DFSpaceIdentifier type () const
250       {
251         return LagrangeSpace_id;
252       }
253 
254       /** \copydoc Dune::Fem::DiscreteFunctionSpaceInterface::basisFunctionSet */
basisFunctionSet(const EntityType & entity) const255       const BasisFunctionSetType basisFunctionSet ( const EntityType &entity ) const
256       {
257         return BasisFunctionSetType( entity, shapeFunctionSet( entity ) );
258       }
259 
260       /** \copydoc Dune::Fem::DiscreteFunctionSpaceInterface::continuous */
continuous() const261       bool continuous () const
262       {
263         return true;
264       }
265 
266       /** \copydoc Dune::Fem::DiscreteFunctionSpaceInterface::continuous */
continuous(const IntersectionType & intersection) const267       bool continuous ( const IntersectionType &intersection ) const
268       {
269         return intersection.conforming();
270       }
271 
272       /** \copydoc Dune::Fem::DiscreteFunctionSpaceInterface::order */
order() const273       int order () const
274       {
275         return polynomialOrder_;
276       }
277 
278       /** \copydoc Dune::Fem::DiscreteFunctionSpaceInterface::blockMapper */
blockMapper() const279       BlockMapperType &blockMapper () const
280       {
281         assert( blockMapper_ );
282         return *blockMapper_;
283       }
284 
285       ///////////////////////////
286       // Non-interface methods //
287       ///////////////////////////
288 
289       /** \brief return interpolation object
290        */
interpolation() const291       InterpolationType interpolation () const
292       {
293         return InterpolationType( *this );
294       }
295 
296       /** \brief return local interpolation for given entity
297        *
298        *  \param[in]  entity  grid part entity
299        */
300       [[deprecated]]
interpolation(const EntityType & entity) const301       LocalInterpolationType interpolation ( const EntityType &entity ) const
302       {
303         return LocalInterpolationType( lagrangePointSetContainer_.compiledLocalKey( entity.type(), polynomialOrder_ ), basisFunctionSet( entity ) );
304       }
305 
306       /** \brief return local interpolation for given entity
307        *
308        *  \param[in]  entity  grid part entity
309        */
localInterpolation(const EntityType & entity) const310       LocalInterpolationType localInterpolation ( const EntityType &entity ) const
311       {
312         return LocalInterpolationType( lagrangePointSetContainer_.compiledLocalKey( entity.type(), polynomialOrder_ ), basisFunctionSet( entity ) );
313       }
314 
315       /** \brief return shape function set for given entity
316        *
317        * \param[in]  entity  entity (of codim 0) for which shape function set
318        *                     is requested
319        *
320        * \returns  ShapeFunctionSetType  shape function set
321        */
shapeFunctionSet(const EntityType & entity) const322       ShapeFunctionSetType shapeFunctionSet ( const EntityType &entity ) const
323       {
324         return shapeFunctionSet( entity.type() );
325       }
326 
327       /** \brief return shape unique function set for geometry type
328        *
329        * \param[in]  type  geometry type (must be a cube) for which
330        *                   shape function set is requested
331        *
332        * \returns  ShapeFunctionSetType  shape function set
333        */
shapeFunctionSet(const GeometryType & type) const334       ShapeFunctionSetType shapeFunctionSet ( const GeometryType &type ) const
335       {
336         return ShapeFunctionSetType( &scalarShapeFunctionSets_[ type ] );
337       }
338 
339       LagrangeDiscreteFunctionSpace ( const ThisType & ) = delete;
340       ThisType &operator= ( const ThisType & ) = delete;
341 
342     protected:
343       std::unique_ptr< BlockMapperType, typename BlockMapperProviderType::Deleter> blockMapper_;
344       ScalarShapeFunctionSetStorageType scalarShapeFunctionSets_;
345       LagrangePointSetContainerType lagrangePointSetContainer_;
346       const int polynomialOrder_;
347     };
348 
349   } // namespace Fem
350 
351 } // namespace Dune
352 
353 #endif // #ifndef DUNE_FEM_SPACE_LAGRANGE_SPACE_HH
354