1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3
4 #ifdef HAVE_CONFIG_H
5 #include "config.h"
6 #endif
7
8 #include <iostream>
9 #include <memory>
10
11 #include <dune/common/exceptions.hh>
12 #include <dune/common/filledarray.hh>
13 #include <dune/common/fvector.hh>
14 #include <dune/common/parallel/mpihelper.hh>
15
16 #include <dune/grid/yaspgrid.hh>
17 #if HAVE_DUNE_ALUGRID
18 #include <dune/alugrid/grid.hh>
19 #endif
20 #include <dune/grid/common/scsgmapper.hh>
21 #include <dune/grid/utility/structuredgridfactory.hh>
22
23 #include <dune/pdelab.hh>
24
25 template<typename GFS>
check_ordering(const GFS & gfs)26 void check_ordering(const GFS& gfs)
27 {
28 const auto& ordering = gfs.ordering();
29
30 Dune::PDELab::LocalFunctionSpace lfs{gfs};
31 Dune::PDELab::LFSIndexCache lfs_cache{lfs};
32
33 for (const auto& element : elements(gfs.entitySet()))
34 {
35 // lfs computes dof indices
36 lfs.bind(element);
37 // cache computes container indices
38 lfs_cache.update();
39
40 for (unsigned i = 0; i < lfs.size(); ++i) {
41 // get i-th local dof index
42 auto di = lfs.dofIndex(i);
43 // map dof index to a container index (from ordering)
44 auto ci_map = ordering.mapIndex(di);
45 // map local index to a container index (from cache)
46 auto ci_cache_i = lfs_cache.containerIndex(i);
47 // map dof index to a container index (from cache)
48 auto ci_cache_di = lfs_cache.containerIndex(di);
49 if (ci_map != ci_cache_i or ci_map != ci_cache_di)
50 DUNE_THROW(Dune::RangeError, "Container index mappings do not match");
51
52 // check that all size suffixes fit the current container index
53 auto size_suffix = ci_map;
54 while (size_suffix.size() != 0) {
55 // get outer container index in the suffix
56 auto block_index = size_suffix.back();
57 // calculate the size for a container that would hold such block
58 size_suffix.pop_back();
59 auto size = ordering.size(size_suffix);
60 // the index should always fit into the size
61 if (size <= block_index)
62 DUNE_THROW(Dune::RangeError, "Size for CI suffix '"
63 << size_suffix
64 << "' cannot fit sub-block '"
65 << block_index << "'");
66 }
67 }
68 }
69
70 using V = Dune::PDELab::Backend::Vector<GFS,double>;
71 V x(gfs);
72 x = 0.0;
73 }
74
75
76 // test function trees
77 template<int dim, bool cube>
78 struct test;
79
80 template<>
81 struct test<2,true> {
82 template<class GV>
testleafgridfunctiontest83 static void testleafgridfunction(const GV& gv)
84 {
85 // instantiate finite element maps
86 auto gt = Dune::GeometryTypes::quadrilateral;
87 typedef Dune::PDELab::P0LocalFiniteElementMap<float,double,GV::dimension> P0FEM;
88 P0FEM p0fem(gt);
89 typedef Dune::PDELab::QkLocalFiniteElementMap<GV,float,double,1> Q12DFEM;
90 Q12DFEM q12dfem(gv);
91 typedef Dune::PDELab::QkLocalFiniteElementMap<GV,float,double,2> Q22DFEM;
92 Q22DFEM q22dfem(gv);
93 typedef Dune::SingleCodimSingleGeomTypeMapper<GV, 0> CellMapper;
94 CellMapper cellmapper(gv);
95 typedef Dune::PDELab::VariableMonomLocalFiniteElementMap<CellMapper,float,double,GV::dimension> MonomFEM;
96 MonomFEM monomfem(cellmapper, gt, 2);
97
98 typedef Dune::PDELab::NoConstraints CON;
99
100 typedef Dune::PDELab::ISTL::VectorBackend<> VBE;
101
102 // make a grid function space
103 typedef Dune::PDELab::GridFunctionSpace<GV,P0FEM,CON,VBE> P0GFS;
104 P0GFS p0gfs(gv,p0fem);
105 typedef Dune::PDELab::GridFunctionSpace<GV,Q12DFEM,CON,VBE> GFS1;
106 {
107 GFS1 gfs1(gv,q12dfem);
108 check_ordering(gfs1);
109 }
110 GFS1 gfs1(gv,q12dfem);
111 typedef Dune::PDELab::GridFunctionSpace<GV,Q22DFEM,CON,VBE> GFS2;
112 GFS2 gfs2(gv,q22dfem);
113 typedef Dune::PDELab::GridFunctionSpace<GV,MonomFEM,CON,VBE> GFS3;
114 {
115 GFS3 gfs3(gv,monomfem);
116 check_ordering(gfs3);
117 }
118 GFS3 gfs3(gv,monomfem);
119
120 typedef Dune::PDELab::PowerGridFunctionSpace<GFS1,3,VBE,Dune::PDELab::InterleavedOrderingTag> P1GFS;
121
122 {
123 // keep this block around to test the vector-based constructor of the ordering tag
124 std::vector<std::size_t> p1_gfs_block_sizes(3);
125 std::fill(p1_gfs_block_sizes.begin(),p1_gfs_block_sizes.end(),1);
126 P1GFS p1gfs(gfs1,gfs1,gfs1,VBE(),p1_gfs_block_sizes);
127 }
128
129 P1GFS p1gfs(gfs1,gfs1,gfs1,VBE(),{{1,1,1}});
130
131 typedef Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed> NVBE;
132
133 using EGFS = Dune::PDELab::GridFunctionSpace<GV,Q12DFEM,CON,VBE>;
134 EGFS egfs(gv,q12dfem);
135 using EPGFS = Dune::PDELab::PowerGridFunctionSpace<EGFS,3,VBE,Dune::PDELab::InterleavedOrderingTag>;
136 {
137 EPGFS epgfs(egfs,egfs,egfs,VBE(),{{1,1,1}});
138 check_ordering(epgfs);
139 }
140
141 EPGFS epgfs(egfs,egfs,egfs,VBE(),{{1,1,1}});
142 typedef Dune::PDELab::PowerGridFunctionSpace<EPGFS,2,NVBE,Dune::PDELab::InterleavedOrderingTag> PGFS;
143 std::vector<std::size_t> p_gfs_block_sizes(2);
144 std::fill(p_gfs_block_sizes.begin(),p_gfs_block_sizes.end(),3);
145 PGFS pgfs(epgfs,epgfs,NVBE(),p_gfs_block_sizes);
146
147
148 typedef Dune::PDELab::PowerGridFunctionSpace<GFS1,3,VBE,Dune::PDELab::EntityBlockedOrderingTag> P2GFS;
149 P2GFS p2gfs(gfs1,gfs1,gfs1);
150
151 typedef Dune::PDELab::PowerGridFunctionSpace<P2GFS,2,NVBE,Dune::PDELab::EntityBlockedOrderingTag> P3GFS;
152 P3GFS p3gfs(p2gfs,p2gfs,NVBE());
153
154 check_ordering(pgfs);
155 //check_ordering(p3gfs);
156
157 typedef Dune::PDELab::PowerGridFunctionSpace<GFS1,3,VBE,Dune::PDELab::EntityBlockedOrderingTag> EBPGFS1;
158 {
159 EBPGFS1 ebpgfs1(gfs1);
160 check_ordering(ebpgfs1);
161 }
162 EBPGFS1 ebpgfs1(gfs1);
163
164 typedef Dune::PDELab::CompositeGridFunctionSpace<
165 VBE,
166 Dune::PDELab::EntityBlockedOrderingTag,
167 GFS1,
168 EBPGFS1,
169 GFS1
170 > EBCGFS1;
171
172 EBCGFS1 ebcgfs1(gfs1,ebpgfs1,gfs1);
173
174 check_ordering(ebcgfs1);
175 }
176 };
177
178 template<>
179 struct test<2,false> {
180 template<class GV>
testleafgridfunctiontest181 static void testleafgridfunction(const GV& gv)
182 {
183 // instantiate finite element maps
184 typedef Dune::SingleCodimSingleGeomTypeMapper<GV, 0> CellMapper;
185 CellMapper cellmapper(gv);
186 typedef Dune::PDELab::VariableMonomLocalFiniteElementMap<CellMapper,float,double,GV::dimension> MonomFEM;
187 MonomFEM monomfem(cellmapper,gv.indexSet().types(0)[0],2);
188
189 typedef Dune::PDELab::NoConstraints CON;
190
191 typedef Dune::PDELab::ISTL::VectorBackend<> VBE;
192
193 // make a grid function space
194 typedef Dune::PDELab::GridFunctionSpace<GV,MonomFEM,CON,VBE> GFS3;
195 GFS3 gfs3(gv,monomfem);
196 check_ordering(gfs3);
197 }
198 };
199
200 template<>
201 struct test<3,true> {
202 template<class GV>
testleafgridfunctiontest203 static void testleafgridfunction(const GV& gv)
204 {
205 // instantiate finite element maps
206 auto gt = Dune::GeometryTypes::hexahedron;
207 typedef Dune::PDELab::P0LocalFiniteElementMap<float,double,GV::dimension> P0FEM;
208 P0FEM p0fem(gt);
209 typedef Dune::PDELab::QkLocalFiniteElementMap<GV,float,double,1> Q1FEM;
210 Q1FEM q1fem(gv);
211
212 // make a grid function space
213 typedef Dune::PDELab::GridFunctionSpace<GV,P0FEM> P0GFS;
214 P0GFS p0gfs(gv,p0fem);
215 check_ordering(p0gfs);
216 // Doesn't work, we need a grid with triangular elemets for that
217 // typedef Dune::PDELab::GridFunctionSpace<GV,P1FEM> P1GFS;
218 // P1GFS p1gfs(gv,p1fem);
219 typedef Dune::PDELab::GridFunctionSpace<GV,Q1FEM> Q1GFS;
220 Q1GFS q1gfs(gv,q1fem);
221 check_ordering(q1gfs);
222 }
223 };
224
225 template<bool cube, class GV>
testleafgridfunction(const GV & gv)226 void testleafgridfunction(const GV& gv)
227 {
228 test<GV::dimension,cube>::testleafgridfunction(gv);
229 }
230
231
main(int argc,char ** argv)232 int main(int argc, char** argv)
233 {
234 //Maybe initialize Mpi
235 Dune::MPIHelper::instance(argc, argv);
236
237 #if HAVE_DUNE_ALUGRID
238 // 2D
239 {
240 std::cout << "2D tests" << std::endl;
241 // need a grid in order to test grid functions
242 // typedef Dune::YaspGrid<2> Grid;
243 typedef Dune::ALUGrid<2,2,Dune::cube,Dune::nonconforming> Grid;
244 Dune::FieldVector<double,2> l(0.0);
245 Dune::FieldVector<double,2> u(1.0);
246 std::array<unsigned int,2> N = {{1,1}};
247 std::shared_ptr<Grid> grid = Dune::StructuredGridFactory<Grid>::createCubeGrid(l,u,N);
248 grid->globalRefine(1);
249
250 std::cout << Dune::GlobalGeometryTypeIndex::index(grid->leafGridView().template begin<0>()->type()) << std::endl;
251 testleafgridfunction<true>(grid->leafGridView());
252 }
253
254 {
255 std::cout << "2D tests" << std::endl;
256 // need a grid in order to test grid functions
257 // typedef Dune::YaspGrid<2> Grid;
258 typedef Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming> Grid;
259 Dune::FieldVector<double,2> l(0.0);
260 Dune::FieldVector<double,2> u(1.0);
261 std::array<unsigned int,2> N = {{1,1}};
262 std::shared_ptr<Grid> grid = Dune::StructuredGridFactory<Grid>::createSimplexGrid(l,u,N);
263 grid->globalRefine(1);
264
265 std::cout << Dune::GlobalGeometryTypeIndex::index(grid->leafGridView().template begin<0>()->type()) << std::endl;
266 testleafgridfunction<false>(grid->leafGridView());
267 }
268 #endif
269
270 // 3D
271 {
272 std::cout << "3D tests" << std::endl;
273 // need a grid in order to test grid functions
274 Dune::FieldVector<double,3> L(1.0);
275 std::array<int,3> N(Dune::filledArray<3,int>(1));
276 Dune::YaspGrid<3> grid(L,N);
277 grid.globalRefine(1);
278
279 testleafgridfunction<true>(grid.leafGridView());
280 }
281
282 // test passed
283 return 0;
284 }
285