dg : support multiple fields equations, implement wave equation and work

on shallow water equations

CMakeLists.txt

@@ -993,3 +993,9 @@ mark_as_advanced(BISON FLEX GMP_LIB GMSH_EXTRA_VERSION HDF5_LIB MAKEINFO
 
 add_executable(gmshdg EXCLUDE_FROM_ALL Solver/dgMainTest.cpp ${GMSH_SRC})
 target_link_libraries(gmshdg ${LINK_LIBRARIES})
+
+add_executable(gmshdgsw EXCLUDE_FROM_ALL Solver/dgMainShallowWater2d.cpp ${GMSH_SRC})
+target_link_libraries(gmshdgsw ${LINK_LIBRARIES})
+
+add_executable(gmshdgwave EXCLUDE_FROM_ALL Solver/dgMainWaveEquation.cpp ${GMSH_SRC})
+target_link_libraries(gmshdgwave ${LINK_LIBRARIES})

Solver/CMakeLists.txt

@@ -15,6 +15,8 @@ set(SRC
   dgAlgorithm.cpp
   dgConservationLaw.cpp
   dgConservationLawAdvection.cpp
+  dgConservationLawShallowWater2d.cpp
+  dgConservationLawWaveEquation.cpp
   function.cpp
 )
 

Solver/dgAlgorithm.cpp

@@ -223,9 +223,10 @@ void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
   fullMatrix<double> residuEl;
   fullMatrix<double> solEl;
   fullMatrix<double> iMassEl;
+  int nFields=sol.size2()/group.getNbElements();
   for(int i=0;i<group.getNbElements();i++) {
-    residuEl.setAsProxy(residu,i,1);
-    solEl.setAsProxy(sol,i,1);
+    residuEl.setAsProxy(residu,i*nFields,nFields);
+    solEl.setAsProxy(sol,i*nFields,nFields);
     iMassEl.setAsProxy(group.getInverseMassMatrix(),i*group.getNbNodes(),group.getNbNodes());
     solEl.gemm(iMassEl,residuEl);
   }
@@ -241,8 +242,6 @@ void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
 			fullMatrix<double> &sol,					
 			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})
 
@@ -256,6 +255,7 @@ void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
   fullMatrix<double> iMassEl;
 
   int nbNodes = eGroups[0]->getNbNodes();
+  int nbFields = sol.size2()/eGroups[0]->getNbElements(); 
 
   for(int j=0; j<orderRK;j++){
     if(j!=0){
@@ -265,8 +265,8 @@ void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
     this->residual(claw,eGroups,fGroups,bGroups,K,residu);
     K.scale(0.);
     for(int i=0;i<eGroups[0]->getNbElements();i++) {
-      residuEl.setAsProxy(residu,i,1);
-      KEl.setAsProxy(K,i,1);
+      residuEl.setAsProxy(residu,i*nbFields,nbFields);
+      KEl.setAsProxy(K,i*nbFields,nbFields);
       iMassEl.setAsProxy(eGroups[0]->getInverseMassMatrix(),i*nbNodes,nbNodes);
       iMassEl.mult(residuEl,KEl);
     }
@@ -372,6 +372,7 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
     fullMatrix<double> &solution, // solution
     fullMatrix<double> &residu) // residual
 {
+  int nbFields=claw.nbFields();
   dgAlgorithm algo;
   residu.scale(0);
   //volume term
@@ -381,19 +382,19 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
   //interface term
   for(size_t i=0;i<fGroups.size() ; i++) {
     dgGroupOfFaces &faces = *fGroups[i];
-    fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2);
-    fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*2);
-    faces.mapToInterface(1, solution, solution, solInterface);
+    fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
+    fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
+    faces.mapToInterface(nbFields, solution, solution, solInterface);
     algo.residualInterface(claw,faces,solInterface,solution,solution,residuInterface);
-    faces.mapFromInterface(1, residuInterface, residu, residu);
+    faces.mapFromInterface(nbFields, residuInterface, residu, residu);
   }
   //boundaries
   for(size_t i=0;i<bGroups.size() ; i++) {
     dgGroupOfFaces &faces = *bGroups[i];
-    fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements());
-    fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements());
-    faces.mapToInterface(1, solution, solution, solInterface);
+    fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
+    fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
+    faces.mapToInterface(nbFields, solution, solution, solInterface);
     algo.residualBoundary(claw,faces,solInterface,solution,residuInterface);
