lua interface can handle explicit time stepping all the way long

Numeric/fullMatrix.cpp

@@ -36,6 +36,18 @@ extern "C" {
                       std::complex<double> *alpha, std::complex<double> *a, int *lda, 
                       std::complex<double> *x, int *incx, std::complex<double> *beta, 
                       std::complex<double> *y, int *incy);
+  void F77NAME(daxpy)(int *n, double *alpha, double *x, int *incx, double *y, int *incy );
+
+
+}
+
+template<> 
+void fullVector<double>::axpy(fullVector<double> &x,double alpha)
+{
+  //  for (int i=0;i<_r;i++) _data[i] += alpha * x._data[i];
+  //  return;
+  int M = _r, INCX = 1, INCY = 1;
+  F77NAME(daxpy)(&M, &alpha, x._data,&INCX, _data, &INCY);
 }
 
 template<> 

Numeric/fullMatrix.h

@@ -71,6 +71,14 @@ class fullVector
     for(int i = 0; i < _r; ++i) s += _data[i] * other._data[i];
     return s;
   }
+  // y <- y + alpha * x
+  void axpy(fullVector<scalar> &x, scalar alpha=1.)
+#if !defined(HAVE_BLAS)
+  {
+    for (int i=0;i<_r;i++) _data[i] += alpha * x._data[i];
+  }
+#endif
+  ;
   void print(const char *name="") const 
   {
     printf("Printing vector %s:\n", name);

Solver/dgAlgorithm.cpp

 #include <set>
 #include "dgAlgorithm.h"
 #include "dgGroupOfElements.h"
+#include "dgSystemOfEquations.h"
 #include "dgConservationLaw.h"
 #include "GEntity.h"
 #include "MElement.h"
@@ -239,45 +240,58 @@ void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
 			      std::vector<dgGroupOfFaces*> &fGroups,		// group of interfacs
 			      std::vector<dgGroupOfFaces*> &bGroups,		// group of boundaries
 			      double h,				        // time-step
-			       std::vector<fullMatrix<double> *> &residu,					
-			       std::vector<fullMatrix<double> *> &sol,					
+			      dgDofContainer &sol,
+			      dgDofContainer &resd,			       
 			      int orderRK)				        // order of RK integrator
  {
+
    // U_{n+1}=U_n+h*(SUM_i a_i*K_i)
    // K_i=M^(-1)R(U_n+b_i*K_{i-1})
    
-   /*
    double a[4] = {h/6.0,h/3.0,h/3.0,h/6.0};
    double b[4] = {0.,h/2.0,h/2.0,h};
    
-     fullMatrix<double> K(sol);
+   // fullMatrix<double> K(sol);
    // Current updated solution
-     fullMatrix<double> Unp(sol);
+   // fullMatrix<double> Unp(sol);
    fullMatrix<double> residuEl, KEl;
    fullMatrix<double> iMassEl;
    
-     int nbNodes = eGroups[0]->getNbNodes();
-     int nbFields = sol.size2()/eGroups[0]->getNbElements(); 
+   int nbFields = claw.nbFields();
+   
+   dgDofContainer K   (eGroups,claw);
+   dgDofContainer Unp (eGroups,claw);
+
+   K._data->scale(0.0);
+   K._data->axpy(*(sol._data));
+   Unp._data->scale(0.0);
+   Unp._data->axpy(*(sol._data));
 
    for(int j=0; j<orderRK;j++){
-     if(j!=0){
-     K.scale(b[j]);
-     K.add(sol);
-     }
-     this->residual(claw,eGroups,fGroups,bGroups,K,residu);
-     K.scale(0.);
-     for(int i=0;i<eGroups[0]->getNbElements();i++) {
-     residuEl.setAsProxy(residu,i*nbFields,nbFields);
-     KEl.setAsProxy(K,i*nbFields,nbFields);
-     iMassEl.setAsProxy(eGroups[0]->getInverseMassMatrix(),i*nbNodes,nbNodes);
+     if(j){
+       K._data->scale(b[j]);
+       K._data->axpy(*(sol._data));
+     }
+     //     printf("iter %d sol size = %d\n",j,sol._dataSize);
+     this->residual(claw,eGroups,fGroups,bGroups,K._dataProxys,resd._dataProxys);
+     K._data->scale(0.);
+     for(int k=0;k < eGroups.size();k++) {
+       int nbNodes = eGroups[k]->getNbNodes();
+       for(int i=0;i<eGroups[k]->getNbElements();i++) {
+	 residuEl.setAsProxy(*(resd._dataProxys[k]),i*nbFields,nbFields);
+	 KEl.setAsProxy(*(K._dataProxys[k]),i*nbFields,nbFields);
+	 iMassEl.setAsProxy(eGroups[k]->getInverseMassMatrix(),i*nbNodes,nbNodes);
 	 iMassEl.mult(residuEl,KEl);
        }
-     Unp.add(K,a[j]);
      }
-     sol=Unp;
+     Unp._data->axpy(*(K._data),a[j]);
    }
    
-   */
+   for (int i=0;i<sol._dataSize;i++){
+     //     printf("tempSol[%d] = %g\n",i,(*tempSol._data)(i));
+     //     memcp
+     (*sol._data)(i)=(*Unp._data)(i);
+   }
  }
 
 void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (maybe useless here)
