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