1 #include <config.h>
2
3 #include <iostream>
4 #include <vector>
5 #include <algorithm>
6 #include <cstdlib>
7
8 #include <dune/common/dynvector.hh>
9
10 #include <dune/fem/test/testgrid.hh>
11 #include <dune/fem/misc/gridwidth.hh>
12
13 #include <dune/fem/gridpart/adaptiveleafgridpart.hh>
14 #include <dune/fem/space/finitevolume.hh>
15 #include <dune/fem/space/discontinuousgalerkin.hh>
16 #include <dune/fem/space/common/adaptationmanager.hh>
17
18 #include <dune/fem/function/blockvectordiscretefunction.hh>
19 #include <dune/fem/function/blockvectorfunction.hh>
20 #include <dune/fem/function/blockvectors/referenceblockvector.hh>
21 #include <dune/fem/function/common/localcontribution.hh>
22 #include <dune/fem/function/petscdiscretefunction.hh>
23 #include <dune/fem/function/vectorfunction.hh>
24 #include <dune/fem/function/vectorfunction/managedvectorfunction.hh>
25
26 #include <dune/fem/function/adaptivefunction.hh>
27 #include <dune/fem/function/combinedfunction.hh>
28
29 #include <dune/fem/misc/mpimanager.hh>
30
31 typedef Dune:: GridSelector::GridType HGridType;
32 typedef Dune::Fem::DGAdaptiveLeafGridPart< HGridType > GridPartType;
33 typedef Dune::Fem::FunctionSpace< double, double, HGridType::dimension, HGridType::dimension+2 > FunctionSpaceType;
34 typedef Dune::Fem::DiscontinuousGalerkinSpace< FunctionSpaceType, GridPartType, 2 > DiscreteFunctionSpaceType;
35
36 template <class DiscreteFunction, class OtherDiscreteFunction>
checkFunction(DiscreteFunction & df,OtherDiscreteFunction & other)37 void checkFunction( DiscreteFunction& df, OtherDiscreteFunction& other )
38 {
39 std::cout << "Checking (" << df.name() << "," << other.name() << "), size = ("<<df.size()<< "," << other.size() << ")....";
40
41 // fill df with zeros
42 df.clear();
43 other.clear();
44
45 // fill df with zeros
46 typedef typename DiscreteFunction :: DofType DofType;
47 std::fill( df.dbegin(), df.dend(), DofType( 0 ) );
48
49 df += other;
50 df -= other;
51
52 // fill df with zeros
53 for( auto&& dof : dofs(df) )
54 dof=static_cast<DofType>(0);
55
56 // fill with increasing values
57 int cont(0);
58 {
59 Dune::Fem::AddLocalContribution< DiscreteFunction > lf( df );
60 for( const auto& entity : entities(df) )
61 {
62 lf.bind( entity );
63 lf.clear();
64 for( int i = 0, n = lf.size(); i < n; ++i,++cont )
65 lf[ i ] = static_cast<DofType>(cont);
66 lf.unbind();
67 }
68 }
69
70 // check block access
71 const std::size_t localBlockSize = DiscreteFunctionSpaceType::localBlockSize;
72 const std::size_t numBlocks = df.blocks();
73 if( df.size() / localBlockSize != numBlocks )
74 DUNE_THROW(Dune::InvalidStateException,"Number of blocks not correct");
75
76 auto dfDofIt(df.dbegin());
77 for(std::size_t i=0;i<numBlocks;++i)
78 for(std::size_t j=0;j<localBlockSize;++j,++dfDofIt)
79 if( std::abs( df.dofVector()[i][j] - *dfDofIt ) > 1e-12 )
80 DUNE_THROW(Dune::InvalidStateException,"Block access did not work");
81
82 // copy to std::vector, sometimes needed for solver interfaces
83 std::vector< DofType > vec( df.size() );
84 std::copy( df.dbegin(), df.dend(), vec.begin() );
85
86 // check copy constructor
87 DiscreteFunction copydf( df );
88 if( !(copydf == df) )
89 DUNE_THROW(Dune::InvalidStateException,"Copy constructor did not work");
90
91 // check move constructor
92 DiscreteFunction movedf( std::move( copydf ) );
93 if( !(movedf == df) )
94 DUNE_THROW(Dune::InvalidStateException,"Move constrcutor did not work");
95
96 df.dofVector()[0] *= 1.0;
97 df.dofVector()[0][0] = 1.0;
98 df.dofVector()[0][HGridType::dimension-1] = 1.0;
99
100 df.assign( other );
101
102 df *= 2.0;
103 df /= 2.0;
104
105 df.axpy( 0.5, other );
106
107 df.enableDofCompression();
108
109 std::stringstream stream;
110 Dune::Fem::StandardOutStream out( stream );
111 df.write( out );
112 Dune::Fem::StandardInStream in( stream );
113 df.read( in );
114
115 df.scalarProductDofs( other );
116 df.scalarProductDofs( df );
117 other.scalarProductDofs( df );
118
119 std::cout << "done!" << std::endl;
120
121 Dune::Fem::RestrictProlongDefault<DiscreteFunction> rp( df );
122 rp.setFatherChildWeight(Dune::DGFGridInfo< HGridType >::refineWeight());
123 //Dune::Fem::AdaptationManager< HGridType, RPDefaultType > adop(grid,rp);
124 }
125
126 // main program
main(int argc,char ** argv)127 int main(int argc, char ** argv)
128 {
129 Dune::Fem::MPIManager :: initialize( argc, argv );
130 try
131 {
132 HGridType &grid = Dune::Fem::TestGrid :: grid();
133
134 GridPartType gridPart( grid );
135 std::cout << "Grid width: " << Dune::Fem::GridWidth :: calcGridWidth( gridPart ) << std::endl;
136
137 DiscreteFunctionSpaceType space( gridPart );
138
139 Dune::Fem::AdaptiveDiscreteFunction< DiscreteFunctionSpaceType > ref ("ref", space);
140
141 Dune::Fem::AdaptiveDiscreteFunction< DiscreteFunctionSpaceType > adf ("adaptive", space);
142 std::cout << "dofs = " << adf.size() << std::endl;
143
144 checkFunction( adf, ref );
145
146 std::vector< double > advec( space.size() );
147 Dune::Fem::AdaptiveDiscreteFunction< DiscreteFunctionSpaceType > adfp ("pointer", space, advec.data() );
148 checkFunction( adfp, ref );
149
150 Dune::DynamicVector< double > vec( space.size() );
151 typedef Dune::Fem::VectorDiscreteFunction< DiscreteFunctionSpaceType, Dune::DynamicVector< double > > VectorDiscreteFunctionType;
152 VectorDiscreteFunctionType vdf ("vector", space, vec);
153 checkFunction( vdf, ref );
154
155 std::vector<double> vec1( space.size() );
156 typedef Dune::Fem::VectorDiscreteFunction< DiscreteFunctionSpaceType, std::vector<double> > StdVectorDiscreteFunctionType;
157 StdVectorDiscreteFunctionType vdf1 ("vector", space, vec1);
158 checkFunction( vdf1, ref );
159
160 Dune::Fem::ManagedDiscreteFunction< VectorDiscreteFunctionType > mdf ("managed", space);
161 checkFunction( mdf, ref );
162
163 typedef Dune::Fem::ReferenceBlockVector< FunctionSpaceType::RangeFieldType, DiscreteFunctionSpaceType::localBlockSize > BlockVectorType;
164 Dune::Fem::BlockVectorDiscreteFunction< DiscreteFunctionSpaceType, BlockVectorType > bdf( "block", space );
165 checkFunction( bdf, ref );
166
167 Dune::Fem::ISTLBlockVectorDiscreteFunction< DiscreteFunctionSpaceType > istldf ("istl", space);
168 checkFunction( istldf, ref );
169
170 checkFunction( adf, istldf );
171 checkFunction( vdf, istldf );
172
173 #if HAVE_PETSC
174 Dune::Fem::PetscDiscreteFunction< DiscreteFunctionSpaceType > petscdf ("petsc", space);
175 checkFunction( petscdf, ref );
176 checkFunction( petscdf, vdf );
177 checkFunction( petscdf, adf );
178 checkFunction( petscdf, istldf );
179 checkFunction( adf, petscdf );
180 checkFunction( istldf, petscdf );
181 #endif
182
183 // refine grid
184 Dune::Fem::GlobalRefine::apply( grid, 1 );
185
186 std::cout << "dofs = " << adf.size() << std::endl;
187 std::cout << "dofs = " << istldf.size() << std::endl;
188 checkFunction( adf, istldf );
189 checkFunction( bdf, istldf );
190 checkFunction( istldf, ref );
191 checkFunction( mdf, ref );
192
193 return 0;
194 }
195 catch( const Dune::Exception& e )
196 {
197 std :: cerr << e.what() << std :: endl;
198 return 1;
199 }
200 }
201