1 #include <config.h>
2
3 #include <algorithm>
4 #include <iostream>
5 #include <vector>
6
7 #include <dune/common/exceptions.hh>
8
9 #include <dune/geometry/referenceelements.hh>
10
11 #include <dune/fem/function/common/common.hh>
12 #include <dune/fem/function/adaptivefunction.hh>
13 #include <dune/fem/function/common/gridfunctionadapter.hh>
14 #include <dune/fem/function/localfunction/temporary.hh>
15 #include <dune/fem/gridpart/leafgridpart.hh>
16 #include <dune/fem/gridpart/adaptiveleafgridpart.hh>
17 #include <dune/fem/gridpart/common/gridpartadapter.hh>
18 #include <dune/fem/gridpart/idgridpart.hh>
19 #include <dune/fem/gridpart/filter/basicfilterwrapper.hh>
20 #include <dune/fem/gridpart/filter/domainfilter.hh>
21 #include <dune/fem/gridpart/filter/radialfilter.hh>
22 #include <dune/fem/gridpart/filteredgridpart.hh>
23 #include <dune/fem/gridpart/geogridpart.hh>
24 #include <dune/fem/misc/gridwidth.hh>
25 #include <dune/fem/space/common/interpolate.hh>
26 #include <dune/fem/space/lagrange.hh>
27
28 #include "failure.hh"
29 #include "checkseed.hh"
30 #include "checkgeometry.hh"
31 #include "checkindexset.hh"
32 #include "checkintersections.hh"
33 #include "checkgridpart.hh"
34 #include <dune/fem/test/testgrid.hh>
35
36
37 template< class GridPart, class LocalFunction >
testExchangeGeometry(const GridPart & gridPart,LocalFunction & localFunction)38 void testExchangeGeometry ( const GridPart &gridPart, LocalFunction &localFunction )
39 {
40 for( const auto& entity : elements( gridPart ) )
41 {
42 const auto refElement = Dune::referenceElement< typename GridPart::ctype, GridPart::dimension >( entity.type() );
43 localFunction.init( entity.impl().hostEntity() );
44 for( int i = 0; i < refElement.size( GridPart::dimension ); ++i )
45 for( int k = 0; k < GridPart::dimensionworld; ++k )
46 localFunction[ i*GridPart::dimensionworld + k ] = refElement.position( i, GridPart::dimension )[ k ];
47
48 auto xchgEntity = gridPart.exchangeGeometry( entity, localFunction );
49 const auto& xchgGeometry = xchgEntity.geometry();
50 if( xchgGeometry.type() != entity.type() )
51 DUNE_THROW( Dune::InvalidStateException, "exchangeGeometry returns wrong geometry type." );
52 if( (xchgGeometry.center() - refElement.position( 0, 0 )).two_norm() > 1e-8 )
53 {
54 std::cerr << "exchangeGeometry returns wrong center: " << xchgGeometry.center() << std::endl;
55 std::cerr << "(real geometry center is: " << entity.geometry().center() << ")." << std::endl;
56 }
57 }
58 };
59
60 template <class GridPartType>
testAll(GridPartType & gridPart)61 void testAll( GridPartType& gridPart )
62 {
63 // run tests
64 std::cout << "Testing entities" << std::endl;
65 testGridPart( gridPart );
66 std::cout << std::endl;
67
68 std::cout << "Testing subentities" << std::endl;
69 testSubEntities< GridPartType::dimension >( gridPart );
70 std::cout << std::endl;
71
72 std::cout << "GridWidth: " << Dune::Fem::GridWidth::calcGridWidth( gridPart ) << std::endl;
73
74 typedef Dune::DefaultFailureHandler FailureHandlerType;
75 FailureHandlerType failureHandler;
76 std::cout << "Testing entity seeds" << std::endl;
77 Dune::Fem::CheckEntitySeed< GridPartType >::check( gridPart );
78 std::cout << "Testing geometries" << std::endl;
79 Dune::Fem::CheckGeometry< GridPartType, FailureHandlerType >::check( gridPart, failureHandler );
80 std::cout << "Testing index set" << std::endl;
81 Dune::Fem::CheckIndexSet< GridPartType, FailureHandlerType >::check( gridPart, failureHandler );
82 std::cout << "Testing intersections" << std::endl;
83 Dune::Fem::CheckIntersections< GridPartType, FailureHandlerType >::check( gridPart, failureHandler );
84
85 std::cout << std::endl << std::endl;
86 }
87
88 template< bool UseConsecutiveIndexSet, class HostGridPartType, class FilterType >
testFilteredGridPart(HostGridPartType & hostGridPart,FilterType & filter)89 void testFilteredGridPart( HostGridPartType& hostGridPart, FilterType& filter )
90 {
91 typedef Dune::Fem::FilteredGridPart< HostGridPartType, FilterType, UseConsecutiveIndexSet > GridPartType;
92 GridPartType gridPart( hostGridPart, filter );
93 testAll( gridPart );
94 }
95
main(int argc,char ** argv)96 int main ( int argc, char ** argv )
97 try
98 {
99 Dune::Fem::MPIManager::initialize( argc, argv );
100
101 // create grid
102 auto& grid = Dune::Fem::TestGrid::grid();
103
104 // refine grid
105 const int step = Dune::Fem::TestGrid::refineStepsForHalf();
106 grid.globalRefine( 2*step );
107 grid.loadBalance();
108
109 typedef typename Dune::Fem::TestGrid::HGridType GridType;
110 // create grid part
111 typedef Dune::Fem::LeafGridPart< GridType > HostGridPartType;
112 HostGridPartType hostGridPart( grid );
113
114 // AdaptiveLeafGridPart
115 {
116 std::cout << "*************************************" << std::endl;
117 std::cout << "*** AdaptiveLeafridPart ***"<< std::endl;
118 std::cout << "*************************************" << std::endl;
119
120 typedef Dune::Fem::AdaptiveLeafGridPart< GridType > GridPartType;
121 typedef typename GridPartType :: GridViewType GridViewType;
122 //GridPartType gridPart( grid );
123 std::unique_ptr< GridPartType > gpPtr;
124 std::unique_ptr< GridViewType > gvPtr;
125
126 gpPtr.reset( new GridPartType( grid ) );
127 GridPartType gp2( grid );
128
129 gvPtr.reset( new GridViewType( gpPtr->gridView() ) );
130
131 gpPtr.reset( new GridPartType( grid ) );
132 gvPtr.reset( new GridViewType( gpPtr->gridView() ) );
133 gpPtr.reset();
134 gvPtr.reset();
135
136 GridPartType& gridPart = gp2;
137
138 testAll( gridPart );
139 }
140
141 // GridPartAdapter
142 {
143 std::cout << "*******************************" << std::endl;
144 std::cout << "*** GridPartAdapter ***"<< std::endl;
145 std::cout << "*******************************" << std::endl;
146
147 // create grid part
148 typedef Dune::GridSelector::GridType GridType;
149 typedef Dune::Fem::GridPartAdapter< typename GridType::LeafGridView > GridPartType;
150
151 GridPartType gridPart( grid.leafGridView() );
152 testAll( gridPart );
153 }
154
155
156 // IdGridPart
157 {
158 std::cout << "*******************************" << std::endl;
159 std::cout << "*** IdGridPart ***"<< std::endl;
160 std::cout << "*******************************" << std::endl;
161
162 typedef Dune::Fem::IdGridPart< HostGridPartType > GridPartType;
163 GridPartType gridPart( hostGridPart );
164
165 testAll( gridPart );
166 }
167
168 // GeoGridPart
169 {
170 std::cout << "*******************************" << std::endl;
171 std::cout << "*** GeoGridPart ***"<< std::endl;
172 std::cout << "*******************************" << std::endl;
173
174 typedef Dune::Fem::FunctionSpace< GridType::ctype, GridType::ctype, GridType::dimensionworld, GridType::dimensionworld > CoordFunctionSpaceType;
175 typedef Dune::Fem::LagrangeDiscreteFunctionSpace< CoordFunctionSpaceType, HostGridPartType, 1 > DiscreteCoordFunctionSpaceType;
176 DiscreteCoordFunctionSpaceType coordFunctionSpace( hostGridPart );
177 typedef Dune::Fem::AdaptiveDiscreteFunction< DiscreteCoordFunctionSpaceType > CoordFunctionType;
178 CoordFunctionType coordFunction( "coordinate function", coordFunctionSpace );
179 typedef Dune::Fem::Identity< CoordFunctionSpaceType > IdentityType;
180 IdentityType identity;
181 Dune::Fem::GridFunctionAdapter< IdentityType, HostGridPartType > identitydDF( "identity", identity, hostGridPart );
182 Dune::Fem::interpolate( identitydDF, coordFunction );
183 typedef Dune::Fem::GeoGridPart< CoordFunctionType > GridPartType;
184 GridPartType gridPart( coordFunction );
185
186 testAll( gridPart );
187
188 std::cout << "Testing exchange geometry" << std::endl;
189 Dune::Fem::TemporaryLocalFunction< DiscreteCoordFunctionSpaceType > tlf( coordFunctionSpace );
190 testExchangeGeometry( gridPart, tlf );
191 std::cout << std::endl;
192 }
193
194 // FilteredGridPart
195 {
196 std::cout << "******************************" << std::endl;
197 std::cout << "*** FilteredGridPart ***"<< std::endl;
198 std::cout << "******************************" << std::endl;
199
200 // create radial filter
201 typedef Dune::Fem::RadialFilter< GridType::ctype, GridType::dimensionworld > RadialFilterType;
202 typedef Dune::Fem::BasicFilterWrapper< HostGridPartType, RadialFilterType > WrapperRadialFilterType;
203 WrapperRadialFilterType wrappedRadialFilter( hostGridPart, RadialFilterType::GlobalCoordinateType( 0 ), 0.25 );
204
205 // test FilteredGridPart with radial filter and allowing non consecutive index set
206 std::cout << std::endl << "Testing FilteredGridPart with radial filter: allow non consecutive index set" << std::endl << std::endl;
207 testFilteredGridPart< false, HostGridPartType, WrapperRadialFilterType >( hostGridPart, wrappedRadialFilter );
208
209 // test FilteredGridPart with radial filter and forcing consecutive index set
210 std::cout << std::endl << "Testing FilteredGridPart with radial filter: force consecutive index set" << std::endl << std::endl;
211 testFilteredGridPart< true, HostGridPartType, WrapperRadialFilterType >( hostGridPart, wrappedRadialFilter );
212 }
213
214
215 {
216 std::cout << "************************************" << std::endl;
217 std::cout << "*** FilteredGridPart (tags) ***"<< std::endl;
218 std::cout << "************************************" << std::endl;
219
220 // create domain filter
221 typedef std::vector< int > DomainArrayType;
222 DomainArrayType tags( grid.size(0), 0 );
223 for( std::size_t i = ( tags.size()/2 ); i < tags.size(); ++i )
224 tags[ i ] = 1;
225 typedef Dune::Fem::DomainFilter< HostGridPartType, DomainArrayType > DomainFilterType;
226 typedef Dune::Fem::BasicFilterWrapper< HostGridPartType, DomainFilterType > WrapperDomainFilterType;
227 WrapperDomainFilterType wrapperDomainFilter( hostGridPart, hostGridPart, tags, 1 );
228
229 // test FilteredGridPart with domain filter and allowing non consecutive index set
230 std::cout << std::endl << "Testing FilteredGridPart with domain filter: allow non consecutive index set" << std::endl << std::endl;
231 testFilteredGridPart< false, HostGridPartType, WrapperDomainFilterType >( hostGridPart, wrapperDomainFilter );
232 }
233
234
235 return 0;
236 }
237 catch( const Dune::Exception &e )
238 {
239 std::cerr << e << std::endl;
240 return 1;
241 }
242 catch( ... )
243 {
244 std::cerr << "Generic exception!" << std::endl;
245 return 2;
246 }
247