-    faces.mapFromInterface(1, residuInterface, residu, residu);
+    faces.mapFromInterface(nbFields, residuInterface, residu, residu);
   }
 }

Solver/dgConservationLaw.h

@@ -52,5 +52,8 @@ class dgConservationLaw {
 };
 
 dgConservationLaw *dgNewConservationLawAdvection(const std::string vname);
+dgConservationLaw *dgNewConservationLawShallowWater2d();
+dgConservationLaw *dgNewConservationLawWaveEquation();
+dgBoundaryCondition *dgNewBoundaryConditionWaveEquationWall();
 
 #endif

Solver/dgConservationLawShallowWater2d.cpp

+#include "dgConservationLaw.h"
+#include "function.h"
+static double g = 9.81;
+class dgConservationLawShallowWater2d : public dgConservationLaw {
+  class advection : public dataCacheDouble {
+    dataCacheDouble /
+    public:
+    advection(dataCacheMap &cacheMap):
+      sol(cacheMap.get("Solution",this))
+      {};
+    void _eval () { 
+      int nQP = sol().size1();
+      if(_value.size1() != nQP)
+        _value=fullMatrix<double>(nQP,9);
+      for(int i=0; i< nQP; i++) {
+        double H = sol(i,0);
+        double U = sol(i,1);
+        double V = sol(i,2);
+        // flux_x
+        _value(i,0) = U;
+        _value(i,1) = U*U/H + g*H*H/2;
+        _value(i,2) = U*V/H;
+        // flux_y
+        _value(i,3) = V;
+        _value(i,4) = V*U/H;
+        _value(i,5) = V*V/H + g*H*H/2;
+        // flux_z
+        _value(i,6) = 0;
+        _value(i,7) = 0;
+        _value(i,8) = 0;
+      }
+      _value.scale(0);
+    }
+  };
+  class source : public dataCacheDouble {
+    dataCacheDouble &xyz, &solution;
+    public :
+    source(dataCacheMap &cacheMap) : 
+      xyz(cacheMap.get("XYZ",this)),
+      solution(cacheMap.get("Solution",this))
+    {}
+    void _eval () {
+      int nQP = xyz().size1();
+      if(_value.size1() != nQP)
+        _value = fullMatrix<double>(nQP,3);
+      for (int i=0; i<nQP; i++) {
+        double H = solution(i,0);
+        double U = solution(i,1);
+        double V = solution(i,2);
+        double wind = 0.1*sin(xyz(i,1)/1e6*M_PI)/H/1000;
+        double f = 1e-4+xyz(i,1)*2e-11;
+        double gamma = 1e-6;
+        f=0;
+        _value (i,0) = 0;
+        _value (i,1) = -gamma*U - f*V + wind;
+        _value (i,2) = -gamma*V + f*U;
+      }
+    }
+  };
+  class riemann : public dataCacheDouble {
+    dataCacheDouble &normals, &solL, &solR;
+    public:
+    riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
+      normals(cacheMapLeft.get("Normals", this)),
+      solL(cacheMapLeft.get("Solution", this)),
+      solR(cacheMapRight.get("Solution", this))
+      {};
+    void _eval () { 
+      int nQP = solL().size1();
+      if(_value.size1() != nQP)
+        _value = fullMatrix<double>(nQP,6);
+      for(int i=0; i< nQP; i++) {
+        double nx = normals(0,i);
+        double ny = normals(1,i);
+        double UnL = nx*solL(i,1) + ny*solL(i,2);
+        double UnR = nx*solR(i,1) + ny*solR(i,2);
+        double UtL = ny*solL(i,1) - nx*solL(i,2);
+        double UtR = ny*solR(i,1) - nx*solR(i,2);
+        double HR = solR(i,0);
+        double HL = solL(i,0);
+        double c = sqrt(g*(HR+HL)/2);
+        double FUn = -(UnL*UnL/HL + g*HL*HL/2 + UnR*UnR/HR + g*HR*HR/2)/2 + c*(UnL-UnR)/2;
+        double FUt = -(UnL*UtL/HL + UnR*UtR/HR)/2 + c*(UtL-UtR)/2;
+        double FH = -(UnL+UnR)/2 + c*(HL-HR)/2;
+        _value(i,0) = FH;
+        _value(i,1) = FUn*nx + FUt*ny;
+        _value(i,2) = FUn*ny - FUt*nx;
+
+      _value(i,1)=0;
+      _value(i,2)=0;
+        _value(i,3) = -_value(i,0);
+        _value(i,4) = -_value(i,1);
+        _value(i,5) = -_value(i,2);
+      }
+    }
+  };
+  public:
+  dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const {
+    return new advection(cacheMap);
+  }
+  dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {
+    return new riemann(cacheMapLeft, cacheMapRight);
+  }
+  dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const {
+    return 0;
+  }
+  dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const {
+    return new source(cacheMap);
+  }
+  dgConservationLawShallowWater2d() 
+  {
+    _nbf = 3; // H U(=Hu) V(=Hv)
+  }
+};
+
+dgConservationLaw *dgNewConservationLawShallowWater2d() {
+  return new dgConservationLawShallowWater2d();
+}

