dg : add functionDerivator

Solver/CMakeLists.txt

@@ -18,6 +18,7 @@ set(SRC
   dgConservationLawShallowWater2d.cpp
   dgConservationLawWaveEquation.cpp
   function.cpp
+  functionDerivator.cpp
 )
 
 file(GLOB HDR RELATIVE ${CMAKE_SOURCE_DIR}/Solver *.h) 

Solver/dgMainWaveEquation.cpp

@@ -6,6 +6,7 @@
 #include "dgConservationLaw.h"
 #include "Gmsh.h"
 #include "function.h"
+#include "functionDerivator.h"
 
 
 #include "MElement.h"
@@ -53,20 +54,31 @@ int main(int argc, char **argv){
   dgConservationLaw *law = dgNewConservationLawWaveEquation();
   fullMatrix<double> sol(nbNodes,elementGroups[0]->getNbElements()*law->nbFields());
   fullMatrix<double> residu(nbNodes,elementGroups[0]->getNbElements()*law->nbFields());
+  /*{ // use this block to test functionDerivator
+    dgConservationLawInitialCondition initLaw;
+    dataCacheMap cacheMap;
+    dataCacheElement &cacheElement = cacheMap.getElement();
+    cacheMap.provideData("UVW").set(elementGroups[0]->getIntegrationPointsMatrix());
+    cacheElement.set(elementGroups[0]->getElement(0));
+    dataCacheDouble &xyz = cacheMap.get("XYZ");
+    dataCacheDouble *source = initLaw.newSourceTerm(cacheMap);
+    functionDerivator fd (*source, xyz, 1e-12);
+    xyz().print();
+    fd.compute();
+    exit(0);
+  }*/
   // initial condition
   {
     dgConservationLawInitialCondition initLaw;
     algo.residualVolume(initLaw,*elementGroups[0],sol,residu);
     algo.multAddInverseMassMatrix(*elementGroups[0],residu,sol);
   }
- 
 
   law->addBoundaryCondition("Left",dgNewBoundaryConditionWaveEquationWall());
   law->addBoundaryCondition("Right",dgNewBoundaryConditionWaveEquationWall());
   law->addBoundaryCondition("Top",dgNewBoundaryConditionWaveEquationWall());
   law->addBoundaryCondition("Bottom",dgNewBoundaryConditionWaveEquationWall());
   
-  
  print("output/p.pos",*elementGroups[0],&sol(0,0),0,3);
  print_vector("output/uv.pos",*elementGroups[0],&sol(0,0),1,3);
   for(int iT=0; iT<10000; iT++) {

Solver/functionDerivator.cpp

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

Solver/functionDerivator.h

+#ifndef _FUNCTION_DERIVATOR_H_
+#define _FUNCTION_DERIVATOR_H_
+
+#include "fullMatrix.h"
+#include "function.h"
+
+class functionDerivator {
+  dataCacheDouble &_f,&_x;
+  fullMatrix<double> _fRef, _xRef,_xDx, _dfdx;
+  double _epsilon;
+  public:
+  functionDerivator (dataCacheDouble &f, dataCacheDouble &x, double epsilon):
+    _f(f),
+    _x(x),
+    _epsilon(epsilon)
+    {}
+  void compute();
+  inline double get(int iQP, int iF, int iX) 
+  {
+    return _dfdx(iQP, iF*_xRef.size2()+iX);
+  }
+};
+
+#endif