1 #include <config.h>
2
3 #include <sstream>
4 #include <string>
5
6 #include <dune/common/exceptions.hh>
7
8 #include <dune/geometry/dimension.hh>
9
10 #include <dune/grid/common/partitionset.hh>
11 #include <dune/grid/common/rangegenerators.hh>
12 #include <dune/grid/io/file/dgfparser/dgfparser.hh>
13
14 #include <dune/fem/common/hybrid.hh>
15 #include <dune/fem/function/hierarchical.hh>
16 #include <dune/fem/function/common/gridfunctionadapter.hh>
17 #include <dune/fem/function/localfunction/const.hh>
18 #include <dune/fem/function/localfunction/temporary.hh>
19 #include <dune/fem/gridpart/leafgridpart.hh>
20 #include <dune/fem/misc/mpimanager.hh>
21 #include <dune/fem/space/discontinuousgalerkin/lagrange.hh>
22 #include <dune/fem/space/discontinuousgalerkin/space.hh>
23 #include <dune/fem/space/discontinuousgalerkin/tuple.hh>
24 #include <dune/fem/space/lagrange.hh>
25 #include <dune/fem/test/exactsolution.hh>
26
27
28 // dgfUnitCube
29 // -----------
30
dgfUnitCube(int dimWorld,int cells,int overlap=1)31 inline static std::string dgfUnitCube ( int dimWorld, int cells, int overlap = 1 )
32 {
33 std::string dgf = "DGF\nINTERVAL\n";
34 for( int i = 0; i < dimWorld; ++i )
35 dgf += " 0";
36 dgf += "\n";
37 for( int i = 0; i < dimWorld; ++i )
38 dgf += " 1";
39 dgf += "\n";
40 for( int i = 0; i < dimWorld; ++i )
41 dgf += (" " + std::to_string( cells ));
42 dgf += "\n#\n\n";
43 dgf += "GRIDPARAMETER\n";
44 dgf += "OVERLAP " + std::to_string( overlap ) + "\n";
45 dgf += "#\n";
46 return dgf;
47 }
48
49
50
51 // FunctionSpace
52 // -------------
53
54 template< class GridPart, int dimRange >
55 using FunctionSpace = Dune::Fem::GridFunctionSpace< GridPart, Dune::Dim< dimRange > >;
56
57
58
59 // TaylorHoodDGSpace
60 // -----------------
61
62 template< class GridPart >
63 using VelocityDGSpace = Dune::Fem::LagrangeDiscontinuousGalerkinSpace< FunctionSpace< GridPart, GridPart::dimensionworld >, GridPart, 2 >;
64
65 template< class GridPart >
66 using PressureDGSpace = Dune::Fem::DiscontinuousGalerkinSpace< FunctionSpace< GridPart, 1 >, GridPart, 1 >;
67
68 template< class GridPart >
69 using TaylorHoodDGSpace = Dune::Fem::TupleDiscontinuousGalerkinSpace< VelocityDGSpace< GridPart >, PressureDGSpace< GridPart > >;
70
71
72
73 // LagrangeSpace
74 // -------------
75
76 template< class GridPart >
77 using LagrangeSpace = Dune::Fem::LagrangeDiscreteFunctionSpace< FunctionSpace< GridPart, GridPart::dimensionworld+1 >, GridPart, 2 >;
78
79
80
81 // equals
82 // ------
83
84 template< class T >
equals(const T & u,const T & v)85 bool equals ( const T &u, const T &v )
86 {
87 typedef std::decay_t< decltype( std::declval< const T & >().two_norm2() ) > real_type;
88
89 T w( u );
90 w -= v;
91 return (w.two_norm2() < 128 * std::numeric_limits< real_type >::epsilon());
92 }
93
94
95
96 // interpolateOnly
97 // ---------------
98
99 template< class GridFunction, class DiscreteFunction >
interpolateOnly(const GridFunction & u,DiscreteFunction & v)100 void interpolateOnly ( const GridFunction &u, DiscreteFunction &v )
101 {
102 v.clear();
103
104 Dune::Fem::ConstLocalFunction< GridFunction > uLocal( u );
105 Dune::Fem::TemporaryLocalFunction< typename DiscreteFunction::DiscreteFunctionSpaceType > vLocal( v.space() );
106 Dune::Fem::LocalInterpolation< typename DiscreteFunction::DiscreteFunctionSpaceType > interpolation( v.space() );
107
108 // iterate over selected partition
109 for( const auto entity : elements( v.gridPart(), Dune::Partitions::all ) )
110 {
111 // initialize u, v to entity
112 auto uGuard = bindGuard( uLocal, entity );
113 auto vGuard = bindGuard( vLocal, entity );
114
115 // bind interpolation to entity
116 auto iGuard = bindGuard( interpolation, entity );
117
118 // perform local interpolation
119 interpolation( uLocal, vLocal.localDofVector() );
120
121 // write interpolation into global DoF vector
122 v.setLocalDofs( entity, vLocal.localDofVector() );
123 }
124 }
125
126
127
128 // performTest
129 // -----------
130
131 template< class DiscreteFunctionSpace >
performTest(const DiscreteFunctionSpace & dfSpace)132 void performTest ( const DiscreteFunctionSpace &dfSpace )
133 {
134 // some stupid type defs
135 typedef Dune::Fem::HierarchicalDiscreteFunction< DiscreteFunctionSpace > DiscreteFunctionType;
136
137 typedef typename DiscreteFunctionSpace::GridPartType GridPartType;
138 typedef typename DiscreteFunctionSpace::LocalBlockIndices BlockIndices;
139
140 typedef typename DiscreteFunctionType::DofType DofType;
141
142 // interpolate a function
143 Dune::Fem::ExactSolution< FunctionSpace< GridPartType, GridPartType::dimensionworld+1 > > uExact;
144 const auto uGridExact = gridFunctionAdapter( "exact solution", uExact, dfSpace.gridPart(), 3 );
145 DiscreteFunctionType u( "solution", dfSpace );
146 interpolateOnly( uGridExact, u );
147
148 // copy discrete function and clear auxiliary dofs
149 DiscreteFunctionType w( u );
150 for( const auto &auxiliaryDof : dfSpace.auxiliaryDofs() )
151 Dune::Hybrid::forEach( BlockIndices(), [ &w, &auxiliaryDof ] ( auto &&j ) { w.dofVector()[ auxiliaryDof ][ j ] = DofType( 0 ); } );
152
153 // make sure u and w differ
154 if( equals( u.dofVector().array(), w.dofVector().array() ) )
155 DUNE_THROW( Dune::InvalidStateException, "Unique representation does not differ from consistent representation." );
156
157 // communicate w
158 w.communicate();
159
160 // now u and w should not differ
161 if( !equals( u.dofVector().array(), w.dofVector().array() ) )
162 DUNE_THROW( Dune::InvalidStateException, "Functions differ after communication." );
163 }
164
165
166
167 // main
168 // ----
169
main(int argc,char ** argv)170 int main ( int argc, char **argv )
171 try
172 {
173 Dune::Fem::MPIManager::initialize( argc, argv );
174
175 // construct unit cube
176 typedef typename Dune::GridSelector::GridType GridType;
177 std::istringstream dgf( dgfUnitCube( GridType::dimensionworld, 16 ) );
178 Dune::GridPtr< GridType > grid( dgf );
179
180 // create leaf grid part
181 typedef Dune::Fem::LeafGridPart< GridType > GridPartType;
182 GridPartType gridPart( *grid );
183
184 // test Taylor-Hood DG space
185 TaylorHoodDGSpace< GridPartType > taylorHoodDGSpace( gridPart );
186 performTest( taylorHoodDGSpace );
187
188 // test Lagrange space
189 LagrangeSpace< GridPartType > lagrangeSpace( gridPart );
190 performTest( lagrangeSpace );
191
192 return 0;
193 }
194 catch( const Dune::Exception &e )
195 {
196 std::cout << "Exception: " << e << std::endl;
197 return 1;
198 }
199