Solver/dgConservationLawWaveEquation.cpp

+#include "dgConservationLaw.h"
+#include "function.h"
+// dp/dt - rho*c^2  div(u,v) = 0
+// du/dt + 1/rho  dp/dx = 0
+// dv/dt + 1/rho  dp/dy = 0
+static double c=1;
+static double rho=1;
+class dgConservationLawWaveEquation : public dgConservationLaw {
+  class hyperbolicFlux : public dataCacheDouble {
+    dataCacheDouble /
+    public:
+    hyperbolicFlux(dataCacheMap &cacheMap):
+      sol(cacheMap.get("Solution",this))
+      {};
+    void _eval () { 
+      int nQP = sol().size1();
+      if(_value.size1() != nQP)
+        _value=fullMatrix<double>(nQP,9);
+      for(int i=0; i< nQP; i++) {
+        double p = sol(i,0);
+        double u = sol(i,1);
+        double v = sol(i,2);
+        // flux_x
+        _value(i,0) = c*c*rho*u;
+        _value(i,1) = p/rho; 
+        _value(i,2) = 0;
+        // flux_y
+        _value(i,3) = c*c*rho*v;
+        _value(i,4) = 0;
+        _value(i,5) = p/rho;
+        // flux_z
+        _value(i,6) = 0;
+        _value(i,7) = 0;
+        _value(i,8) = 0;
+      }
+    }
+  };
+  class riemann : public dataCacheDouble {
+    dataCacheDouble &normals, &solL, &solR;
+    public:
+    riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
+      normals(cacheMapLeft.get("Normals", this)),
+      solL(cacheMapLeft.get("Solution", this)),
+      solR(cacheMapRight.get("Solution", this))
+      {};
+    void _eval () { 
+      int nQP = solL().size1();
+      if(_value.size1() != nQP)
+        _value = fullMatrix<double>(nQP,6);
+      for(int i=0; i< nQP; i++) {
+        double nx = normals(0,i);
+        double ny = normals(1,i);
+        double unL = nx*solL(i,1) + ny*solL(i,2);
+        double unR = nx*solR(i,1) + ny*solR(i,2);
+        double pR = solR(i,0);
+        double pL = solL(i,0);
+
+        double pRiemann =  0.5 * (pL+pR - (unR-unL)*(rho*c));
+        double unRiemann =  0.5 * (unL+unR - (pR-pL)/(rho*c));
+
+        double Fp = -rho*c*c*unRiemann;
+        double Fun = -pRiemann/rho;
+
+        _value(i,0) = Fp;
+        _value(i,1) = Fun*nx;
+        _value(i,2) = Fun*ny;
+
+        _value(i,3) = -_value(i,0);
+        _value(i,4) = -_value(i,1);
+        _value(i,5) = -_value(i,2);
+      }
+    }
+  };
+  public:
+  dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const {
+    return new hyperbolicFlux(cacheMap);
+  }
+  dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {
+    return new riemann(cacheMapLeft, cacheMapRight);
+  }
+  dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const {
+    return 0;
+  }
+  dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const {
+    return 0;
+  }
+  dgConservationLawWaveEquation() 
+  {
+    _nbf = 3; // H U(=Hu) V(=Hv)
+  }
+};
+
+class dgBoundaryConditionWaveEquationWall : public dgBoundaryCondition {
+  class term : public dataCacheDouble {
+    dataCacheDouble &sol,&normals;
+    public:
+    term(dataCacheMap &cacheMap):
+    sol(cacheMap.get("Solution",this)),
+    normals(cacheMap.get("Normals",this)){}
+    void _eval () { 
+      int nQP = sol().size1();
+      if(_value.size1() != nQP)
+        _value = fullMatrix<double>(nQP,3);
+      for(int i=0; i< nQP; i++) {
+        double nx = normals(0,i);
+        double ny = normals(1,i);
+        double p = sol(i,0);
+
+        _value(i,0) = 0;
+        _value(i,1) = -p/rho*nx;
+        _value(i,2) = -p/rho*ny;
+      }
+    }
+  };
+  public:
+  dgBoundaryConditionWaveEquationWall(){}
+  dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
+    return new term(cacheMapLeft);
+  }
+};
+
+dgConservationLaw *dgNewConservationLawWaveEquation() {
+  return new dgConservationLawWaveEquation();
+}
+dgBoundaryCondition *dgNewBoundaryConditionWaveEquationWall() {
+  return new dgBoundaryConditionWaveEquationWall();
+}

