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