1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 
4 #include <algorithm>
5 #include <limits>
6 #include <iostream>
7 #include <type_traits>
8 
9 #include <config.h>
10 
11 #include <dune/common/quadmath.hh>
12 #include <dune/geometry/referenceelements.hh>
13 #include <dune/geometry/quadraturerules.hh>
14 #include <dune/geometry/quadraturerules/compositequadraturerule.hh>
15 #include <dune/geometry/refinement.hh>
16 
17 bool success = true;
18 
19 template <class ctype, std::enable_if_t<std::numeric_limits<ctype>::is_specialized, int> = 0>
20 ctype eps ()
21 {
22   return std::numeric_limits<ctype>::epsilon();
23 }
24 
25 template <class ctype, std::enable_if_t<!std::numeric_limits<ctype>::is_specialized, int> = 0>
26 ctype eps ()
27 {
28   return std::numeric_limits<float>::epsilon();
29 }
30 
31 /*
32    This is a simple accuracy test on the reference element. It integrates
33    x^p and y^p with the quadrature rule of order p, which should give
34    an exact result.
35  */
36 
37 /*
38    Exact (analytical) solution on different reference elements.
39  */
40 
41 template <class ctype, int dim>
42 ctype analyticalSolution (Dune::GeometryType t, int p, int direction )
43 {
44   using Dune::GeometryType;
45   ctype exact=0;
46 
47   if (t.isCube())
48   {
49     exact = ctype( 1 )/(p+1);
50     return exact;
51   }
52 
53   if (t.isSimplex())
54   {
55     /* 1/(prod(k=1..dim,(p+k)) */
56     exact = ctype( 1 );
57     for( int k = 1; k <= dim; ++k )
58       exact *= p+k;
59     exact = ctype( 1 ) / exact;
60     return exact;
61   }
62 
63   if (t.isPrism())
64   {
65     const int pdim = (dim > 0 ? dim-1 : 0);
66     if( direction < dim-1 )
67     {
68       Dune::GeometryType nt = Dune::GeometryTypes::simplex( dim-1 );
69       if( dim > 0 )
70         exact = analyticalSolution< ctype, pdim >( nt, p, direction );
71       else
72         exact = ctype( 1 );
73     }
74     else
75       exact = ctype( 1 ) / ctype( Dune::Factorial< pdim >::factorial * (p+1));
76     return exact;
77   }
78 
79   if (t.isPyramid())
80   {
81     switch( direction )
82     {
83     case 0 :
84     case 1 :
85       exact = ctype( 1 )/((p+3)*(p+1));
86       break;
87     case 2 :
88       exact = ctype( 2 )/((p+1)*(p+2)*(p+3));
89       break;
90     };
91     return exact;
92   }
93 
94   DUNE_THROW(Dune::NotImplemented, __func__ << " for " << t);
95   return exact;
96 }
97 
98 template<class QuadratureRule>
checkQuadrature(const QuadratureRule & quad)99 void checkQuadrature(const QuadratureRule &quad)
100 {
101   using std::pow; using std::abs;
102   using namespace Dune;
103   typedef typename QuadratureRule::CoordType ctype;
104   const unsigned int dim = QuadratureRule::d;
105   const unsigned int p = quad.order();
106   const Dune::GeometryType& t = quad.type();
107   FieldVector<ctype,dim> integral(0);
108   for (typename QuadratureRule::iterator qp=quad.begin(); qp!=quad.end(); ++qp)
109   {
110     // pos of integration point
111     const FieldVector< ctype, dim > &x = qp->position();
112     const ctype weight = qp->weight();
113 
114     for (unsigned int d=0; d<dim; d++)
115       integral[d] += weight*pow(x[d], p);
116   }
117 
118   ctype maxRelativeError = 0;
119   for(unsigned int d=0; d<dim; d++)
120   {
121     ctype exact = analyticalSolution<ctype,dim>(t,p,d);
122     ctype relativeError = abs(integral[d] - exact) / (abs(integral[d]) + abs(exact));
123     if (relativeError > maxRelativeError)
124       maxRelativeError = relativeError;
125   }
126   ctype epsilon = pow(2.0,p)*p*eps<ctype>();
127   if (p==0)
128     epsilon = 2.0*eps<ctype>();
129   if (maxRelativeError > epsilon) {
130     std::cerr << "Error: Quadrature for " << t << " and order=" << p << " failed" << std::endl;
131     for (unsigned int d=0; d<dim; d++)
132     {
133       ctype exact = analyticalSolution<ctype,dim>(t,p,d);
134       ctype relativeError = abs(integral[d] - exact) / (abs(integral[d]) + abs(exact));
135       std::cerr << "       relative error " << relativeError << " in direction " << d << " (exact = " << exact << " numerical = " << integral[d] << ")" << std::endl;
136     }
137     success = false;
138   }
139 }
140 
141 template<class QuadratureRule>
checkWeights(const QuadratureRule & quad)142 void checkWeights(const QuadratureRule &quad)
143 {
144   using std::abs;
145   typedef typename QuadratureRule::CoordType ctype;
146   const unsigned int dim = QuadratureRule::d;
147   const unsigned int p = quad.order();
148   const Dune::GeometryType& t = quad.type();
149   typedef typename QuadratureRule::iterator QuadIterator;
150   ctype volume = 0;
151   QuadIterator qp = quad.begin();
152   QuadIterator qend = quad.end();
153   for (; qp!=qend; ++qp)
154   {
155     volume += qp->weight();
156   }
157   if (abs(volume - Dune::ReferenceElements<ctype, dim>::general(t).volume())
158       > quad.size()*eps<ctype>())
159   {
160     std::cerr << "Error: Quadrature for " << t << " and order=" << p
161               << " does not sum to volume of RefElem" << std::endl;
162     std::cerr << "\tSums to " << volume << "( RefElem.volume() = "
163               << Dune::ReferenceElements<ctype, dim>::general(t).volume()
164               << ")" << "(difference " << volume -
165     Dune::ReferenceElements<ctype, dim>::general(t).volume()
166               << ")" << std::endl;
167     success = false;
168   }
169 }
170 
171 template<class ctype, int dim>
check(Dune::GeometryType type,unsigned int maxOrder,Dune::QuadratureType::Enum qt=Dune::QuadratureType::GaussLegendre)172 void check(Dune::GeometryType type,
173            unsigned int maxOrder,
174            Dune::QuadratureType::Enum qt = Dune::QuadratureType::GaussLegendre)
175 {
176   typedef Dune::QuadratureRule<ctype, dim> Quad;
177   std::string qt_str = "";
178   switch (qt) {
179     case Dune::QuadratureType::GaussLegendre: qt_str = "GaussLegendre"; break;
180     case Dune::QuadratureType::GaussJacobi_1_0: qt_str = "GaussJacobi_1_0"; break;
181     case Dune::QuadratureType::GaussJacobi_2_0: qt_str = "GaussJacobi_2_0"; break;
182     case Dune::QuadratureType::GaussJacobi_n_0: qt_str = "GaussJacobi_n_0"; break;
183     case Dune::QuadratureType::GaussLobatto: qt_str = "GaussLobatto"; break;
184     case Dune::QuadratureType::GaussRadauLeft: qt_str = "GaussRadauLeft"; break;
185     case Dune::QuadratureType::GaussRadauRight: qt_str = "GaussRadauRight"; break;
186     default: qt_str = "unknown";
187   }
188   std::cout << "check(Quadrature of type " << qt_str << ")" << std::endl;
189   for (unsigned int p=0; p<=maxOrder; ++p)
190   {
191     const Quad & quad = Dune::QuadratureRules<ctype,dim>::rule(type, p, qt);
192     if (quad.type() != type || unsigned(quad.order()) < p) {
193       std::cerr << "Error: Type mismatch! Requested Quadrature for " << type
194                 << " and order=" << p << "." << std::endl
195                 << "\tGot Quadrature for " << quad.type() << " and order="
196                 << quad.order() << std::endl;
197       success = false;
198       return;
199     }
200     checkWeights(quad);
201     checkQuadrature(quad);
202   }
203   if (dim>0 && (dim>3 || type.isCube() || type.isSimplex()))
204   {
205     type = type.isCube() ? Dune::GeometryTypes::cube(dim-1) : Dune::GeometryTypes::simplex(dim-1);
206     check<ctype,((dim==0) ? 0 : dim-1)>(type , maxOrder, qt);
207   }
208 }
209 
210 template<class ctype, int dim>
checkCompositeRule(Dune::GeometryType type,unsigned int maxOrder,unsigned int maxRefinement,Dune::QuadratureType::Enum qt=Dune::QuadratureType::GaussLegendre)211 void checkCompositeRule(Dune::GeometryType type,
212                         unsigned int maxOrder,
213                         unsigned int maxRefinement,
214                         Dune::QuadratureType::Enum qt = Dune::QuadratureType::GaussLegendre)
215 {
216   typedef Dune::QuadratureRule<ctype, dim> BaseQuad;
217   typedef Dune::CompositeQuadratureRule<ctype, dim> Quad;
218 
219   for (unsigned int p=0; p<=maxOrder; ++p)
220   {
221     const BaseQuad& baseQuad = Dune::QuadratureRules<ctype,dim>::rule(type, p, qt);
222     Quad quad = Quad(baseQuad, Dune::refinementLevels(maxRefinement));
223 
224     checkWeights(quad);
225     checkQuadrature(quad);
226   }
227   if (dim>0 && (dim>3 || type.isCube() || type.isSimplex()))
228   {
229     type = type.isCube() ? Dune::GeometryTypes::cube(dim-1) : Dune::GeometryTypes::simplex(dim-1);
230     check<ctype,((dim==0) ? 0 : dim-1)>(type, maxOrder, qt);
231   }
232 }
233 
main(int argc,char ** argv)234 int main (int argc, char** argv)
235 {
236   unsigned int maxOrder = 45;
237   if (argc > 1)
238   {
239     maxOrder = std::atoi(argv[1]);
240     std::cout << "maxOrder = " << maxOrder << std::endl;
241   }
242   try {
243     check<double,4>(Dune::GeometryTypes::cube(4), maxOrder);
244     check<double,4>(Dune::GeometryTypes::cube(4), std::min(maxOrder, unsigned(31)),
245                     Dune::QuadratureType::GaussLobatto);
246     check<double,4>(Dune::GeometryTypes::cube(4), std::min(maxOrder, unsigned(30)),
247                     Dune::QuadratureType::GaussRadauLeft);
248     check<double,4>(Dune::GeometryTypes::cube(4), std::min(maxOrder, unsigned(30)),
249                     Dune::QuadratureType::GaussRadauRight);
250     check<double,4>(Dune::GeometryTypes::simplex(4), maxOrder);
251 
252 #if HAVE_LAPACK
253     check<double,4>(Dune::GeometryTypes::simplex(4), maxOrder, Dune::QuadratureType::GaussJacobi_n_0);
254 #else
255     std::cout << "Skip GaussJacobi_n_0 tests as LAPACK has not been found" << std::endl;
256 #endif
257 
258     check<double,3>(Dune::GeometryTypes::prism, maxOrder);
259     check<double,3>(Dune::GeometryTypes::pyramid, maxOrder);
260 
261     unsigned int maxRefinement = 4;
262 
263     checkCompositeRule<double,2>(Dune::GeometryTypes::triangle, maxOrder, maxRefinement);
264 
265 #if HAVE_QUADMATH
266     check<Dune::Float128,4>(Dune::GeometryTypes::cube(4), maxOrder);
267     check<Dune::Float128,4>(Dune::GeometryTypes::cube(4), std::min(maxOrder, unsigned(31)),
268                     Dune::QuadratureType::GaussLobatto);
269     check<Dune::Float128,4>(Dune::GeometryTypes::cube(4), std::min(maxOrder, unsigned(30)),
270                     Dune::QuadratureType::GaussRadauLeft);
271     check<Dune::Float128,4>(Dune::GeometryTypes::cube(4), std::min(maxOrder, unsigned(30)),
272                     Dune::QuadratureType::GaussRadauRight);
273 #else
274     std::cout << "Skip Float128 tests as Quadmath is not supported" << std::endl;
275 #endif
276   }
277   catch( const Dune::Exception &e )
278   {
279     std::cerr << e << std::endl;
280     throw;
281   }
282   catch (...) {
283     std::cerr << "Generic exception!" << std::endl;
284     throw;
285   }
286 
287   return success ? 0 : 1;
288 }
289