Travaux d'Ulrich à reporter

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

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()

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()

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))

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

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__":

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):
         """