dg : faster integrals on boundary

Solver/functionDerivator.cpp

 #include "functionDerivator.h"
 #include "function.h"
 
-void functionDerivator::compute() {
+const fullMatrix<double> &functionDerivator::compute() {
   _xRef = _x.get();
   _fRef = _f.get();
-  _dfdx = fullMatrix<double>(_fRef.size1(),_fRef.size2()*_xRef.size2());
+  _dfdx.resize(_fRef.size1(),_fRef.size2()*_xRef.size2(), false);
   for (int j=0;j<_xRef.size2();j++) {
-    _xDx = _xRef;
+    fullMatrix<double> &x =_x.set();
     for (int i=0;i<_fRef.size1();i++)
-      _xDx(i,j) += _epsilon;
-    _x.set()=_xDx;
+      x(i,j) += _epsilon;
     const fullMatrix<double> &f = _f.get();
-    for (int i=0; i<_fRef.size1(); i++) 
-      for (int k=0; k<_fRef.size2(); k++) 
+    for (int k=0; k<_fRef.size2(); k++) 
+      for (int i=0; i<_fRef.size1(); i++) 
         _dfdx(i, k*_xRef.size2()+j) = (f(i,k)-_fRef(i,k))/_epsilon;
+    for (int i=0;i<_fRef.size1();i++)
+      x(i,j) = _xRef(i,j);
   }
-  _x.set()=_xRef;
+  _x.set();
+  return _dfdx;
 };

Solver/functionDerivator.h

@@ -6,7 +6,7 @@
 
 class functionDerivator {
   dataCacheDouble &_f,&_x;
-  fullMatrix<double> _fRef, _xRef,_xDx, _dfdx;
+  fullMatrix<double> _fRef, _xRef, _dfdx;
   double _epsilon;
   public:
   functionDerivator (dataCacheDouble &f, dataCacheDouble &x, double epsilon):
@@ -14,7 +14,7 @@ class functionDerivator {
     _x(x),
     _epsilon(epsilon)
     {}
-  void compute();
+  const fullMatrix<double> &compute();
   inline double get(int iQP, int iF, int iX) 
   {
     return _dfdx(iQP, iF*_xRef.size2()+iX);