add the fd check for the first derivative

[[Imported from SVN: r9278]]

test/finitedifferencecheck.hh

+#ifndef FINITE_DIFFERENCE_CHECK_HH
+#define FINITE_DIFFERENCE_CHECK_HH
+
+
+template <class AssemblerType, class VectorType>
+void finiteDifferenceCheck(AssemblerType& assembler, const VectorType& x)
+{
+    VectorType rhs(x.size());
+    rhs = 0;
+    double eps = 1e-6;
+
+    const int blocksize = VectorType::block_type::size;
+
+    // ///////////////////////////////////////////////////////////
+    //   Assemble gradient and hessian at the current position
+    // ///////////////////////////////////////////////////////////
+
+    assembler.assembleMatrix(x, rhs);
+
+
+    // ///////////////////////////////////////////////////////////
+    //   Compare rhs with a finite-difference approximation
+    // ///////////////////////////////////////////////////////////
+    for (int i=0; i<rhs.size(); i++) {
+
+        for (int j=0; j<blocksize; j++) {
+
+            VectorType value0 = x;
+            VectorType value1 = x;
+
+            value0[i][j] -= eps;
+            value1[i][j] += eps;
+
+            double fd = (assembler.computeEnergy(value1) - assembler.computeEnergy(value0)) / (2*eps);
+
+            // relative error. NOTE: assembleMatrix currently assembles the NEGATIVE gradient!
+            double error = (-rhs[i][j] - fd)/(-rhs[i][j] + fd);
+            if (std::fabs(-rhs[i][j] - fd) > 1) {
+                std::cout << "ERROR! gradient: " << rhs[i][j] << ",   fd: " << fd << std::endl;
+                exit(1);
+            }
+                
+            //std::cout << "gradient: " << rhs[i][j] << ",   fd: " << fd << std::endl;
+
+        }
+
+    }
+    //exit(0);
+    // ///////////////////////////////////////////////////////////
+    //   Compare hessian with a finite-difference approximation
+    // ///////////////////////////////////////////////////////////
+
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,blocksize,blocksize> > MatrixType;
+
+    const MatrixType& M = *assembler.getMatrix();
+
+    for (int i=0; i<M.N(); i++) {
+
+        typename MatrixType::row_type::ConstIterator cIt    = M[i].begin();
+        typename MatrixType::row_type::ConstIterator cEndIt = M[i].end();
+
+        for (; cIt!=cEndIt; ++cIt) {
+
+            // Loop over all entries of the matrix block
+            for (int j=0; j<blocksize; j++) {
+
+                for (int k=0; k<blocksize; k++) {
+
+                    // Is the matrix symmetric?
+                    if ( std::abs((*cIt)[j][k] - M[cIt.index()][i][k][j]) > 1e-6) {
+
+                        std::cout << "Matrix asymmetry: " 
+                                  << (*cIt)[j][k] << " -- " << M[cIt.index()][i][k][j] << std::endl;
+                        //exit(1);
+                    }
+                    //assert( std::abs((*cIt)[j][k] - M[cIt.index()][i][k][j]) < 1e-6);
+
+                }
+
+            }
+
+        }
+
+    }        
+}
+
+#endif

test/neohooketest.cc

@@ -6,6 +6,7 @@
 
 #include "makesymmgrids.hh"
 #include "checksymmetry.hh"
+#include "finitedifferencecheck.hh"
 
 using namespace Dune;
 
@@ -33,7 +34,16 @@ int main() try
 
     // Check rhs for symmetry
     checkSymmetry(rhs, perm_tire2d);
-    
+
+    grid.globalRefine(2);
+    x.resize(grid.size(grid.maxLevel(),2));
+    x = 0;
+    for (int i=0; i<10; i++) {
+        x[i%x.size()][i%2] += i;
+        // Compare derivatives with a finite-difference approximation
+        finiteDifferenceCheck(assembler, x);
+    }
+
     return 0;
 } catch (Exception e) {