From: michael <michael@localhost.localdomain>
Date: Thu, 17 Dec 2020 10:52:44 +0000 (+0100)
Subject: Renamed test source file
X-Git-Tag: V9_7_0~81
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=ee4ef4fc3319f5e0c202e935468109ff6ddcb643;p=tools%2Fsolverlab.git

Renamed test source file
---

diff --git a/CoreFlows/examples/C/CMakeLists.txt b/CoreFlows/examples/C/CMakeLists.txt
index b966a6c..587e1d3 100755
--- a/CoreFlows/examples/C/CMakeLists.txt
+++ b/CoreFlows/examples/C/CMakeLists.txt
@@ -98,14 +98,14 @@ CreateTestExecAndInstall(StationaryDiffusionEquation_3DFV_StructuredTetrahedra.c
 CreateTestExecAndInstall(testEOS.cxx  "${libs_for_tests}" )
 
  
-  add_executable(WaveSystem_FV_SphericalExplosion_CDMATH.exe WaveSystem_2DFV_SphericalExplosion_CDMATH.cxx)
+  add_executable(WaveSystem_FV_SphericalExplosion_CDMATH.exe WaveSystem_FV_SphericalExplosion_CDMATH.cxx)
   target_link_libraries(WaveSystem_FV_SphericalExplosion_CDMATH.exe CoreFlowsLibs )
   install(TARGETS WaveSystem_FV_SphericalExplosion_CDMATH.exe DESTINATION share/examples)
   add_test(NAME WaveSystem_2DFV_SphericalExplosion_CDMATH_SQUARE  COMMAND "./WaveSystem_FV_SphericalExplosion_CDMATH.exe")     
   add_test(NAME WaveSystem_3DFV_SphericalExplosion_CDMATH_CUBE    COMMAND "./WaveSystem_FV_SphericalExplosion_CDMATH.exe" resources/meshCube.med)     
 
 if( SOLVERLAB_WITH_MPI )
-  add_executable(WaveSystem_2DFV_SphericalExplosion_MPI.exe WaveSystem_2DFV_SphericalExplosion_MPI.cxx)                    # compilation of the testxxx.exe 
+  add_executable(WaveSystem_2DFV_SphericalExplosion_MPI.exe WaveSystem_FV_SphericalExplosion_MPI.cxx)                    # compilation of the testxxx.exe 
   target_link_libraries(WaveSystem_2DFV_SphericalExplosion_MPI.exe CoreFlowsLibs  ${MPI_LIBRARY})              # provide required lib for testxxx.exe 
   install(TARGETS WaveSystem_2DFV_SphericalExplosion_MPI.exe DESTINATION share/examples)
   add_test(NAME WaveSystem_2DFV_SphericalExplosion_MPI_SEQ_SQUARE  COMMAND "${MPIEXEC}" "-n" "1" "./WaveSystem_2DFV_SphericalExplosion_MPI.exe")     
diff --git a/CoreFlows/examples/C/WaveSystem_2DFV_SphericalExplosion_CDMATH.cxx b/CoreFlows/examples/C/WaveSystem_2DFV_SphericalExplosion_CDMATH.cxx
deleted file mode 100755
index 0299ef2..0000000
--- a/CoreFlows/examples/C/WaveSystem_2DFV_SphericalExplosion_CDMATH.cxx
+++ /dev/null
@@ -1,298 +0,0 @@
-//============================================================================
-// Author      : Michael NDJINGA
-// Date        : November 2020
-// Description : multiD linear wave system
-//============================================================================
-
-#include <iostream>
-#include <string>
-#include <cmath>
-
-#include "Mesh.hxx"
-#include "Cell.hxx"
-#include "Face.hxx"
-#include "Field.hxx"
-#include "SparseMatrixPetsc.hxx"
-#include "CdmathException.hxx"
-
-using namespace std;
-
-double p0  =155e5;   //reference pressure in a pressurised nuclear vessel
-double c0  =700.;    //reference sound speed for water at 155 bars, 600K
-double rho0=p0/c0*c0;//reference density
-double precision=1e-5;
-
-void initial_conditions_shock(Mesh my_mesh,Field& pressure_field,Field& velocity_field)
-{
-    double rayon=0.35;
-    double xcentre=0.;
-    double ycentre=0;
-    double zcentre=0;
-	
-	double x, y, z;
-	double val, valX, valY, valZ;
-	
-    int dim    =my_mesh.getMeshDimension();
-    int nbCells=my_mesh.getNumberOfCells();
-
-    for (int j=0 ; j<nbCells ; j++)
-    {
-        x = my_mesh.getCell(j).x() ;
-		if(dim>1)
-		{
-			y = my_mesh.getCell(j).y() ;
-			if(dim==3)
-				z = my_mesh.getCell(j).z() ;
-		}
-
-        valX=(x-xcentre)*(x-xcentre);
-		if(dim==1)
-			val=sqrt(valX);
-		else if(dim==2)
-		{
-			valY=(y-ycentre)*(y-ycentre);
-			val=sqrt(valX+valY);		
-		}
-		else if(dim==3)
-		{
-			valY=(y-ycentre)*(y-ycentre);
-			valZ=(z-zcentre)*(z-zcentre);
-			val=sqrt(valX+valY+valZ);		
-		}
-		
-		for(int idim=0; idim<dim; idim++)
-			velocity_field[j,idim]=0;
-			
-        if (val<rayon)
-            pressure_field(j) = 155e5;
-        else
-            pressure_field(j) = 70e5;
-    }
-}
-
-Matrix jacobianMatrices(Vector normal, double coeff)
-{
-    int dim=normal.size();
-    Matrix    A(dim+1,dim+1);
-    Matrix absA(dim+1,dim+1);
-
-    absA(0,0)=c0*coeff;
-    for(int i=0 ; i<dim; i++)
-    {
-        A(i+1,0)=      normal[i]*coeff;
-        A(0,i+1)=c0*c0*normal[i]*coeff;
-        for( int j=0 ; j<dim; j++)
-            absA(i+1,j+1)=c0*normal[i]*normal[j]*coeff;
-    }
-    return (A - absA)*(1./2);
-}
-    
-SparseMatrixPetsc computeDivergenceMatrix(Mesh my_mesh, int nbVoisinsMax, double dt)
-{
-    int nbCells = my_mesh.getNumberOfCells();
-    int dim=my_mesh.getMeshDimension();
-    int nbComp=dim+1;
-    Vector normal(dim);
-
-    SparseMatrixPetsc implMat(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp);
-
-    Matrix idMoinsJacCL(nbComp);
-    
-    for(int j=0; j<nbCells; j++)//On parcourt les cellules
-    {
-        Cell Cj = my_mesh.getCell(j);
-        int nbFaces = Cj.getNumberOfFaces();
-
-        for(int k=0; k<nbFaces; k++)
-        {
-            int indexFace = Cj.getFacesId()[k];
-            Face Fk = my_mesh.getFace(indexFace);
-            for(int i =0; i<dim ; i++)
-                normal[i] = Cj.getNormalVector(k, i);//normale sortante
-
-            Matrix Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure());
-
-            int cellAutre =-1;
-            if ( not Fk.isBorder())
-            {
-                /* hypothese: La cellule d'index indexC1 est la cellule courante index j */
-                if (Fk.getCellsId()[0] == j) 
-                    // hypothese verifiée 
-                    cellAutre = Fk.getCellsId()[1];
-                else if(Fk.getCellsId()[1] == j) 
-                    // hypothese non verifiée 
-                    cellAutre = Fk.getCellsId()[0];
-                else
-                    throw CdmathException("computeDivergenceMatrix: problem with mesh, unknown cell number");
-                    
-                implMat.addValue(j*nbComp,cellAutre*nbComp,Am);
-                implMat.addValue(j*nbComp,        j*nbComp,Am*(-1.));
-            }
-            else 
-            {    
-				if( Fk.getGroupName() != "Periodic" && Fk.getGroupName() != "Neumann")//Wall boundary condition unless Periodic/Neumann specified explicitly
-                {
-					Vector v(dim+1);
-                    for(int i=0; i<dim; i++)
-                        v[i+1]=normal[i];
-                    Matrix idMoinsJacCL=v.tensProduct(v)*2;
-                    
-                    implMat.addValue(j*nbComp,j*nbComp,Am*(-1.)*idMoinsJacCL);
-				}
-                else if( Fk.getGroupName() == "Periodic")//Periodic boundary condition
-                {
-                    int indexFP=my_mesh.getIndexFacePeriodic(indexFace);
-                    Face Fp = my_mesh.getFace(indexFP);
-                    cellAutre = Fp.getCellsId()[0];
-                    
-                    implMat.addValue(j*nbComp,cellAutre*nbComp,Am);
-                    implMat.addValue(j*nbComp,        j*nbComp,Am*(-1.));
-				}
-                else if(Fk.getGroupName() != "Neumann")//Nothing to do for Neumann boundary condition
-                {
-                    cout<< Fk.getGroupName() <<endl;
-                    throw CdmathException("computeDivergenceMatrix: Unknown boundary condition name");
-				}
-            }
-        }   
-    }     
-    return implMat;
-}
-
-void WaveSystem(double tmax, int ntmax, double cfl, int output_freq, const Mesh& my_mesh, const string file)
-{
-	/* Retrieve mesh data */
-    int dim=my_mesh.getMeshDimension();
-    int nbCells = my_mesh.getNumberOfCells();
-    std::string meshName=my_mesh.getName();
-    int nbVoisinsMax=my_mesh.getMaxNbNeighbours(CELLS);
-    double dx_min=my_mesh.minRatioVolSurf();
-
-
-    /* Initial conditions */
-    cout << "Construction of the initial condition" << endl;
-    
-    Field pressure_field("Pressure",CELLS,my_mesh,1) ;
-    Field velocity_field("Velocity",CELLS,my_mesh,dim) ;
-    initial_conditions_shock(my_mesh,pressure_field, velocity_field);
-
-    /* iteration vectors */
-    Vector  Un(nbCells*(dim+1));
-    Vector dUn(nbCells*(dim+1));
-    
-    for(int k =0; k<nbCells; k++){
-        Un[k*(dim+1)+0] =     pressure_field[k];
-        Un[k*(dim+1)+1] =rho0*velocity_field[k,0];
-        if(dim>=2)
-            Un[k + 2*nbCells] = rho0*velocity_field[k,1] ;
-        if(dim==3)
-            Un[k + 3*nbCells] = rho0*velocity_field[k,2];
-		}
-
-    /*
-     * MED output of the initial condition at t=0 and iter = 0
-     */
-    int it=0;
-    bool isStationary=false;
-    double time=0.;
-    double dt = cfl * dx_min / c0;
-    
-    cout << "Saving the solution at T=" << time << endl;
-    pressure_field.setTime(time,it);
-    pressure_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_pressure");
-    velocity_field.setTime(time,it);
-    velocity_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_velocity");
-    /* --------------------------------------------- */
-
-    SparseMatrixPetsc divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt);
-
-    /* Time loop */
-    cout<< "Starting computation of the linear wave system with an explicit UPWIND scheme" << endl;
-    while (it<ntmax && time <= tmax && ! isStationary)
-    {
-        dUn=divMat*Un;
-        Un-=dUn;
-        
-        time=time+dt;
-        it=it+1;
- 
-         isStationary = dUn.norm()<precision;
-        /* Sauvegardes */
-        if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax)
-        {
-            cout<<"-- Iteration: " << it << ", Time: " << time << ", dt: " << dt<<endl;
-
-            for(int k=0; k<nbCells; k++)
-            {
-                pressure_field[k]  =Un[k*(dim+1)+0];
-                velocity_field[k,0]=Un[k*(dim+1)+1]/rho0;
-                if(dim>1)
-                {
-                    velocity_field[k,1]=Un[k*(dim+1)+2]/rho0;
-                    if(dim>2)
-                        velocity_field[k,2]=Un[k*(dim+1)+3]/rho0;
-				}
-			}
-            pressure_field.setTime(time,it);
-            pressure_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_pressure",false);
-            velocity_field.setTime(time,it);
-            velocity_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_velocity",false);
-		}
-    }
-    cout<<"End of calculation -- Iteration: " << it << ", Time: "<< time<< ", dt: " << dt<<endl;
-
-    if(it>=ntmax)
-        cout<< "Nombre de pas de temps maximum ntmax= "<< ntmax<< " atteint"<<endl;
-    else if(isStationary)
-        cout<< "Régime stationnaire atteint au pas de temps "<< it<< ", t= "<< time<<endl;       
-    else
-        cout<< "Temps maximum Tmax= "<< tmax<< " atteint"<<endl;
-}
- 
-int main(int argc, char *argv[])
-{
-	cout << "-- Starting the RESOLUTION OF THE 2D WAVE SYSTEM"<<endl;
-	cout << "- Numerical scheme : Upwind explicit scheme" << endl;
-	cout << "- Boundary conditions : WALL" << endl;
-
-    // Problem data
-    double cfl=0.49;
-    double tmax=1.;
-    double ntmax=3;//20000;
-    int freqSortie=100;
-    string fileOutPut="SphericalWave";
-
-	Mesh myMesh;
-	
-	if(argc<2)
-	{
-		    cout << "- DOMAIN : SQUARE" << endl;
-		    cout << "- MESH : CARTESIAN, GENERATED INTERNALLY WITH CDMATH" << endl<< endl;
-		    cout << "Construction of a cartesian mesh" << endl;
-	    double xinf=-0.5;
-	    double xsup= 0.5;
-	    double yinf=-0.5;
-	    double ysup= 0.5;
-	    int nx=50;
-	    int ny=50;
-	    myMesh=Mesh(xinf,xsup,nx,yinf,ysup,ny);
-	    double eps=precision;
-	    myMesh.setGroupAtPlan(xsup,0,eps,"RightEdge");
-	    myMesh.setGroupAtPlan(xinf,0,eps,"LeftEdge");
-	    myMesh.setGroupAtPlan(yinf,1,eps,"BottomEdge");
-	    myMesh.setGroupAtPlan(ysup,1,eps,"TopEdge");
-	}
-	else
-	{
-	    cout << "- MESH:  GENERATED EXTERNALLY WITH SALOME" << endl;
-	    cout << "Loading of a mesh named "<<argv[1] << endl;
-	    string filename = argv[1];
-	    myMesh=Mesh(filename);
-	}
-
-	WaveSystem(tmax,ntmax,cfl,freqSortie,myMesh,fileOutPut);
-	
-    cout << "Simulation complete." << endl;
-
-    return 0;
-}
diff --git a/CoreFlows/examples/C/WaveSystem_2DFV_SphericalExplosion_MPI.cxx b/CoreFlows/examples/C/WaveSystem_2DFV_SphericalExplosion_MPI.cxx
deleted file mode 100755
index 0356da9..0000000
--- a/CoreFlows/examples/C/WaveSystem_2DFV_SphericalExplosion_MPI.cxx
+++ /dev/null
@@ -1,404 +0,0 @@
-//============================================================================
-// Author      : Michael NDJINGA
-// Date        : November 2020
-// Description : multiD linear wave system run in parallel with calls to MPI
-//               Test used in the MPI version of the Salome platform
-//============================================================================
-
-#include <iostream>
-#include <string>
-#include <cmath>
-
-#include "Mesh.hxx"
-#include "Cell.hxx"
-#include "Face.hxx"
-#include "Field.hxx"
-#include "CdmathException.hxx"
-
-#include <petscksp.h>
-
-using namespace std;
-
-double p0  =155e5;   //reference pressure in a pressurised nuclear vessel
-double c0  =700.;    //reference sound speed for water at 155 bars, 600K
-double rho0=p0/c0*c0;//reference density
-double precision=1e-5;
-
-void initial_conditions_shock(Mesh my_mesh,Field& pressure_field,Field& velocity_field)
-{
-    double rayon=0.35;
-    double xcentre=0.;
-    double ycentre=0;
-    double zcentre=0;
-	
-	double x, y, z;
-	double val, valX, valY, valZ;
-	
-    int dim    =my_mesh.getMeshDimension();
-    int nbCells=my_mesh.getNumberOfCells();
-
-    for (int j=0 ; j<nbCells ; j++)
-    {
-        x = my_mesh.getCell(j).x() ;
-		if(dim>1)
-		{
-			y = my_mesh.getCell(j).y() ;
-			if(dim==3)
-				z = my_mesh.getCell(j).z() ;
-		}
-
-        valX=(x-xcentre)*(x-xcentre);
-		if(dim==1)
-			val=sqrt(valX);
-		else if(dim==2)
-		{
-			valY=(y-ycentre)*(y-ycentre);
-			val=sqrt(valX+valY);		
-		}
-		else if(dim==3)
-		{
-			valY=(y-ycentre)*(y-ycentre);
-			valZ=(z-zcentre)*(z-zcentre);
-			val=sqrt(valX+valY+valZ);		
-		}
-		
-		for(int idim=0; idim<dim; idim++)
-			velocity_field[j,idim]=0;
-			
-        if (val<rayon)
-            pressure_field(j) = 155e5;
-        else
-            pressure_field(j) = 70e5;
-    }
-}
-
-void addValue( int i, int j, Matrix M, Mat * mat )
-{
-    int I,J;
-    for (int k=0; k<M.getNumberOfRows(); k++)
-        for (int l=0; l<M.getNumberOfColumns(); l++)
-        {
-            I=i+k;
-            J=j+l;
-            MatSetValues( *mat,1, &I, 1, &J, &M(k,l), ADD_VALUES);
-        }
-}
-
-Matrix jacobianMatrices(Vector normal, double coeff)
-{
-    int dim=normal.size();
-    Matrix    A(dim+1,dim+1);
-    Matrix absA(dim+1,dim+1);
-
-    absA(0,0)=c0*coeff;
-    for(int i=0 ; i<dim; i++)
-    {
-        A(i+1,0)=      normal[i]*coeff;
-        A(0,i+1)=c0*c0*normal[i]*coeff;
-        for( int j=0 ; j<dim; j++)
-            absA(i+1,j+1)=c0*normal[i]*normal[j]*coeff;
-    }
-    return (A - absA)*(1./2);
-}
-    
-void computeDivergenceMatrix(Mesh my_mesh, Mat * implMat, double dt)
-{
-    int nbCells = my_mesh.getNumberOfCells();
-    int dim=my_mesh.getMeshDimension();
-    int nbComp=dim+1;
-    Vector normal(dim);
-
-    Matrix idMoinsJacCL(nbComp);
-    
-    for(int j=0; j<nbCells; j++)//On parcourt les cellules
-    {
-        Cell Cj = my_mesh.getCell(j);
-        int nbFaces = Cj.getNumberOfFaces();
-
-        for(int k=0; k<nbFaces; k++)
-        {
-            int indexFace = Cj.getFacesId()[k];
-            Face Fk = my_mesh.getFace(indexFace);
-            for(int i =0; i<dim ; i++)
-                normal[i] = Cj.getNormalVector(k, i);//normale sortante
-
-            Matrix Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure());
-
-            int cellAutre =-1;
-            if ( not Fk.isBorder())
-            {
-                /* hypothese: La cellule d'index indexC1 est la cellule courante index j */
-                if (Fk.getCellsId()[0] == j) 
-                    // hypothese verifiée 
-                    cellAutre = Fk.getCellsId()[1];
-                else if(Fk.getCellsId()[1] == j) 
-                    // hypothese non verifiée 
-                    cellAutre = Fk.getCellsId()[0];
-                else
-                    throw CdmathException("computeDivergenceMatrix: problem with mesh, unknown cell number");
-                    
-                addValue(j*nbComp,cellAutre*nbComp,Am      ,implMat);
-                addValue(j*nbComp,        j*nbComp,Am*(-1.),implMat);
-            }
-            else 
-            {    
-				if( Fk.getGroupName() != "Periodic" && Fk.getGroupName() != "Neumann")//Wall boundary condition unless Periodic/Neumann specified explicitly
-                {
-					Vector v(dim+1);
-                    for(int i=0; i<dim; i++)
-                        v[i+1]=normal[i];
-                    Matrix idMoinsJacCL=v.tensProduct(v)*2;
-                    
-                    addValue(j*nbComp,j*nbComp,Am*(-1.)*idMoinsJacCL,implMat);
-				}
-                else if( Fk.getGroupName() == "Periodic")//Periodic boundary condition
-                {
-                    int indexFP=my_mesh.getIndexFacePeriodic(indexFace);
-                    Face Fp = my_mesh.getFace(indexFP);
-                    cellAutre = Fp.getCellsId()[0];
-                    
-                    addValue(j*nbComp,cellAutre*nbComp,Am      ,implMat);
-                    addValue(j*nbComp,        j*nbComp,Am*(-1.),implMat);
-				}
-                else if(Fk.getGroupName() != "Neumann")//Nothing to do for Neumann boundary condition
-                {
-                    cout<< Fk.getGroupName() <<endl;
-                    throw CdmathException("computeDivergenceMatrix: Unknown boundary condition name");
-				}
-            }
-        }   
-    }     
-}
-
-void WaveSystem(double tmax, int ntmax, double cfl, int output_freq, const Mesh& my_mesh, const string file, int rank, int size, string resultDirectory)
-{
-	/* Time iteration variables */
-    int it=0;
-    bool isStationary=false;
-    double time=0.;
-    double dt;
-    double norm;
-	    
-	/* PETSc variables */
-	int globalNbUnknowns;
-	int  localNbUnknowns;
-	int d_nnz, o_nnz;
-	Vec Un, dUn, Un_seq;
-    Mat divMat;
-	VecScatter scat;
-
-	
-	/* Mesh parameters managed only by proc 0 */
-    int nbCells;
-    int dim;
-    int nbComp;	    
-	int idx;//Index where to add the vector values
-	double value;//value to add in the vector	
-    Field pressure_field, velocity_field;
-    std::string meshName;
-
-	if(rank == 0)
-	    globalNbUnknowns=my_mesh.getNumberOfCells()*(my_mesh.getMeshDimension()+1);//nbCells*nbComp
-	
-	MPI_Bcast(&globalNbUnknowns, 1, MPI_INT, 0, MPI_COMM_WORLD);
- 
-    /* iteration vectors */
-	VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE    ,globalNbUnknowns,&Un);
-	VecDuplicate (Un,&dUn);
-	if(rank == 0)
-		VecCreateSeq(PETSC_COMM_SELF,globalNbUnknowns,&Un_seq);//For saving results on proc 0
-
-	VecScatterCreateToZero(Un,&scat,&Un_seq);
-
-	/* System matrix */
-	VecGetLocalSize(Un, &localNbUnknowns);
-
-	if(rank == 0)
-	{
-		int nbVoisinsMax = my_mesh.getMaxNbNeighbours(CELLS);
-	    d_nnz=(nbVoisinsMax+1)*(my_mesh.getMeshDimension()+1);//(nbVoisinsMax+1)*nbComp
-	    o_nnz= nbVoisinsMax   *(my_mesh.getMeshDimension()+1);//                 nbComp
-	}
-	MPI_Bcast(&d_nnz, 1, MPI_INT, 0, MPI_COMM_WORLD);
-	MPI_Bcast(&o_nnz, 1, MPI_INT, 0, MPI_COMM_WORLD);
-
-   	MatCreateAIJ(PETSC_COMM_WORLD,localNbUnknowns,localNbUnknowns,globalNbUnknowns,globalNbUnknowns,d_nnz,NULL,o_nnz,NULL,&divMat);
-
-	if(rank == 0)
-		{
-		/* Retrieve mesh data */
-	    nbCells = my_mesh.getNumberOfCells();
-	    dim=my_mesh.getMeshDimension();
-	    nbComp=dim+1;	    
-	    meshName=my_mesh.getName();
-	    double dx_min=my_mesh.minRatioVolSurf();
-		dt = cfl * dx_min / c0;
-		
-	    /* Initial conditions */
-	    cout<<"Building the initial condition on processor 0" << endl;
-	    
-	    pressure_field=Field("Pressure",CELLS,my_mesh,1) ;
-	    velocity_field=Field("Velocity",CELLS,my_mesh,dim) ;
-	    initial_conditions_shock(my_mesh,pressure_field, velocity_field);
-	
-	    cout << "Saving the solution at T=" << time <<"  on processor 0"<<endl;
-	    pressure_field.setTime(time,it);
-	    pressure_field.writeMED(resultDirectory+"/WaveSystem"+to_string(dim)+"DUpwind_"+to_string(size)+"Procs_"+meshName+"_pressure");
-	    velocity_field.setTime(time,it);
-	    velocity_field.writeMED(resultDirectory+"/WaveSystem"+to_string(dim)+"DUpwind_"+to_string(size)+"Procs_"+meshName+"_velocity");
-	    /* --------------------------------------------- */
-	
-
-	    for(int k =0; k<nbCells; k++)
-	    {
-			idx = k*nbComp;
-			value=pressure_field[k];//vale to add in the vector
-			VecSetValues(Un,1,&idx,&value,INSERT_VALUES);
-			for(int idim =0; idim<dim; idim++)
-			{
-				idx = k*nbComp+1+idim;
-				value =rho0*velocity_field[k,idim];
-				VecSetValues(Un,1,&idx,&value,INSERT_VALUES);
-			}
-		}
-	    computeDivergenceMatrix(my_mesh,&divMat,dt);
-	}		
-
-	VecAssemblyBegin(Un);
-	VecAssemblyEnd(Un);
-	
-	//VecView(Un, PETSC_VIEWER_STDOUT_WORLD );
-
-	MatAssemblyBegin(divMat, MAT_FINAL_ASSEMBLY);
-	MatAssemblyEnd(  divMat, MAT_FINAL_ASSEMBLY);
-
-	//MatView(divMat,	PETSC_VIEWER_STDOUT_WORLD );
-
-	MPI_Bcast(&dt, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
-
-    /* Time loop */
-	PetscPrintf(PETSC_COMM_WORLD,"Starting computation of the linear wave system on all processors : \n\n");
-    while (it<ntmax && time <= tmax && !isStationary)
-    {
-		MatMult(divMat,Un,dUn);
-        VecAXPY(Un,-1,dUn);
-        
-        time=time+dt;
-        it=it+1;
- 
-		VecNorm(dUn,NORM_2,&norm);
-        isStationary = norm<precision;
-        /* Sauvegardes */
-        if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax )
-        {
-			PetscPrintf(PETSC_COMM_WORLD,"-- Iteration: %d, Time: %f, dt: %f, saving results on processor 0 \n", it, time, dt);
-			VecScatterBegin(scat,Un,Un_seq,INSERT_VALUES,SCATTER_FORWARD);
-			VecScatterEnd(  scat,Un,Un_seq,INSERT_VALUES,SCATTER_FORWARD);
-
-			if(rank == 0)
-			{
-	            for(int k=0; k<nbCells; k++)
-	            {
-					idx = k*(dim+1)+0;
-					VecGetValues(Un_seq,1,&idx,&value);
-	                pressure_field[k]  =value;
-					for(int idim =0; idim<dim; idim++)
-					{
-						idx = k*nbComp+1+idim;
-						VecGetValues(Un_seq,1,&idx,&value);
-						velocity_field[k,idim] = value/rho0;
-					}
-				}
-	            pressure_field.setTime(time,it);
-	            pressure_field.writeMED(resultDirectory+"/WaveSystem"+to_string(dim)+"DUpwind_"+to_string(size)+"Procs_"+meshName+"_pressure",false);
-	            velocity_field.setTime(time,it);
-	            velocity_field.writeMED(resultDirectory+"/WaveSystem"+to_string(dim)+"DUpwind_"+to_string(size)+"Procs_"+meshName+"_velocity",false);
-			}
-		}
-    }
-
-	PetscPrintf(PETSC_COMM_WORLD,"\n End of calculations at iteration: %d, and time: %f\n", it, time);
-    if(it>=ntmax)
-        PetscPrintf(PETSC_COMM_WORLD, "Nombre de pas de temps maximum ntmax= %d atteint\n", ntmax);
-    else if(isStationary)
-        PetscPrintf(PETSC_COMM_WORLD, "Régime stationnaire atteint au pas de temps %d, t= %f\n", it, time);       
-    else
-        PetscPrintf(PETSC_COMM_WORLD, "Temps maximum tmax= %f atteint\n", tmax);
-
-	VecDestroy(&Un);
-	VecDestroy(&Un_seq);
-	VecDestroy(&dUn);
-	MatDestroy(&divMat);
-}
- 
-int main(int argc, char *argv[])
-{
-	/* PETSc initialisation */
-	PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
-	PetscMPIInt    size;        /* size of communicator */
-	PetscMPIInt    rank;        /* processor rank */
-	MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
-	MPI_Comm_size(PETSC_COMM_WORLD,&size);
-	
-    // Problem data
-    double cfl=0.49;
-    double tmax=1.;
-    int ntmax=2000;
-    int freqSortie=10;
-    string fileOutPut="SphericalWave";
-    Mesh myMesh;
-	string resultDirectory="./";
-	
-	if(size>1)
-		PetscPrintf(PETSC_COMM_WORLD,"---- More than one processor detected : running a parallel simulation ----\n");
-		PetscPrintf(PETSC_COMM_WORLD,"---- Limited parallelism : input and output remain sequential ----\n");
-		PetscPrintf(PETSC_COMM_WORLD,"---- Only the matrix-vector products are done in parallel ----\n");
-		PetscPrintf(PETSC_COMM_WORLD,"---- Processor 0 is in charge of building the mesh, saving the results, filling and then distributing the matrix to other processors.\n\n");
-		
-	if(rank == 0)
-	{
-	    cout << "-- Starting the RESOLUTION OF THE 2D WAVE SYSTEM on "<< size <<" processors"<<endl;
-	    cout << "- Numerical scheme : Upwind explicit scheme" << endl;
-	    cout << "- Boundary conditions : WALL" << endl;
-	
-		/* Read or create mesh */
-		if(argc<2)
-		{
-		    cout << "- DOMAIN : SQUARE" << endl;
-		    cout << "- MESH : CARTESIAN, GENERATED INTERNALLY WITH CDMATH" << endl<< endl;
-		    cout << "Construction of a cartesian mesh on processor 0" << endl;
-		    double xinf=-0.5;
-		    double xsup= 0.5;
-		    double yinf=-0.5;
-		    double ysup= 0.5;
-		    int nx=50;
-		    int ny=50;
-		    myMesh=Mesh(xinf,xsup,nx,yinf,ysup,ny);
-		    
-		    double eps=precision;
-		    myMesh.setGroupAtPlan(xsup,0,eps,"RightEdge");
-		    myMesh.setGroupAtPlan(xinf,0,eps,"LeftEdge");
-		    myMesh.setGroupAtPlan(yinf,1,eps,"BottomEdge");
-		    myMesh.setGroupAtPlan(ysup,1,eps,"TopEdge");
-		}
-		else
-		{
-		    cout << "- MESH:  GENERATED EXTERNALLY WITH SALOME" << endl<< endl;
-		    cout << "Loading of a mesh named "<<argv[1]<<" on processor 0" << endl;
-		    string filename = argv[1];
-		    myMesh=Mesh(filename);
-		}
-
-		/*Detect directory where to same results*/
-		if(argc>2)
-		    resultDirectory = argv[2];
-	}
-	
-	WaveSystem(tmax,ntmax,cfl,freqSortie,myMesh,fileOutPut, rank, size, resultDirectory);
-
-	if(rank == 0)
-		cout << "Simulation complete." << endl;
-
-	PetscFinalize();
-    return 0;
-}
diff --git a/CoreFlows/examples/C/WaveSystem_FV_SphericalExplosion_CDMATH.cxx b/CoreFlows/examples/C/WaveSystem_FV_SphericalExplosion_CDMATH.cxx
new file mode 100755
index 0000000..0299ef2
--- /dev/null
+++ b/CoreFlows/examples/C/WaveSystem_FV_SphericalExplosion_CDMATH.cxx
@@ -0,0 +1,298 @@
+//============================================================================
+// Author      : Michael NDJINGA
+// Date        : November 2020
+// Description : multiD linear wave system
+//============================================================================
+
+#include <iostream>
+#include <string>
+#include <cmath>
+
+#include "Mesh.hxx"
+#include "Cell.hxx"
+#include "Face.hxx"
+#include "Field.hxx"
+#include "SparseMatrixPetsc.hxx"
+#include "CdmathException.hxx"
+
+using namespace std;
+
+double p0  =155e5;   //reference pressure in a pressurised nuclear vessel
+double c0  =700.;    //reference sound speed for water at 155 bars, 600K
+double rho0=p0/c0*c0;//reference density
+double precision=1e-5;
+
+void initial_conditions_shock(Mesh my_mesh,Field& pressure_field,Field& velocity_field)
+{
+    double rayon=0.35;
+    double xcentre=0.;
+    double ycentre=0;
+    double zcentre=0;
+	
+	double x, y, z;
+	double val, valX, valY, valZ;
+	
+    int dim    =my_mesh.getMeshDimension();
+    int nbCells=my_mesh.getNumberOfCells();
+
+    for (int j=0 ; j<nbCells ; j++)
+    {
+        x = my_mesh.getCell(j).x() ;
+		if(dim>1)
+		{
+			y = my_mesh.getCell(j).y() ;
+			if(dim==3)
+				z = my_mesh.getCell(j).z() ;
+		}
+
+        valX=(x-xcentre)*(x-xcentre);
+		if(dim==1)
+			val=sqrt(valX);
+		else if(dim==2)
+		{
+			valY=(y-ycentre)*(y-ycentre);
+			val=sqrt(valX+valY);		
+		}
+		else if(dim==3)
+		{
+			valY=(y-ycentre)*(y-ycentre);
+			valZ=(z-zcentre)*(z-zcentre);
+			val=sqrt(valX+valY+valZ);		
+		}
+		
+		for(int idim=0; idim<dim; idim++)
+			velocity_field[j,idim]=0;
+			
+        if (val<rayon)
+            pressure_field(j) = 155e5;
+        else
+            pressure_field(j) = 70e5;
+    }
+}
+
+Matrix jacobianMatrices(Vector normal, double coeff)
+{
+    int dim=normal.size();
+    Matrix    A(dim+1,dim+1);
+    Matrix absA(dim+1,dim+1);
+
+    absA(0,0)=c0*coeff;
+    for(int i=0 ; i<dim; i++)
+    {
+        A(i+1,0)=      normal[i]*coeff;
+        A(0,i+1)=c0*c0*normal[i]*coeff;
+        for( int j=0 ; j<dim; j++)
+            absA(i+1,j+1)=c0*normal[i]*normal[j]*coeff;
+    }
+    return (A - absA)*(1./2);
+}
+    
+SparseMatrixPetsc computeDivergenceMatrix(Mesh my_mesh, int nbVoisinsMax, double dt)
+{
+    int nbCells = my_mesh.getNumberOfCells();
+    int dim=my_mesh.getMeshDimension();
+    int nbComp=dim+1;
+    Vector normal(dim);
+
+    SparseMatrixPetsc implMat(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp);
+
+    Matrix idMoinsJacCL(nbComp);
+    
+    for(int j=0; j<nbCells; j++)//On parcourt les cellules
+    {
+        Cell Cj = my_mesh.getCell(j);
+        int nbFaces = Cj.getNumberOfFaces();
+
+        for(int k=0; k<nbFaces; k++)
+        {
+            int indexFace = Cj.getFacesId()[k];
+            Face Fk = my_mesh.getFace(indexFace);
+            for(int i =0; i<dim ; i++)
+                normal[i] = Cj.getNormalVector(k, i);//normale sortante
+
+            Matrix Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure());
+
+            int cellAutre =-1;
+            if ( not Fk.isBorder())
+            {
+                /* hypothese: La cellule d'index indexC1 est la cellule courante index j */
+                if (Fk.getCellsId()[0] == j) 
+                    // hypothese verifiée 
+                    cellAutre = Fk.getCellsId()[1];
+                else if(Fk.getCellsId()[1] == j) 
+                    // hypothese non verifiée 
+                    cellAutre = Fk.getCellsId()[0];
+                else
+                    throw CdmathException("computeDivergenceMatrix: problem with mesh, unknown cell number");
+                    
+                implMat.addValue(j*nbComp,cellAutre*nbComp,Am);
+                implMat.addValue(j*nbComp,        j*nbComp,Am*(-1.));
+            }
+            else 
+            {    
+				if( Fk.getGroupName() != "Periodic" && Fk.getGroupName() != "Neumann")//Wall boundary condition unless Periodic/Neumann specified explicitly
+                {
+					Vector v(dim+1);
+                    for(int i=0; i<dim; i++)
+                        v[i+1]=normal[i];
+                    Matrix idMoinsJacCL=v.tensProduct(v)*2;
+                    
+                    implMat.addValue(j*nbComp,j*nbComp,Am*(-1.)*idMoinsJacCL);
+				}
+                else if( Fk.getGroupName() == "Periodic")//Periodic boundary condition
+                {
+                    int indexFP=my_mesh.getIndexFacePeriodic(indexFace);
+                    Face Fp = my_mesh.getFace(indexFP);
+                    cellAutre = Fp.getCellsId()[0];
+                    
+                    implMat.addValue(j*nbComp,cellAutre*nbComp,Am);
+                    implMat.addValue(j*nbComp,        j*nbComp,Am*(-1.));
+				}
+                else if(Fk.getGroupName() != "Neumann")//Nothing to do for Neumann boundary condition
+                {
+                    cout<< Fk.getGroupName() <<endl;
+                    throw CdmathException("computeDivergenceMatrix: Unknown boundary condition name");
+				}
+            }
+        }   
+    }     
+    return implMat;
+}
+
+void WaveSystem(double tmax, int ntmax, double cfl, int output_freq, const Mesh& my_mesh, const string file)
+{
+	/* Retrieve mesh data */
+    int dim=my_mesh.getMeshDimension();
+    int nbCells = my_mesh.getNumberOfCells();
+    std::string meshName=my_mesh.getName();
+    int nbVoisinsMax=my_mesh.getMaxNbNeighbours(CELLS);
+    double dx_min=my_mesh.minRatioVolSurf();
+
+
+    /* Initial conditions */
+    cout << "Construction of the initial condition" << endl;
+    
+    Field pressure_field("Pressure",CELLS,my_mesh,1) ;
+    Field velocity_field("Velocity",CELLS,my_mesh,dim) ;
+    initial_conditions_shock(my_mesh,pressure_field, velocity_field);
+
+    /* iteration vectors */
+    Vector  Un(nbCells*(dim+1));
+    Vector dUn(nbCells*(dim+1));
+    
+    for(int k =0; k<nbCells; k++){
+        Un[k*(dim+1)+0] =     pressure_field[k];
+        Un[k*(dim+1)+1] =rho0*velocity_field[k,0];
+        if(dim>=2)
+            Un[k + 2*nbCells] = rho0*velocity_field[k,1] ;
+        if(dim==3)
+            Un[k + 3*nbCells] = rho0*velocity_field[k,2];
+		}
+
+    /*
+     * MED output of the initial condition at t=0 and iter = 0
+     */
+    int it=0;
+    bool isStationary=false;
+    double time=0.;
+    double dt = cfl * dx_min / c0;
+    
+    cout << "Saving the solution at T=" << time << endl;
+    pressure_field.setTime(time,it);
+    pressure_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_pressure");
+    velocity_field.setTime(time,it);
+    velocity_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_velocity");
+    /* --------------------------------------------- */
+
+    SparseMatrixPetsc divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt);
+
+    /* Time loop */
+    cout<< "Starting computation of the linear wave system with an explicit UPWIND scheme" << endl;
+    while (it<ntmax && time <= tmax && ! isStationary)
+    {
+        dUn=divMat*Un;
+        Un-=dUn;
+        
+        time=time+dt;
+        it=it+1;
+ 
+         isStationary = dUn.norm()<precision;
+        /* Sauvegardes */
+        if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax)
+        {
+            cout<<"-- Iteration: " << it << ", Time: " << time << ", dt: " << dt<<endl;
+
+            for(int k=0; k<nbCells; k++)
+            {
+                pressure_field[k]  =Un[k*(dim+1)+0];
+                velocity_field[k,0]=Un[k*(dim+1)+1]/rho0;
+                if(dim>1)
+                {
+                    velocity_field[k,1]=Un[k*(dim+1)+2]/rho0;
+                    if(dim>2)
+                        velocity_field[k,2]=Un[k*(dim+1)+3]/rho0;
+				}
+			}
+            pressure_field.setTime(time,it);
+            pressure_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_pressure",false);
+            velocity_field.setTime(time,it);
+            velocity_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_velocity",false);
+		}
+    }
+    cout<<"End of calculation -- Iteration: " << it << ", Time: "<< time<< ", dt: " << dt<<endl;
+
+    if(it>=ntmax)
+        cout<< "Nombre de pas de temps maximum ntmax= "<< ntmax<< " atteint"<<endl;
+    else if(isStationary)
+        cout<< "Régime stationnaire atteint au pas de temps "<< it<< ", t= "<< time<<endl;       
+    else
+        cout<< "Temps maximum Tmax= "<< tmax<< " atteint"<<endl;
+}
+ 
+int main(int argc, char *argv[])
+{
+	cout << "-- Starting the RESOLUTION OF THE 2D WAVE SYSTEM"<<endl;
+	cout << "- Numerical scheme : Upwind explicit scheme" << endl;
+	cout << "- Boundary conditions : WALL" << endl;
+
+    // Problem data
+    double cfl=0.49;
+    double tmax=1.;
+    double ntmax=3;//20000;
+    int freqSortie=100;
+    string fileOutPut="SphericalWave";
+
+	Mesh myMesh;
+	
+	if(argc<2)
+	{
+		    cout << "- DOMAIN : SQUARE" << endl;
+		    cout << "- MESH : CARTESIAN, GENERATED INTERNALLY WITH CDMATH" << endl<< endl;
+		    cout << "Construction of a cartesian mesh" << endl;
+	    double xinf=-0.5;
+	    double xsup= 0.5;
+	    double yinf=-0.5;
+	    double ysup= 0.5;
+	    int nx=50;
+	    int ny=50;
+	    myMesh=Mesh(xinf,xsup,nx,yinf,ysup,ny);
+	    double eps=precision;
+	    myMesh.setGroupAtPlan(xsup,0,eps,"RightEdge");
+	    myMesh.setGroupAtPlan(xinf,0,eps,"LeftEdge");
+	    myMesh.setGroupAtPlan(yinf,1,eps,"BottomEdge");
+	    myMesh.setGroupAtPlan(ysup,1,eps,"TopEdge");
+	}
+	else
+	{
+	    cout << "- MESH:  GENERATED EXTERNALLY WITH SALOME" << endl;
+	    cout << "Loading of a mesh named "<<argv[1] << endl;
+	    string filename = argv[1];
+	    myMesh=Mesh(filename);
+	}
+
+	WaveSystem(tmax,ntmax,cfl,freqSortie,myMesh,fileOutPut);
+	
+    cout << "Simulation complete." << endl;
+
+    return 0;
+}
diff --git a/CoreFlows/examples/C/WaveSystem_FV_SphericalExplosion_MPI.cxx b/CoreFlows/examples/C/WaveSystem_FV_SphericalExplosion_MPI.cxx
new file mode 100755
index 0000000..0356da9
--- /dev/null
+++ b/CoreFlows/examples/C/WaveSystem_FV_SphericalExplosion_MPI.cxx
@@ -0,0 +1,404 @@
+//============================================================================
+// Author      : Michael NDJINGA
+// Date        : November 2020
+// Description : multiD linear wave system run in parallel with calls to MPI
+//               Test used in the MPI version of the Salome platform
+//============================================================================
+
+#include <iostream>
+#include <string>
+#include <cmath>
+
+#include "Mesh.hxx"
+#include "Cell.hxx"
+#include "Face.hxx"
+#include "Field.hxx"
+#include "CdmathException.hxx"
+
+#include <petscksp.h>
+
+using namespace std;
+
+double p0  =155e5;   //reference pressure in a pressurised nuclear vessel
+double c0  =700.;    //reference sound speed for water at 155 bars, 600K
+double rho0=p0/c0*c0;//reference density
+double precision=1e-5;
+
+void initial_conditions_shock(Mesh my_mesh,Field& pressure_field,Field& velocity_field)
+{
+    double rayon=0.35;
+    double xcentre=0.;
+    double ycentre=0;
+    double zcentre=0;
+	
+	double x, y, z;
+	double val, valX, valY, valZ;
+	
+    int dim    =my_mesh.getMeshDimension();
+    int nbCells=my_mesh.getNumberOfCells();
+
+    for (int j=0 ; j<nbCells ; j++)
+    {
+        x = my_mesh.getCell(j).x() ;
+		if(dim>1)
+		{
+			y = my_mesh.getCell(j).y() ;
+			if(dim==3)
+				z = my_mesh.getCell(j).z() ;
+		}
+
+        valX=(x-xcentre)*(x-xcentre);
+		if(dim==1)
+			val=sqrt(valX);
+		else if(dim==2)
+		{
+			valY=(y-ycentre)*(y-ycentre);
+			val=sqrt(valX+valY);		
+		}
+		else if(dim==3)
+		{
+			valY=(y-ycentre)*(y-ycentre);
+			valZ=(z-zcentre)*(z-zcentre);
+			val=sqrt(valX+valY+valZ);		
+		}
+		
+		for(int idim=0; idim<dim; idim++)
+			velocity_field[j,idim]=0;
+			
+        if (val<rayon)
+            pressure_field(j) = 155e5;
+        else
+            pressure_field(j) = 70e5;
+    }
+}
+
+void addValue( int i, int j, Matrix M, Mat * mat )
+{
+    int I,J;
+    for (int k=0; k<M.getNumberOfRows(); k++)
+        for (int l=0; l<M.getNumberOfColumns(); l++)
+        {
+            I=i+k;
+            J=j+l;
+            MatSetValues( *mat,1, &I, 1, &J, &M(k,l), ADD_VALUES);
+        }
+}
+
+Matrix jacobianMatrices(Vector normal, double coeff)
+{
+    int dim=normal.size();
+    Matrix    A(dim+1,dim+1);
+    Matrix absA(dim+1,dim+1);
+
+    absA(0,0)=c0*coeff;
+    for(int i=0 ; i<dim; i++)
+    {
+        A(i+1,0)=      normal[i]*coeff;
+        A(0,i+1)=c0*c0*normal[i]*coeff;
+        for( int j=0 ; j<dim; j++)
+            absA(i+1,j+1)=c0*normal[i]*normal[j]*coeff;
+    }
+    return (A - absA)*(1./2);
+}
+    
+void computeDivergenceMatrix(Mesh my_mesh, Mat * implMat, double dt)
+{
+    int nbCells = my_mesh.getNumberOfCells();
+    int dim=my_mesh.getMeshDimension();
+    int nbComp=dim+1;
+    Vector normal(dim);
+
+    Matrix idMoinsJacCL(nbComp);
+    
+    for(int j=0; j<nbCells; j++)//On parcourt les cellules
+    {
+        Cell Cj = my_mesh.getCell(j);
+        int nbFaces = Cj.getNumberOfFaces();
+
+        for(int k=0; k<nbFaces; k++)
+        {
+            int indexFace = Cj.getFacesId()[k];
+            Face Fk = my_mesh.getFace(indexFace);
+            for(int i =0; i<dim ; i++)
+                normal[i] = Cj.getNormalVector(k, i);//normale sortante
+
+            Matrix Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure());
+
+            int cellAutre =-1;
+            if ( not Fk.isBorder())
+            {
+                /* hypothese: La cellule d'index indexC1 est la cellule courante index j */
+                if (Fk.getCellsId()[0] == j) 
+                    // hypothese verifiée 
+                    cellAutre = Fk.getCellsId()[1];
+                else if(Fk.getCellsId()[1] == j) 
+                    // hypothese non verifiée 
+                    cellAutre = Fk.getCellsId()[0];
+                else
+                    throw CdmathException("computeDivergenceMatrix: problem with mesh, unknown cell number");
+                    
+                addValue(j*nbComp,cellAutre*nbComp,Am      ,implMat);
+                addValue(j*nbComp,        j*nbComp,Am*(-1.),implMat);
+            }
+            else 
+            {    
+				if( Fk.getGroupName() != "Periodic" && Fk.getGroupName() != "Neumann")//Wall boundary condition unless Periodic/Neumann specified explicitly
+                {
+					Vector v(dim+1);
+                    for(int i=0; i<dim; i++)
+                        v[i+1]=normal[i];
+                    Matrix idMoinsJacCL=v.tensProduct(v)*2;
+                    
+                    addValue(j*nbComp,j*nbComp,Am*(-1.)*idMoinsJacCL,implMat);
+				}
+                else if( Fk.getGroupName() == "Periodic")//Periodic boundary condition
+                {
+                    int indexFP=my_mesh.getIndexFacePeriodic(indexFace);
+                    Face Fp = my_mesh.getFace(indexFP);
+                    cellAutre = Fp.getCellsId()[0];
+                    
+                    addValue(j*nbComp,cellAutre*nbComp,Am      ,implMat);
+                    addValue(j*nbComp,        j*nbComp,Am*(-1.),implMat);
+				}
+                else if(Fk.getGroupName() != "Neumann")//Nothing to do for Neumann boundary condition
+                {
+                    cout<< Fk.getGroupName() <<endl;
+                    throw CdmathException("computeDivergenceMatrix: Unknown boundary condition name");
+				}
+            }
+        }   
+    }     
+}
+
+void WaveSystem(double tmax, int ntmax, double cfl, int output_freq, const Mesh& my_mesh, const string file, int rank, int size, string resultDirectory)
+{
+	/* Time iteration variables */
+    int it=0;
+    bool isStationary=false;
+    double time=0.;
+    double dt;
+    double norm;
+	    
+	/* PETSc variables */
+	int globalNbUnknowns;
+	int  localNbUnknowns;
+	int d_nnz, o_nnz;
+	Vec Un, dUn, Un_seq;
+    Mat divMat;
+	VecScatter scat;
+
+	
+	/* Mesh parameters managed only by proc 0 */
+    int nbCells;
+    int dim;
+    int nbComp;	    
+	int idx;//Index where to add the vector values
+	double value;//value to add in the vector	
+    Field pressure_field, velocity_field;
+    std::string meshName;
+
+	if(rank == 0)
+	    globalNbUnknowns=my_mesh.getNumberOfCells()*(my_mesh.getMeshDimension()+1);//nbCells*nbComp
+	
+	MPI_Bcast(&globalNbUnknowns, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ 
+    /* iteration vectors */
+	VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE    ,globalNbUnknowns,&Un);
+	VecDuplicate (Un,&dUn);
+	if(rank == 0)
+		VecCreateSeq(PETSC_COMM_SELF,globalNbUnknowns,&Un_seq);//For saving results on proc 0
+
+	VecScatterCreateToZero(Un,&scat,&Un_seq);
+
+	/* System matrix */
+	VecGetLocalSize(Un, &localNbUnknowns);
+
+	if(rank == 0)
+	{
+		int nbVoisinsMax = my_mesh.getMaxNbNeighbours(CELLS);
+	    d_nnz=(nbVoisinsMax+1)*(my_mesh.getMeshDimension()+1);//(nbVoisinsMax+1)*nbComp
+	    o_nnz= nbVoisinsMax   *(my_mesh.getMeshDimension()+1);//                 nbComp
+	}
+	MPI_Bcast(&d_nnz, 1, MPI_INT, 0, MPI_COMM_WORLD);
+	MPI_Bcast(&o_nnz, 1, MPI_INT, 0, MPI_COMM_WORLD);
+
+   	MatCreateAIJ(PETSC_COMM_WORLD,localNbUnknowns,localNbUnknowns,globalNbUnknowns,globalNbUnknowns,d_nnz,NULL,o_nnz,NULL,&divMat);
+
+	if(rank == 0)
+		{
+		/* Retrieve mesh data */
+	    nbCells = my_mesh.getNumberOfCells();
+	    dim=my_mesh.getMeshDimension();
+	    nbComp=dim+1;	    
+	    meshName=my_mesh.getName();
+	    double dx_min=my_mesh.minRatioVolSurf();
+		dt = cfl * dx_min / c0;
+		
+	    /* Initial conditions */
+	    cout<<"Building the initial condition on processor 0" << endl;
+	    
+	    pressure_field=Field("Pressure",CELLS,my_mesh,1) ;
+	    velocity_field=Field("Velocity",CELLS,my_mesh,dim) ;
+	    initial_conditions_shock(my_mesh,pressure_field, velocity_field);
+	
+	    cout << "Saving the solution at T=" << time <<"  on processor 0"<<endl;
+	    pressure_field.setTime(time,it);
+	    pressure_field.writeMED(resultDirectory+"/WaveSystem"+to_string(dim)+"DUpwind_"+to_string(size)+"Procs_"+meshName+"_pressure");
+	    velocity_field.setTime(time,it);
+	    velocity_field.writeMED(resultDirectory+"/WaveSystem"+to_string(dim)+"DUpwind_"+to_string(size)+"Procs_"+meshName+"_velocity");
+	    /* --------------------------------------------- */
+	
+
+	    for(int k =0; k<nbCells; k++)
+	    {
+			idx = k*nbComp;
+			value=pressure_field[k];//vale to add in the vector
+			VecSetValues(Un,1,&idx,&value,INSERT_VALUES);
+			for(int idim =0; idim<dim; idim++)
+			{
+				idx = k*nbComp+1+idim;
+				value =rho0*velocity_field[k,idim];
+				VecSetValues(Un,1,&idx,&value,INSERT_VALUES);
+			}
+		}
+	    computeDivergenceMatrix(my_mesh,&divMat,dt);
+	}		
+
+	VecAssemblyBegin(Un);
+	VecAssemblyEnd(Un);
+	
+	//VecView(Un, PETSC_VIEWER_STDOUT_WORLD );
+
+	MatAssemblyBegin(divMat, MAT_FINAL_ASSEMBLY);
+	MatAssemblyEnd(  divMat, MAT_FINAL_ASSEMBLY);
+
+	//MatView(divMat,	PETSC_VIEWER_STDOUT_WORLD );
+
+	MPI_Bcast(&dt, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+
+    /* Time loop */
+	PetscPrintf(PETSC_COMM_WORLD,"Starting computation of the linear wave system on all processors : \n\n");
+    while (it<ntmax && time <= tmax && !isStationary)
+    {
+		MatMult(divMat,Un,dUn);
+        VecAXPY(Un,-1,dUn);
+        
+        time=time+dt;
+        it=it+1;
+ 
+		VecNorm(dUn,NORM_2,&norm);
+        isStationary = norm<precision;
+        /* Sauvegardes */
+        if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax )
+        {
+			PetscPrintf(PETSC_COMM_WORLD,"-- Iteration: %d, Time: %f, dt: %f, saving results on processor 0 \n", it, time, dt);
+			VecScatterBegin(scat,Un,Un_seq,INSERT_VALUES,SCATTER_FORWARD);
+			VecScatterEnd(  scat,Un,Un_seq,INSERT_VALUES,SCATTER_FORWARD);
+
+			if(rank == 0)
+			{
+	            for(int k=0; k<nbCells; k++)
+	            {
+					idx = k*(dim+1)+0;
+					VecGetValues(Un_seq,1,&idx,&value);
+	                pressure_field[k]  =value;
+					for(int idim =0; idim<dim; idim++)
+					{
+						idx = k*nbComp+1+idim;
+						VecGetValues(Un_seq,1,&idx,&value);
+						velocity_field[k,idim] = value/rho0;
+					}
+				}
+	            pressure_field.setTime(time,it);
+	            pressure_field.writeMED(resultDirectory+"/WaveSystem"+to_string(dim)+"DUpwind_"+to_string(size)+"Procs_"+meshName+"_pressure",false);
+	            velocity_field.setTime(time,it);
+	            velocity_field.writeMED(resultDirectory+"/WaveSystem"+to_string(dim)+"DUpwind_"+to_string(size)+"Procs_"+meshName+"_velocity",false);
+			}
+		}
+    }
+
+	PetscPrintf(PETSC_COMM_WORLD,"\n End of calculations at iteration: %d, and time: %f\n", it, time);
+    if(it>=ntmax)
+        PetscPrintf(PETSC_COMM_WORLD, "Nombre de pas de temps maximum ntmax= %d atteint\n", ntmax);
+    else if(isStationary)
+        PetscPrintf(PETSC_COMM_WORLD, "Régime stationnaire atteint au pas de temps %d, t= %f\n", it, time);       
+    else
+        PetscPrintf(PETSC_COMM_WORLD, "Temps maximum tmax= %f atteint\n", tmax);
+
+	VecDestroy(&Un);
+	VecDestroy(&Un_seq);
+	VecDestroy(&dUn);
+	MatDestroy(&divMat);
+}
+ 
+int main(int argc, char *argv[])
+{
+	/* PETSc initialisation */
+	PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
+	PetscMPIInt    size;        /* size of communicator */
+	PetscMPIInt    rank;        /* processor rank */
+	MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
+	MPI_Comm_size(PETSC_COMM_WORLD,&size);
+	
+    // Problem data
+    double cfl=0.49;
+    double tmax=1.;
+    int ntmax=2000;
+    int freqSortie=10;
+    string fileOutPut="SphericalWave";
+    Mesh myMesh;
+	string resultDirectory="./";
+	
+	if(size>1)
+		PetscPrintf(PETSC_COMM_WORLD,"---- More than one processor detected : running a parallel simulation ----\n");
+		PetscPrintf(PETSC_COMM_WORLD,"---- Limited parallelism : input and output remain sequential ----\n");
+		PetscPrintf(PETSC_COMM_WORLD,"---- Only the matrix-vector products are done in parallel ----\n");
+		PetscPrintf(PETSC_COMM_WORLD,"---- Processor 0 is in charge of building the mesh, saving the results, filling and then distributing the matrix to other processors.\n\n");
+		
+	if(rank == 0)
+	{
+	    cout << "-- Starting the RESOLUTION OF THE 2D WAVE SYSTEM on "<< size <<" processors"<<endl;
+	    cout << "- Numerical scheme : Upwind explicit scheme" << endl;
+	    cout << "- Boundary conditions : WALL" << endl;
+	
+		/* Read or create mesh */
+		if(argc<2)
+		{
+		    cout << "- DOMAIN : SQUARE" << endl;
+		    cout << "- MESH : CARTESIAN, GENERATED INTERNALLY WITH CDMATH" << endl<< endl;
+		    cout << "Construction of a cartesian mesh on processor 0" << endl;
+		    double xinf=-0.5;
+		    double xsup= 0.5;
+		    double yinf=-0.5;
+		    double ysup= 0.5;
+		    int nx=50;
+		    int ny=50;
+		    myMesh=Mesh(xinf,xsup,nx,yinf,ysup,ny);
+		    
+		    double eps=precision;
+		    myMesh.setGroupAtPlan(xsup,0,eps,"RightEdge");
+		    myMesh.setGroupAtPlan(xinf,0,eps,"LeftEdge");
+		    myMesh.setGroupAtPlan(yinf,1,eps,"BottomEdge");
+		    myMesh.setGroupAtPlan(ysup,1,eps,"TopEdge");
+		}
+		else
+		{
+		    cout << "- MESH:  GENERATED EXTERNALLY WITH SALOME" << endl<< endl;
+		    cout << "Loading of a mesh named "<<argv[1]<<" on processor 0" << endl;
+		    string filename = argv[1];
+		    myMesh=Mesh(filename);
+		}
+
+		/*Detect directory where to same results*/
+		if(argc>2)
+		    resultDirectory = argv[2];
+	}
+	
+	WaveSystem(tmax,ntmax,cfl,freqSortie,myMesh,fileOutPut, rank, size, resultDirectory);
+
+	if(rank == 0)
+		cout << "Simulation complete." << endl;
+
+	PetscFinalize();
+    return 0;
+}