Solver/dgMainShallowWater2d.cpp

+#include <stdio.h>
+#include <vector>
+#include "GModel.h"
+#include "dgGroupOfElements.h"
+#include "dgAlgorithm.h"
+#include "dgConservationLaw.h"
+#include "Gmsh.h"
+#include "function.h"
+
+
+#include "MElement.h"
+void print (const char *filename,const dgGroupOfElements &els, double *v,int iField=1, int nbFields=1);
+
+class dgConservationLawInitialCondition : public dgConservationLaw {
+  class gaussian : public dataCacheDouble {
+    dataCacheDouble &xyz;
+    double _xc,_yc;
+    public:
+    gaussian(dataCacheMap &cacheMap,double xc, double yc):xyz(cacheMap.get("XYZ",this)),_xc(xc),_yc(yc){};
+    void _eval () { 
+      if(_value.size1() != xyz().size1())
+        _value=fullMatrix<double>(xyz().size1(),1);
+      for(int i=0; i< _value.size1(); i++) {
+        _value(i,0)=exp(-(pow(xyz(i,0)-_xc,2)+pow(xyz(i,1)-_yc,2))*100);
+      }
+    }
+  };
+  public:
+  dgConservationLawInitialCondition() {
+    _nbf = 1;
+  }
+  dataCacheDouble *newSourceTerm(dataCacheMap &cacheMap)const {
+    return new gaussian(cacheMap,0.2,0.3);
+  }
+};
+
+/*int main(int argc, char **argv) {
+  GmshMergeFile("input/square1.msh");
+  //we probably need a class to group those three ones
+  std::vector<dgGroupOfElements*> elementGroups;
+  std::vector<dgGroupOfFaces*> faceGroups;
+  std::vector<dgGroupOfFaces*> boundaryGroups;
+  int order=1;
+  int dimension=2;
+  dgAlgorithm algo;
+  function::registerDefaultFunctions();
+  algo.buildGroups(GModel::current(),dimension,order,elementGroups,faceGroups,boundaryGroups);
+
+  //for now, we suppose there is only one group of elements
+  int nbNodes = elementGroups[0]->getNbNodes();
+  fullMatrix<double> sol(nbNodes,elementGroups[0]->getNbElements());
+  fullMatrix<double> residu(nbNodes,elementGroups[0]->getNbElements());
+  // initial condition
+  {
+    dgConservationLawInitialCondition initLaw;
+    algo.residualVolume(initLaw,*elementGroups[0],sol,residu);
+    algo.multAddInverseMassMatrix(*elementGroups[0],residu,sol);
+  }
+
+  fullMatrix<double> advectionSpeed(3,1);
+  advectionSpeed(0,0)=0.15;
+  advectionSpeed(1,0)=0.05;
+  advectionSpeed(2,0)=0.;
+  function::add("advectionSpeed",function::newFunctionConstant(advectionSpeed));
+
+  dgConservationLaw *law = dgNewConservationLawAdvection("advectionSpeed");
+
+  fullMatrix<double> outsideValue(1,1);
+  outsideValue(0,0)=0;
+  function::add("outsideValue",function::newFunctionConstant(outsideValue));
+  law->addBoundaryCondition("Left",dgBoundaryCondition::newOutsideValueCondition(*law,"outsideValue"));
+  law->addBoundaryCondition("Right",dgBoundaryCondition::newOutsideValueCondition(*law,"outsideValue"));
+
+  law->addBoundaryCondition("Top",dgBoundaryCondition::new0FluxCondition(*law));
+  law->addBoundaryCondition("Bottom",dgBoundaryCondition::new0FluxCondition(*law));
+
+  
+  
+ print("output/init.pos",*elementGroups[0],&sol(0,0));
+  for(int iT=0; iT<100; iT++) {
+    algo.rungeKutta(*law,elementGroups,faceGroups,boundaryGroups,0.01,residu,sol);
+    if(iT%10==0){
+      char name[100];
+      sprintf(name,"output/test_%05i.pos",iT/10);
+      printf("%i\n",iT);
+      print(name,*elementGroups[0],&sol(0,0));
+    }
+  }
+  delete law;
+}*/
+class dgConservationLawStommel : public dgConservationLaw {
+  class initbath : public dataCacheDouble {
+    dataCacheDouble &uvw;
+    public:
+    initbath(dataCacheMap &cacheMap):
+      uvw(cacheMap.get("UVW"))
+    {}
+    void _eval () { 
+      if(_value.size1() != uvw().size1())
+        _value=fullMatrix<double>(uvw().size1(),3);
+      for(int i=0; i< _value.size1(); i++) {
+        _value(i,0)=1000;
+        _value(i,1)=0;
+        _value(i,2)=0;
+      }
+    }
+  };
+  public:
+  dgConservationLawStommel() {
+    _nbf = 3;
+  }
+  dataCacheDouble *newSourceTerm(dataCacheMap &cacheMap)const {
+    return new initbath(cacheMap);
+  }
+};
+int main(int argc, char **argv) {
+  GmshMergeFile("input/stommel.msh");
+  //we probably need a class to group those three ones
+  std::vector<dgGroupOfElements*> elementGroups;
+  std::vector<dgGroupOfFaces*> faceGroups;
+  std::vector<dgGroupOfFaces*> boundaryGroups;
+  int order=1;
+  int dimension=2;
+  dgAlgorithm algo;
+  function::registerDefaultFunctions();
+  algo.buildGroups(GModel::current(),dimension,order,elementGroups,faceGroups,boundaryGroups);
+
+  //for now, we suppose there is only one group of elements
+  int nbNodes = elementGroups[0]->getNbNodes();
+  fullMatrix<double> sol(nbNodes,elementGroups[0]->getNbElements()*3);
+  fullMatrix<double> residu(nbNodes,elementGroups[0]->getNbElements()*3);
+  // initial condition
+  {
+    dgConservationLawStommel initLaw;
+    algo.residualVolume(initLaw,*elementGroups[0],sol,residu);
+    algo.multAddInverseMassMatrix(*elementGroups[0],residu,sol);
+  }
+
+  dgConservationLaw *law = dgNewConservationLawShallowWater2d();
+
+  law->addBoundaryCondition("Top",dgBoundaryCondition::new0FluxCondition(*law));
+  law->addBoundaryCondition("Bottom",dgBoundaryCondition::new0FluxCondition(*law));
+  law->addBoundaryCondition("Left",dgBoundaryCondition::new0FluxCondition(*law));
+  law->addBoundaryCondition("Right",dgBoundaryCondition::new0FluxCondition(*law));
+  
+  for(int iT=0; iT<100; iT++) {
+    if(iT%10==0){
+      printf("%i\n",iT);
+      char name[100];
+      sprintf(name,"output/H_%05i.pos",iT/10);
+      print(name,*elementGroups[0],&sol(0,0),0,3);
+      sprintf(name,"output/U_%05i.pos",iT/10);
+      print(name,*elementGroups[0],&sol(0,0),1,3);
+      sprintf(name,"output/V_%05i.pos",iT/10);
+      print(name,*elementGroups[0],&sol(0,0),2,3);
+    }
+    algo.rungeKutta(*law,elementGroups,faceGroups,boundaryGroups,1e-8,residu,sol);
+  }
+  delete law;
+}
+
+void print (const char *filename,const dgGroupOfElements &els, double *v,int iField, int nbFields) {
+  FILE *file = fopen(filename,"w");
+  fprintf(file,"View \"%s\" {\n", filename);
+  for(int iel=0;iel<els.getNbElements();iel++){
+    MElement *el = els.getElement(iel);
+    fprintf(file,"ST (");
+    for (int iv=0; iv<el->getNumVertices(); iv++) {
+      MVertex *vertex = el->getVertex(iv);
+      SPoint3 coord = vertex->point();
+      fprintf(file,"%e, %e, %e%c ",coord.x(),coord.y(),coord.z(),iv==el->getNumVertices()-1?')':',');
+    }
+    fprintf(file,"{");
+    v+=(iField)*el->getNumVertices();
+    for (int iv=0; iv<el->getNumVertices(); iv++){
+      fprintf(file,"%e%c ",*(v++),iv==el->getNumVertices()-1?'}':',');
+    }
+    v+=(nbFields-1-iField)*el->getNumVertices();
+    fprintf(file,";\n");
+  }
+  fprintf(file,"};");
+  fclose(file);
+}