@@ -409,8 +423,8 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
     dgGroupOfFaces &faces = *fGroups[i];
     int iGroupLeft = -1, iGroupRight = -1;
     for(size_t j=0;j<eGroups.size() ; j++) {
-      if (eGroups[i] == &faces.getGroupLeft())iGroupLeft = j;
-      if (eGroups[i] == &faces.getGroupRight())iGroupRight= j;
+      if (eGroups[j] == &faces.getGroupLeft())iGroupLeft = j;
+      if (eGroups[j] == &faces.getGroupRight())iGroupRight= j;
     }
     
     fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
@@ -424,8 +438,8 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
     dgGroupOfFaces &faces = *bGroups[i];
     int iGroupLeft = -1, iGroupRight = -1;
     for(size_t j=0;j<eGroups.size() ; j++) {
-      if (eGroups[i] == &faces.getGroupLeft())iGroupLeft = j;
-      if (eGroups[i] == &faces.getGroupRight())iGroupRight= j;
+      if (eGroups[j] == &faces.getGroupLeft())iGroupLeft = j;
+      if (eGroups[j] == &faces.getGroupRight())iGroupRight= j;
     }
     fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
     fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);

Solver/dgAlgorithm.h

@@ -7,6 +7,7 @@ class GModel;
 class dgGroupOfElements;
 class dgGroupOfFaces;
 class dgConservationLaw;
