diff --git a/xfv/XtendedFiniteVolume.py b/xfv/XtendedFiniteVolume.py
old mode 100755
new mode 100644
index 246274a..efc38aa
--- a/xfv/XtendedFiniteVolume.py
+++ b/xfv/XtendedFiniteVolume.py
@@ -31,11 +31,11 @@ def __create_mesh(meshfile: Path) -> Mesh1dEnriched:
     :param meshfile: path to the mesh file
     :param enrichment_type: type of enrichment desired
     """
-    coord_mesh = np.loadtxt(meshfile, dtype=np.float64, skiprows=2, usecols=(1,))
-    nodes_number = coord_mesh.shape[0]
+    coord_mesh = np.loadtxt(meshfile, dtype=np.float64, skiprows=2, usecols=(1,)) 
+    nodes_number = coord_mesh.shape[0] # donne le nombre des noeuds
     coord_init = np.zeros([nodes_number, 1], dtype=np.float64, order='C')
-    coord_init[:, 0] = coord_mesh
-    vit_init = np.zeros([nodes_number, 1], dtype=np.float64, order='C')
+    coord_init[:, 0] = coord_mesh  
+    vit_init = np.zeros([nodes_number, 1], dtype=np.float64, order='C') # fixe les vitesses initales des noeuds à zéro
     return Mesh1dEnriched(initial_coordinates=coord_init, initial_velocities=vit_init)
 
 
@@ -58,7 +58,7 @@ def __init_velocity(nodes, data):
         vitesse_cible = 0
 
     print("Projectile velocity initialized to : {:} m/s".format(vitesse_projectile))
-    nodes.upundemi[nodes.nodes_in_projectile, 0] = vitesse_projectile
+    nodes.upundemi[nodes.nodes_in_projectile, 0] = vitesse_projectile  # demander signification
     nodes.umundemi[nodes.nodes_in_projectile, 0] = vitesse_projectile
     print("Target velocity initialized to : {:} m/s".format(vitesse_cible))
     nodes.upundemi[nodes.nodes_in_target, 0] = vitesse_cible
@@ -253,7 +253,7 @@ def main(directory: Path) -> None:
         projectile_model = data.material_projectile.constitutive_model
         (projectile_elasticity, projectile_plasticity, projectile_shear_modulus,
          projectile_yield_stress, projectile_plasticity_criterion) = \
-            _build_material_constitutive_model(projectile_model)
+         _build_material_constitutive_model(projectile_model)
     else:
         projectile_elasticity, projectile_plasticity = False, False
         projectile_shear_modulus = None
diff --git a/xfv/post_processing/cohesive_dissipated_energy.py b/xfv/post_processing/cohesive_dissipated_energy.py
new file mode 100644
index 0000000..89cde5c
--- /dev/null
+++ b/xfv/post_processing/cohesive_dissipated_energy.py
@@ -0,0 +1,107 @@
+# -*- coding: utf-8 -*-
+"""
+Plot the dissipated energy by the cohesive law eventually with experimental data
+"""
+
+import argparse
+import pathlib
+import numpy as np
+import matplotlib.pyplot as plt
+from collections import namedtuple
+
+from xfv.src.output_manager.outputdatabaseexploit import OutputDatabaseExploit
+
+CellInfo = namedtuple("CellInfo", ["ouverture_min", "ouverture_max", "temps_apparition"])
+
+def run():
+    """
+    Run post processing program
+    """
+    # ----------------------------------------------------------
+    # Read instructions
+    # ----------------------------------------------------------
+    parser = argparse.ArgumentParser()
+    parser.add_argument("-case", action='append', nargs='+',
+                        help="the path to the output repository")
+    parser.add_argument("-item_ids", type=int, action='append', nargs='+', help="the id of the item to look at")
+    parser.add_argument("--output_filename", default="all_fields.hdf5",
+                        help="the name of the output hdf5 band (default = all_fields.hdf5)")
+    parser.add_argument("--write_data", action="store_true",
+                        help="Write a file with time and dissipated energy")
+    args = parser.parse_args()
+
+    if args.case is None:
+        raise ValueError("At least one case is needed. Use -case to specify cases to plot")
+
+    # ----------------------------------------------------------
+    # Prepare figure
+    # ----------------------------------------------------------
+    plt.figure(1)
+    plt.title("Evolution of the dissipated energy by the cohesive model", fontweight='bold', fontsize=18)
+    plt.xlabel("Time [mus]", fontsize=16)
+    plt.ylabel("Dissipated energy [J ou J/kg selon les cas]", fontsize=16)
+
+
+    # ----------------------------------------------------------
+    # Read dissipated energy for each cell
+    # ----------------------------------------------------------
+    for case in args.case[0]:
+        dissipated_energy_dict = {}
+        path_to_db = pathlib.Path.cwd().joinpath(case, args.output_filename)
+        # Read database :
+        hd_band = OutputDatabaseExploit(path_to_db)
+        for time in hd_band.saved_times:
+            cell_status = hd_band.extract_field_at_time("CellStatus", time)[:]
+            enriched_cells = np.where(cell_status)[0]
+            if len(enriched_cells) > 0:
+                dissipated_energy = hd_band.extract_field_at_time("AdditionalDissipatedEnergy", time)[:]
+                #print('time=',time)
+
+                for i in range(len(dissipated_energy[:, 0])):
+                    cell_id = dissipated_energy[i, 0]
+                    dis_energy = dissipated_energy[i, 1]
+                    try:
+                        dissipated_energy_dict[cell_id].append([time, dis_energy])
+                    except KeyError:  # 1ere fois que la maille est enrichie => init list opening
+                        dissipated_energy_dict[cell_id] = [[time, dis_energy]]
+                    #print(dissipated_energy_dict[cell_id][0,1])
+
+        # ----------------------------------------------------------
+        # Permet d'aficher les energies dissipées pour toutes les cellules
+        # ----------------------------------------------------------         
+        if args.item_ids[0][0] == 1000:
+            args.item_ids[0] = []
+            args.item_ids[0] = dissipated_energy[:,0]
+        
+
+        #Transformation en array
+        for key in dissipated_energy_dict:
+            dissipated_energy_dict[key] = np.array(dissipated_energy_dict[key])
+
+        # import ipdb; ipdb.set_trace()
+
+        for item_id in args.item_ids[0]:
+            # Read database :
+            j = 0.
+            for ane in dissipated_energy[:, 0]:
+                if ane == item_id:
+                    j += 1.
+            if j < 0.99:
+                print(dissipated_energy[:, 0])
+                exit(f'item_ids selectionné ne fait pas parti de la liste des cellules enrichies. Choississez item_ids (ie. le numéro de cellules) dans la liste ci-dessus ou tapez 1000 pour selectionner toutes les cellules')
+            # Plot field :
+            plt.plot(dissipated_energy_dict[item_id][:,0],dissipated_energy_dict[item_id][:,1],label= 'cell n°'+str(item_id))
+            if args.write_data:
+                data_path = f"Field_evolution_cohesive_dissipated_energy_at_cell_{item_id}.dat"
+                with open(data_path, "w") as file_object:
+                    for x_data, y_data in zip(dissipated_energy_dict[item_id][:,0], dissipated_energy_dict[item_id][:,1]):
+                        file_object.write("{:20.18g}\t{:20.18g}\n".format(x_data, y_data))
+                print("Data written in {:s}".format(data_path))
+
+if __name__ == "__main__":
+    run()
+    # ----------------------------------------------------------
+    # Show figure
+    # ----------------------------------------------------------
+    plt.legend(loc="right")
+    plt.show()
diff --git a/xfv/post_processing/cohesive_law.py b/xfv/post_processing/cohesive_law.py
new file mode 100755
index 0000000..8eeb4de
--- /dev/null
+++ b/xfv/post_processing/cohesive_law.py
@@ -0,0 +1,108 @@
+# -*- coding: utf-8 -*-
+"""
+Plot the free surface velocity eventually with experimental data
+"""
+
+import argparse
+import pathlib
+import numpy as np
+import matplotlib.pyplot as plt
+from collections import namedtuple
+
+from xfv.src.output_manager.outputdatabaseexploit import OutputDatabaseExploit
+
+CellInfo = namedtuple("CellInfo", ["ouverture_min", "ouverture_max", "temps_apparition"])
+
+def run():
+    """
+    Run post processing program
+    """
+    # ----------------------------------------------------------
+    # Read instructions
+    # ----------------------------------------------------------
+    parser = argparse.ArgumentParser()
+    parser.add_argument("-case", action='append', nargs='+',
+                        help="the path to the output repository")
+    parser.add_argument("-item_ids", type=int, action='append', nargs='+', help="the id of the item to look at")
+    parser.add_argument("--output_filename", default="all_fields.hdf5",
+                        help="the name of the output hdf5 band (default = all_fields.hdf5)")
+    parser.add_argument("--write_data", action="store_true", help="To write data in an output file")
+    args = parser.parse_args()
+
+    if args.case is None:
+        raise ValueError("At least one case is needed. Use -case to specify cases to plot")
+
+    # ----------------------------------------------------------
+    # Prepare figure
+    # ----------------------------------------------------------
+    plt.figure(1)
+    plt.title("Evolution of the cohesive law", fontweight='bold', fontsize=18)
+    plt.xlabel("discontinuity opening [m]", fontsize=16)
+    plt.ylabel("cohesive force []", fontsize=16)
+
+
+    # ----------------------------------------------------------
+    # Read discontinuity opening for each cell
+    # ----------------------------------------------------------
+    for case in args.case[0]:
+        opening_dict = {}
+        force_dict = {}
+        path_to_db = pathlib.Path.cwd().joinpath(case, args.output_filename)
+        # Read database :
+        hd_band = OutputDatabaseExploit(path_to_db)
+        for time in hd_band.saved_times:
+            cell_status = hd_band.extract_field_at_time("CellStatus", time)[:]
+            enriched_cells = np.where(cell_status)[0]
+            if len(enriched_cells) > 0:
+                opening = hd_band.extract_field_at_time("AdditionalDiscontinuityOpening", time)[:]
+                cohesive_force = hd_band.extract_field_at_time("AdditionalCohesiveForce", time)[:]
+                for i in range(len(opening[:, 0])):
+                    cell_id = opening[i, 0]
+                    op = opening[i, 1]
+                    force = cohesive_force[i, 1]
+                    try:
+                        opening_dict[cell_id].append([time, op])
+                        force_dict[cell_id].append([time, force])
+                    except KeyError:  # 1ere fois que la maille est enrichie => init list opening
+                        opening_dict[cell_id] = [[time, op]]
+                        force_dict[cell_id] = [[time, force]]
+        
+    # Permet de selectionner les item_ids pour toutes les cellules enrichies
+    if args.item_ids[0][0] == 1000:
+        args.item_ids[0] = []
+        args.item_ids[0] = opening[:,0]
+        print(args.item_ids[0])
+
+    #Transformation en array
+    for key in opening_dict:
+        opening_dict[key] = np.array(opening_dict[key])
+    for key in force_dict:
+        force_dict[key] = np.array(force_dict[key])
+        # import ipdb; ipdb.set_trace()
+
+    for item_id in args.item_ids[0]:
+        # Read database :
+        j = 0.
+        for ane in opening[:, 0]:
+            if ane == item_id:
+                j += 1.
+                if j < 0.99:
+                    print(opening[:, 0])
+                    exit(f'Choissisez item_ids dans la liste ci-dessus')
+            # Plot field :
+        plt.plot(opening_dict[item_id][:,1],force_dict[item_id][:,1], '+', label= 'cell n°'+str(item_id))
+
+        if args.write_data:
+            data_path = f"Field_evolution_cohesive_law_{item_id}.dat"
+            with open(data_path, "w") as file_object:
+                for x_data, y_data in zip(opening_dict[item_id][:,1], force_dict[item_id][:,1]):
+                    file_object.write("{:20.18g}\t{:20.18g}\n".format(x_data, y_data))
+            print("Data written in {:s}".format(data_path))
+
+if __name__ == "__main__":
+    run()
+    # ----------------------------------------------------------
+    # Show figure
+    # ----------------------------------------------------------
+    plt.legend(loc="best")
+    plt.show()
diff --git a/xfv/post_processing/field_evolution_in_time.py b/xfv/post_processing/field_evolution_in_time.py
index b98669d..8b791f0 100755
--- a/xfv/post_processing/field_evolution_in_time.py
+++ b/xfv/post_processing/field_evolution_in_time.py
@@ -6,6 +6,7 @@
 import argparse
 import pathlib
 import matplotlib.pyplot as plt
+import numpy as np
 
 from xfv.post_processing.tools.hdf5_postprocessing_tools import get_field_evolution_in_time_for_item
 
@@ -68,6 +69,10 @@ def run():
         if args.verbose:
             print("Path to database : {:}".format(path_to_db))
             print("Read field " + field + " in database... ")
+        # Permet de selectionner les items ids de toutes les cellules
+        if args.item_ids[0][0] == 1000:
+            args.item_ids[0] = []
+            args.item_ids[0] = [i for i in range(350,734)]
 
         for item_id in args.item_ids[0]:
             # Read database :
@@ -76,9 +81,10 @@ def run():
                 print("Done !")
                 print("~~~~~~~~~~~~~")
             # Plot field :
-            plt.plot(item_history[:, 0] * 1.e+6, item_history[:, 1], '.-', label=case)
+            plt.plot(item_history[:, 0] * 1.e+6, item_history[:, 1], '.-', label= 'cell n°'+ str(item_id))
             if args.write_data:
-                data_path = f"{case}Field_evolution_{field}_{item_id}.dat"
+                data_path = f"Field_evolution_{field}_at_cell_{item_id}.dat"
+                #data_path = pathlib.Path.cwd().joinpath(case, f"Field_evolution_{field}_{item_id}.txt")
                 with open(data_path, "w") as file_object:
                     for x_data, y_data in zip(item_history[:, 0], item_history[:, 1]):
                         file_object.write("{:20.18g}\t{:20.18g}\n".format(x_data, y_data))
diff --git a/xfv/post_processing/free_surface_velocity.py b/xfv/post_processing/free_surface_velocity.py
old mode 100755
new mode 100644
index e958125..567cc9b
--- a/xfv/post_processing/free_surface_velocity.py
+++ b/xfv/post_processing/free_surface_velocity.py
@@ -7,6 +7,7 @@
 import pathlib
 import numpy as np
 import matplotlib.pyplot as plt
+import tikzplotlib
 
 from xfv.post_processing.tools.hdf5_postprocessing_tools import get_field_evolution_in_time_for_item
 
@@ -82,11 +83,12 @@ def run():
                 print("New t0 is : " + str(time_0))
         plt.plot((time - time_0) * 1.e+6, velocity, label=case)
         if args.write_data:
-            data_path = pathlib.Path.cwd().joinpath(case, args.file_write_data)
+            data_path = f"Free_surface_velocity_evolution.dat"
+            #data_path = pathlib.Path.cwd().joinpath(case, args.file_write_data)
             with open(data_path, "w") as file_object:
                 for x_data, y_data in zip(time, velocity):
                     file_object.write("{:20.18g}\t{:20.18g}\n".format(x_data, y_data))
-            print("Data written in {:s}".format(data_path))
+            print("Data written in {:s}."+str(data_path))
 
     # ----------------------------------------------------------
     # Plot experimental data
diff --git a/xfv/post_processing/profile_figure_in_space.py b/xfv/post_processing/profile_figure_in_space.py
index 2a2ddf9..40e0482 100755
--- a/xfv/post_processing/profile_figure_in_space.py
+++ b/xfv/post_processing/profile_figure_in_space.py
@@ -39,14 +39,13 @@ def run():
                 print("Done !")
                 print("~~~~~~~~~~~~~")
             # Plot field :
-            plt.plot(coord * 1.e+03, field_value, label=case)
-
-        if (ARGS.write_data):
-            data_path="{:s}Profile_{:s}_{:3.1e}.dat".format(case, ARGS.field, current_time)
-            with open(data_path, "w") as file_object:
-                for x_data, y_data in zip(coord, field_value):
-                    file_object.write("{:20.18g}\t{:20.18g}\n".format(x_data, y_data))
-            print("Data written in {:s}".format(data_path))
+            plt.plot(coord * 1.e+03, field_value, label='temps = '+str(current_time)+' s')
+            if (ARGS.write_data):
+                data_path="Profile_{:s}_at_time_{:3.1e}.dat".format(ARGS.field, current_time)
+                with open(data_path, "w") as file_object:
+                    for x_data, y_data in zip(coord, field_value):
+                        file_object.write("{:20.18g}\t{:20.18g}\n".format(x_data, y_data))
+                print("Data written in {:s}".format(data_path))
 
 
 if __name__ == "__main__":
diff --git a/xfv/src/cell/cell.py b/xfv/src/cell/cell.py
old mode 100755
new mode 100644
index 21feece..a9fc05d
--- a/xfv/src/cell/cell.py
+++ b/xfv/src/cell/cell.py
@@ -85,6 +85,8 @@ def __init__(self, nbr_of_cells: int):
             self._nbr_of_cells, material_data.shear_modulus_init, material_data.shear_modulus_init)
         self._fields_manager["YieldStress"] = Field(
             self._nbr_of_cells, material_data.yield_stress_init, material_data.yield_stress_init)
+        # Plasticity
+        self._fields_manager["EquivalentPlasticStrain"] = Field(self._nbr_of_cells, 0., 0.)
         # Porosity
         self._fields_manager["Porosity"] = Field(
             self._nbr_of_cells, material_data.porosity_init, material_data.porosity_init)
@@ -227,6 +229,13 @@ def yield_stress(self):
         """
         return self._fields_manager["YieldStress"]
 
