Modified integrate_one_over_r test for new features in the

QGaussOneOverR class.

git-svn-id: https://svn.dealii.org/trunk@18550 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/tests/codim_one/integrate_one_over_r.cc b/tests/codim_one/integrate_one_over_r.cc
index 48f261602bb700a9e899cbb8f5dab1b5255a79ac..97bda0d1913826fe42d4d4a869d7c340ce7f1996 100644 (file)
--- a/tests/codim_one/integrate_one_over_r.cc
+++ b/tests/codim_one/integrate_one_over_r.cc
@@ -22,6 +22,7 @@
 
 // all include files needed for the program
 #include <base/quadrature_lib.h>
+#include <base/geometry_info.h>
 
 
 #include <fstream>
@@ -41,7 +42,7 @@ ofstream logfile("integrate_one_over_r/output");
 int main()
 {
     deallog.attach(logfile);
-    deallog.depth_console(0);
+    deallog.depth_console(1);
     deallog<< std::fixed;
     
     deallog << endl 
@@ -50,6 +51,8 @@ int main()
            << "the distance from (x,y) to four vertices of the square." << endl 
            << endl;
 
+    double eps = 1e-10;
+
     for(unsigned int m=1; m<7; ++m) {
        deallog << " =========Quadrature Order: " << m 
                << " =============================== " << endl;
@@ -58,11 +61,13 @@ int main()
            deallog << " ===============Vertex Index: " << index 
                    << " ============================= " << endl;
            QGaussOneOverR<2> quad(m, index);
+           QGaussOneOverR<2> quad_de(m, index, true);
            for (unsigned int i=0; i<6; ++i) {
                for (unsigned int j=0; j<6; ++j) {
                   
                    double exact_integral  = exact_integral_one_over_r(index, i,j);
                    double approx_integral = 0;
+                   double approx_integral_2 = 0;
                    
                    for (unsigned int q=0; q< quad.size(); ++q) {
                        double x = quad.point(q)[0];
@@ -70,12 +75,21 @@ int main()
                        approx_integral += ( pow(x, (double)i) * 
                                             pow(y, (double)j) * 
                                             quad.weight(q) );
+                       approx_integral_2 += ( pow(x, (double)i) * 
+                                              pow(y, (double)j) * 
+                                              quad_de.weight(q) /
+                                              (quad_de.point(q)-GeometryInfo<2>::unit_cell_vertex(index)).norm());
                    }
 
                    deallog << "f(x,y) = x^" << i
                            << " y^" << j 
                            << ", Error = " 
-                           << approx_integral - exact_integral << endl;
+                           << approx_integral - exact_integral;
+                   if (std::abs(approx_integral-approx_integral_2) < eps)
+                       deallog << endl;
+                   else 
+                       deallog <<  ", desing: "  
+                               << approx_integral_2 - exact_integral << endl;
                }
            }
        }