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
10 #include <dune/common/filledarray.hh>
11 #include <dune/common/parallel/mpihelper.hh>
12
13 #include <dune/grid/yaspgrid.hh>
14
15 #include <dune/pdelab/finiteelementmap/qkdg.hh>
16 #include <dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
17
18 // Wrapped Qk local basis which can have a zero size on some elements
19 template <class T, int k, int d>
20 struct WrappedBasis {
21 using Basis = Dune::QkStuff::QkLocalBasis<T, T, k, d>;
22
23 Basis basis;
24 bool setZero;
25
WrappedBasisWrappedBasis26 WrappedBasis()
27 : setZero(false)
28 {
29 }
30
sizeWrappedBasis31 unsigned int size() const
32 {
33 return setZero ? 0 : basis.size();
34 }
35
evaluateFunctionWrappedBasis36 inline void evaluateFunction(const typename Basis::Traits::DomainType& in,
37 std::vector<typename Basis::Traits::RangeType>& out) const
38 {
39 if (setZero)
40 out.clear();
41 else
42 basis.evaluateFunction(in, out);
43 }
44
evaluateJacobianWrappedBasis45 void evaluateJacobian(const typename Basis::Traits::DomainType& in,
46 std::vector<typename Basis::Traits::JacobianType>& out) const
47 {
48 if (setZero)
49 out.clear();
50 else
51 basis.evaluateFunction(in, out);
52 }
53
orderWrappedBasis54 unsigned int order() const
55 {
56 return setZero ? 0 : basis.order();
57 }
58 };
59
60 // wrapped QkDG local coefficients which can have a zero size on some elements
61 template <int k, int d>
62 struct WrappedCoefficients {
63 using Coefficients = Dune::QkStuff::QkDGLocalCoefficients<k, d>;
64
65 Coefficients coeffs;
66 bool setZero;
67
WrappedCoefficientsWrappedCoefficients68 WrappedCoefficients()
69 : setZero(false)
70 {
71 }
72
sizeWrappedCoefficients73 std::size_t size() const
74 {
75 return setZero ? 0 : coeffs.size();
76 }
77
localKeyWrappedCoefficients78 const Dune::LocalKey& localKey(std::size_t i) const
79 {
80 if (setZero) {
81 DUNE_THROW(Dune::Exception, "set zero");
82 } else {
83 return coeffs.localKey(i);
84 }
85 }
86 };
87
88 // wrapped Qk local interpolation which can have a zero size on some elements
89 template <class T, int k, int d>
90 struct WrappedInterpolation {
91 using Interpolation
92 = Dune::QkStuff::QkLocalInterpolation<k, d, Dune::QkStuff::QkLocalBasis<T, T, k, d> >;
93
94 Interpolation interpolation;
95 bool setZero;
96
WrappedInterpolationWrappedInterpolation97 WrappedInterpolation()
98 : setZero(false)
99 {
100 }
101
102 template <typename F, typename C>
interpolateWrappedInterpolation103 void interpolate(const F& f, std::vector<C>& out) const
104 {
105 if (setZero)
106 out.clear();
107 else
108 interpolation.interpolate(f, out);
109 }
110 };
111
112 // wrapped Qk local finite element which can have a zero size on some elements
113 template <class T, int k, int d>
114 class WrappedFiniteElement
115 {
116 using Basis = WrappedBasis<T, k, d>;
117 using Coefficients = WrappedCoefficients<k, d>;
118 using Interpolation = WrappedInterpolation<T, k, d>;
119
120 public:
121 typedef Dune::LocalFiniteElementTraits<Basis, Coefficients, Interpolation> Traits;
122
WrappedFiniteElement()123 WrappedFiniteElement()
124 : gt(Dune::GeometryTypes::cube(d))
125 {}
126
localBasis() const127 const typename Traits::LocalBasisType& localBasis() const
128 {
129 return basis;
130 }
131
localCoefficients() const132 const typename Traits::LocalCoefficientsType& localCoefficients() const
133 {
134 return coefficients;
135 }
136
localInterpolation() const137 const typename Traits::LocalInterpolationType& localInterpolation() const
138 {
139 return interpolation;
140 }
141
type() const142 Dune::GeometryType type() const
143 {
144 return gt;
145 }
146
clone() const147 WrappedFiniteElement* clone() const
148 {
149 return new WrappedFiniteElement(*this);
150 }
151
setZero(bool value)152 void setZero(bool value)
153 {
154 basis.setZero = value;
155 coefficients.setZero = value;
156 interpolation.setZero = value;
157 }
158
159 private:
160 Basis basis;
161 Coefficients coefficients;
162 Interpolation interpolation;
163 Dune::GeometryType gt;
164 };
165
166 // wrapped QkDG local finite element map which has no local basis functions on elements with an even
167 // index
168 template <class GV, class T, int k, int d>
169 class WrappedFiniteElementMap
170 : public Dune::PDELab::
171 LocalFiniteElementMapInterface<Dune::PDELab::
172 LocalFiniteElementMapTraits<WrappedFiniteElement<T, k,
173 d> >,
174 WrappedFiniteElementMap<GV, T, k, d> >
175 {
176
177 public:
178 using Traits = Dune::PDELab::LocalFiniteElementMapTraits<WrappedFiniteElement<T, k, d> >;
179
WrappedFiniteElementMap(const GV & gv_,std::function<bool (std::size_t)> zeroPredicate_)180 WrappedFiniteElementMap(const GV& gv_, std::function<bool(std::size_t)> zeroPredicate_)
181 : gv(gv_), zeroPredicate(zeroPredicate_)
182 {
183 }
184
fixedSize() const185 bool fixedSize() const
186 {
187 return false;
188 }
189
hasDOFs(int codim) const190 bool hasDOFs(int codim) const
191 {
192 return codim == 0;
193 }
194
size(Dune::GeometryType gt) const195 std::size_t size(Dune::GeometryType gt) const
196 {
197 DUNE_THROW(Dune::PDELab::VariableElementSize, "the fem has variable element size");
198 }
199
maxLocalSize() const200 std::size_t maxLocalSize() const
201 {
202 return Dune::QkStuff::QkSize<k, d>::value;
203 }
204
205 template <class EntityType>
find(const EntityType & e) const206 const typename Traits::FiniteElementType& find(const EntityType& e) const
207 {
208 fe.setZero(zeroPredicate(gv.indexSet().index(e)));
209 return fe;
210 }
211
212 private:
213 mutable WrappedFiniteElement<T, k, d> fe;
214 GV gv;
215 std::function<bool(std::size_t)> zeroPredicate;
216 };
217
218 /**
219 * check if the size of the ordering corresponds to the number of basis functions reported by the
220 * finite element map
221 */
222 template <class GFS, class FEM>
testSize(const GFS & gfs,const FEM & fem)223 bool testSize(const GFS& gfs, const FEM& fem)
224 {
225 std::cout << "ordering fixedSize: " << gfs.ordering().fixedSize() << "\n";
226 std::cout << "ordering size: " << gfs.ordering().size() << "\n";
227 std::cout << "ordering blockCount: " << gfs.ordering().blockCount() << "\n";
228
229 std::size_t size = 0;
230 for (const auto& e : Dune::elements(gfs.gridView())) {
231 size += fem.find(e).localBasis().size();
232 }
233 std::cout << "size by summing: " << size << "\n";
234
235 return size == gfs.ordering().size();
236 }
237
238 /**
239 * create a grid function space on a structured grid and test wether the size of the ordering is
240 * correct
241 */
242 template <int dim>
testPredicate(std::size_t cells,bool should_be_fixed_size,std::function<bool (std::size_t)> predicate)243 bool testPredicate(std::size_t cells, bool should_be_fixed_size, std::function<bool(std::size_t)> predicate)
244 {
245 Dune::FieldVector<double, dim> L(1.0);
246 std::array<int, dim> N(Dune::filledArray<dim, int>(cells));
247 Dune::YaspGrid<dim> grid(L, N);
248
249 using GV = typename Dune::YaspGrid<dim>::LeafGridView;
250 using FEM = WrappedFiniteElementMap<GV, double, 1, dim>;
251 using GFS = Dune::PDELab::GridFunctionSpace<GV, FEM>;
252
253 GV gv = grid.leafGridView();
254 FEM fem(gv, predicate);
255 GFS gfs(gv, fem);
256 gfs.update();
257 return testSize(gfs, fem) and gfs.ordering().fixedSize() == should_be_fixed_size;
258 }
259
main(int argc,char ** argv)260 int main(int argc, char** argv)
261 {
262 try {
263 // Maybe initialize Mpi
264 Dune::MPIHelper::instance(argc, argv);
265
266 const std::size_t cells = 15;
267
268 // In these two basic scenarios, PDELab should detect that the ordering actually is fixed
269 // size and switch to fixed size mode after collecting the size information
270 std::cout << "testing 2D all cells active" << std::endl;
271 if (!testPredicate<2>(cells, true, [](auto i) { return true; })) {
272 return -1;
273 }
274 std::cout << "testing 2D all cells disabled" << std::endl;
275 if (!testPredicate<2>(cells, true, [](auto i) { return false; })) {
276 return -1;
277 }
278
279
280 // in the following we test different scenarios of a grid function space with a finite element
281 // that has 0 or c number of DOFs on different entities, where c is a constant. Wether an entity
282 // has 0 DOFs is determined by the predicates provided below.
283
284 std::cout << "testing 1D alternating sizes" << std::endl;
285 if (!testPredicate<1>(cells, false, [](auto i) { return i % 2 == 0; })) {
286 return -1;
287 }
288 std::cout << "testing 2D alternating sizes" << std::endl;
289 if (!testPredicate<2>(cells, false, [](auto i) { return i % 2 == 0; })) {
290 return -1;
291 }
292 std::cout << "testing 1D empty fem at the end" << std::endl;
293 if (!testPredicate<1>(cells, false, [cells](auto i) { return i > cells / 2; })) {
294 return -1;
295 }
296 std::cout << "testing 1D empty fem at the beginning" << std::endl;
297 if (!testPredicate<1>(cells, false, [cells](auto i) { return i <= cells / 2; })) {
298 return -1;
299 }
300 return 0;
301
302 } catch (Dune::Exception& e) {
303 std::cerr << "Dune reported error: " << e << std::endl;
304 return 1;
305 } catch (...) {
306 std::cerr << "Unknown exception thrown!" << std::endl;
307 return 1;
308 }
309 }
310