+    @property
+    def equivalent_plastic_strain(self):
+        """
+        equivalent plastic strain
+        """
+        return self._fields_manager["EquivalentPlasticStrain"]
+
     @property
     def stress(self):
         """
diff --git a/xfv/src/cell/one_dimension_cell.py b/xfv/src/cell/one_dimension_cell.py
old mode 100755
new mode 100644
index a019f14..47869d1
--- a/xfv/src/cell/one_dimension_cell.py
+++ b/xfv/src/cell/one_dimension_cell.py
@@ -96,7 +96,7 @@ def apply_equation_of_state(cls, cell: Cell, eos, density: np.array, density_new
             if np.isnan(sound_velocity_new_value).any():
                 negative_vson = np.where(np.isnan(sound_velocity_new_value))
                 msg = "Sound speed square < 0 in cells {}\n".format(np.where(negative_vson))
-                msg += "density = {}\n".format(my_variables['NewDensity'][negative_vson])
+                msg += "density = {}\n".format(density[negative_vson])
                 msg += "energy = {}\n".format(energy_new_value[negative_vson])
                 msg += "pressure = {}\n".format(pressure_new_value[negative_vson])
                 raise ValueError(msg)
@@ -237,7 +237,7 @@ def __init__(self, number_of_elements: int):
                                                         order='C')
         self._plastic_strain_rate = np.zeros([number_of_elements, 3], dtype=np.float64,
                                              order='C')  # Dp_xx, Dp_yy, Dp_zz
-        # EOS :
+	# EOS :
         self._target_eos = self.data.material_target.constitutive_model.eos.build_eos_obj()
         self._projectile_eos = None
         if self.data.data_contains_a_projectile:
@@ -545,7 +545,7 @@ def compute_shear_modulus(self, shear_modulus_model, mask):
         :type mask: np.array([nbr_cells, 1], dtype=bool)
         """
         self.shear_modulus.new_value[mask] = shear_modulus_model.compute(
-            self.density.new_value[mask])
+            self.density.new_value[mask], self.pressure.current_value[mask])
 
         if self.data.material_target.porosity_model is not None:
             self.shear_modulus.new_value[mask] = self._compute_micro_to_macro_shear_modulus(
@@ -579,7 +579,8 @@ def compute_yield_stress(self, yield_stress_model, mask):
 
         :type mask: np.array([nbr_cells, 1], dtype=bool)
         """
-        self.yield_stress.new_value[mask] = yield_stress_model.compute(self.density.new_value[mask])
+        self.yield_stress.new_value[mask] = yield_stress_model.compute(self.density.new_value[mask],
+                                                                       self.equivalent_plastic_strain.current_value[mask], self.shear_modulus.new_value[mask])
 
         if self.data.material_target.porosity_model is not None:
             self.yield_stress.new_value[mask] = self._compute_micro_to_macro_yield_stress(
@@ -707,6 +708,13 @@ def _compute_equivalent_plastic_strain_rate(j2, shear_modulus, yield_stress, dt)
         """
         return  (j2 - yield_stress) / (3. * shear_modulus * dt)
 
+    def _compute_equivalent_plastic_strain(self,mask,delta_t):
+        """
+        Compute the equivalent plastic strain
+        param t_time : temps t
+        """
+        self.equivalent_plastic_strain.new_value[:] = self._equivalent_plastic_strain_rate[:]*delta_t + self.equivalent_plastic_strain.current_value[:]
+    
     def apply_plasticity(self, mask, delta_t):
         """
         Apply plasticity treatment if criterion is activated :
@@ -717,6 +725,7 @@ def apply_plasticity(self, mask, delta_t):
         :param mask: boolean array to identify cells to be computed
         :param dt: time step
         """
+
         invariant_j2_el = np.sqrt(compute_second_invariant(self.deviatoric_stress_new[mask, :]))
         yield_stress = self.yield_stress.new_value[mask]
         shear_modulus = self.shear_modulus.new_value[mask]
@@ -727,7 +736,8 @@ def apply_plasticity(self, mask, delta_t):
         self._equivalent_plastic_strain_rate[mask] = self._compute_equivalent_plastic_strain_rate(
             invariant_j2_el, shear_modulus, yield_stress, delta_t)
         self._deviatoric_stress_new[mask] *= radial_return[np.newaxis].T
-
+        self._compute_equivalent_plastic_strain(mask,delta_t)
+        
     def impose_pressure(self, ind_cell: int, pressure: float):
         """
         Pressure imposition
@@ -745,6 +755,7 @@ def increment_variables(self):
         """
         super().increment_variables()
         self._deviatoric_stress_current[:, :] = self._deviatoric_stress_new[:, :]
+        self.equivalent_plastic_strain.current_value = self.equivalent_plastic_strain.new_value
 
     def compute_new_coordinates(self, topology, x_coord):
         """
diff --git a/xfv/src/cell/one_dimension_enriched_cell_hansbo.py b/xfv/src/cell/one_dimension_enriched_cell_hansbo.py
old mode 100755
new mode 100644
index 01239e6..e62e4b0
--- a/xfv/src/cell/one_dimension_enriched_cell_hansbo.py
+++ b/xfv/src/cell/one_dimension_enriched_cell_hansbo.py
@@ -90,11 +90,15 @@ def __init__(self, n_cells: int):
         self._enr_deviatoric_stress_new = np.zeros([n_cells, 3])
         self._enr_deviatoric_strain_rate = np.zeros([n_cells, 3])
         self._enr_equivalent_plastic_strain_rate = np.zeros([n_cells])
+        self._enr_equivalent_plastic_strain = Field(n_cells, 0, 0)
         self._enr_plastic_strain_rate = np.zeros([n_cells, 3])
 
         # Indicator right part of enr cell is plastic
         self.plastic_enr_cells = np.zeros([n_cells], dtype=bool)
 
+        # Porosity
+        self._enr_porosity = Field(n_cells, np.zeros([n_cells]), np.zeros([n_cells]))
+
     def initialize_additional_cell_dof(self, disc: Discontinuity):
         """
         Values to initialize the right part fields when discontinuity disc is created
@@ -120,6 +124,9 @@ def initialize_additional_cell_dof(self, disc: Discontinuity):
         self.enr_yield_stress.current_value[enr_cell] = \
             np.copy(self.yield_stress.current_value[enr_cell])
         self._enr_coordinates_x[enr_cell] = np.copy(self._coordinates_x[enr_cell])
+        self._enr_equivalent_plastic_strain.current_value[enr_cell] = \
+            np.copy(self.equivalent_plastic_strain.current_value[enr_cell])
+        self._enr_porosity.current_value[enr_cell] = np.copy(self.porosity.current_value[enr_cell])
 
         # Initialization of new value field
         # (so that the current value is not erased if the field is not updated in current step)
@@ -137,6 +144,10 @@ def initialize_additional_cell_dof(self, disc: Discontinuity):
             np.copy(self.yield_stress.new_value[enr_cell])
         self._enr_deviatoric_stress_new[enr_cell] = \
             np.copy(self._deviatoric_stress_new[enr_cell])
+        self._enr_equivalent_plastic_strain.new_value[enr_cell] = \
+            np.copy(self.equivalent_plastic_strain.new_value[enr_cell])
+        self._enr_porosity.new_value[enr_cell] = np.copy(self.porosity.new_value[enr_cell])
+
         # Other quantities initialization
         self._enr_deviatoric_strain_rate[enr_cell] = \
             np.copy(self._deviatoric_strain_rate[enr_cell])
@@ -294,6 +305,13 @@ def enr_equivalent_plastic_strain_rate(self):
         """
         return self._enr_equivalent_plastic_strain_rate
 
+    @property
+    def enr_equivalent_plastic_strain(self):
+        """
+        Accessor on the right part of cracked cell equivalent plastic strain at time t
+        """
+        return self._enr_equivalent_plastic_strain
+
     @property
     def enr_plastic_strain_rate(self):
         """
@@ -301,6 +319,13 @@ def enr_plastic_strain_rate(self):
         """
         return self._enr_plastic_strain_rate
 
+    @property
+    def enr_porosity(self):
+        """
+        Accessor on the right part of cracked cell porosity at time t
+        """
+        return self._enr_porosity
+
     @property
     def pressure_field(self):
         """
@@ -455,6 +480,20 @@ def print_infos(self):
                 format(self.enr_artificial_viscosity.current_value[cell_i])
         print(message)
 
+    def compute_enriched_elements_new_porosity(self, delta_t, porosity_model):
+        """
+        Compute the new porosity according to the porosity model in XDATA
+
+        :param delta_t: model to compute the shear modulus
+        :param porosity_model: porosity model to compute
+        """
+        mask = self.enriched
+
+        self._enr_porosity.new_value[mask] = porosity_model.compute_porosity(
+            delta_t,
+            self.enr_porosity.current_value[mask],
+            self.enr_pressure.current_value[mask])
+
     def compute_enriched_elements_new_pressure(self, delta_t):
         """
         Compute pressure, internal energy and sound velocity in left and right parts of
@@ -469,7 +508,18 @@ def compute_enriched_elements_new_pressure(self, delta_t):
         plasticity_activated = (target_model.plasticity_model is not None)
 
         mask = self.enriched
+
+        if self.data.material_target.porosity_model is not None:
+                (self.enr_density.current_value[mask], self.enr_density.new_value[mask],
+                 self.enr_pressure.current_value[mask]) = self._compute_macro_to_micro_hydro(
+                     self.enr_density.current_value[mask],
+                     self.enr_density.new_value[mask],
+                     self.enr_pressure.current_value[mask],
+                     self.enr_porosity.new_value[mask])
+
         if elasticity_activated or plasticity_activated:
+            # Not sure there is cells in the intersection mask_classic / target
+            # => test it before starting Newton
             self.enr_energy.current_value[mask] += \
                 OneDimensionCell.add_elastic_energy_method(
                     delta_t, self.enr_density.current_value[mask],
@@ -495,6 +545,17 @@ def compute_enriched_elements_new_pressure(self, delta_t):
                 pressure_right_new, energy_right, energy_right_new,
                 pseudo_right, cson_right_new)
 
+        if self.data.material_target.porosity_model is not None:
+            self.enr_density.current_value[mask],\
+            self.enr_density.new_value[mask],\
+            self.enr_pressure.current_value[mask],\
+            pressure_new_right_value = self._compute_micro_to_macro_hydro(
+                density_right,
+                density_right_new,
+                pressure_right,
+                pressure_new_right_value,
+                self.enr_porosity.new_value[mask])
+
         # Save results :
         self.enr_pressure.new_value[mask] = pressure_new_right_value
         self.enr_energy.new_value[mask] = energy_new_right_value
@@ -681,7 +742,14 @@ def compute_enriched_shear_modulus(self, shear_modulus_model):
         if not mask.any():
             return
         self.enr_shear_modulus.new_value[mask] = \
-            shear_modulus_model.compute(self.enr_density.new_value[mask])
+            shear_modulus_model.compute(self.enr_density.new_value[mask],self.enr_pressure.current_value[mask])
+
+        if self.data.material_target.porosity_model is not None:
+            self._enr_shear_modulus.new_value[mask] = self._compute_micro_to_macro_shear_modulus(
+                self.enr_shear_modulus.new_value[mask],
+                self.enr_density.new_value[mask],
+                self.enr_sound_velocity.current_value[mask],
+                self.enr_porosity.new_value[mask])
 
     def apply_plasticity_enr(self, mask_mesh, delta_t):
         """
@@ -707,6 +775,16 @@ def apply_plasticity_enr(self, mask_mesh, delta_t):
             self._compute_equivalent_plastic_strain_rate(invariant_j2_el, shear_modulus,
                                                          yield_stress, delta_t)
         self._enr_deviatoric_stress_new[mask] *= radial_return[np.newaxis].T
+        self._enr_equivalent_plastic_strain.new_value[:] = \
+            self._compute_enr_equivalent_plastic_strain(mask,delta_t)
+
+    def _compute_enr_equivalent_plastic_strain(self,mask,delta_t):
+        """
+        Compute the equivalent plastic strain
+        param t_time : temps t
+        """
+        enr_equivalent_plastic_strain = delta_t* self._enr_equivalent_plastic_strain_rate[:] + self._enr_equivalent_plastic_strain.current_value[:]
+        return enr_equivalent_plastic_strain
 
     def compute_enriched_yield_stress(self, yield_stress_model):
         """
@@ -716,7 +794,13 @@ def compute_enriched_yield_stress(self, yield_stress_model):
         """
         mask = self.enriched
         self.enr_yield_stress.new_value[mask] = \
-            yield_stress_model.compute(self.enr_density.new_value[mask])
+            yield_stress_model.compute(self.enr_density.new_value[mask],
+                                       self._enr_equivalent_plastic_strain.current_value[mask],self.enr_shear_modulus.new_value[mask])
+
+        if self.data.material_target.porosity_model is not None:
+            self._enr_yield_stress.new_value[mask] = self._compute_micro_to_macro_yield_stress(
+                self.enr_yield_stress.new_value[mask],
+                self.enr_porosity.new_value[mask])
 
     def compute_enriched_elements_new_time_step(self):
         """
@@ -766,6 +850,13 @@ def cell_enr_increment(self):
         self._right_part_size.increment_values()
         # Elasticity
         self._enr_deviatoric_stress_current[:] = self._enr_deviatoric_stress_new[:]
+        # Plasticity
+        self._enr_equivalent_plastic_strain.increment_values()
+        # Rheology
+        self._enr_shear_modulus.increment_values()
+        self._enr_yield_stress.increment_values()
+        # Porosity
+        self._enr_porosity.increment_values()
 
     def compute_new_coordinates(self, topology, node_coord):
         """
diff --git a/xfv/src/cohesive_calculation/__init__.py b/xfv/src/cohesive_calculation/__init__.py
new file mode 100644
index 0000000..ff96786
--- /dev/null
+++ b/xfv/src/cohesive_calculation/__init__.py
@@ -0,0 +1,7 @@
+"""
+Package for cohesive module
+"""
+from xfv.src.cohesive_calculation.lineardata import LinearData
+from xfv.src.cohesive_calculation.linearenergy import LinearEnergy
+from xfv.src.cohesive_calculation.linearpercent import LinearPercent
+from xfv.src.cohesive_calculation.cohesivecalculationmodel import CohesiveCalculationModel
\ No newline at end of file
diff --git a/xfv/src/cohesive_calculation/cohesivecalculationmodel.py b/xfv/src/cohesive_calculation/cohesivecalculationmodel.py
new file mode 100644
index 0000000..b067094
--- /dev/null
+++ b/xfv/src/cohesive_calculation/cohesivecalculationmodel.py
@@ -0,0 +1,26 @@
+# -*- coding: utf-8 -*-
+"""
+Interface for the cohesive calculation model
+"""
+from abc import abstractmethod
+import numpy as np
+
+
+class CohesiveCalculationModel:  # pylint: disable=too-few-public-methods
+    """
+    Abstract class for the cohesive calculation model
+    """
+
+
+    @abstractmethod
+    def get_values(energy, stress, mass, section, ind, cohesive_model)-> float:
+        """
+        Compute and return critical strength and critical separation
+
+        :param energy: array for internal energy of cells [J/kg]
+        :param stress: array for stress of cells [Pa]
+        :param mass: array for mass of cells [kg]
+        :param section: float for section of material [m^2]
+        :param ind: integer for indice of cell
+        :cohesive_model: type of cohesive model 
+        """
diff --git a/xfv/src/cohesive_calculation/lineardata.py b/xfv/src/cohesive_calculation/lineardata.py
new file mode 100644
index 0000000..44b1484
--- /dev/null
+++ b/xfv/src/cohesive_calculation/lineardata.py
@@ -0,0 +1,31 @@
+# -*- coding: utf-8 -*-
+"""
+Implementation of the LinearData class
+"""
+
+import numpy as np
+from xfv.src.cohesive_calculation.cohesivecalculationmodel import CohesiveCalculationModel
+
+class LinearData(CohesiveCalculationModel):  # pylint: disable=too-few-public-methods
+    """
+    Class for cohesive linear data
+    """
+
+    def get_values(energy, stress, mass, section, ind, cohesive_model) -> float : 
+        """
+        Compute and return critical strength and critical separation
+
+        :param energy: array for internal energy of cells [J/kg]
+        :param stress: array for stress of cells [Pa]
+        :param mass: array for mass of cells [kg]
+        :param section: float for section of material [m^2]
+        :param ind: integer for indice of cell
+        :cohesive_model: type of cohesive model 
+        """
+        critical_strength = abs(stress[ind,0])
+        critical_separation = cohesive_model.critical_separation
+        dissipated_energy = critical_strength*critical_separation*section/2.
+        print('critical strength =', critical_strength)
+        print('critical separation =', critical_separation)
+        print('cohesive dissipated energy =', dissipated_energy)
+        return critical_strength, critical_separation
diff --git a/xfv/src/cohesive_calculation/linearenergy.py b/xfv/src/cohesive_calculation/linearenergy.py
new file mode 100644
index 0000000..fbd12f9
--- /dev/null
+++ b/xfv/src/cohesive_calculation/linearenergy.py
@@ -0,0 +1,33 @@
+# -*- coding: utf-8 -*-
+"""
+Implementation of the LinearEnergy class
+"""
+
+import numpy as np
+from xfv.src.cohesive_calculation.cohesivecalculationmodel import CohesiveCalculationModel
+
+class LinearEnergy(CohesiveCalculationModel) :  # pylint: disable=too-few-public-methods
+    """
+    Class for cohesive linear data
+    """
+
+    def get_values(energy, stress, mass, section, ind, cohesive_model) -> float :
+        """
+        Compute and return critical strength and critical separation
+
+        :param energy: array for internal energy of cells [J/kg]
+        :param stress: array for stress of cells [Pa]
+        :param mass: array for mass of cells [kg]
+        :param section: float for section of material [m^2]
+        :param ind: integer for indice of cell
+        :cohesive_model: type of cohesive model 
+        """
+
+        #if cohesive_model.dissipated_energy > energy[ind]*mass[ind]:
+        #    raise ValueError("L energie à dissiper par le modèle cohésif est superieure à l'energie interne de la cellule n°"+str(ind))
+        critical_strength = abs(stress[ind, 0])
+        critical_separation = 2.*cohesive_model.dissipated_energy/(critical_strength*section)
+        print('critical strength =', critical_strength)
+        print('critical separation =', critical_separation)
+        print('cohesive dissipated energy =', cohesive_model.dissipated_energy)
+        return critical_strength, critical_separation
diff --git a/xfv/src/cohesive_calculation/linearpercent.py b/xfv/src/cohesive_calculation/linearpercent.py
new file mode 100644
index 0000000..0d5c41e
--- /dev/null
+++ b/xfv/src/cohesive_calculation/linearpercent.py
@@ -0,0 +1,33 @@
+# -*- coding: utf-8 -*-
+"""
+Implementation of the LinearPercent class
+"""
+
+from distutils.errors import DistutilsInternalError
+import numpy as np
+from xfv.src.cohesive_calculation.cohesivecalculationmodel import CohesiveCalculationModel
+
+class LinearPercent(CohesiveCalculationModel) :  # pylint: disable=too-few-public-methods
+    """
+    Class for cohesive linear Percent 
+    """
+
+    def get_values(energy, stress, mass, section, ind, cohesive_model) -> float :
+        """
+        Compute and return critical strength and critical separation
+
+        :param energy: array for internal energy of cells [J/kg]
+        :param stress: array for stress of cells [Pa]
+        :param mass: array for mass of cells [kg]
+        :param section: float for section of material [m^2]
+        :param ind: integer for indice of cell
+        :cohesive_model: type of cohesive model 
+        """
+
+        critical_strength = abs(stress[ind, 0])
+        dissipated_energy = energy[ind]*mass[ind]*cohesive_model.purcentage
+        critical_separation = 2.*dissipated_energy/(critical_strength*section)
+        print('critical strength =', critical_strength)
+        print('critical separation =', critical_separation)
+        print('cohesive dissipated energy =', dissipated_energy)
+        return critical_strength, critical_separation
diff --git a/xfv/src/cohesive_model/cohesive_zone_model.py b/xfv/src/cohesive_model/cohesive_zone_model.py
index 7fd1680..49bb2d9 100644
--- a/xfv/src/cohesive_model/cohesive_zone_model.py
+++ b/xfv/src/cohesive_model/cohesive_zone_model.py
@@ -12,17 +12,67 @@ class CohesiveZoneModel:
     """
     A class for the computation of the cohesive force
     """
-    def __init__(self, cohesive_law_points: np.array, unloading_model: UnloadingModelBase):
+    def __init__(self, cohesive_law_points: np.array, unloading_model: UnloadingModelBase, cohesive_zone_model_name, dissipated_energy, purcentage, cohesive_calculation_model):
         """
         Construction d'un modèle cohésif
 
         :param cohesive_law_points: array describing the stress - opening curve of the
         cohesive model
         :param unloading_model:
+        :param cohesive_zone_model_name: string for the name of calculation cohesive model
+        :param dissipated_energy: float for energy to dissipate by cohesive model if critical separation is reached (required for LinearPercent model)
+        :param purcentage: float for percentage of internal energy to dissipate (required for LinearPercent model)
+        :param cohesive_calculation_model: class for calculation of cohesive model
         """
         self._critical_separation = cohesive_law_points[-1, 0]
-        self._cohesive_law = CohesiveLaw(cohesive_law_points)
+        self._critical_strength = cohesive_law_points[0,-1]
         self._unloading_model = unloading_model
+        self._cohesive_zone_model_name = cohesive_zone_model_name
+        self._dissipated_energy = dissipated_energy
+        self._purcentage = purcentage
+        self._cohesive_calculation_model = cohesive_calculation_model
+
+    @property
+    def cohesive_zone_model_name(self):  # pylint: disable=invalid-name
+        """
+        Cohesive zone model name
+        """
+        return self._cohesive_zone_model_name
+
+    @property
+    def cohesive_calculation_model(self):  # pylint: disable=invalid-name
+        """
+        _cohesive calculation model
+        """
+        return self._cohesive_calculation_model
+
+    @property
+    def critical_separation(self):  # pylint: disable=invalid-name
+        """
+        critical separation
+        """
+        return self._critical_separation
+
+    @property
+    def critical_strength(self):  # pylint: disable=invalid-name
+        """
+        critical strength
+        """
+        return self._critical_strength
+
+    @property
+    def dissipated_energy(self):  # pylint: disable=invalid-name
+        """
+        dissipated energy
+        """
+        return self._dissipated_energy
+
+    @property
+    def purcentage(self):  # pylint: disable=invalid-name
+        """
+        percentage of internal energy to be dissipated
+        """
+        return self._purcentage
 
     def compute_cohesive_stress(self, disc):
         """
@@ -30,29 +80,34 @@ def compute_cohesive_stress(self, disc):
         current discontinuity opening
 
         :param disc: discontinuity
+
+	    Attention : formule de la dissipation d'energie du modèle cohesif seulement dans le cas de la décharge à zéro
         """
         cohesive_force = 0.
         new_opening = disc.discontinuity_opening.new_value[0]
-
-        if new_opening <= 0.:
+        if new_opening <= 0. and disc.history_max_opening < disc.critical_separation :
             cohesive_force = 0.
+            disc.dissipated_energy.new_value = disc.critical_strength*disc.history_max_opening/2.
 
-        elif 0. < new_opening < disc.history_max_opening:
+        elif 0. < new_opening < disc.history_max_opening and disc.history_max_opening < disc.critical_separation :
             cohesive_force = \
                 self._unloading_model.compute_unloading_reloading_condition(disc, new_opening)
+            disc.dissipated_energy.new_value = disc.critical_strength*disc.history_max_opening/2.
 
-        elif disc.history_max_opening <= new_opening < self._critical_separation:
-            cohesive_force = self._cohesive_law.compute_cohesive_force(new_opening)
+        elif disc.history_max_opening <= new_opening < disc.critical_separation :
+            cohesive_force = disc._cohesive_law.compute_cohesive_force(new_opening)
             # Update the discontinuity indicators
             disc.history_max_opening = max(abs(disc.history_max_opening), abs(new_opening))
             disc.history_min_cohesive_force = \
-                self._cohesive_law.compute_cohesive_force(disc.history_max_opening)
-            disc.damage_variable.new_value = new_opening / self._critical_separation
+                disc._cohesive_law.compute_cohesive_force(disc.history_max_opening)
+            disc.damage_variable.new_value = new_opening / disc.critical_separation
+            disc.dissipated_energy.new_value = disc.critical_strength*disc.history_max_opening/2.
 
-        elif new_opening >= self._critical_separation:
+        else :
             disc.damage_variable.new_value = 1.
             cohesive_force = 0.
             disc.history_max_opening = max(abs(disc.history_max_opening), abs(new_opening))
             disc.history_min_cohesive_force = 0.
-
+            disc.dissipated_energy.new_value = disc.critical_strength*disc.critical_separation/2.
+            
         return cohesive_force
diff --git a/xfv/src/data/cohesive_model_props.py b/xfv/src/data/cohesive_model_props.py
index df6f6c1..06a4e2e 100644
--- a/xfv/src/data/cohesive_model_props.py
+++ b/xfv/src/data/cohesive_model_props.py
@@ -9,13 +9,19 @@
 from xfv.src.data.type_checked_dataclass import TypeCheckedDataClass
 from xfv.src.data.unloading_model_props import UnloadingModelProps
 from xfv.src.cohesive_model.cohesive_zone_model import CohesiveZoneModel
-
+from xfv.src.cohesive_calculation.lineardata import LinearData
+from xfv.src.cohesive_calculation.linearenergy import LinearEnergy
+from xfv.src.cohesive_calculation.linearpercent import LinearPercent
+from xfv.src.cohesive_calculation.cohesivecalculationmodel import CohesiveCalculationModel
 
 @dataclass  # pylint: disable=missing-class-docstring
 class CohesiveZoneModelProps(TypeCheckedDataClass):
     cohesive_strength: float
     critical_opening: float
+    _cohesive_zone_model_name: str
     unloading_model: UnloadingModelProps
+    _dissipated_energy: float
+    _purcentage: float
     _cohesive_zone_model_class = CohesiveZoneModel
 
     def _build_cohesive_law(self):
@@ -24,21 +30,67 @@ def _build_cohesive_law(self):
         """
         raise NotImplementedError("This is an abstract method!")
 
+
+@dataclass  # pylint: disable=missing-class-docstring  
+class LinearDataCohesiveZoneModelProps(CohesiveZoneModelProps):
+
+    _cohesive_calculation_model = LinearData
+
+    def _build_cohesive_law(self):
+        """
+        Build and return the CohesiveLaw
+        """
+        return np.array([
+            [0, self.cohesive_strength],
+            [self.critical_opening, 0]])   
+
     def build_cohesive_model_obj(self):
         """
         Build and return the CohesiveModel object
         """
         return self._cohesive_zone_model_class(self._build_cohesive_law(),
-                                               self.unloading_model.build_unloading_model_obj())
+                                               self.unloading_model.build_unloading_model_obj(),
+                                               self._cohesive_zone_model_name, self._dissipated_energy,
+                                               self._purcentage, self._cohesive_calculation_model)
+
+@dataclass  # pylint: disable=missing-class-docstring
+class LinearPercentCohesiveZoneModelProps(CohesiveZoneModelProps):
+
+    _cohesive_calculation_model = LinearPercent
+
+    def _build_cohesive_law(self):
+        return np.array([
+            [0, self.cohesive_strength],
+            [self.critical_opening, 0]])
+
+    def build_cohesive_model_obj(self):
+        """
+        Build and return the CohesiveModel object
+        """
+        return self._cohesive_zone_model_class(self._build_cohesive_law(),
+                                               self.unloading_model.build_unloading_model_obj(),
+                                               self._cohesive_zone_model_name, self._dissipated_energy,
+                                               self._purcentage, self._cohesive_calculation_model)
 
 
 @dataclass  # pylint: disable=missing-class-docstring
-class LinearCohesiveZoneModelProps(CohesiveZoneModelProps):
+class LinearEnergyCohesiveZoneModelProps(CohesiveZoneModelProps):
+
+    _cohesive_calculation_model = LinearEnergy
 
     def _build_cohesive_law(self):
         return np.array([
             [0, self.cohesive_strength],
             [self.critical_opening, 0]])
+    
+    def build_cohesive_model_obj(self):
+        """
+        Build and return the CohesiveModel object
+        """
+        return self._cohesive_zone_model_class(self._build_cohesive_law(),
+                                               self.unloading_model.build_unloading_model_obj(),
+                                               self._cohesive_zone_model_name, self._dissipated_energy,
+                                               self._purcentage,self._cohesive_calculation_model)
 
 
 @dataclass  # pylint: disable=missing-class-docstring
@@ -59,7 +111,7 @@ class TrilinearCohesiveZoneModelProps(CohesiveZoneModelProps):
     stress_1: float
     opening_2: float
     stress_2: float
-
+  
     def _build_cohesive_law(self):
         return np.array([
             [0, self.cohesive_strength],
diff --git a/xfv/src/data/data_container.py b/xfv/src/data/data_container.py
old mode 100755
new mode 100644
index 603a855..b6ddb3a
--- a/xfv/src/data/data_container.py
+++ b/xfv/src/data/data_container.py
@@ -16,7 +16,9 @@
                                                        MarchTableFunctionProps,
                                                        SuccessiveRampFunctionProps)
 from xfv.src.data.cohesive_model_props import (CohesiveZoneModelProps,
-                                               LinearCohesiveZoneModelProps,
+                                               LinearDataCohesiveZoneModelProps,
+                                               LinearPercentCohesiveZoneModelProps,
+                                               LinearEnergyCohesiveZoneModelProps,
                                                BilinearCohesiveZoneModelProps,
                                                TrilinearCohesiveZoneModelProps)
 from xfv.src.data.unloading_model_props import (UnloadingModelProps,
@@ -25,11 +27,11 @@
 from xfv.src.data.contact_props import (ContactProps, PenaltyContactProps,
                                         LagrangianMultiplierProps)
 from xfv.src.data.equation_of_state_props import (EquationOfStateProps, MieGruneisenProps)
-from xfv.src.data.yield_stress_props import (YieldStressProps, ConstantYieldStressProps)
-from xfv.src.data.shear_modulus_props import (ShearModulusProps, ConstantShearModulusProps)
+from xfv.src.data.yield_stress_props import (YieldStressProps, ConstantYieldStressProps, SCGYieldStressProps)
+from xfv.src.data.shear_modulus_props import (ShearModulusProps, ConstantShearModulusProps, SCGShearModulusProps)
 from xfv.src.data.plasticity_criterion_props import (PlasticityCriterionProps,
                                                      VonMisesCriterionProps)
-from xfv.src.data.rupture_criterion_props import (RuptureCriterionProps,
+from xfv.src.data.rupture_criterion_props import (DoubleCriterionProps, RuptureCriterionProps,
                                                   DamageCriterionProps,
                                                   PorosityCriterionProps,
                                                   HalfRodComparisonCriterionProps,
@@ -99,7 +101,7 @@ def __post_init__(self):
 ALL_VARIABLES = ["NodeVelocity", "NodeCoordinates", "CellSize", "Pressure", "Density",
                  "InternalEnergy", "SoundVelocity", "ArtificialViscosity", "Stress",
                  "DeviatoricStress", "EquivalentPlasticStrainRate", "PlasticStrainRate",
-                 "Porosity", "CohesiveForce", "DiscontinuityOpening", "ShearModulus", "YieldStress"]
+                 "Porosity", "CohesiveForce", "DissipatedEnergy","DiscontinuityOpening", "ShearModulus", "YieldStress"]
 
 
 @dataclass  # pylint: disable=missing-class-docstring
@@ -385,8 +387,6 @@ def __get_cohesive_model_props(matter) -> Optional[CohesiveZoneModelProps]:
 
         cohesive_model_name = params['name'].lower()
 
-        cohesive_strength = params['coefficients']['cohesive-strength']
-        critical_separation = params['coefficients']['critical-separation']
 
         unloading_model_name = params['unloading-model']['name'].lower()
         unloading_model_slope: Optional[float] = params['unloading-model'].get('slope')
@@ -400,25 +400,61 @@ def __get_cohesive_model_props(matter) -> Optional[CohesiveZoneModelProps]:
             raise ValueError(f"Unknwown unloading model name: {unloading_model_name} "
                              "Please choose among (progressiveunloading, "
                              " lossofstiffnessunloading)")
-
-        if cohesive_model_name == "linear":
-            cohesive_model_props: CohesiveZoneModelProps = LinearCohesiveZoneModelProps(
-                cohesive_strength, critical_separation, unloading_model_props)
+        print('cohesive model is', cohesive_model_name)
+         
+        if cohesive_model_name == "lineardata":
+            cohesive_strength = params['coefficients']['cohesive-strength']
+            critical_separation = params['coefficients']['critical-separation']
+            purcentage_internal_energy = 0.
+            dissipated_energy = 0.
+            cohesive_model_props = LinearDataCohesiveZoneModelProps(cohesive_strength, critical_separation,
+                                                                cohesive_model_name, unloading_model_props, dissipated_energy, purcentage_internal_energy)
+        elif cohesive_model_name == "linearenergy":
+            cohesive_strength = 0.
+            critical_separation = 0.
+            dissipated_energy = params['coefficients']['dissipated-energy']
+            purcentage_internal_energy = 0.
+
+            # permet d'affecter la contrainte cohesive lorsque le critère de rupture est 'maximalstress'
+            failure_data = matter.get('failure')
+            failure_criterion_data = failure_data['failure-criterion']
+            fail_crit_name: str = failure_criterion_data['name']
+            fail_crit_value: Optional[float] = failure_criterion_data.get('value')
+            if fail_crit_name =='MaximalStress':
+                cohesive_strength = fail_crit_value
+            cohesive_model_props = LinearEnergyCohesiveZoneModelProps(cohesive_strength, critical_separation,
+                                                                cohesive_model_name, unloading_model_props,
+                                                                     dissipated_energy, purcentage_internal_energy)
+        elif cohesive_model_name == "linearpercent":
+            dissipated_energy = 0.
+            purcentage_internal_energy = params['coefficients']['purcentage-internal-energy']
+            cohesive_model_props = LinearPercentCohesiveZoneModelProps(0., 0.,
+                                                                cohesive_model_name, unloading_model_props,
+                                                                     dissipated_energy, purcentage_internal_energy)
         elif cohesive_model_name == "bilinear":
+            cohesive_strength = params['coefficients']['cohesive-strength']
+            critical_separation = params['coefficients']['critical-separation']
+            purcentage_internal_energy = 0.
+            dissipated_energy = 0.
             cohesive_model_props = BilinearCohesiveZoneModelProps(
-                cohesive_strength, critical_separation, unloading_model_props,
+                cohesive_strength, critical_separation, cohesive_model_name, unloading_model_props,
+                dissipated_energy, purcentage_internal_energy, purcentage_internal_energy, purcentage_internal_energy,
                 params['coefficients']['separation-at-point-1'],
                 params['coefficients']['stress-at-point-1'])
+            
         elif cohesive_model_name == "trilinear":
+            cohesive_strength = params['coefficients']['cohesive-strength']
+            critical_separation = params['coefficients']['critical-separation']
             cohesive_model_props = TrilinearCohesiveZoneModelProps(
-                cohesive_strength, critical_separation, unloading_model_props,
+                cohesive_strength, critical_separation, cohesive_model_name, unloading_model_props,
+                dissipated_energy, purcentage_internal_energy, purcentage_internal_energy, purcentage_internal_energy,
                 params['coefficients']['separation-at-point-1'],
                 params['coefficients']['stress-at-point-1'],
                 params['coefficients']['separation-at-point-2'],
                 params['coefficients']['stress-at-point-2'])
         else:
             raise ValueError(f"Unknwon cohesive model: {cohesive_model_name} ."
-                             "Please choose among (linear, bilinear, trilinear)")
+                             "Please choose among (LinearData, LinearPercent, LinearEnergy, bilinear, trilinear)")
 
         return cohesive_model_props
 
@@ -440,10 +476,24 @@ def __get_porosity_model_props(matter) -> Optional[Tuple[PorosityModelProps, str
                 initial_porosity_for_johnson = params['coefficients']['initial-porosity']
                 effective_strength_for_johnson = params['coefficients']['effective-strength']
                 viscosity_for_johnson = params['coefficients']['viscosity']
+                maximal_porosity_for_johnson = 1.e30
+
+                # Lorsque le critère de rupture est la porosité
+                failure_data = matter.get('failure')
+                failure_criterion_data = failure_data['failure-criterion']
+                fail_crit_name: str = failure_criterion_data['name']
+                fail_crit_value: Optional[float] = failure_criterion_data.get('value')
+                if fail_crit_name == 'Porosity':
+                    maximal_porosity_for_johnson = fail_crit_value
+                    print('maximal porosity = ', maximal_porosity_for_johnson)
+                elif fail_crit_name == 'DoubleCriterion':
+                    maximal_porosity_for_johnson = failure_criterion_data.get('value-one')
+                    print('maximal porosity = ', maximal_porosity_for_johnson)
+
                 porosity_model_props: PorosityModelProps = JohnsonModelProps(
                     initial_porosity_for_johnson,
                     effective_strength_for_johnson,
-                    viscosity_for_johnson)
+                    viscosity_for_johnson, maximal_porosity_for_johnson)
         except:
             raise ValueError(f"No keyword 'name' for porosity model name: {porosity_model_name}."
                              "Please choose among (JohnsonModel)."
@@ -510,6 +560,8 @@ def __fill_in_material_props(self, material) -> Tuple[InitialValues,
                     isinstance(failure_criterion, MaximalStressCriterionProps)):
                 print("Failure criterion value and cohesive strength have different value. "
                       "This may result in errors in the future")
+                print('failure criterion =', failure_criterion)
+                print('cohesive strength =', cohesive_model.cohesive_strength)
 
         # Porosity model
         porosity_model_props = self.__get_porosity_model_props(material)
@@ -553,7 +605,7 @@ def __get_initial_values(self, matter) -> Tuple[float, float, float, float,
         return (velocity, pressure, temperature, density, internal_energy,
                 yield_stress, shear_modulus, porosity)
 
-    def __get_yield_stress_and_shear_modulus(self, matter) -> Tuple[float, float]:
+    def __get_yield_stress_and_shear_modulus(self, matter) -> Tuple[float, float,Optional[float],Optional[float],Optional[float],Optional[float]]:
         """
         Returns the yield stress and shear modulus read from the json file specified
         under the coefficients key
@@ -568,8 +620,8 @@ def __get_yield_stress_and_shear_modulus(self, matter) -> Tuple[float, float]:
 
         with open(self._datafile_dir / coefficients_file, 'r') as json_fid:
             coef = json.load(json_fid)
-            yield_stress = float(coef['EPP']["yield_stress"])
-            shear_modulus = float(coef['EPP']["shear_modulus"])
+            yield_stress = float(coef[params.get('plasticity-model')]["yield_stress"])
+            shear_modulus = float(coef[params.get('plasticity-model')]["shear_modulus"])
         return yield_stress, shear_modulus
 
     def __get_initial_porosity(self, matter)->float:
@@ -612,23 +664,47 @@ def __get_rheology_props(self, matter) -> Tuple[Optional[ShearModulusProps],
                  plasticity_model : plasticity model
                  plasticity_criterion : plasticity criterion
         """
+        yield_stress, shear_modulus = self.__get_yield_stress_and_shear_modulus(matter)
+        params = matter['initialization']
+        json_path = params['init-thermo']
+        
+        with open(self._datafile_dir / json_path, 'r') as json_fid:
+            coef = json.load(json_fid)
+            coef = coef["InitThermo"]
+            rho_0 = float(coef["initial_density"])
         params = matter.get('rheology')
-        if not params:
+        if params is None:
             return None, None, None
+        coefficients_file = params.get('coefficients')
+        with open(self._datafile_dir / coefficients_file, 'r') as json_fid:
+            coef = json.load(json_fid)
 
-        yield_stress, shear_modulus = self.__get_yield_stress_and_shear_modulus(matter)
         elasticity_model_name = params['elasticity-model']
-        if elasticity_model_name != "Linear":
-            raise ValueError(f"Model {elasticity_model_name} not implemented. Choose Linear")
-        elasticity_model = ConstantShearModulusProps(shear_modulus)
+        if elasticity_model_name == "Linear":
+            Gp_prime = 0.
+            elasticity_model = ConstantShearModulusProps(shear_modulus, rho_0, Gp_prime)
+        elif elasticity_model_name == "SCG":
+            Gp_prime = float(coef[params.get('elasticity-model')]["Gp_prime"])
+            elasticity_model = SCGShearModulusProps(shear_modulus, rho_0, Gp_prime)
+        else:
+            raise ValueError(f"Model {elasticity_model_name} not implemented. Choose Linear or SCG")
 
         plasticity_model_name = params.get('plasticity-model')
         if not plasticity_model_name:
             return elasticity_model, None, None
 
-        if plasticity_model_name != "EPP":
-            raise ValueError(f"Model {plasticity_model_name} not implemented. Choose EPP")
-        plasticity_model = ConstantYieldStressProps(yield_stress)
+        if plasticity_model_name == "EPP" :
+            beta = 0.
+            m = 0.
+            Y_max = 0.
+            plasticity_model = ConstantYieldStressProps(yield_stress, shear_modulus, Y_max, beta, m)
+        elif plasticity_model_name == "SCG":
+            beta = float(coef[params.get('plasticity-model')]["beta_ecrouissage"])
+            m = float(coef[params.get('plasticity-model')]["m_ecrouissage"])
+            Y_max = float(coef[params.get('plasticity-model')]["yield_stress_max"])
+            plasticity_model = SCGYieldStressProps(yield_stress, shear_modulus, Y_max, beta, m)
+        else:
+            raise ValueError(f"Model {plasticity_model_name} not implemented. Choose EPP or SCG")
 
         plasticity_criterion_name = params['plasticity-criterion']
         if plasticity_criterion_name == "VonMises":
@@ -714,6 +790,10 @@ def __get_failure_criterion_props(matter) -> Tuple[
         elif fail_crit_name == "NonLocalStress":
             radius: Optional[int] = failure_criterion_data.get('radius')
             failure_criterion = NonLocalStressCriterionProps(fail_crit_value, radius)
+        elif fail_crit_name == 'DoubleCriterion':
+            value_one : Optional[float] = failure_criterion_data.get('value-one')
+            value_two : Optional[int] = failure_criterion_data.get('value-two')
+            failure_criterion = DoubleCriterionProps(value_one, value_two)
         else:
             raise ValueError(f"Unknown failure criterion {fail_crit_name}. "
                              "Please choose among (MinimumPressure, Damage, "
diff --git a/xfv/src/data/porosity_model_props.py b/xfv/src/data/porosity_model_props.py
index ea46d5e..8254d2c 100644
--- a/xfv/src/data/porosity_model_props.py
+++ b/xfv/src/data/porosity_model_props.py
@@ -44,6 +44,7 @@ class JohnsonModelProps(PorosityModelProps):
     initial_porosity_for_johnson: float
     effective_strength_for_johnson: float
     viscosity_for_johnson: float  # pylint: disable=invalid-name
+    maximal_porosity_for_Johnson: float
     _porosity_model_class = JohnsonModel
 
     def __post_init__(self):
diff --git a/xfv/src/data/rupture_criterion_props.py b/xfv/src/data/rupture_criterion_props.py
index c890520..694767e 100644
--- a/xfv/src/data/rupture_criterion_props.py
+++ b/xfv/src/data/rupture_criterion_props.py
@@ -13,7 +13,7 @@
 from xfv.src.rupturecriterion.maximalstress import MaximalStressCriterion
 from xfv.src.rupturecriterion.minimumpressure import MinimumPressureCriterion
 from xfv.src.rupturecriterion.nonlocalstress import NonLocalStressCriterion
-
+from xfv.src.rupturecriterion.doublecriterion import DoubleCriterion
 
 @dataclass  # pylint: disable=missing-class-docstring
 class RuptureCriterionProps(TypeCheckedDataClass):
@@ -101,4 +101,14 @@ class NonLocalStressCriterionProps(RuptureCriterionProps):
     def __post_init__(self):
         super().__post_init__()  # typecheck first
         self._ensure_defined('value', 'NonLocalStressCriterionProps',
-                             'failure/failure-criterion/value')  # ensures that exists
\ No newline at end of file
+                             'failure/failure-criterion/value')  # ensures that exists
+
+@dataclass  # pylint: disable=missing-class-docstring
+class DoubleCriterionProps(RuptureCriterionProps):
+    p_limit: Optional[float]
+    number_cell: Optional[int]  # optional to personalize the error message
+
+    _rupture_criterion_class = DoubleCriterion
+
+    def __post_init__(self):
+        super().__post_init__()  # typecheck first
diff --git a/xfv/src/data/shear_modulus_props.py b/xfv/src/data/shear_modulus_props.py
index d35951f..160cadd 100644
--- a/xfv/src/data/shear_modulus_props.py
+++ b/xfv/src/data/shear_modulus_props.py
@@ -8,6 +8,7 @@
 from xfv.src.data.type_checked_dataclass import TypeCheckedDataClass
 from xfv.src.rheology.shearmodulus import ShearModulus
 from xfv.src.rheology.constantshearmodulus import ConstantShearModulus
+from xfv.src.rheology.scgshearmodulus import SCGShearModulus
 
 @dataclass  # pylint: disable=missing-class-docstring
 class ShearModulusProps(TypeCheckedDataClass):
@@ -34,9 +35,14 @@ def build_shear_modulus_obj(self):
 @dataclass  # pylint: disable=missing-class-docstring, too-many-instance-attributes
 class ConstantShearModulusProps(ShearModulusProps):
     init_value: float
+    init_density: float
+    Gp_prime: float
     _shear_modulus_class = ConstantShearModulus
 
     def __post_init__(self):
         super().__post_init__()  # typecheck first
         self._ensure_positivity('init_value')
 
+
+class SCGShearModulusProps(ConstantShearModulusProps):
+    _shear_modulus_class = SCGShearModulus
diff --git a/xfv/src/data/yield_stress_props.py b/xfv/src/data/yield_stress_props.py
index 87217a4..ec7363e 100644
--- a/xfv/src/data/yield_stress_props.py
+++ b/xfv/src/data/yield_stress_props.py
@@ -8,6 +8,7 @@
 from xfv.src.data.type_checked_dataclass import TypeCheckedDataClass
 from xfv.src.rheology.yieldstress import YieldStress
 from xfv.src.rheology.constantyieldstress import ConstantYieldStress
+from xfv.src.rheology.scgyieldstress import SCGYieldStress
 
 
 @dataclass  # pylint: disable=missing-class-docstring
@@ -35,4 +36,11 @@ def build_yield_stress_obj(self):
 @dataclass  # pylint: disable=missing-class-docstring, too-many-instance-attributes
 class ConstantYieldStressProps(YieldStressProps):
     init_value: float
+    init_shear_modulus: float
+    Y_max: float
+    beta: float
+    m: float
     _yield_stress_class = ConstantYieldStress
+
+class SCGYieldStressProps(ConstantYieldStressProps):
+    _yield_stress_class = SCGYieldStress
diff --git a/xfv/src/discontinuity/discontinuity.py b/xfv/src/discontinuity/discontinuity.py
old mode 100755
new mode 100644
index 3d11332..69759d0
--- a/xfv/src/discontinuity/discontinuity.py
+++ b/xfv/src/discontinuity/discontinuity.py
@@ -5,6 +5,10 @@
 import numpy as np
 from xfv.src.fields.field import Field
 from xfv.src.data.enriched_mass_matrix_props import EnrichedMassMatrixProps
+from xfv.src.cohesive_calculation.lineardata import LinearData
+from xfv.src.cohesive_calculation.linearenergy import LinearEnergy
+from xfv.src.cohesive_calculation.linearpercent import LinearPercent
+from xfv.src.cohesive_model.cohesive_law import CohesiveLaw
 
 
 class Discontinuity:
@@ -28,6 +32,8 @@ class Discontinuity:
     in_nodes = np.zeros([], dtype=int)
     out_nodes = np.zeros([], dtype=int)
 
+
+
     def __init__(self, cell_id: int, mask_in_nodes: np.array, mask_out_nodes: np.array,
                  discontinuity_position_in_ruptured_element: float,
                  enriched_mass_matrix_props: EnrichedMassMatrixProps):
@@ -81,6 +87,7 @@ def __init__(self, cell_id: int, mask_in_nodes: np.array, mask_out_nodes: np.arr
                 Discontinuity.in_nodes, np.zeros((1, 1), dtype=int), axis=0)
             Discontinuity.out_nodes = np.append(
                 Discontinuity.out_nodes, np.zeros((1, 1), dtype=int), axis=0)
+
         for ind, disc in enumerate(Discontinuity.discontinuity_list()):
             disc.enr_velocity_current = Discontinuity.enr_velocity_current[ind]
             disc.enr_velocity_new = Discontinuity.enr_velocity_new[ind]
@@ -91,6 +98,8 @@ def __init__(self, cell_id: int, mask_in_nodes: np.array, mask_out_nodes: np.arr
             disc.ruptured_cell_id = Discontinuity.ruptured_cell_id[ind]
             disc.in_nodes = Discontinuity.in_nodes[ind]
             disc.out_nodes = Discontinuity.out_nodes[ind]
+            
+
 
         self.__mask_in_nodes = mask_in_nodes
         self.in_nodes[:] = np.where(self.__mask_in_nodes)[0]
@@ -111,9 +120,13 @@ def __init__(self, cell_id: int, mask_in_nodes: np.array, mask_out_nodes: np.arr
         # (Always created but null if no damage data in the XDATA.json file...)
         self.cohesive_force = Field(1, current_value=0., new_value=0.)
         self.discontinuity_opening = Field(1, current_value=0., new_value=0.)
+        self.dissipated_energy = Field(1, current_value=0., new_value=0.)
         self.damage_variable = Field(1, current_value=0., new_value=0.)
         self.history_max_opening = 0.
         self.history_min_cohesive_force = 1.e+30
+        self._cohesive_law = []
+        self.critical_separation = None
+        self.critical_strength = None
 
         # Creation of the enriched mass matrix
         self.mass_matrix_enriched = enriched_mass_matrix_props.build_enriched_mass_matrix_obj()
@@ -231,6 +244,23 @@ def compute_discontinuity_new_opening(self, node_position: np.array):
         xd_new = (1 - epsilon) * enr_coord_d + epsilon * coord_d
         self.discontinuity_opening.new_value = (xd_new - xg_new)[0][0]
 
+    def create_cohesive_law(self, cells, section, cohesive_model):
+        if cohesive_model.cohesive_zone_model_name is None :
+            return
+        ind = self.ruptured_cell_id
+        self.critical_strength, self.critical_separation = cohesive_model.cohesive_calculation_model.get_values(cells.energy.new_value, cells.stress, cells.mass, section, ind, cohesive_model)
+        points = np.array([
+            [0, self.critical_strength],
+            [self.critical_separation, 0]])
+        self._cohesive_law = CohesiveLaw(points)
+
+    def compute_critical_value(self, stress, energy):
+        for ind in range(len(Discontinuity.critical_strength)):
+            if abs(Discontinuity.critical_strength[ind]) < 1.e-16:
+                Discontinuity.critical_strength[ind] = abs(stress[Discontinuity.ruptured_cell_id[ind]])
+            if abs(Discontinuity.critical_separation[ind]) < 1.e-16:
+                Discontinuity.critical_separation[ind] = 2*energy[Discontinuity.ruptured_cell_id[ind]]/Discontinuity.critical_strength[ind]
+        
     def reinitialize_kinematics_after_contact(self):
         """
         Set the new velocity to the old one to cancel the increment that has lead to contact
@@ -245,7 +275,9 @@ def enr_increment(self):
         # Kinematics
         self.enr_velocity_current[:] = self.enr_velocity_new[:]
         self.enr_coordinates_current[:] = self.enr_coordinates_new[:]
+        
         # Cohesive model
         self.cohesive_force.increment_values()
         self.damage_variable.increment_values()
         self.discontinuity_opening.increment_values()
+        self.dissipated_energy.increment_values()
diff --git a/xfv/src/mesh/mesh1denriched.py b/xfv/src/mesh/mesh1denriched.py
old mode 100755
new mode 100644
index 4e4e478..7b1c6e8
--- a/xfv/src/mesh/mesh1denriched.py
+++ b/xfv/src/mesh/mesh1denriched.py
@@ -299,6 +299,10 @@ def compute_new_cells_porosity(self, delta_t: float, porosity_model):
         :param porosity_model: model to compute the porosity model
         """
         self.cells.compute_new_porosity(delta_t, porosity_model, self.cells.classical)
+        # right part of enriched cells
+        if self.cells.enriched.any():
+            # test if any enriched cell in order to avoid error in the Newton initialization
+            self.cells.compute_enriched_elements_new_porosity(delta_t, porosity_model)
 
     def compute_new_cells_pseudo_viscosity(self, delta_t: float):
         """
@@ -322,7 +326,9 @@ def compute_new_cohesive_forces(self):
         Computation of cohesive forces at t+dt
         """
         if self.cohesive_zone_model is not None:
-            self.nodes.compute_enriched_nodes_cohesive_forces(self.cohesive_zone_model)
+            self.nodes.compute_enriched_nodes_cohesive_forces(self.cohesive_zone_model,
+                                                              self.cells.stress[:,0],
+                                                self.cells.energy.new_value)
 
     def increment(self):
         """
@@ -453,7 +459,8 @@ def apply_rupture_treatment(self, treatment, time: float):
         :type treatment: RuptureTreatment
         """
         treatment.apply_treatment(self.cells, self.__ruptured_cells,
-                                  self.nodes, self.__topology, time)
+                                  self.nodes, self.__topology, time,
+                                   self.cohesive_zone_model, self.nodes.section)
 
     def apply_plasticity(self, delta_t: float, yield_stress_model, plasticity_criterion,
                          mask_mesh: np.array):
@@ -473,7 +480,7 @@ def apply_plasticity(self, delta_t: float, yield_stress_model, plasticity_criter
         # la variable dev_stress_new et doit donc être appelée à la fin de
         # l'étape du calcul de plasticité pour conserver la prédiction élastique dans
         # le calcul du taux de déformation plastique, plasticité cumulée, ...
-
+        
         # 1) Compute yield stress
         self.cells.compute_yield_stress(yield_stress_model, mask_mesh)
         if mask_mesh.any() and self.cells.enriched.any():
@@ -592,6 +599,13 @@ def stress_xx_field(self) -> np.array:
         """
         return self.cells.stress_xx_field
 
+    @property
+    def equivalent_plastic_strain_field(self) -> np.array:
+        """
+        equivalent plastic strain field
+        """
+        return self.cells.equivalent_plastic_strain_field
+
     @property
     def porosity_field(self) -> np.array:
         """
diff --git a/xfv/src/node/one_dimension_enriched_node_hansbo.py b/xfv/src/node/one_dimension_enriched_node_hansbo.py
old mode 100755
new mode 100644
index e51deaf..4f64fa5
--- a/xfv/src/node/one_dimension_enriched_node_hansbo.py
+++ b/xfv/src/node/one_dimension_enriched_node_hansbo.py
@@ -192,7 +192,7 @@ def compute_enriched_nodes_new_force(self, contrainte_xx: np.array, enr_contrain
             self._force[disc.mask_in_nodes] += f_node_left_minus * self.section  # F1-
             self._force[disc.mask_out_nodes] += f_node_right_plus * self.section  # F2+
 
-    def compute_enriched_nodes_cohesive_forces(self, cohesive_model):
+    def compute_enriched_nodes_cohesive_forces(self, cohesive_model, stress, energy):
         """
         Compute the cohesive forces for the enriched nodes
 
@@ -204,10 +204,10 @@ def compute_enriched_nodes_cohesive_forces(self, cohesive_model):
         nb_disc = len(disc_list)
 
         applied_force_arr = np.ndarray((nb_disc,))
-        for ind, disc in enumerate(disc_list):
-            # Compute cohesive stress
+        for ind, disc in enumerate(disc_list):  
             cohesive_stress = cohesive_model.compute_cohesive_stress(disc)
             disc.cohesive_force.new_value = cohesive_stress
+            disc.dissipated_energy.new_value *= self.section
             applied_force_arr[ind] = self.section * cohesive_stress
 
         self.apply_force_on_discontinuity_boundaries_arr(applied_force_arr)
diff --git a/xfv/src/node/one_dimension_node.py b/xfv/src/node/one_dimension_node.py
old mode 100755
new mode 100644
index 3162c46..97d3091
--- a/xfv/src/node/one_dimension_node.py
+++ b/xfv/src/node/one_dimension_node.py
@@ -78,7 +78,7 @@ def compute_new_force(self, topologie, contrainte, classical_cell_mask: np.array
         Calcul des forces agissant sur les noeuds
 
         :param topologie: topologie du calcul
-        :param contrainte: tenseur des contriante de cauchy sigma xx
+        :param contrainte: tenseur des contrainte de cauchy sigma xx
         :param classical_cell_mask: masks of the classical cells
         :type topologie: Topology
         :type contrainte: numpy.array([nbr_of_node-1, 1], dtype=np.float64, order='C')
diff --git a/xfv/src/output_manager/outputdatabaseexploit.py b/xfv/src/output_manager/outputdatabaseexploit.py
old mode 100755
new mode 100644
index c993ab5..b54e3a7
--- a/xfv/src/output_manager/outputdatabaseexploit.py
+++ b/xfv/src/output_manager/outputdatabaseexploit.py
@@ -30,10 +30,11 @@ class OutputDatabaseExploit:
     field_type_converter["NodeVelocity"] = ("ClassicalNodeVelocity", "AdditionalNodeVelocity")
     field_type_converter["NodeCoordinates"] = ("NodeCoordinates", "None")
     field_type_converter["CohesiveForce"] = ("AdditionalCohesiveForce", "None")
-    field_type_converter["Porosity"] = ("ClassicalPorosity", "None")
-    field_type_converter["ShearModulus"] = ("ClassicalShearModulus", "None")
-    field_type_converter["YieldStress"] = ("ClassicalYieldStress", "None")
+    field_type_converter["Porosity"] = ("ClassicalPorosity", "AdditionalPorosity")
+    field_type_converter["ShearModulus"] = ("ClassicalShearModulus", "AdditionalShearModulus")
+    field_type_converter["YieldStress"] = ("ClassicalYieldStress", "AdditionalYieldStress")
     field_type_converter["DiscontinuityOpening"] = ("AdditionalDiscontinuityOpening", "None")
+    field_type_converter["DissipatedEnergy"] = ("AdditionalDissipatedEnergy", "None")
     field_type_converter["NodeStatus"] = ("NodeStatus", "None")
     field_type_converter["CellStatus"] = ("CellStatus", "None")
 
diff --git a/xfv/src/output_manager/outputmanager.py b/xfv/src/output_manager/outputmanager.py
old mode 100755
new mode 100644
index 68a576b..5212c4d
--- a/xfv/src/output_manager/outputmanager.py
+++ b/xfv/src/output_manager/outputmanager.py
@@ -53,6 +53,9 @@
 enr_field_list["DiscontinuityOpening"] = FieldConstruction("AdditionalDiscontinuityOpening",
                                                            None, ("discontinuity_opening",
                                                                   "current_value"))
+enr_field_list["DissipatedEnergy"] = FieldConstruction("AdditionalDissipatedEnergy",
+                                                           None, ("dissipated_energy",
+                                                                  "current_value"))
 enr_field_list["Pressure"] = FieldConstruction("AdditionalPressure", "cells",
                                                ("enr_pressure", "current_value"))
 enr_field_list["Density"] = FieldConstruction("AdditionalDensity", "cells",
@@ -71,7 +74,12 @@
     "AdditionalEquivalentPlasticStrainRate", "cells", ("enr_equivalent_plastic_strain_rate", ))
 enr_field_list["PlasticStrainRate"] = FieldConstruction(
     "AdditionalPlasticStrainRate", "cells", ("enr_plastic_strain_rate", ))
-
+enr_field_list["Porosity"] = FieldConstruction("AdditionalPorosity", "cells",
+                                               ("enr_porosity", "current_value"))
+field_list["ShearModulus"] = FieldConstruction("AdditionalShearModulus",
+                                               "cells", ("enr_shear_modulus", "current_value"))
+field_list["YieldStress"] = FieldConstruction("AdditionalYieldStress",
+                                              "cells", ("enr_yield_stress", "current_value"))
 
 class OutputManager(metaclass=Singleton):
     """
diff --git a/xfv/src/porosity_model/johnsonmodel.py b/xfv/src/porosity_model/johnsonmodel.py
index 2a0bd86..a2d8fb2 100644
--- a/xfv/src/porosity_model/johnsonmodel.py
+++ b/xfv/src/porosity_model/johnsonmodel.py
@@ -12,10 +12,11 @@ class JohnsonModel(PorosityModelBase):  # pylint: disable=too-few-public-methods
     def __init__(self,
                  initial_porosity_for_johnson,
                  effective_strength_for_johnson,
-                 viscosity_for_johnson):
+                 viscosity_for_johnson, maximal_porosity_for_Johnson):
         self.initial_porosity_for_johnson = initial_porosity_for_johnson
         self.effective_strength_for_johnson = effective_strength_for_johnson
         self.viscosity_for_johnson = viscosity_for_johnson
+        self.maximal_porosity_for_Johnson = maximal_porosity_for_Johnson
 
     def compute_porosity(self, delta_t: float,
                          porosity: np.array,
@@ -52,6 +53,7 @@ def _compute_johnson_porosity(self, delta_t: float,
         Compute the new value of porosity
         """
         initial_porosity_for_johnson = self.initial_porosity_for_johnson
+        maximal_porosity_for_Johnson = self.maximal_porosity_for_Johnson
         viscosity_for_johnson = self.viscosity_for_johnson
 
         power = 1.0/3.0
@@ -62,5 +64,6 @@ def _compute_johnson_porosity(self, delta_t: float,
         dalphadt = dalphadt/viscosity_for_johnson
         porosity_new = porosity + dalphadt*delta_t
         porosity_new = np.maximum(porosity_new, initial_porosity_for_johnson)
-        porosity_new = np.where(porosity > 1., porosity_new, 1.)
+        porosity_new = np.minimum(porosity_new, maximal_porosity_for_Johnson)
+        porosity_new = np.where(porosity_new > 1., porosity_new, 1.)
         return porosity_new
diff --git a/xfv/src/rheology/__init__.py b/xfv/src/rheology/__init__.py
old mode 100755
new mode 100644
index 2623d93..18f16a2
--- a/xfv/src/rheology/__init__.py
+++ b/xfv/src/rheology/__init__.py
@@ -3,5 +3,7 @@
 """
 from xfv.src.rheology.shearmodulus import ShearModulus
 from xfv.src.rheology.constantshearmodulus import ConstantShearModulus
+from xfv.src.rheology.scgshearmodulus import SCGShearModulus
 from xfv.src.rheology.yieldstress import YieldStress
 from xfv.src.rheology.constantyieldstress import ConstantYieldStress
+from xfv.src.rheology.scgyieldstress import SCGYieldStress
diff --git a/xfv/src/rheology/constantshearmodulus.py b/xfv/src/rheology/constantshearmodulus.py
old mode 100755
new mode 100644
index 9cddf2b..df0b935
--- a/xfv/src/rheology/constantshearmodulus.py
+++ b/xfv/src/rheology/constantshearmodulus.py
@@ -12,20 +12,22 @@ class ConstantShearModulus(ShearModulus):  # pylint: disable=too-few-public-meth
     Class for constant shear modulus
     """
 
-    def __init__(self, init_value):
+    def __init__(self, init_value, init_density, Gp_prime):
         """
         Init of the class
 
         :param init_value: Value of the shear modulus
         """
-        super().__init__(init_value)
+        super().__init__(init_value, init_density, Gp_prime)
         self.init_value = init_value
+        self.init_density = init_density
 
-    def compute(self, density: np.array) -> np.array:
+    def compute(self, density: np.array, pressure) -> np.array:
         """
         Compute the shear modulus => returns constant value of shear modulus
 
         :param density: the current density
         :return: the computed shear modulus
         """
-        return np.ones_like(density) * self.init_value
\ No newline at end of file
+        
+        return np.ones_like(density) * self.init_value
diff --git a/xfv/src/rheology/constantyieldstress.py b/xfv/src/rheology/constantyieldstress.py
old mode 100755
new mode 100644
index 6a8189f..64d2f95
--- a/xfv/src/rheology/constantyieldstress.py
+++ b/xfv/src/rheology/constantyieldstress.py
@@ -12,16 +12,20 @@ class ConstantYieldStress(YieldStress):  # pylint: disable=too-few-public-method
     A class for constant yield stress calculation
     """
 
-    def __init__(self, init_value):
+    def __init__(self, init_value, init_shear_modulus, Y_max, beta, m):
         """
         Initialization of the constant yield stress class
 
         :param init_value: initial yield stress
         """
-        super().__init__(init_value)
+        super().__init__(init_value, init_shear_modulus, Y_max, beta, m)
         self.init_value = init_value
+        self.init_shear_modulus = init_shear_modulus
+        self.Y_max = Y_max
+        self.beta = beta
+        self.m = m 
 
-    def compute(self, density: np.array) -> np.array:
+    def compute(self, density: np.array, strain_plastic_eq, G) -> np.array:
         """
         Compute the value of the yield stress
 
diff --git a/xfv/src/rheology/scgshearmodulus.py b/xfv/src/rheology/scgshearmodulus.py
new file mode 100644
index 0000000..5c9428b
--- /dev/null
+++ b/xfv/src/rheology/scgshearmodulus.py
@@ -0,0 +1,41 @@
+# -*- coding: utf-8 -*-
+"""
+Implementation of the SCGShearModulus class
+"""
+
+import numpy as np
+from xfv.src.rheology.shearmodulus import ShearModulus
+
+
+class SCGShearModulus(ShearModulus):  # pylint: disable=too-few-public-methods
+    """
+    Class for scg shear modulus
+    """
+
+    def __init__(self, init_value, init_density, Gp_prime):
+        """
+        Init of the class
+
+        :param init_value: Value of the shear modulus
+        """
+        super().__init__(init_value, init_density, Gp_prime)
+        self.init_value = init_value
+        self.init_density = init_density
+        self.Gp_prime = Gp_prime
+	
+
+    def compute(self, density: np.array, pressure) -> np.array:
+        """
+        Compute the shear modulus => returns SCG value of shear modulus
+
+        :param density: the current density
+        :return: the computed shear modulus
+        """
+        is_pressure_below_zero = pressure < 0.
+        #is_pressure_ok = np.logical(is_pressure_below_zero)
+        pressure[is_pressure_below_zero] = 0.
+        density_part = np.power(density/self.init_density,(1./3.))
+        G = np.ones_like(density) * self.init_value *(1. + self.Gp_prime/self.init_value * pressure / density_part)
+        #print('G_max=',max(G))
+        #print('G_min=',min(G))
+        return G
diff --git a/xfv/src/rheology/scgyieldstress.py b/xfv/src/rheology/scgyieldstress.py
new file mode 100644
index 0000000..18aad3e
--- /dev/null
+++ b/xfv/src/rheology/scgyieldstress.py
@@ -0,0 +1,43 @@
+# -*- coding: utf-8 -*-
+"""
+Implementation of the ConstantYieldStress class
+"""
+
+import numpy as np
+from xfv.src.rheology.yieldstress import YieldStress
+
+
+class SCGYieldStress(YieldStress):  # pylint: disable=too-few-public-methods
+    """
+    A class for constant yield stress calculation
+    """
+
+    def __init__(self, init_value, init_shear_modulus, Y_max, beta, m):
+        """
+        Initialization of the constant yield stress class
+
+        :param init_value: initial yield stress
+        """
+        super().__init__(init_value, init_shear_modulus, Y_max, beta, m)
+        self.init_value = init_value
+        self.init_shear_modulus = init_shear_modulus
+        self.Y_max = Y_max
+        self.beta = beta
+        self.m = m 
+
+    def compute(self, density: np.array, strain_plastic_eq, G) -> np.array:
+        """
+        Compute the value of the yield stress
+
+        :param density: the current density
+        :return: the computed yield stress
+        """
+        #print('max epsilon_eq=',max(strain_plastic_eq))
+        part_Y = np.ones_like(density) * self.init_value*(1.+self.beta*np.power(strain_plastic_eq,self.m))
+        is_Ymax_below_partY = self.Y_max < part_Y
+        #is_Y_ok = np.logical(is_Ymax_below_partY)
+        part_Y[is_Ymax_below_partY] = self.Y_max
+        Y = part_Y*G/self.init_shear_modulus
+        #print('Y_max',max(Y))
+        #print('Y_min',min(Y))
+        return Y
diff --git a/xfv/src/rheology/shearmodulus.py b/xfv/src/rheology/shearmodulus.py
old mode 100755
new mode 100644
index 5bfdf1f..6539f6e
--- a/xfv/src/rheology/shearmodulus.py
+++ b/xfv/src/rheology/shearmodulus.py
@@ -10,11 +10,13 @@ class ShearModulus:  # pylint: disable=too-few-public-methods
     """
     Abstract class for the shear modulus computation
     """
-    def __init__(self, initial_value):
+    def __init__(self, initial_value, init_density, Gp_prime):
         self.shear_modulus = initial_value
+        self.initial_density = init_density
+        self.Gp_prime = Gp_prime
 
     @abstractmethod
-    def compute(self, density: np.array) -> np.array:
+    def compute(self, density: np.array, pressure) -> np.array:
         """
         Compute the new value of shear modulus
 
diff --git a/xfv/src/rheology/yieldstress.py b/xfv/src/rheology/yieldstress.py
old mode 100755
new mode 100644
index 36aeeb0..d768fda
--- a/xfv/src/rheology/yieldstress.py
+++ b/xfv/src/rheology/yieldstress.py
@@ -10,15 +10,20 @@ class YieldStress:  # pylint: disable=too-few-public-methods
     """
     Interface for yield stress computation
     """
-    def __init__(self, initial_value):
+    def __init__(self, init_value, init_shear_modulus, Y_max, beta, m):
         """
-        Initialization of the class
-        :param initial_value: yield_stress initial value
+        Initialization of the constant yield stress class
+
+        :param init_value: initial yield stress
         """
-        self.yield_stress = initial_value
+        self.init_value = init_value
+        self.init_shear_modulus = init_shear_modulus
+        self.Y_max = Y_max
+        self.beta = beta
+        self.m = m 
 
     @abstractmethod
-    def compute(self, density: np.array) -> np.array:
+    def compute(self, density: np.array, strain_plastic_eq, G) -> np.array:
         """
         Compute the new value of shear modulus
 
diff --git a/xfv/src/rupturecriterion/__init__.py b/xfv/src/rupturecriterion/__init__.py
index c66f7ef..1ed1e5b 100644
--- a/xfv/src/rupturecriterion/__init__.py
+++ b/xfv/src/rupturecriterion/__init__.py
@@ -7,3 +7,4 @@
 from xfv.src.rupturecriterion.porosity_criterion import PorosityCriterion
 from xfv.src.rupturecriterion.minimumpressure import MinimumPressureCriterion
 from xfv.src.rupturecriterion.maximalstress import MaximalStressCriterion
+from xfv.src.rupturecriterion.doublecriterion import DoubleCriterion
\ No newline at end of file
diff --git a/xfv/src/rupturecriterion/doublecriterion.py b/xfv/src/rupturecriterion/doublecriterion.py
new file mode 100644
index 0000000..96d1681
--- /dev/null
+++ b/xfv/src/rupturecriterion/doublecriterion.py
@@ -0,0 +1,35 @@
+# -*- coding: utf-8 -*-
+"""
+Implementation of PorosityCriterion class
+"""
+
+import numpy as np
+
+from xfv.src.rupturecriterion.rupturecriterion import RuptureCriterion
+from xfv.src.rupturecriterion.halfrodcomparison import HalfRodComparisonCriterion
+from xfv.src.rupturecriterion.porosity_criterion import PorosityCriterion
+
+class DoubleCriterion(RuptureCriterion):   # pylint: disable=too-few-public-methods
+    """
+    A rupture criterion based on porosity value
+    """
+    def __init__(self, p_limit, number_cell):
+        super().__init__()
+        self.__limit_porosity = p_limit
+        self._number_cell = number_cell
+        self.criterion_one = HalfRodComparisonCriterion(number_cell)
+        self.criterion_two = PorosityCriterion(p_limit)
+
+
+    def check_criterion(self, cells, *args, **kwargs):
+        """
+        Return the mask of the cells where porosity is above the threshold porosity
+
+        :param cells: cells on which to check the criterion
+        :return: the mask of the cells where porosity is above the threshold porosity
+        """
+        
+        array_criterion_one = self.criterion_one.check_criterion(cells)
+        array_criterion_two = self.criterion_two.check_criterion(cells)
+    
+        return np.logical_and(array_criterion_one, array_criterion_two)
diff --git a/xfv/src/rupturecriterion/porosity_criterion.py b/xfv/src/rupturecriterion/porosity_criterion.py
index 66af121..729978f 100644
--- a/xfv/src/rupturecriterion/porosity_criterion.py
+++ b/xfv/src/rupturecriterion/porosity_criterion.py
@@ -20,4 +20,4 @@ def check_criterion(self, cells, *args, **kwargs):
         :param cells: cells on which to check the criterion
         :return: the mask of the cells where porosity is above the threshold porosity
         """
-        return cells.porosity.new_value > self.__limit_porosity
+        return cells.porosity.new_value >= self.__limit_porosity
diff --git a/xfv/src/rupturetreatment/enrichelement.py b/xfv/src/rupturetreatment/enrichelement.py
index 1b29802..989c0e7 100755
--- a/xfv/src/rupturetreatment/enrichelement.py
+++ b/xfv/src/rupturetreatment/enrichelement.py
@@ -35,7 +35,7 @@ def lump_style(self) -> EnrichedMassMatrixProps:
         """
         return self.__lump
 
-    def apply_treatment(self, cells, ruptured_cells, nodes, topology, time):
+    def apply_treatment(self, cells, ruptured_cells, nodes, topology, time, cohesive_model, section):
         """
         Apply the rupture treatment by enriching one of the cells that is marked as ruptured cells
 
@@ -43,6 +43,9 @@ def apply_treatment(self, cells, ruptured_cells, nodes, topology, time):
         :param ruptured_cells: boolean array marking the ruptured cells
         :param nodes: array of all nodes
         :param topology: topology of the problem
+        :param time:
+        :param cohesive_model: class for the cohesive model
+        :param section: float for section of material
         """
         if ruptured_cells.any():  # Enrichment is made once for all
             cells_to_be_enr = np.logical_and(ruptured_cells, ~cells.enriched)
@@ -68,7 +71,7 @@ def apply_treatment(self, cells, ruptured_cells, nodes, topology, time):
                     print("==> In nodes : ", np.nonzero(in_nodes))
                     print("==> Out nodes : ", np.nonzero(out_nodes))
                     nodes.classical[nodes_to_be_enr] = False
-                    # Calcul coordonnées de la disc
+                    # Calcul coordonn�es de la disc
                     x_left = nodes.xt[nodes_to_be_enr[0]]
                     x_right = nodes.xt[nodes_to_be_enr[1]]
                     assert x_left < x_right
@@ -78,6 +81,7 @@ def apply_treatment(self, cells, ruptured_cells, nodes, topology, time):
                     disc = Discontinuity(cell_tb_enr, in_nodes, out_nodes,
                                          self.__position_rupture, self.__lump)
                     cells.classical[cell_tb_enr] = False
+                    disc.create_cohesive_law(cells, section, cohesive_model)
                     # Initialisation de la partie droite des champs + cell size
                     self.initialize_cracked_cell_size(cells, cell_tb_enr)
                     cells.initialize_additional_cell_dof(disc)
@@ -102,11 +106,11 @@ def initialize_cracked_cell_size(self, cells, cell_tb_enr):
         cells.left_part_size.new_value = \
             self.__position_rupture * cells.size_t_plus_dt[cell_tb_enr]
         # L'initialisation des tailles gauches et droites courantes n'est
-        # pas nécessaire. On initialise simplement avec des tailles fictives de
-        # sorte qu'on peut gérer le calcul des forces cohésives à l'itération où la
-        # discontinuité est créée. Cette taille ficitive permet simplement
-        # d'obtenir une ouverture nulle de la fissure à l'itération de création de
-        # la discontinuité.  Elle sera écrasée après ce calcul lors de l'appel
+        # pas n�cessaire. On initialise simplement avec des tailles fictives de
+        # sorte qu'on peut g�rer le calcul des forces coh�sives � l'it�ration o� la
+        # discontinuit� est cr��e. Cette taille ficitive permet simplement
+        # d'obtenir une ouverture nulle de la fissure � l'it�ration de cr�ation de
+        # la discontinuit�.  Elle sera �cras�e apr�s ce calcul lors de l'appel
         # de mesh.increment().
         cells.right_part_size.current_value = \
             (1. - self.__position_rupture) * cells.size_t[cell_tb_enr]