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