+class dgDofContainer;
 class dgTerm;
 
 class dgAlgorithm {
@@ -44,8 +45,8 @@ class dgAlgorithm {
 		   std::vector<dgGroupOfFaces*> &fGroups,  // group of interfacs
 		   std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
 		   double h,	
-		   std::vector<fullMatrix<double> *> &residu,
-		   std::vector<fullMatrix<double> *> &sol,
+		   dgDofContainer &residu,
+		   dgDofContainer &sol,
 		   int orderRK=4);
   static void multAddInverseMassMatrix ( /*dofManager &dof,*/
 					const dgGroupOfElements & group,

Solver/dgSystemOfEquations.cpp

@@ -13,9 +13,7 @@ const char dgSystemOfEquations::className[] = "dgSystemOfEquations";
 // Define the methods we will expose to Lua
 #define _method(class, name) {#name, &class::name}
 Luna<dgSystemOfEquations>::RegType dgSystemOfEquations::methods[] = {
-   _method(dgSystemOfEquations, getRHS),
-   _method(dgSystemOfEquations, computeRHS),
-   _method(dgSystemOfEquations, getSolution),
+   _method(dgSystemOfEquations, RK44),
    _method(dgSystemOfEquations, exportSolution),
    _method(dgSystemOfEquations, setConservationLaw),
    _method(dgSystemOfEquations, setOrder),
@@ -56,6 +54,7 @@ dgSystemOfEquations::dgSystemOfEquations(lua_State *L){
   if (!obj)throw;
   _gm = obj->getGModel();
   _dimension = _gm->getNumRegions() ? 3 : 2;
+  _solution = 0;
 }
 
 // set the conservation law as a string (for now)
@@ -98,8 +97,8 @@ int dgSystemOfEquations::addBoundaryCondition (lua_State *L){
 int dgSystemOfEquations::L2Projection (lua_State *L){
   dgConservationLawL2Projection Law(std::string(luaL_checkstring(L, 1)),*_claw);
   for (int i=0;i<_elementGroups.size();i++){
-    _algo->residualVolume(Law,*_elementGroups[i],*_solutionProxys[i],*_rightHandSideProxys[i]);
-    _algo->multAddInverseMassMatrix(*_elementGroups[i],*_rightHandSideProxys[i],*_solutionProxys[i]);
+    _algo->residualVolume(Law,*_elementGroups[i],*_solution->_dataProxys[i],*_rightHandSide->_dataProxys[i]);
+    _algo->multAddInverseMassMatrix(*_elementGroups[i],*_rightHandSide->_dataProxys[i],*_solution->_dataProxys[i]);
   }
 }
 
@@ -113,60 +112,35 @@ int dgSystemOfEquations::setup(lua_State *L){
 		     _faceGroups,
 		     _boundaryGroups);
   // compute the total size of the soution
-  _dataSize = 0;
-  int nbFields = _claw->nbFields();
-  for (int i=0;i<_elementGroups.size();i++){
-    int nbNodes    = _elementGroups[i]->getNbNodes();
-    int nbElements = _elementGroups[i]->getNbElements();
-    _dataSize += nbNodes*nbFields*nbElements;
-  }
-  // allocate the big vectors
-  _solution      = new double [_dataSize];
-  _rightHandSide = new double [_dataSize];
-  // create proxys for each group
-  int offset = 0;
-  _solutionProxys.resize(_elementGroups.size());
-  _rightHandSideProxys.resize(_elementGroups.size());
-  for (int i=0;i<_elementGroups.size();i++){    
-    int nbNodes    = _elementGroups[i]->getNbNodes();
-    int nbElements = _elementGroups[i]->getNbElements();
-    _solutionProxys[i] = new fullMatrix<double> (&_solution[offset*sizeof(double)],nbNodes, nbFields*nbElements);
-    _rightHandSideProxys[i] = new fullMatrix<double> (&_rightHandSide[offset*sizeof(double)],nbNodes, nbFields*nbElements);
-    offset += nbNodes*nbFields*nbElements;
-  }
+  _solution = new dgDofContainer(_elementGroups,*_claw);
+  _rightHandSide = new dgDofContainer(_elementGroups,*_claw);
+  //  printf("aaaaaaaaaaaaa\n");
 }
 
 
-int dgSystemOfEquations::getSolution(lua_State *L){
-  lua_pushlightuserdata (L, _solution);
-  return 1;
-}
-
-int dgSystemOfEquations::computeRHS(lua_State *L){
-  _algo->residual(*_claw, _elementGroups, _faceGroups, _boundaryGroups, _solutionProxys, _rightHandSideProxys);
-  lua_pushlightuserdata (L, _solution);
-  return 1;
-}
-
-int dgSystemOfEquations::getRHS(lua_State *L){
-  lua_pushlightuserdata (L, _rightHandSide);
+int dgSystemOfEquations::RK44(lua_State *L){
+  double dt = luaL_checknumber(L, 1);
+  _algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt,  *_solution, *_rightHandSide);
+  double normSolution = _solution->_data->norm();
+  lua_pushnumber (L, normSolution);
   return 1;
 }
 
 int dgSystemOfEquations::exportSolution(lua_State *L){
   std::string outputFile(luaL_checkstring(L, 1));
   export_solution_as_is(outputFile);
+  return 0;
 }
 #endif // HAVE_LUA
 
 dgSystemOfEquations::~dgSystemOfEquations(){
   for (int i=0;i<_elementGroups.size();i++){
-    delete _solutionProxys[i];
-    delete _rightHandSideProxys[i];
     delete _elementGroups[i];
   }
-  delete [] _solution;
-  delete [] _rightHandSide;
+  if (_solution){
+    delete _solution;
+    delete _rightHandSide;
+  }
 }
 
 void dgSystemOfEquations::export_solution_as_is (const std::string &name){
@@ -229,3 +203,31 @@ void dgSystemOfEquations::export_solution_as_is (const std::string &name){
   pv->getData()->writeMSH(name+".msh", false); 
   delete pv;
 }
+
+dgDofContainer::dgDofContainer (std::vector<dgGroupOfElements*> &elementGroups, const dgConservationLaw &claw){
+  _dataSize = 0;
+  int nbFields = claw.nbFields();
+  for (int i=0;i<elementGroups.size();i++){
+    int nbNodes    = elementGroups[i]->getNbNodes();
+    int nbElements = elementGroups[i]->getNbElements();
+    _dataSize += nbNodes*nbFields*nbElements;
+  }
+  // allocate the big vectors
+  _data      = new fullVector<double>(_dataSize);
+  // create proxys for each group
+  int offset = 0;
+  _dataProxys.resize(elementGroups.size());
+  for (int i=0;i<elementGroups.size();i++){    
+    int nbNodes    = elementGroups[i]->getNbNodes();
+    int nbElements = elementGroups[i]->getNbElements();
+    _dataProxys[i] = new fullMatrix<double> (&(*_data)(offset*sizeof(double)),nbNodes, nbFields*nbElements);
+    offset += nbNodes*nbFields*nbElements;
+  }  
+  //  printf("datasize = %d\n",_dataSize);
+}
+
+dgDofContainer::~dgDofContainer (){
+  if (!_dataSize)return;
+  for (int i=0;i<_dataProxys.size();++i) delete _dataProxys[i];
+  delete _data;
+}

Solver/dgSystemOfEquations.h

@@ -19,6 +19,18 @@ extern "C" {
 #include "luna.h"
 #endif // HAVE LUA
 
+struct dgDofContainer {
+private:
+  dgDofContainer (const dgDofContainer&);  
+public:
+  int _dataSize; // the full data size i.e. concerning all groups
+  fullVector<double> * _data; // the full data itself
+  std::vector<fullMatrix<double> *> _dataProxys; // proxys 
+  dgDofContainer (std::vector<dgGroupOfElements*> &groups, const dgConservationLaw &claw);
+  ~dgDofContainer ();  
+};
+
+
 class dgSystemOfEquations {
   // the mesh and the model
   GModel *_gm;
@@ -31,13 +43,9 @@ class dgSystemOfEquations {
   int _order;
   // dimension of the problem
   int _dimension;
-  // solution and righ hand sides as a large arrays of doubles
-  int _dataSize;
-  double * _solution;
-  double * _rightHandSide;
-  // proxis of the solution and of the right hand side for each group;
-  std::vector<fullMatrix<double> *> _solutionProxys;
-  std::vector<fullMatrix<double> *> _rightHandSideProxys;
+  // solution and righ hand sides
+  dgDofContainer *_solution;
+  dgDofContainer *_rightHandSide;
   // groups of elements (volume terms)
   std::vector<dgGroupOfElements*> _elementGroups;
   // groups of faces (interface terms)
@@ -51,19 +59,17 @@ public:
   static const char className[];
   static Luna<dgSystemOfEquations>::RegType methods[];
   static void Register(lua_State *L);
-  int computeRHS      (lua_State *L); // get the Right Hand
-  int getRHS      (lua_State *L); // get the Right Hand
-  int getSolution (lua_State *L); // get the Solution
   int exportSolution (lua_State *L); // export the solution
   int setOrder (lua_State *L); // set the polynomial order
   int setConservationLaw (lua_State *L); // set the conservationLaw
   int addBoundaryCondition (lua_State *L); // add a boundary condition : "physical name", "type", [options]
   int setup (lua_State *L); // setup the groups and allocate
   int L2Projection (lua_State *L); // assign the solution to a given function
+  int RK44 (lua_State *L); // assign the solution to a given function
   dgSystemOfEquations(lua_State *L);
 #endif // HAVE LUA
   inline const fullMatrix<double> getSolutionProxy (int iGroup, int iElement){
-    return fullMatrix<double> ( *_solutionProxys [iGroup] ,
+    return fullMatrix<double> ( *_solution->_dataProxys [iGroup] ,
 				iElement * _claw->nbFields(),_claw->nbFields());
   }
   void export_solution_as_is (const std::string &fileName);

contrib/DiscreteIntegration/Integration3D.h

@@ -43,11 +43,11 @@ class DI_Point
   DI_Point (double x, double y, double z, const double ls) : x_(x), y_(y), z_(z) {addLs(ls);}
   DI_Point (double x, double y, double z, const DI_Element *e,
             const std::vector<const gLevelset *> &RPNi) : x_(x), y_(y), z_(z) {computeLs(e, RPNi);}
-  DI_Point(const DI_Point &p) : x_(p.x()), y_(p.y()), z_(p.z()) {Ls.clear(); Ls = p.Ls;}
+  //  DI_Point(const DI_Point &p) : x_(p.x()), y_(p.y()), z_(p.z()) {Ls.clear(); Ls = p.Ls;}
   // assignment
   DI_Point & operator=(const DI_Point & rhs);
   // destructor
-  virtual ~DI_Point () {}
+  //  virtual ~DI_Point () {}
   // add a levelset value (adjusted to 0 if ls<ZERO_LS_TOL) into the vector Ls
   void addLs (const double ls);
   // add a levelset value evaluated into the element e