1 #include <config.h>
2
3 #include <iostream>
4
5 using namespace Dune;
6
7 #include <dune/fem/function/adaptivefunction.hh>
8 #include <dune/fem/function/petscdiscretefunction.hh>
9 #include <dune/fem/function/vectorfunction.hh>
10 #include <dune/fem/function/blockvectorfunction.hh>
11 #include <dune/fem/function/blockvectordiscretefunction.hh>
12 #include <dune/fem/space/discontinuousgalerkin.hh>
13 #include <dune/fem/space/padaptivespace/padaptivespace.hh>
14 #include <dune/fem/space/hpdg/orthogonal.hh>
15
16 #if USE_COMBINED_SPACE
17 #include <dune/fem/space/combinedspace.hh>
18 #endif
19
20 #include <dune/fem/quadrature/cachingquadrature.hh>
21 #include <dune/fem/gridpart/adaptiveleafgridpart.hh>
22
23 #include <dune/fem/misc/l2norm.hh>
24 #include <dune/fem/misc/capabilities.hh>
25
26 #if HAVE_GRAPE && WANT_GRAPE && GRIDDIM > 1
27 #define USE_GRAPE 1
28 #else
29 #define USE_GRAPE 0
30 #endif
31
32 #if USE_GRAPE && GRIDDIM > 1
33 #include <dune/grid/io/visual/grapedatadisplay.hh>
34 #endif
35
36 #include <dune/fem/io/parameter.hh>
37
38 #include <dune/fem/test/testgrid.hh>
39
40 // polynom approximation order of quadratures,
41 // at least poolynom order of basis functions
42 static const int polOrd = POLORDER;
43 static const int dimRange = 3 ;
44
45 using namespace Fem;
46
47 //***********************************************************************
48 /*! L2 Projection of a function f:
49 */
50 //***********************************************************************
51
52 //! the index set we are using
53 typedef GridSelector::GridType MyGridType;
54 typedef DGAdaptiveLeafGridPart< MyGridType > GridPartType;
55
56 //! define the function space, \f[ \R^2 \rightarrow \R \f]
57 // see dune/common/functionspace.hh
58 // typedef MatrixFunctionSpace < double , double, MyGridType::dimensionworld, 2,5 > FuncSpace;
59
60 #ifdef USE_COMBINED_SPACE
61 typedef DiscontinuousGalerkinSpace< FunctionSpace < double , double, MyGridType::dimensionworld, 1 >,
62 GridPartType, polOrd, CachingStorage> ContainedDiscreteFunctionSpaceType;
63 #ifdef POINTBASED
64 static const std::string usingSpaceName("Using CombinedSpace< dimRange, PointBased >");
65 typedef CombinedSpace< ContainedDiscreteFunctionSpaceType, dimRange, PointBased > DiscreteFunctionSpaceType ;
66 #else
67 static const std::string usingSpaceName("Using CombinedSpace< dimRange, VariableBased >");
68 typedef CombinedSpace< ContainedDiscreteFunctionSpaceType, dimRange, VariableBased > DiscreteFunctionSpaceType ;
69 #endif
70 #else
71 static const std::string usingSpaceName("Using DiscontinuousGalerkinSpace< dimRange >");
72 typedef DiscontinuousGalerkinSpace< FunctionSpace < double , double, MyGridType::dimensionworld, dimRange >,
73 GridPartType, polOrd, CachingStorage> DiscreteFunctionSpaceType;
74 //typedef hpDG::OrthogonalDiscontinuousGalerkinSpace< FunctionSpace < double , double, MyGridType::dimensionworld, dimRange >,
75 // GridPartType, polOrd, true > DiscreteFunctionSpaceType;
76 //typedef PAdaptiveDGSpace< FunctionSpace < double , double, MyGridType::dimensionworld, dimRange >,
77 // GridPartType, polOrd, CachingStorage > DiscreteFunctionSpaceType;
78 #endif
79
80
81 ////typedef CombinedSpace< ContainedDiscreteFunctionSpaceType, dimRange, PointBased > DiscreteFunctionSpaceType ;
82
83 typedef typename DiscreteFunctionSpaceType :: FunctionSpaceType FunctionSpaceType;
84
85 //! define the type of discrete function we are using , see
86
87 #if HAVE_PETSC
88 // PetscDiscreteFunction uses AdaptiveDiscreteFunction for dof prolongation and
89 // resttriction
90 typedef PetscDiscreteFunction< DiscreteFunctionSpaceType > DiscreteFunctionType;
91 #else
92 typedef AdaptiveDiscreteFunction< DiscreteFunctionSpaceType > DiscreteFunctionType;
93 #endif
94
95
96 //typedef Dune::Fem::ManagedDiscreteFunction< VectorDiscreteFunction< DiscreteFunctionSpaceType, std::vector< double > > > DiscreteFunctionType;
97 //typedef ISTLBlockVectorDiscreteFunction< DiscreteFunctionSpaceType > DiscreteFunctionType;
98 //typedef Dune::Fem::ReferenceBlockVector< FunctionSpaceType::RangeFieldType,
99 // DiscreteFunctionSpaceType::localBlockSize > BlockVectorType;
100 //typedef Dune::Fem::BlockVectorDiscreteFunction< DiscreteFunctionSpaceType, BlockVectorType > DiscreteFunctionType;
101
102 typedef DofManager< MyGridType > DofManagerType;
103
104 typedef AdaptationManager< MyGridType, RestrictProlongDefault< DiscreteFunctionType > > AdaptationManagerType;
105
106 // ***********************************************************
107 // the exact solution to the problem for EOC calculation
108 struct ExactSolution
109 : public Fem::Function< FunctionSpaceType, ExactSolution >
110 {
111 typedef FunctionSpaceType::RangeType RangeType;
112 typedef FunctionSpaceType::RangeFieldType RangeFieldType;
113 typedef FunctionSpaceType::DomainType DomainType;
114
115 //! f(x,y) = x*(1-x)*y*(1-y)
evaluateExactSolution116 void evaluate ( const DomainType &x, RangeType &ret ) const
117 {
118 ret = 2.; // maximum of function is 2
119 for( int i = 0; i < DomainType::dimension; ++i )
120 ret *= sin( x[ i ]*(1.0 -x[ i ])*4.);
121 }
122
evaluateExactSolution123 void evaluate ( const DomainType &x, RangeFieldType time, RangeType &ret ) const
124 {
125 evaluate( x, ret );
126 }
127 };
128
129 // ********************************************************************
adapt(MyGridType & grid,DiscreteFunctionType & solution,int step)130 void adapt( MyGridType &grid, DiscreteFunctionType &solution, int step )
131 {
132 const DiscreteFunctionType::DiscreteFunctionSpaceType &space = solution.space();
133 RestrictProlongDefault<DiscreteFunctionType> rp(solution);
134 rp.setFatherChildWeight(DGFGridInfo< MyGridType >::refineWeight());
135
136 AdaptationManagerType adop(grid,rp);
137
138 std::string message;
139
140 int mark = 1;
141 int count = std::abs(step);
142
143 if(step < 0)
144 {
145 message += "Coarsening...";
146 mark = -1;
147 }
148 else
149 message += "Refining...";
150
151 if( Parameter :: verbose() )
152 {
153 std::cout << "Grid leaf size: " << grid.size( 0 ) << std::endl;
154 std::cout << "AdaptiveLeafIndexSet.size: " << space.indexSet().size( 0 ) << std::endl;
155 }
156
157 for( int i = 0; i < count; ++i )
158 {
159 for( const auto& entity : space)
160 grid.mark( mark, entity );
161 adop.adapt();
162 if( Parameter::verbose() )
163 std::cout << message << std::endl;
164 }
165
166 if( Parameter :: verbose() )
167 {
168 std::cout << "Grid leaf size: " << grid.size( 0 ) << std::endl;
169 std::cout << "AdaptiveLeafIndexSet.size: " << space.indexSet().size( 0 ) << std::endl;
170 }
171
172 }
173 // ********************************************************************
algorithm(MyGridType & grid,DiscreteFunctionType & solution,int step,int turn)174 double algorithm ( MyGridType &grid, DiscreteFunctionType &solution, int step, int turn )
175 {
176 const unsigned int order = solution.space().order();
177
178 DiscreteFunctionType tmp ( solution );
179
180 ExactSolution f;
181 auto gridFunc = gridFunctionAdapter(f, solution.space().gridPart(), 2);
182 {
183 Dune::Fem::interpolate( gridFunc, solution );
184 Dune :: Fem :: L2Norm< GridPartType > l2norm ( solution.space().gridPart(), 2*order+2 ) ;
185 double new_error = l2norm.distance( f ,solution );
186 if( Parameter::verbose() )
187 std::cout << "before ref." << new_error << "\n\n";
188 }
189
190 adapt(grid,solution,step);
191
192 // tmp solution should be zero after adapt
193 double tmpErr = tmp.normSquaredDofs();
194 if( tmpErr > 0 )
195 {
196 // return big error
197 return 1e80;
198 }
199
200 #if USE_GRAPE
201 // if Grape was found, then display last solution
202 if(turn > 0)
203 {
204 std::cerr << "GRAPE 1" << std::endl;
205 GrapeDataDisplay < MyGridType > grape(grid);
206 grape.dataDisplay( solution );
207 }
208 #endif
209
210 // calculation L2 error on refined grid
211 // pol ord for calculation the error chould by higher than
212 // pol for evaluation the basefunctions
213 Dune :: Fem :: L2Norm< GridPartType > l2norm ( solution.space().gridPart(), 2*order+2 ) ;
214 double error = l2norm.distance( gridFunc, solution );
215
216 #if USE_GRAPE
217 // if Grape was found, then display last solution
218 if(turn > 0)
219 {
220 std::cerr << "GRAPE 2" << std::endl;
221 GrapeDataDisplay< MyGridType > grape(grid);
222 grape.dataDisplay( solution );
223 }
224 #endif
225
226 //! perform l2-projection to refined grid
227 Dune::Fem::interpolate( gridFunc, solution );
228 double new_error = l2norm.distance( gridFunc, solution );
229 if( Parameter::verbose() )
230 {
231 std::cout << "\nL2 Error : " << error << " on new grid " << new_error << "\n\n";
232 }
233 #if USE_GRAPE
234 // if Grape was found, then display last solution
235 if(turn > 0)
236 {
237 std::cerr << "SIZE: " << solution.space().size()
238 << " GRID: " << grid.size(0) << std::endl;
239 std::cerr << "GRAPE 3" << std::endl;
240 GrapeDataDisplay< MyGridType > grape(grid);
241 grape.dataDisplay( solution );
242 }
243 #endif
244 return error;
245 }
246
247 //**************************************************
248 //
249 // main programm, run algorithm twice to calc EOC
250 //
251 //**************************************************
main(int argc,char * argv[])252 int main( int argc, char *argv[] )
253 try {
254 MPIManager :: initialize( argc, argv );
255
256 const char* paramName = "parameter";
257 if( argc < 2 )
258 {
259 std :: cerr << "Usage: " << argv[ 0 ] << "<parameter>" << std :: endl;
260 }
261 else
262 paramName = argv[1];
263
264 std::string paramFile( paramName );
265 std::cout << usingSpaceName << std::endl;
266
267 // append parameter
268 Parameter :: append( argc , argv );
269 Parameter :: append( paramFile );
270
271 int ml = 2 ; // default value = 2
272 ml = Parameter :: getValue ("lagrangeadapt.maxlevel", ml);
273
274 std::vector<double> error(ml);
275
276 MyGridType &grid = Dune::Fem::TestGrid::grid();
277 const int step = Dune::Fem::TestGrid::refineStepsForHalf();
278
279 GridPartType part ( grid );
280 DiscreteFunctionSpaceType space( part );
281
282 // threshold for EOC difference to predicted value
283 const double eocThreshold = Parameter :: getValue("adapt.eocthreshold", double(0.2) );
284
285 const bool isLocallyAdaptive = Dune::Fem::Capabilities::isLocallyAdaptive< MyGridType > :: v ;
286
287 DiscreteFunctionType solution ( "sol", space );
288 solution.clear();
289 if( Parameter::verbose() )
290 std::cout << "------------ Refining:" << std::endl;
291 for(int i=0; i<ml; i+=1)
292 {
293 error[i] = algorithm ( grid , solution, step, (i==ml-1));
294 if (i>0)
295 {
296 if ( isLocallyAdaptive )
297 {
298 double eoc = log( error[i-1]/error[i]) / M_LN2;
299 if( Parameter::verbose() )
300 std::cout << "EOC = " << eoc << std::endl;
301 if( std::abs( eoc - (space.order()+eocThreshold) ) < 0 )
302 {
303 DUNE_THROW(InvalidStateException,"EOC check of refinement failed");
304 }
305 }
306 else
307 std::cout << "no EOC for non-adaptive grid" << std::endl;
308 }
309 }
310 if( Parameter::verbose() )
311 std::cout << "------------ Coarsening:" << std::endl;
312 for(int i=ml-1; i>=0; i-=1)
313 {
314 error[i] = algorithm ( grid , solution,-step, 1);
315 if (i<ml-1)
316 {
317 if( isLocallyAdaptive )
318 {
319 double eoc = log( error[i+1]/error[i]) / M_LN2;
320 if( Parameter::verbose() )
321 std::cout << "EOC = " << eoc << std::endl;
322 if( std::abs( eoc + (space.order()+eocThreshold) ) < 0 )
323 {
324 DUNE_THROW(InvalidStateException,"EOC check of coarsening failed");
325 }
326 }
327 else
328 std::cout << "no EOC for non-adaptive grid" << std::endl;
329 }
330 }
331 return 0;
332 }
333 catch( const Dune :: Exception &exception )
334 {
335 std :: cerr << exception << std :: endl;
336 return 1;
337 }
338
339