1 // undef NDEBUG so we can always use assert.
2 #undef NDEBUG
3
4 #ifdef YASPGRID
5 #define LEAFGRIDPART
6 #endif
7
8 #include <config.h>
9 #include <iostream>
10
11 #include <dune/fem/gridpart/adaptiveleafgridpart.hh>
12 #include <dune/fem/gridpart/leafgridpart.hh>
13
14 #include <dune/fem/misc/double.hh>
15 #include <dune/fem/space/lagrange.hh>
16 #include <dune/fem/function/adaptivefunction.hh>
17 #include <dune/fem/function/common/gridfunctionadapter.hh>
18 #include <dune/fem/space/common/interpolate.hh>
19 #include <dune/fem/io/streams/binarystreams.hh>
20 #include <dune/fem/io/streams/asciistreams.hh>
21 #include <dune/fem/io/streams/virtualstreams.hh>
22 #include <dune/fem/io/file/vtkio.hh>
23
24 #include "testgrid.hh"
25 #include "exactsolution.hh"
26
27 using namespace Dune;
28 using namespace Fem;
29
30 typedef GridSelector::GridType MyGridType;
31
32 #ifdef LEAFGRIDPART
33 // AdaptiveLeafGridPart does not work with YaspGrid for p > 1, since not entities for
34 // codim > 0 and codim < dim are available, which are needed to build the index set
35 typedef LeafGridPart< MyGridType > GridPartType;
36 #else
37 typedef AdaptiveLeafGridPart< MyGridType > GridPartType;
38 #endif
39
40 typedef FunctionSpace< double, double, GridPartType::dimensionworld, 1 >
41 FunctionSpaceType;
42
43 typedef ExactSolution< FunctionSpaceType > ExactSolutionType;
44
45
46 template <class DiscreteFunctionType>
writeOut(VirtualOutStream out,const DiscreteFunctionType & solution)47 void writeOut ( VirtualOutStream out, const DiscreteFunctionType &solution )
48 {
49 out << solution;
50 out.flush();
51 }
52
53 template <class DiscreteFunctionType>
readBack(VirtualInStream in,DiscreteFunctionType & solution)54 void readBack ( VirtualInStream in, DiscreteFunctionType &solution )
55 {
56 solution.clear();
57 in >> solution;
58 }
59
60 template <class DiscreteFunctionSpaceType>
algorithm(MyGridType & grid,const int polOrder)61 int algorithm( MyGridType& grid, const int polOrder )
62 {
63 const int step = TestGrid :: refineStepsForHalf();
64 grid.globalRefine( 2*step );
65 const int rank = grid.comm().rank();
66
67 ////////////////////////////////////////////////////////
68 // create data structures (after grid has been refined)
69 ////////////////////////////////////////////////////////
70
71 GridPartType gridPart( grid );
72
73 DiscreteFunctionSpaceType discreteFunctionSpace( gridPart, polOrder );
74
75 typedef AdaptiveDiscreteFunction< DiscreteFunctionSpaceType > DiscreteFunctionType;
76 DiscreteFunctionType solution( "solution", discreteFunctionSpace );
77 solution.clear();
78
79 // interpolate
80 interpolate( gridFunctionAdapter( ExactSolutionType(), gridPart, discreteFunctionSpace.order() + 2 ), solution );
81
82 // let's check on IO
83 DiscreteFunctionType readback( "readback", discreteFunctionSpace );
84
85 std::string casename( "lagrangeinterpolation_r" + std::to_string(rank) + "_p" + std::to_string(polOrder) + "_" );
86
87 BinaryFileOutStream bout( casename + "solution-xdr.tmp" );
88 writeOut( virtualize( bout ), solution );
89
90 BinaryFileInStream bin( casename + "solution-xdr.tmp" );
91 readBack( virtualize( bin ), readback );
92 if( readback != solution )
93 {
94 std :: cerr << "binary read/write gives different function." << std :: endl;
95 assert( readback == solution );
96 return 1;
97 }
98
99 ASCIIOutStream aout( casename + "solution-ascii.tmp" );
100 writeOut( virtualize( aout ), solution );
101
102 ASCIIInStream ain( casename + "solution-ascii.tmp" );
103 readBack( virtualize( ain ), readback );
104 assert( readback == solution );
105 if( readback != solution )
106 {
107 std :: cerr << "ascii read/write gives different function." << std :: endl;
108 assert( readback == solution );
109 return 1;
110 }
111
112 // output to vtk file
113 VTKIO<GridPartType> vtkWriter(gridPart);
114 vtkWriter.addVertexData(solution);
115 vtkWriter.pwrite(casename+"lagrangeinterpol",
116 Parameter::commonOutputPath().c_str(),"",
117 Dune::VTK::ascii);
118
119 // reset grid
120 grid.globalRefine( -2*step );
121 return 0;
122 }
123
run(MyGridType & grid,const int polOrder)124 void run( MyGridType& grid, const int polOrder )
125 {
126 // test dynamic Lagrange space
127 {
128 typedef DynamicLagrangeDiscreteFunctionSpace< FunctionSpaceType, GridPartType > DiscreteFunctionSpaceType;
129 algorithm< DiscreteFunctionSpaceType >(grid, polOrder);
130 }
131
132 #if HAVE_DUNE_LOCALFUNCTIONS
133 // test Lagrange space using LocalFiniteElementSpace
134 {
135 typedef LagrangeSpace< FunctionSpaceType, GridPartType > DiscreteFunctionSpaceType;
136 algorithm< DiscreteFunctionSpaceType >(grid, polOrder);
137 }
138 #endif
139 }
140
141
main(int argc,char ** argv)142 int main(int argc, char ** argv)
143 {
144 MPIManager :: initialize( argc, argv );
145 try
146 {
147 // add command line parameters to global parameter table
148 Dune::Fem::Parameter::append( argc, argv );
149 // append parameters from the parameter file
150 Dune::Fem::Parameter::append( (argc < 2) ? "parameter" : argv[ 1 ] );
151
152 MyGridType &grid = TestGrid :: grid();
153
154 // if parameter was specified then only run test for this pol order
155 // this is primarily for debugging
156 if( Dune::Fem::Parameter::exists("fem.lagrange.polynomialOrder") )
157 {
158 int p = Dune::Fem::Parameter::getValue< int >( "fem.lagrange.polynomialOrder");
159 run( grid, p );
160 }
161 else
162 {
163 for( int p=1; p <= 3; ++p )
164 {
165 run( grid, p );
166 }
167 }
168
169 return 0;
170 }
171 catch( const Exception& e )
172 {
173 std :: cerr << e.what() << std :: endl;
174 return 1;
175 }
176 }
177
178