1 #ifdef HAVE_CONFIG_H
2 #include "config.h"
3 #endif
4 
5 #include <iostream>
6 #include <cassert>
7 
8 #include <dune/common/filledarray.hh>
9 #include <dune/common/math.hh>
10 #include <dune/common/parallel/mpihelper.hh>
11 
12 #include <dune/grid/yaspgrid.hh>
13 
14 #include <dune/pdelab.hh>
15 
16 // an analytic scalar function
17 template<typename T>
18 class F : public Dune::PDELab::FunctionInterface<
19   Dune::PDELab::FunctionTraits<T,2,Dune::FieldVector<T,2>,
20                                T,1,Dune::FieldVector<T,1> >,
21   F<T> >
22 {
23 public:
evaluate(const Dune::FieldVector<T,2> & x,Dune::FieldVector<T,1> & y) const24   inline void evaluate (const Dune::FieldVector<T,2>& x,
25                         Dune::FieldVector<T,1>& y) const
26   {
27     y = sin(3*Dune::StandardMathematicalConstants<T>::pi()*x[0])
28       * cos(7*Dune::StandardMathematicalConstants<T>::pi()*x[1]);
29   }
30 };
31 
32 // an analytic vector-valued function
33 template<typename T>
34 class G : public Dune::PDELab::FunctionInterface<
35   Dune::PDELab::FunctionTraits<T,2,Dune::FieldVector<T,2>,
36                                T,2,Dune::FieldVector<T,2> >,
37   G<T> >
38 {
39 public:
evaluate(const Dune::FieldVector<T,2> & x,Dune::FieldVector<T,2> & y) const40   inline void evaluate (const Dune::FieldVector<T,2>& x,
41                         Dune::FieldVector<T,2>& y) const
42   {
43     y[0] =  x.two_norm()*x[1];
44     y[1] = -x.two_norm()*x[0];
45   }
46 };
47 
48 
49 // iterate over grid view and use analytic function as grid function through adapter
50 template<class GV, class T>
testgridfunction(const GV & gv,const T & t)51 void testgridfunction (const GV& gv, const T& t)
52 {
53   // make a grid function from the analytic function
54   typedef Dune::PDELab::FunctionToGridFunctionAdapter<GV,T> GF;
55   GF gf(gv,t);
56 
57   typedef typename GV::Traits::template Codim<0>::Iterator ElementIterator;
58   for (ElementIterator it = gv.template begin<0>();
59        it!=gv.template end<0>(); ++it)
60     {
61       typename T::Traits::DomainType x(0.0);
62       typename T::Traits::RangeType y;
63       gf.evaluate(*it,x,y);
64 
65       // make a function in local coordinates from the grid function
66       Dune::PDELab::GridFunctionToLocalFunctionAdapter<GF>
67         lf(gf,*it);
68       lf.evaluate(x,y);
69 
70       // make a function in local coordinates from the global function
71       Dune::PDELab::GlobalFunctionToLocalFunctionAdapter<T,typename ElementIterator::Entity>
72         lg(t,*it);
73       lg.evaluate(x,y);
74     }
75 }
76 
77 // iterate over grid view and use analytic function as grid function through adapter
78 template<class GV, class T>
testvtkexport(const GV & gv,const T & t)79 void testvtkexport (const GV& gv, const T& t)
80 {
81   // make a grid function from the analytic function
82   typedef Dune::PDELab::FunctionToGridFunctionAdapter<GV,T> GF;
83   GF gf(gv,t);
84 
85   // make a VTKFunction from grid function
86   Dune::PDELab::VTKGridFunctionAdapter<GF> vtkf(gf,"blub");
87 
88   Dune::VTKWriter<GV> vtkwriter(gv,Dune::VTK::conforming);
89   vtkwriter.addVertexData(std::make_shared< Dune::PDELab::VTKGridFunctionAdapter<GF> >(gf,"blub")); // VTKWriter takes control
90   vtkwriter.write("single",Dune::VTK::ascii);
91 }
92 
93 // a grid function
94 template<typename G, typename T>
95 class L : public Dune::PDELab::GridFunctionBase<
96   Dune::PDELab::GridFunctionTraits<G,T,1,Dune::FieldVector<T,1> >,
97   L<G,T> >
98 {
99   typedef typename G::Traits::template Codim<0>::Entity ElementType;
100 public:
L(const G & g_)101   L (const G& g_) : g(g_) {}
102 
evaluate(const Dune::FieldVector<T,2> & x,Dune::FieldVector<T,1> & y) const103   inline void evaluate (const Dune::FieldVector<T,2>& x,
104                         Dune::FieldVector<T,1>& y) const
105   {
106     y = sin(3*3.1415*x[0])*cos(7*3.1415*x[1]);
107   }
108 
evaluate(const ElementType & e,const Dune::FieldVector<T,2> & x,Dune::FieldVector<T,1> & y) const109   inline void evaluate (const ElementType& e,
110                         const Dune::FieldVector<T,2>& x,
111                         Dune::FieldVector<T,1>& y) const
112   {
113     evaluate(e.geometry().global(x),y);
114   }
115 
getGridView() const116   inline const G& getGridView () const
117   {
118     return g;
119   }
120 
121 private:
122   const G& g;
123 };
124 
125 
126 // test function trees
127 template<class GV>
testfunctiontree(const GV & gv)128 void testfunctiontree (const GV& gv)
129 {
130   // a leaf function
131   typedef L<GV,typename GV::Grid::ctype> A;
132   A a(gv);
133 
134   // test power
135   typedef Dune::PDELab::PowerGridFunction<A,10> B10;
136   B10 b10(a);
137   typedef Dune::PDELab::PowerGridFunction<A,2> B2;
138   B2 b2(a,a);
139   typedef Dune::PDELab::PowerGridFunction<A,3> B3;
140   B3 b3(a,a,a);
141   typedef Dune::PDELab::PowerGridFunction<A,4> B4;
142   B4 b4(a,a,a,a);
143   typedef Dune::PDELab::PowerGridFunction<A,5> B5;
144   B5 b5(a,a,a,a,a);
145   typedef Dune::PDELab::PowerGridFunction<A,6> B6;
146   B6 b6(a,a,a,a,a,a);
147   typedef Dune::PDELab::PowerGridFunction<A,7> B7;
148   B7 b7(a,a,a,a,a,a,a);
149   typedef Dune::PDELab::PowerGridFunction<A,8> B8;
150   B8 b8(a,a,a,a,a,a,a,a);
151   typedef Dune::PDELab::PowerGridFunction<A,9> B9;
152   B9 b9(a,a,a,a,a,a,a,a,a);
153 
154   // test composite
155   typedef Dune::PDELab::CompositeGridFunction<A,A,A,A,A,A,A,A,A> C9;
156   C9 c9(a,a,a,a,a,a,a,a,a);
157   typedef Dune::PDELab::CompositeGridFunction<B10,B2> C2;
158   C2 c2(b10,b2);
159   typedef Dune::PDELab::CompositeGridFunction<B10,B2,B3> C3;
160   C3 c3(b10,b2,b3);
161   typedef Dune::PDELab::CompositeGridFunction<B10,B2,B3,B4> C4;
162   C4 c4(b10,b2,b3,b4);
163   typedef Dune::PDELab::CompositeGridFunction<B10,B2,B3,B4,B5> C5;
164   C5 c5(b10,b2,b3,b4,b5);
165   typedef Dune::PDELab::CompositeGridFunction<B10,B2,B3,B4,B5,B6> C6;
166   C6 c6(b10,b2,b3,b4,b5,b6);
167   typedef Dune::PDELab::CompositeGridFunction<B10,B2,B3,B4,B5,B6,B7> C7;
168   C7 c7(b10,b2,b3,b4,b5,b6,b7);
169   typedef Dune::PDELab::CompositeGridFunction<B10,B2,B3,B4,B5,B6,B7,B8> C8;
170   C8 c8(b10,b2,b3,b4,b5,b6,b7,b8);
171 
172   typedef Dune::PDELab::CompositeGridFunction<C2,C9> T;
173   T t(c2,c9);
174 
175   std::cout << "depth of T is " << Dune::TypeTree::TreeInfo<T>::depth << std::endl;
176   std::cout << "number of nodes in T is " << Dune::TypeTree::TreeInfo<T>::nodeCount << std::endl;
177   std::cout << "number of leaves in T is " << Dune::TypeTree::TreeInfo<T>::leafCount << std::endl;
178 
179   Dune::VTKWriter<GV> vtkwriter(gv,Dune::VTK::conforming);
180   Dune::PDELab::vtkwriter_tree_addvertexdata(vtkwriter,t);
181   vtkwriter.write("multi",Dune::VTK::ascii);
182 }
183 
184 template<int k, class GV>
testgridviewfunction(const GV & gv)185 void testgridviewfunction (const GV& gv)
186 {
187     enum { dim = GV::dimension };
188     using FEM = Dune::PDELab::QkLocalFiniteElementMap<GV,float,double,k>;
189     FEM fem(gv);
190     // make a grid function space
191     using GFS = Dune::PDELab::GridFunctionSpace<GV,FEM>;
192     GFS gfs(gv,fem);
193     // make vector
194     using Vector = Dune::PDELab::Backend::Vector<GFS, double>;
195     Vector x(gfs);
196     // make functions
197     typedef Dune::PDELab::DiscreteGridViewFunction<GFS,Vector> DiscreteFunction;
198     DiscreteFunction dgvf(gfs,x);
199     // interpolate linear function
200     using Domain = Dune::FieldVector<typename GV::ctype, GV::dimension>;
201     using std::pow;
202     auto f = [](const Domain& x) {return pow(x[0],k);};
203     Dune::PDELab::interpolate(f,gfs,x);
204     // make local functions
205     auto localf = localFunction(dgvf);
206     // iterate grid and evaluate local function
207     static const int maxDiffOrder = 1;
208     std::cout << "max diff order: " << maxDiffOrder << std::endl;
209     std::cout << "checking for:\n";
210     std::cout << "\tevaluate\n";
211     if (maxDiffOrder >= 1)
212         std::cout << "\tjacobian\n";
213     if (maxDiffOrder >= 2)
214         std::cout << "\thessian\n";
215     if (maxDiffOrder >= 3)
216         std::cout << "\tdiff(3)\n";
217     for (auto it=gv.template begin<0>(); it!=gv.template end<0>(); ++it)
218     {
219         localf.bind(*it);
220         Dune::FieldVector<double,1> value;
221         Dune::FieldMatrix<double,1,dim> jacobian;
222         Dune::FieldMatrix<double,dim,dim> hessian;
223         Dune::FieldVector<double,dim> pos(0.0);
224         auto gpos = it->geometry().global(pos);
225         value = localf(pos);
226         assert(std::abs(value - pow(gpos[0],k)) < 1e-6);
227         if (maxDiffOrder >= 1) {
228             jacobian = derivative(localf)(pos);
229             assert(std::abs(jacobian[0][0] - k*pow(gpos[0],k-1)) < 1e-6);
230             assert(std::abs(jacobian[0][1]) < 1e-6);
231         }
232 /* The following block is currently disabled as the tested QkLocalFiniteElementMap
233  * does not implement it and the maxDiffOrder based magic to not have these trigger
234  * compile errors had to be removed because of upstream changes in dune-localfunctions.
235         if (maxDiffOrder >= 2) {
236             hessian = derivative(derivative(localf))(pos);
237             assert(std::abs(hessian[0][0] - k*(k-1)*pow(gpos[0],k-2)) < 1e-6);
238             assert(std::abs(hessian[0][1]) < 1e-6);
239             assert(std::abs(hessian[1][1]) < 1e-6);
240             assert(std::abs(hessian[1][1]) < 1e-6);
241         }
242         if (maxDiffOrder >= 3)
243             derivative(
244                 derivative(
245                     derivative(localf)));
246 */
247         localf.unbind();
248     }
249 }
250 
main(int argc,char ** argv)251 int main(int argc, char** argv)
252 {
253   try{
254     //Maybe initialize Mpi
255     Dune::MPIHelper::instance(argc, argv);
256 
257     std::cout << "instantiate and evaluate some functions" << std::endl;
258 
259     // instantiate F and evaluate
260     F<float> f;
261     Dune::FieldVector<float,2> x;
262     x[0] = 1.0; x[1] = 2.0;
263     Dune::FieldVector<float,1> y;
264     f.evaluate(x,y);
265 
266     // instantiate G and evaluate
267     G<float> g;
268     Dune::FieldVector<float,2> u;
269     g.evaluate(x,u);
270 
271     // need a grid in order to test grid functions
272     Dune::FieldVector<double,2> L(1.0);
273     std::array<int,2> N(Dune::filledArray<2,int>(1));
274     Dune::YaspGrid<2> grid(L,N);
275     grid.globalRefine(6);
276 
277     // run algorithm on a grid
278     std::cout << "instantiate grid functions on a grid" << std::endl;
279     testgridfunction(grid.leafGridView(),F<Dune::YaspGrid<2>::ctype>());
280 
281     // run algorithm on a grid
282     std::cout << "testing vtk output" << std::endl;
283     testvtkexport(grid.leafGridView(),F<Dune::YaspGrid<2>::ctype>());
284     testfunctiontree(grid.leafGridView());
285 
286     testgridviewfunction<1>(grid.leafGridView());
287     testgridviewfunction<2>(grid.leafGridView());
288     // testgridviewfunction<3>(grid.leafGridView());
289 
290     // test passed
291     return 0;
292 
293   }
294   catch (Dune::Exception &e){
295     std::cerr << "Dune reported error: " << e << std::endl;
296     return 1;
297   }
298   catch (...){
299     std::cerr << "Unknown exception thrown!" << std::endl;
300     return 1;
301   }
302 }
303