Solver/dgMainWaveEquation.cpp

+#include <stdio.h>
+#include <vector>
+#include "GModel.h"
+#include "dgGroupOfElements.h"
+#include "dgAlgorithm.h"
+#include "dgConservationLaw.h"
+#include "Gmsh.h"
+#include "function.h"
+
+
+#include "MElement.h"
+void print (const char *filename,const dgGroupOfElements &els, double *v,int iField=1, int nbFields=1);
+void print_vector (const char *filename,const dgGroupOfElements &els, double *v,int iField, int nbFields);
+
+class dgConservationLawInitialCondition : public dgConservationLaw {
+  class gaussian : public dataCacheDouble {
+    dataCacheDouble &xyz;
+    double _xc,_yc;
+    public:
+    gaussian(dataCacheMap &cacheMap,double xc, double yc):xyz(cacheMap.get("XYZ",this)),_xc(xc),_yc(yc){};
+    void _eval () { 
+      if(_value.size1() != xyz().size1())
+        _value=fullMatrix<double>(xyz().size1(),3);
+      for(int i=0; i< _value.size1(); i++) {
+        _value(i,0)=exp(-(pow(xyz(i,0)-_xc,2)+pow(xyz(i,1)-_yc,2))*100)+1;
+      }
+    }
+  };
+  public:
+  dgConservationLawInitialCondition() {
+    _nbf = 3;
+  }
+  dataCacheDouble *newSourceTerm(dataCacheMap &cacheMap)const {
+    //return new gaussian(cacheMap,0.2,0.3);
+    return new gaussian(cacheMap,0.5,0.5);
+  }
+};
+
+int main(int argc, char **argv){
+  GmshMergeFile("input/square1.msh");
+  //we probably need a class to group those three ones
+  std::vector<dgGroupOfElements*> elementGroups;
+  std::vector<dgGroupOfFaces*> faceGroups;
+  std::vector<dgGroupOfFaces*> boundaryGroups;
+  int order=1;
+  int dimension=2;
+  dgAlgorithm algo;
+  function::registerDefaultFunctions();
+  algo.buildGroups(GModel::current(),dimension,order,elementGroups,faceGroups,boundaryGroups);
+
+  //for now, we suppose there is only one group of elements
+  int nbNodes = elementGroups[0]->getNbNodes();
+  dgConservationLaw *law = dgNewConservationLawWaveEquation();
+  fullMatrix<double> sol(nbNodes,elementGroups[0]->getNbElements()*law->nbFields());
+  fullMatrix<double> residu(nbNodes,elementGroups[0]->getNbElements()*law->nbFields());
+  // 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++) {
+    algo.rungeKutta(*law,elementGroups,faceGroups,boundaryGroups,1e-3,residu,sol);
+    if(iT%10==0){
+      char name[100];
+      printf("%i\n",iT);
+      sprintf(name,"output/p_%05i.pos",iT);
+      print(name,*elementGroups[0],&sol(0,0),0,3);
+      sprintf(name,"output/uv_%05i.pos",iT);
+      print_vector(name,*elementGroups[0],&sol(0,0),1,3);
+      sprintf(name,"output/v_%05i.pos",iT);
+      print(name,*elementGroups[0],&sol(0,0),2,3);
+      sprintf(name,"output/u_%05i.pos",iT);
+      print(name,*elementGroups[0],&sol(0,0),1,3);
+    }
+  }
+  delete law;
+}
+
+void print (const char *filename,const dgGroupOfElements &els, double *v,int iField, int nbFields) {
+  FILE *file = fopen(filename,"w");
+  fprintf(file,"View \"%s\" {\n", filename);
+  for(int iel=0;iel<els.getNbElements();iel++){
+    MElement *el = els.getElement(iel);
+    fprintf(file,"ST (");
+    for (int iv=0; iv<el->getNumVertices(); iv++) {
+      MVertex *vertex = el->getVertex(iv);
+      SPoint3 coord = vertex->point();
+      fprintf(file,"%e, %e, %e%c ",coord.x(),coord.y(),coord.z(),iv==el->getNumVertices()-1?')':',');
+    }
+    fprintf(file,"{");
+    v+=(iField)*el->getNumVertices();
+    for (int iv=0; iv<el->getNumVertices(); iv++){
+      fprintf(file,"%e%c ",*(v++),iv==el->getNumVertices()-1?'}':',');
+    }
+    v+=(nbFields-1-iField)*el->getNumVertices();
+    fprintf(file,";\n");
+  }
+  fprintf(file,"};");
+  fclose(file);
+}
+void print_vector (const char *filename,const dgGroupOfElements &els, double *v,int iField, int nbFields) {
+  FILE *file = fopen(filename,"w");
+  fprintf(file,"View \"%s\" {\n", filename);
+  for(int iel=0;iel<els.getNbElements();iel++){
+    MElement *el = els.getElement(iel);
+    fprintf(file,"VT (");
+    for (int iv=0; iv<el->getNumVertices(); iv++) {
+      MVertex *vertex = el->getVertex(iv);
+      SPoint3 coord = vertex->point();
+      fprintf(file,"%e, %e, %e%c ",coord.x(),coord.y(),coord.z(),iv==el->getNumVertices()-1?')':',');
+    }
+    v+=(iField)*el->getNumVertices();
+    fprintf(file,"{%e, %e, 0, %e, %e, 0, %e, %e, 0};\n", v[0],v[3],v[1],v[4],v[2],v[5]);
+    v+=(nbFields-iField)*el->getNumVertices();
+  }
+  fprintf(file,"};");
+  fclose(file);
+}