1 #include <config.h>
2 
3 #include <iostream>
4 #include <sstream>
5 #include <string>
6 #include <tuple>
7 #include <vector>
8 
9 static const int dimw = Dune::GridSelector::dimworld;
10 
11 #include <dune/fem/quadrature/cachingquadrature.hh>
12 #include <dune/fem/space/combinedspace.hh>
13 #include <dune/fem/space/discontinuousgalerkin.hh>
14 #include <dune/fem/space/lagrange.hh>
15 
16 #include <dune/fem/gridpart/adaptiveleafgridpart.hh>
17 #include <dune/fem/gridpart/leafgridpart.hh>
18 
19 #include <dune/fem/misc/double.hh>
20 
21 #include <dune/fem/test/testgrid.hh>
22 
23 
24 // polynom approximation order of quadratures,
25 // at least polynom order of basis functions
26 const int polOrd = POLORDER;
27 
28 //***********************************************************************
29 /*! L2 Projection of a function f:
30 
31    This is an example how to solve the equation on
32    \f[\Omega = (0,1)^2 \f]
33 
34    \f[ \int_{\Omega} u \phi = \int_{\Omega} f \phi  \ \ \ in \Omega \f]
35    \f[ f(x,y) = x ( 1 - x) y ( 1 - y ) \f]
36 
37    Here u is the L_2 projection of f.
38 
39    The Projection should converge to the given function f.
40    with the finite element method using lagrangian elements of polynom order +1.
41  */
42 //***********************************************************************
43 
44 //! the index set we are using
45 typedef Dune::GridSelector::GridType MyGridType;
46 typedef Dune::Fem::AdaptiveLeafGridPart< MyGridType > GridPartType;
47 
48 // see dune/common/functionspace.hh
49 typedef Dune::Fem::FunctionSpace< MyGridType::ctype, double, dimw, 2 > FuncSpace1;
50 typedef Dune::Fem::FunctionSpace< MyGridType::ctype, double, dimw, 1 > FuncSpace2;
51 
52 typedef FuncSpace1::RangeType RangeType1;
53 typedef FuncSpace2::RangeType RangeType2;
54 typedef FuncSpace1::JacobianRangeType JacobianRangeType1;
55 typedef FuncSpace2::JacobianRangeType JacobianRangeType2;
56 typedef FuncSpace1::HessianRangeType HessianRangeType1;
57 typedef FuncSpace2::HessianRangeType HessianRangeType2;
58 
59 //! define the function space our unkown belong to
60 //! see dune/fem/lagrangebase.hh
61 typedef Dune::Fem::LagrangeDiscontinuousGalerkinSpace< FuncSpace1, GridPartType,
62                                             polOrd+1, Dune::Fem::CachingStorage > DiscreteFunctionSpaceType1;
63 
64 typedef Dune::Fem::LagrangeDiscontinuousGalerkinSpace< FuncSpace2, GridPartType,
65                                             polOrd, Dune::Fem::CachingStorage > DiscreteFunctionSpaceType2;
66 
67 typedef Dune::Fem::TupleDiscreteFunctionSpace< DiscreteFunctionSpaceType1, DiscreteFunctionSpaceType2 > DiscreteFunctionSpaceType;
68 
69 typedef DiscreteFunctionSpaceType::IteratorType::Entity EntityType;
70 
71 typedef DiscreteFunctionSpaceType::FunctionSpaceType FuncSpace;
72 typedef FuncSpace::RangeType RangeType;
73 typedef FuncSpace::JacobianRangeType JacobianRangeType;
74 typedef FuncSpace::HessianRangeType HessianRangeType;
75 
76 
77 const int dimRange = FuncSpace::dimRange;
78 const int dimRange1 = FuncSpace1::dimRange;
79 const int dimRange2 = FuncSpace2::dimRange;
80 const int dimDomain = FuncSpace::dimDomain;
81 
82 
checkBasisSet(const EntityType & entity,const DiscreteFunctionSpaceType & space)83 void checkBasisSet ( const EntityType &entity, const DiscreteFunctionSpaceType &space )
84 {
85   auto bSet = space.basisFunctionSet( entity );
86 
87   auto bSet1 = space.subDiscreteFunctionSpace< 0 >().basisFunctionSet( entity );
88   auto bSet2 = space.subDiscreteFunctionSpace< 1 >().basisFunctionSet( entity );
89 
90   const int nrBasis = bSet.size();
91   const int nrBasis1 = bSet1.size();
92   const int nrBasis2 = bSet2.size();
93 
94   std::vector< RangeType > ranges( nrBasis, RangeType( 0. ) );
95   std::vector< RangeType1 > ranges1( nrBasis1, RangeType1( 0. ) );
96   std::vector< RangeType2 > ranges2( nrBasis2, RangeType2( 0. ) );
97 
98   std::vector< JacobianRangeType > jacs( nrBasis, JacobianRangeType( 0. ) );
99   std::vector< JacobianRangeType1 > jacs1( nrBasis1, JacobianRangeType1( 0. ) );
100   std::vector< JacobianRangeType2 > jacs2( nrBasis2, JacobianRangeType2( 0. ) );
101 
102   HessianRangeType null;
103   HessianRangeType1 null1;
104   HessianRangeType2 null2;
105 
106   for( int i = 0; i < dimDomain; ++i )
107     for( int j = 0; j < dimDomain; ++j )
108     {
109       for( int r = 0; r < dimRange; ++r )
110         null[ r ][ i ][ j ] = .0;
111       for( int r = 0; r < dimRange1; ++r )
112         null1[ r ][ i ][ j ] = .0;
113       for( int r = 0; r < dimRange2; ++r )
114         null2[ r ][ i ][ j ] = .0;
115     }
116 
117   std::vector< HessianRangeType > hessians( nrBasis, null );
118   std::vector< HessianRangeType1 > hessians1( nrBasis1, null1 );
119   std::vector< HessianRangeType2 > hessians2( nrBasis2, null2 );
120 
121   Dune::Fem::CachingQuadrature< GridPartType, 0 > quadrature( entity, space.order() + 1 );
122   for( const auto& qp : quadrature )
123   {
124 
125     bSet.evaluateAll( qp, ranges );
126     bSet1.evaluateAll( qp, ranges1 );
127     bSet2.evaluateAll( qp, ranges2 );
128 
129     bSet.jacobianAll( qp, jacs );
130     bSet1.jacobianAll( qp, jacs1 );
131     bSet2.jacobianAll( qp, jacs2 );
132 
133     bSet.hessianAll( qp, hessians );
134     bSet1.hessianAll( qp, hessians1 );
135     bSet2.hessianAll( qp, hessians2 );
136 
137     RangeType value( 0 );
138     JacobianRangeType jac( 0 );
139     HessianRangeType hessian( null );
140 
141     for( int i = 0; i < nrBasis; ++i )
142     {
143       for( int r = 0; r < dimRange; ++r )
144       {
145         if( r < dimRange1 )
146         {
147           for( int j = 0; j < dimDomain; ++j )
148           {
149             jac[ r ][ j ] = ( i < nrBasis1 ) ? jacs1[ i ][ r ][ j ] : 0.;
150 
151             for( int k = 0; k < dimDomain; ++k )
152               hessian[ r ][ j ][ k ] = ( i < nrBasis1 ) ? hessians1[ i ][ r ][ j ][ k ] : 0.;
153           }
154 
155           value[ r ] = ( i < nrBasis1 ) ? ranges1[ i ][ r ] : 0.;
156         }
157         else
158         {
159           value[ r ] = ( i < nrBasis1 ) ? 0. : ranges2[ i - nrBasis1 ][ r - dimRange1 ];
160 
161           for( int j = 0; j < dimDomain; ++j )
162           {
163             jac[ r ][ j ] = ( i < nrBasis1 ) ? 0. : jacs2[ i-nrBasis1 ][ r - dimRange1 ][ j ];
164 
165             for( int k = 0; k < dimDomain; ++k )
166               hessian[ r ][ j ][ k ] = ( i < nrBasis1 ) ? 0. : hessians2[ i - nrBasis1 ][ r - dimRange1 ][ j ][ k ];
167           }
168         }
169       }
170 
171       // check function value
172       value -= ranges[ i ];
173       if( value.two_norm() > 1e-8 )
174         DUNE_THROW( Dune::InvalidStateException, "Basisfunction::evaluate returns wrong value." );
175 
176       JacobianRangeType jacHelp( jac );
177       // check jac
178       jacHelp -= jacs[ i ];
179       if( jacHelp.frobenius_norm() > 1e-8 )
180         DUNE_THROW( Dune::InvalidStateException, "Basisfunction::jacobian returns wrong value." );
181 
182       for( int r = 0; r < dimRange; ++r )
183         hessian[ r ] -= hessians[ i ][ r ];
184 
185       double val = 0;
186       for( int r = 0; r < dimRange; ++r )
187         val += hessian[ r ].frobenius_norm2();
188 
189       if( std::sqrt( val ) > 1e-8 )
190         DUNE_THROW( Dune::InvalidStateException, "Basisfunction::hessian returns wrong value." );
191     }
192   }
193 }
194 
195 
196 struct DummyLocalFunction
197 {
198   typedef FuncSpace FunctionSpaceType;
199 
200   typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
201   typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
202 
203   typedef typename FunctionSpaceType::DomainType DomainType;
204   typedef typename FunctionSpaceType::RangeType RangeType;
205   typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
206   typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
207 
208   static const int dimRange = FunctionSpaceType::dimRange;
209   static const int dimDomain = FunctionSpaceType::dimDomain;
210 
211   typedef typename DiscreteFunctionSpaceType::EntityType EntityType;
212 
DummyLocalFunctionDummyLocalFunction213   DummyLocalFunction ( const EntityType &entity ) : entity_( entity ){}
214 
215   template< class Point >
evaluateDummyLocalFunction216   void evaluate ( const Point &arg, RangeType &dest ) const
217   {
218     dest = 1.;
219   }
220 
221   template< class Point >
jacobianDummyLocalFunction222   void jacobian ( const Point &arg, JacobianRangeType &jac ) const
223   {
224     jac = 0;
225   }
226 
227   template< class Point >
hessianDummyLocalFunction228   void hessian ( const Point &arg, HessianRangeType &hes ) const
229   {}
230 
orderDummyLocalFunction231   int order () const { return 1; }
entityDummyLocalFunction232   const EntityType &entity () const { return entity_; }
233 
initDummyLocalFunction234   void init ( const EntityType &entity ) {}
235 
236 protected:
237   const EntityType &entity_;
238 };
239 
240 
checkInterpolation(const EntityType & entity,const DiscreteFunctionSpaceType & space)241 void checkInterpolation ( const EntityType &entity, const DiscreteFunctionSpaceType &space )
242 {
243   std::vector< double > ldv;
244   ldv.resize( space.basisFunctionSet( entity ).size() );
245 
246   DummyLocalFunction lf( entity );
247   space.interpolation( entity ) ( lf, ldv );
248 
249   //TODO check wether interpolation is successfull
250 }
251 
252 
main(int argc,char ** argv)253 int main ( int argc, char **argv )
254 {
255   Dune::Fem::MPIManager::initialize( argc, argv );
256   try
257   {
258     MyGridType &grid = Dune::Fem::TestGrid::grid();
259 
260     GridPartType part( grid );
261     DiscreteFunctionSpaceType space( part );
262 
263     EntityType entity = *(++space.begin());
264 
265     checkBasisSet( entity, space );
266     checkInterpolation( entity, space );
267 
268     return 0;
269   }
270   catch( const Dune::Exception &exception )
271   {
272     std::cerr << exception << std::endl;
273     return 1;
274   }
275 }
276 
277