From 40ea8dd678bad72511c611706d3ef3364a7f52e8 Mon Sep 17 00:00:00 2001 From: Ulrich Houire Date: Wed, 27 Apr 2022 10:43:28 +0200 Subject: [PATCH 01/16] ajout saut de ligne --- xfv/src/discontinuity/discontinuity.py | 1 + 1 file changed, 1 insertion(+) diff --git a/xfv/src/discontinuity/discontinuity.py b/xfv/src/discontinuity/discontinuity.py index 3d11332..40b0c62 100755 --- a/xfv/src/discontinuity/discontinuity.py +++ b/xfv/src/discontinuity/discontinuity.py @@ -81,6 +81,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] From 32903df4384393997db2b80e493442cde2373117 Mon Sep 17 00:00:00 2001 From: Ulrich Houire Date: Wed, 27 Apr 2022 10:48:01 +0200 Subject: [PATCH 02/16] =?UTF-8?q?D=C3=A9mo=20saut=20de=20ligne?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- xfv/src/discontinuity/discontinuity.py | 1 + 1 file changed, 1 insertion(+) diff --git a/xfv/src/discontinuity/discontinuity.py b/xfv/src/discontinuity/discontinuity.py index 40b0c62..b10913d 100755 --- a/xfv/src/discontinuity/discontinuity.py +++ b/xfv/src/discontinuity/discontinuity.py @@ -246,6 +246,7 @@ 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() From c63f9872a647bfdd55650656d068cb96dae82c67 Mon Sep 17 00:00:00 2001 From: Ulrich Houire Date: Wed, 27 Apr 2022 10:57:10 +0200 Subject: [PATCH 03/16] commentaires --- xfv/XtendedFiniteVolume.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) mode change 100755 => 100644 xfv/XtendedFiniteVolume.py diff --git a/xfv/XtendedFiniteVolume.py b/xfv/XtendedFiniteVolume.py old mode 100755 new mode 100644 index 246274a..2cbe8ed --- 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 @@ -187,7 +187,7 @@ def main(directory: Path) -> None: # ---- # TIME MANAGEMENT final_time = data.time.final_time initial_time_step = data.time.initial_time_step - +# data.material_target.initial_values.rho_init # ---- # LOADING left_bc = data.boundary_condition.left_BC left_boundary_condition = _build_boundary_function(left_bc) @@ -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 @@ -348,7 +348,7 @@ def main(directory: Path) -> None: # CELLS DEVIATOR STRESSES COMPUTATION # # ---------------------------------------------# if projectile_elasticity: - my_mesh.apply_elasticity(dt, projectile_shear_modulus, my_mesh.cells.cell_in_projectile) + my_mesh.apply_elasticity(dt, projectile_shear_modulus, my_mesh.cells.cell_in_projectile,) if target_elasticity: my_mesh.apply_elasticity(dt, target_shear_modulus, my_mesh.cells.cell_in_target) # ---------------------------------------------# @@ -382,6 +382,7 @@ def main(directory: Path) -> None: # NODES FORCES COMPUTATION # # ---------------------------------------------# my_mesh.compute_new_nodes_forces() + #my_mesh.compute_new_cohesive_forces(data.material_target.cohesive_model.cohesive_zone_model_name) my_mesh.compute_new_cohesive_forces() # ---------------------------------------------# # LOADING # From ebd195a94f80b7fec8b84a83e2004aa78d3ea50c Mon Sep 17 00:00:00 2001 From: Ulrich Houire Date: Wed, 27 Apr 2022 10:59:01 +0200 Subject: [PATCH 04/16] Commentaires --- xfv/XtendedFiniteVolume.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/xfv/XtendedFiniteVolume.py b/xfv/XtendedFiniteVolume.py index 2cbe8ed..efc38aa 100644 --- a/xfv/XtendedFiniteVolume.py +++ b/xfv/XtendedFiniteVolume.py @@ -187,7 +187,7 @@ def main(directory: Path) -> None: # ---- # TIME MANAGEMENT final_time = data.time.final_time initial_time_step = data.time.initial_time_step -# data.material_target.initial_values.rho_init + # ---- # LOADING left_bc = data.boundary_condition.left_BC left_boundary_condition = _build_boundary_function(left_bc) @@ -348,7 +348,7 @@ def main(directory: Path) -> None: # CELLS DEVIATOR STRESSES COMPUTATION # # ---------------------------------------------# if projectile_elasticity: - my_mesh.apply_elasticity(dt, projectile_shear_modulus, my_mesh.cells.cell_in_projectile,) + my_mesh.apply_elasticity(dt, projectile_shear_modulus, my_mesh.cells.cell_in_projectile) if target_elasticity: my_mesh.apply_elasticity(dt, target_shear_modulus, my_mesh.cells.cell_in_target) # ---------------------------------------------# @@ -382,7 +382,6 @@ def main(directory: Path) -> None: # NODES FORCES COMPUTATION # # ---------------------------------------------# my_mesh.compute_new_nodes_forces() - #my_mesh.compute_new_cohesive_forces(data.material_target.cohesive_model.cohesive_zone_model_name) my_mesh.compute_new_cohesive_forces() # ---------------------------------------------# # LOADING # From 39d0b4570d49364c8149b291eb221fed084edc5b Mon Sep 17 00:00:00 2001 From: Ulrich Houire Date: Wed, 27 Apr 2022 11:23:35 +0200 Subject: [PATCH 05/16] =?UTF-8?q?mod=C3=A8le=20SCG?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- xfv/src/cell/cell.py | 9 ++++ xfv/src/cell/one_dimension_cell.py | 19 +++++-- .../one_dimension_enriched_cell_hansbo.py | 27 +++++++++- xfv/src/data/data_container.py | 50 +++++++++++++------ xfv/src/data/shear_modulus_props.py | 6 +++ xfv/src/data/yield_stress_props.py | 8 +++ xfv/src/mesh/mesh1denriched.py | 7 +++ xfv/src/rheology/__init__.py | 2 + xfv/src/rheology/constantshearmodulus.py | 10 ++-- xfv/src/rheology/constantyieldstress.py | 10 ++-- xfv/src/rheology/scgshearmodulus.py | 41 +++++++++++++++ xfv/src/rheology/scgyieldstress.py | 43 ++++++++++++++++ xfv/src/rheology/shearmodulus.py | 6 ++- xfv/src/rheology/yieldstress.py | 15 ++++-- 14 files changed, 219 insertions(+), 34 deletions(-) mode change 100755 => 100644 xfv/src/cell/cell.py mode change 100755 => 100644 xfv/src/cell/one_dimension_cell.py mode change 100755 => 100644 xfv/src/cell/one_dimension_enriched_cell_hansbo.py mode change 100755 => 100644 xfv/src/rheology/__init__.py mode change 100755 => 100644 xfv/src/rheology/constantshearmodulus.py mode change 100755 => 100644 xfv/src/rheology/constantyieldstress.py create mode 100644 xfv/src/rheology/scgshearmodulus.py create mode 100644 xfv/src/rheology/scgyieldstress.py mode change 100755 => 100644 xfv/src/rheology/shearmodulus.py mode change 100755 => 100644 xfv/src/rheology/yieldstress.py 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..a2539c1 --- a/xfv/src/cell/one_dimension_cell.py +++ b/xfv/src/cell/one_dimension_cell.py @@ -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..4374f9f --- a/xfv/src/cell/one_dimension_enriched_cell_hansbo.py +++ b/xfv/src/cell/one_dimension_enriched_cell_hansbo.py @@ -90,6 +90,7 @@ 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 @@ -143,6 +144,8 @@ def initialize_additional_cell_dof(self, disc: Discontinuity): self._enr_stress[enr_cell] = np.copy(self._stress[enr_cell]) self._enr_equivalent_plastic_strain_rate[enr_cell] = \ np.copy(self._equivalent_plastic_strain_rate[enr_cell]) + self._enr_equivalent_plastic_strain.current_value[enr_cell] = \ + np.copy(self.equivalent_plastic_strain.current_value[enr_cell]) def reconstruct_enriched_hydro_field(self, classical_field: Field, enriched_field_name: str): """ @@ -294,6 +297,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): """ @@ -681,7 +691,7 @@ 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]) def apply_plasticity_enr(self, mask_mesh, delta_t): """ @@ -707,6 +717,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 +736,8 @@ 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]) def compute_enriched_elements_new_time_step(self): """ @@ -766,6 +787,8 @@ 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() def compute_new_coordinates(self, topology, node_coord): """ diff --git a/xfv/src/data/data_container.py b/xfv/src/data/data_container.py index 603a855..4fc98f2 100755 --- a/xfv/src/data/data_container.py +++ b/xfv/src/data/data_container.py @@ -25,8 +25,8 @@ 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, @@ -553,7 +553,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 +568,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 +612,45 @@ 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: - 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": 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/mesh/mesh1denriched.py b/xfv/src/mesh/mesh1denriched.py index 4e4e478..68bb955 100755 --- a/xfv/src/mesh/mesh1denriched.py +++ b/xfv/src/mesh/mesh1denriched.py @@ -592,6 +592,13 @@ def stress_xx_field(self) -> np.array: """ return self.cells.stress_xx_field + @property + def equivalent_plastic_strain_field(self) -> np.array: + """ + Porosity field + """ + return self.cells.equivalent_plastic_strain_field + @property def porosity_field(self) -> np.array: """ 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 From 48f4a3364b51dff0cd98d3385c922e1cb93fabd5 Mon Sep 17 00:00:00 2001 From: Ulrich Houire Date: Wed, 27 Apr 2022 11:32:07 +0200 Subject: [PATCH 06/16] correction faute d'orthographe --- xfv/src/node/one_dimension_node.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) mode change 100755 => 100644 xfv/src/node/one_dimension_node.py 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') From cd76a1b80fb0573675a26d5509cc521929a69e7a Mon Sep 17 00:00:00 2001 From: Ulrich Houire Date: Wed, 27 Apr 2022 14:10:15 +0200 Subject: [PATCH 07/16] =?UTF-8?q?Impl=C3=A9mentation=20endommagement,=20co?= =?UTF-8?q?hesive,=20xfv?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- xfv/src/cohesive_model/cohesive_zone_model.py | 26 +++++++++----- xfv/src/data/cohesive_model_props.py | 28 +++++++++++++-- xfv/src/data/data_container.py | 35 +++++++++++++++---- xfv/src/discontinuity/discontinuity.py | 21 ++++++++++- xfv/src/mesh/mesh1denriched.py | 6 ++-- .../one_dimension_enriched_node_hansbo.py | 18 ++++++++-- .../rupturecriterion/porosity_criterion.py | 4 ++- 7 files changed, 113 insertions(+), 25 deletions(-) mode change 100755 => 100644 xfv/src/data/data_container.py mode change 100755 => 100644 xfv/src/discontinuity/discontinuity.py mode change 100755 => 100644 xfv/src/mesh/mesh1denriched.py mode change 100755 => 100644 xfv/src/node/one_dimension_enriched_node_hansbo.py diff --git a/xfv/src/cohesive_model/cohesive_zone_model.py b/xfv/src/cohesive_model/cohesive_zone_model.py index 7fd1680..d002b0d 100644 --- a/xfv/src/cohesive_model/cohesive_zone_model.py +++ b/xfv/src/cohesive_model/cohesive_zone_model.py @@ -12,7 +12,7 @@ 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): """ Construction d'un modèle cohésif @@ -21,38 +21,48 @@ def __init__(self, cohesive_law_points: np.array, unloading_model: UnloadingMode :param unloading_model: """ self._critical_separation = cohesive_law_points[-1, 0] + self._critical_strength = cohesive_law_points[0,-1] self._cohesive_law = CohesiveLaw(cohesive_law_points) self._unloading_model = unloading_model - + self._cohesive_zone_model_name = cohesive_zone_model_name + self._dissipated_energy = dissipated_energy + self._purcentage = purcentage + def compute_cohesive_stress(self, disc): """ Compute the cohesive force for the current opening of discontinuity according to the current discontinuity opening :param disc: discontinuity + + attention : formule de la dissipiation 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] + self._cohesive_law = CohesiveLaw(np.array([[0, disc.critical_strength], [disc.critical_separation, 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: + elif disc.history_max_opening <= new_opening < disc.critical_separation : cohesive_force = self._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 - - elif new_opening >= self._critical_separation: + disc.damage_variable.new_value = new_opening / disc.critical_separation + disc.dissipated_energy.new_value = disc.critical_strength*disc.history_max_opening/2. + 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..8f7e0f3 100644 --- a/xfv/src/data/cohesive_model_props.py +++ b/xfv/src/data/cohesive_model_props.py @@ -15,24 +15,46 @@ class CohesiveZoneModelProps(TypeCheckedDataClass): cohesive_strength: float critical_opening: float + cohesive_zone_model_name: str unloading_model: UnloadingModelProps + dissipated_energy: float + purcentage_internal_energy: float _cohesive_zone_model_class = CohesiveZoneModel - + def _build_cohesive_law(self): """ Build the cohesive law that is needed by the CohesiveModel """ raise NotImplementedError("This is an abstract method!") + + def cohesive_model_name(self): + return self.cohesive_zone_model_name + + def purcentage(self): + return self.purcentage_internal_energy + + def dissip_energy(self): + return self.dissipated_energy 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_model_name(), self.dissip_energy(), + self.purcentage()) @dataclass # pylint: disable=missing-class-docstring +class LinearMixedCohesiveZoneModelProps(CohesiveZoneModelProps): + + def _build_cohesive_law(self): + return np.array([ + [0, self.cohesive_strength], + [self.critical_opening, 0]]) + + +@dataclass # pylint: disable=missing-class-docstring class LinearCohesiveZoneModelProps(CohesiveZoneModelProps): def _build_cohesive_law(self): diff --git a/xfv/src/data/data_container.py b/xfv/src/data/data_container.py old mode 100755 new mode 100644 index 4fc98f2..a78e3c7 --- a/xfv/src/data/data_container.py +++ b/xfv/src/data/data_container.py @@ -17,6 +17,7 @@ SuccessiveRampFunctionProps) from xfv.src.data.cohesive_model_props import (CohesiveZoneModelProps, LinearCohesiveZoneModelProps, + LinearMixedCohesiveZoneModelProps, BilinearCohesiveZoneModelProps, TrilinearCohesiveZoneModelProps) from xfv.src.data.unloading_model_props import (UnloadingModelProps, @@ -99,7 +100,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 @@ -400,25 +401,45 @@ def __get_cohesive_model_props(matter) -> Optional[CohesiveZoneModelProps]: raise ValueError(f"Unknwown unloading model name: {unloading_model_name} " "Please choose among (progressiveunloading, " " lossofstiffnessunloading)") - + print('cohesive model is', cohesive_model_name) + if cohesive_model_name == "linear": - cohesive_model_props: CohesiveZoneModelProps = LinearCohesiveZoneModelProps( - cohesive_strength, critical_separation, unloading_model_props) + purcentage_internal_energy = 0. + dissipated_energy = 0. + cohesive_model_props = LinearCohesiveZoneModelProps(cohesive_strength,critical_separation, + cohesive_model_name,unloading_model_props,dissipated_energy,purcentage_internal_energy) + elif cohesive_model_name == "linear_mixed_purcentage": + dissipated_energy = 0. + purcentage_internal_energy = params['coefficients']['purcentage-internal-energy'] + cohesive_model_props = LinearMixedCohesiveZoneModelProps(0.,0., + cohesive_model_name,unloading_model_props, + dissipated_energy, purcentage_internal_energy) + elif cohesive_model_name == "linear_mixed_fixed": + dissipated_energy = params['coefficients']['dissipated-energy'] + purcentage_internal_energy = 0. + cohesive_model_props = LinearMixedCohesiveZoneModelProps(0.,0., + cohesive_model_name,unloading_model_props, + dissipated_energy, purcentage_internal_energy) elif cohesive_model_name == "bilinear": + 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_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 (linear, linear_mixed_purcentage, linear_mixed_fixed, bilinear, trilinear)") return cohesive_model_props diff --git a/xfv/src/discontinuity/discontinuity.py b/xfv/src/discontinuity/discontinuity.py old mode 100755 new mode 100644 index b10913d..6a72323 --- a/xfv/src/discontinuity/discontinuity.py +++ b/xfv/src/discontinuity/discontinuity.py @@ -63,6 +63,8 @@ def __init__(self, cell_id: int, mask_in_nodes: np.array, mask_out_nodes: np.arr Discontinuity.ruptured_cell_id = np.zeros((1, 1), dtype=int) Discontinuity.in_nodes = np.zeros((1, 1), dtype=int) Discontinuity.out_nodes = np.zeros((1, 1), dtype=int) + Discontinuity.critical_strength = np.zeros((1,1), dtype = float) + Discontinuity.critical_separation = np.zeros((1,1), dtype = float) else: Discontinuity.enr_velocity_current = np.append( Discontinuity.enr_velocity_current, init, axis=0) @@ -81,7 +83,11 @@ 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) - + Discontinuity.critical_strength = np.append( + Discontinuity.critical_strength, np.zeros((1, 1), dtype=float), axis=0) + Discontinuity.critical_separation = np.append( + Discontinuity.critical_separation, np.zeros((1, 1), dtype=float), 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] @@ -92,6 +98,10 @@ 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] + disc.critical_strength = Discontinuity.critical_strength[ind] + disc.critical_separation = Discontinuity.critical_separation[ind] + + self.__mask_in_nodes = mask_in_nodes self.in_nodes[:] = np.where(self.__mask_in_nodes)[0] @@ -112,6 +122,7 @@ 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 @@ -232,6 +243,13 @@ 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 compute_critical_value(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 @@ -251,3 +269,4 @@ def enr_increment(self): 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 68bb955..f94228b --- a/xfv/src/mesh/mesh1denriched.py +++ b/xfv/src/mesh/mesh1denriched.py @@ -322,7 +322,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): """ @@ -473,7 +475,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(): 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..b280a92 --- 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,22 @@ 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 + if cohesive_model._cohesive_zone_model_name == "linear_mixed_purcentage" : + Discontinuity.compute_critical_value(stress, cohesive_model._purcentage*energy/self.section) + elif cohesive_model._cohesive_zone_model_name == "linear_mixed_fixed" : + dissipated_energy_array = cohesive_model._dissipated_energy*np.ones_like(energy) + Discontinuity.compute_critical_value(stress, dissipated_energy_array/self.section) + for ind, disc in enumerate(disc_list): + if dissipated_energy_array[0] > energy[Discontinuity.ruptured_cell_id[ind]]: + raise ValueError(f" Dissipated energy bigger than internal energy") + else: + for ind, disc in enumerate(disc_list): + disc.critical_separation = cohesive_model._critical_separation + disc.critical_strength = cohesive_model._critical_strength + 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/rupturecriterion/porosity_criterion.py b/xfv/src/rupturecriterion/porosity_criterion.py index 66af121..0fdee88 100644 --- a/xfv/src/rupturecriterion/porosity_criterion.py +++ b/xfv/src/rupturecriterion/porosity_criterion.py @@ -20,4 +20,6 @@ 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 + is_porosity_over_limit_porosity = cells.porosity.new_value >= self.__limit_porosity + cells.porosity.new_value[is_porosity_over_limit_porosity] = self.__limit_porosity + return cells.porosity.new_value >= self.__limit_porosity From 95774b42a50a5cc1c72b0b24dcbfd18338a30b32 Mon Sep 17 00:00:00 2001 From: Ulrich Houire Date: Wed, 27 Apr 2022 14:28:04 +0200 Subject: [PATCH 08/16] Correction respect des donnees locales --- xfv/src/cohesive_model/cohesive_zone_model.py | 39 ++++++++++++++++++- xfv/src/data/cohesive_model_props.py | 2 +- .../one_dimension_enriched_node_hansbo.py | 12 +++--- 3 files changed, 44 insertions(+), 9 deletions(-) diff --git a/xfv/src/cohesive_model/cohesive_zone_model.py b/xfv/src/cohesive_model/cohesive_zone_model.py index d002b0d..13037f4 100644 --- a/xfv/src/cohesive_model/cohesive_zone_model.py +++ b/xfv/src/cohesive_model/cohesive_zone_model.py @@ -27,7 +27,42 @@ def __init__(self, cohesive_law_points: np.array, unloading_model: UnloadingMode self._cohesive_zone_model_name = cohesive_zone_model_name self._dissipated_energy = dissipated_energy self._purcentage = purcentage - + + @property + def cohesive_zone_model_name(self): # pylint: disable=invalid-name + """ + Cohesive zone model name + """ + return self._cohesive_zone_model_name + + @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): """ Compute the cohesive force for the current opening of discontinuity according to the @@ -35,7 +70,7 @@ def compute_cohesive_stress(self, disc): :param disc: discontinuity - attention : formule de la dissipiation d'energie du modèle cohesif seulement dans le cas de la décharge à zéro + 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] diff --git a/xfv/src/data/cohesive_model_props.py b/xfv/src/data/cohesive_model_props.py index 8f7e0f3..bc7d4cf 100644 --- a/xfv/src/data/cohesive_model_props.py +++ b/xfv/src/data/cohesive_model_props.py @@ -20,7 +20,7 @@ class CohesiveZoneModelProps(TypeCheckedDataClass): dissipated_energy: float purcentage_internal_energy: float _cohesive_zone_model_class = CohesiveZoneModel - + def _build_cohesive_law(self): """ Build the cohesive law that is needed by the CohesiveModel diff --git a/xfv/src/node/one_dimension_enriched_node_hansbo.py b/xfv/src/node/one_dimension_enriched_node_hansbo.py index b280a92..123ad21 100644 --- a/xfv/src/node/one_dimension_enriched_node_hansbo.py +++ b/xfv/src/node/one_dimension_enriched_node_hansbo.py @@ -204,18 +204,18 @@ def compute_enriched_nodes_cohesive_forces(self, cohesive_model, stress, energy) nb_disc = len(disc_list) applied_force_arr = np.ndarray((nb_disc,)) - if cohesive_model._cohesive_zone_model_name == "linear_mixed_purcentage" : - Discontinuity.compute_critical_value(stress, cohesive_model._purcentage*energy/self.section) - elif cohesive_model._cohesive_zone_model_name == "linear_mixed_fixed" : - dissipated_energy_array = cohesive_model._dissipated_energy*np.ones_like(energy) + if cohesive_model.cohesive_zone_model_name == "linear_mixed_purcentage" : + Discontinuity.compute_critical_value(stress, cohesive_model.purcentage*energy/self.section) + elif cohesive_model.cohesive_zone_model_name == "linear_mixed_fixed" : + dissipated_energy_array = cohesive_model.dissipated_energy*np.ones_like(energy) Discontinuity.compute_critical_value(stress, dissipated_energy_array/self.section) for ind, disc in enumerate(disc_list): if dissipated_energy_array[0] > energy[Discontinuity.ruptured_cell_id[ind]]: raise ValueError(f" Dissipated energy bigger than internal energy") else: for ind, disc in enumerate(disc_list): - disc.critical_separation = cohesive_model._critical_separation - disc.critical_strength = cohesive_model._critical_strength + disc.critical_separation = cohesive_model.critical_separation + disc.critical_strength = cohesive_model.critical_strength for ind, disc in enumerate(disc_list): cohesive_stress = cohesive_model.compute_cohesive_stress(disc) disc.cohesive_force.new_value = cohesive_stress From a66c0b0a8bf282434e10b003bebd39217dbe00cb Mon Sep 17 00:00:00 2001 From: Ulrich Houire Date: Thu, 5 May 2022 09:11:20 +0200 Subject: [PATCH 09/16] Post-processing changes --- .../cohesive_dissipated_energy.py | 107 +++++++++++++++++ xfv/post_processing/cohesive_law.py | 108 ++++++++++++++++++ .../field_evolution_in_time.py | 10 +- xfv/post_processing/free_surface_velocity.py | 6 +- .../profile_figure_in_space.py | 4 +- 5 files changed, 229 insertions(+), 6 deletions(-) create mode 100644 xfv/post_processing/cohesive_dissipated_energy.py create mode 100755 xfv/post_processing/cohesive_law.py mode change 100755 => 100644 xfv/post_processing/free_surface_velocity.py 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..8557d6d 100755 --- a/xfv/post_processing/profile_figure_in_space.py +++ b/xfv/post_processing/profile_figure_in_space.py @@ -39,10 +39,10 @@ def run(): print("Done !") print("~~~~~~~~~~~~~") # Plot field : - plt.plot(coord * 1.e+03, field_value, label=case) + plt.plot(coord * 1.e+03, field_value, label='temps '+str(current_time)) if (ARGS.write_data): - data_path="{:s}Profile_{:s}_{:3.1e}.dat".format(case, ARGS.field, current_time) + 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)) From 4567f867fcba436474e13f05d7133c97c3aaae6c Mon Sep 17 00:00:00 2001 From: Ulrich Houire Date: Thu, 5 May 2022 09:15:21 +0200 Subject: [PATCH 10/16] Post processing changes (additionnal) --- xfv/src/output_manager/outputdatabaseexploit.py | 1 + xfv/src/output_manager/outputmanager.py | 3 +++ 2 files changed, 4 insertions(+) mode change 100755 => 100644 xfv/src/output_manager/outputdatabaseexploit.py mode change 100755 => 100644 xfv/src/output_manager/outputmanager.py diff --git a/xfv/src/output_manager/outputdatabaseexploit.py b/xfv/src/output_manager/outputdatabaseexploit.py old mode 100755 new mode 100644 index c993ab5..61d2b47 --- a/xfv/src/output_manager/outputdatabaseexploit.py +++ b/xfv/src/output_manager/outputdatabaseexploit.py @@ -34,6 +34,7 @@ class OutputDatabaseExploit: field_type_converter["ShearModulus"] = ("ClassicalShearModulus", "None") field_type_converter["YieldStress"] = ("ClassicalYieldStress", "None") 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..8f0c982 --- 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", From 640996e0b0c681f92853e96ec463abdc1f6e3080 Mon Sep 17 00:00:00 2001 From: Ulrich Houire Date: Thu, 5 May 2022 09:55:38 +0200 Subject: [PATCH 11/16] Improvement of coupling method (xfv, czm, endo) --- xfv/src/cohesive_calculation/__init__.py | 7 ++ .../cohesivecalculationmodel.py | 26 ++++++++ xfv/src/cohesive_calculation/lineardata.py | 31 +++++++++ xfv/src/cohesive_calculation/linearenergy.py | 33 ++++++++++ xfv/src/cohesive_calculation/linearpercent.py | 33 ++++++++++ xfv/src/cohesive_model/cohesive_zone_model.py | 24 +++++-- xfv/src/data/cohesive_model_props.py | 66 ++++++++++++++----- xfv/src/data/data_container.py | 55 ++++++++++------ xfv/src/discontinuity/discontinuity.py | 29 +++++--- xfv/src/mesh/mesh1denriched.py | 3 +- .../one_dimension_enriched_node_hansbo.py | 12 ---- xfv/src/rupturetreatment/enrichelement.py | 18 +++-- 12 files changed, 264 insertions(+), 73 deletions(-) create mode 100644 xfv/src/cohesive_calculation/__init__.py create mode 100644 xfv/src/cohesive_calculation/cohesivecalculationmodel.py create mode 100644 xfv/src/cohesive_calculation/lineardata.py create mode 100644 xfv/src/cohesive_calculation/linearenergy.py create mode 100644 xfv/src/cohesive_calculation/linearpercent.py 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 13037f4..49bb2d9 100644 --- a/xfv/src/cohesive_model/cohesive_zone_model.py +++ b/xfv/src/cohesive_model/cohesive_zone_model.py @@ -12,21 +12,25 @@ class CohesiveZoneModel: """ A class for the computation of the cohesive force """ - def __init__(self, cohesive_law_points: np.array, unloading_model: UnloadingModelBase, cohesive_zone_model_name, dissipated_energy, purcentage): + 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._critical_strength = cohesive_law_points[0,-1] - self._cohesive_law = CohesiveLaw(cohesive_law_points) 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 @@ -35,6 +39,13 @@ def cohesive_zone_model_name(self): # pylint: disable=invalid-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 """ @@ -74,8 +85,6 @@ def compute_cohesive_stress(self, disc): """ cohesive_force = 0. new_opening = disc.discontinuity_opening.new_value[0] - self._cohesive_law = CohesiveLaw(np.array([[0, disc.critical_strength], [disc.critical_separation, 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. @@ -86,18 +95,19 @@ def compute_cohesive_stress(self, disc): disc.dissipated_energy.new_value = disc.critical_strength*disc.history_max_opening/2. elif disc.history_max_opening <= new_opening < disc.critical_separation : - cohesive_force = self._cohesive_law.compute_cohesive_force(new_opening) + 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._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. + 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 bc7d4cf..06a4e2e 100644 --- a/xfv/src/data/cohesive_model_props.py +++ b/xfv/src/data/cohesive_model_props.py @@ -9,16 +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 + _cohesive_zone_model_name: str unloading_model: UnloadingModelProps - dissipated_energy: float - purcentage_internal_energy: float + _dissipated_energy: float + _purcentage: float _cohesive_zone_model_class = CohesiveZoneModel def _build_cohesive_law(self): @@ -26,15 +29,20 @@ def _build_cohesive_law(self): Build the cohesive law that is needed by the CohesiveModel """ raise NotImplementedError("This is an abstract method!") - - def cohesive_model_name(self): - return self.cohesive_zone_model_name - def purcentage(self): - return self.purcentage_internal_energy - - def dissip_energy(self): - return self.dissipated_energy + +@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): """ @@ -42,25 +50,47 @@ def build_cohesive_model_obj(self): """ return self._cohesive_zone_model_class(self._build_cohesive_law(), self.unloading_model.build_unloading_model_obj(), - self.cohesive_model_name(), self.dissip_energy(), - self.purcentage()) + self._cohesive_zone_model_name, self._dissipated_energy, + self._purcentage, self._cohesive_calculation_model) @dataclass # pylint: disable=missing-class-docstring -class LinearMixedCohesiveZoneModelProps(CohesiveZoneModelProps): +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): + +@dataclass # pylint: disable=missing-class-docstring +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 @@ -81,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 index a78e3c7..f74fd69 100644 --- a/xfv/src/data/data_container.py +++ b/xfv/src/data/data_container.py @@ -16,8 +16,9 @@ MarchTableFunctionProps, SuccessiveRampFunctionProps) from xfv.src.data.cohesive_model_props import (CohesiveZoneModelProps, - LinearCohesiveZoneModelProps, - LinearMixedCohesiveZoneModelProps, + LinearDataCohesiveZoneModelProps, + LinearPercentCohesiveZoneModelProps, + LinearEnergyCohesiveZoneModelProps, BilinearCohesiveZoneModelProps, TrilinearCohesiveZoneModelProps) from xfv.src.data.unloading_model_props import (UnloadingModelProps, @@ -386,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') @@ -403,43 +402,59 @@ def __get_cohesive_model_props(matter) -> Optional[CohesiveZoneModelProps]: " lossofstiffnessunloading)") print('cohesive model is', cohesive_model_name) - if cohesive_model_name == "linear": + 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 = LinearCohesiveZoneModelProps(cohesive_strength,critical_separation, - cohesive_model_name,unloading_model_props,dissipated_energy,purcentage_internal_energy) - elif cohesive_model_name == "linear_mixed_purcentage": - dissipated_energy = 0. - purcentage_internal_energy = params['coefficients']['purcentage-internal-energy'] - cohesive_model_props = LinearMixedCohesiveZoneModelProps(0.,0., - cohesive_model_name,unloading_model_props, - dissipated_energy, purcentage_internal_energy) - elif cohesive_model_name == "linear_mixed_fixed": + 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. - cohesive_model_props = LinearMixedCohesiveZoneModelProps(0.,0., - cohesive_model_name,unloading_model_props, + + # 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, cohesive_model_name, unloading_model_props, - dissipated_energy, purcentage_internal_energy,purcentage_internal_energy,purcentage_internal_energy, + 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, cohesive_model_name, unloading_model_props, - dissipated_energy, purcentage_internal_energy,purcentage_internal_energy,purcentage_internal_energy, + 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, linear_mixed_purcentage, linear_mixed_fixed, bilinear, trilinear)") + "Please choose among (LinearData, LinearPercent, LinearEnergy, bilinear, trilinear)") return cohesive_model_props @@ -531,6 +546,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) diff --git a/xfv/src/discontinuity/discontinuity.py b/xfv/src/discontinuity/discontinuity.py index 6a72323..69759d0 100644 --- 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): @@ -63,8 +69,6 @@ def __init__(self, cell_id: int, mask_in_nodes: np.array, mask_out_nodes: np.arr Discontinuity.ruptured_cell_id = np.zeros((1, 1), dtype=int) Discontinuity.in_nodes = np.zeros((1, 1), dtype=int) Discontinuity.out_nodes = np.zeros((1, 1), dtype=int) - Discontinuity.critical_strength = np.zeros((1,1), dtype = float) - Discontinuity.critical_separation = np.zeros((1,1), dtype = float) else: Discontinuity.enr_velocity_current = np.append( Discontinuity.enr_velocity_current, init, axis=0) @@ -83,10 +87,6 @@ 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) - Discontinuity.critical_strength = np.append( - Discontinuity.critical_strength, np.zeros((1, 1), dtype=float), axis=0) - Discontinuity.critical_separation = np.append( - Discontinuity.critical_separation, np.zeros((1, 1), dtype=float), axis=0) for ind, disc in enumerate(Discontinuity.discontinuity_list()): disc.enr_velocity_current = Discontinuity.enr_velocity_current[ind] @@ -98,8 +98,6 @@ 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] - disc.critical_strength = Discontinuity.critical_strength[ind] - disc.critical_separation = Discontinuity.critical_separation[ind] @@ -126,6 +124,9 @@ def __init__(self, cell_id: int, mask_in_nodes: np.array, mask_out_nodes: np.arr 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() @@ -243,7 +244,17 @@ 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 compute_critical_value(stress, energy): + 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]]) diff --git a/xfv/src/mesh/mesh1denriched.py b/xfv/src/mesh/mesh1denriched.py index f94228b..1a260b6 100644 --- a/xfv/src/mesh/mesh1denriched.py +++ b/xfv/src/mesh/mesh1denriched.py @@ -455,7 +455,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): diff --git a/xfv/src/node/one_dimension_enriched_node_hansbo.py b/xfv/src/node/one_dimension_enriched_node_hansbo.py index 123ad21..4f64fa5 100644 --- a/xfv/src/node/one_dimension_enriched_node_hansbo.py +++ b/xfv/src/node/one_dimension_enriched_node_hansbo.py @@ -204,18 +204,6 @@ def compute_enriched_nodes_cohesive_forces(self, cohesive_model, stress, energy) nb_disc = len(disc_list) applied_force_arr = np.ndarray((nb_disc,)) - if cohesive_model.cohesive_zone_model_name == "linear_mixed_purcentage" : - Discontinuity.compute_critical_value(stress, cohesive_model.purcentage*energy/self.section) - elif cohesive_model.cohesive_zone_model_name == "linear_mixed_fixed" : - dissipated_energy_array = cohesive_model.dissipated_energy*np.ones_like(energy) - Discontinuity.compute_critical_value(stress, dissipated_energy_array/self.section) - for ind, disc in enumerate(disc_list): - if dissipated_energy_array[0] > energy[Discontinuity.ruptured_cell_id[ind]]: - raise ValueError(f" Dissipated energy bigger than internal energy") - else: - for ind, disc in enumerate(disc_list): - disc.critical_separation = cohesive_model.critical_separation - disc.critical_strength = cohesive_model.critical_strength for ind, disc in enumerate(disc_list): cohesive_stress = cohesive_model.compute_cohesive_stress(disc) disc.cohesive_force.new_value = cohesive_stress 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] From b750d2fdec254c545cc627e630a0842eff2f3cb5 Mon Sep 17 00:00:00 2001 From: Ulrich Houire Date: Tue, 10 May 2022 09:39:58 +0200 Subject: [PATCH 12/16] =?UTF-8?q?Compl=C3=A9ment=20sur=20la=20rh=C3=A9olog?= =?UTF-8?q?ie=20des=20elements=20enrichis?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../one_dimension_enriched_cell_hansbo.py | 74 ++++++++++++++++++- 1 file changed, 71 insertions(+), 3 deletions(-) diff --git a/xfv/src/cell/one_dimension_enriched_cell_hansbo.py b/xfv/src/cell/one_dimension_enriched_cell_hansbo.py index 4374f9f..e62e4b0 100644 --- a/xfv/src/cell/one_dimension_enriched_cell_hansbo.py +++ b/xfv/src/cell/one_dimension_enriched_cell_hansbo.py @@ -96,6 +96,9 @@ def __init__(self, n_cells: int): # 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 @@ -121,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) @@ -138,14 +144,16 @@ 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]) self._enr_stress[enr_cell] = np.copy(self._stress[enr_cell]) self._enr_equivalent_plastic_strain_rate[enr_cell] = \ np.copy(self._equivalent_plastic_strain_rate[enr_cell]) - self._enr_equivalent_plastic_strain.current_value[enr_cell] = \ - np.copy(self.equivalent_plastic_strain.current_value[enr_cell]) def reconstruct_enriched_hydro_field(self, classical_field: Field, enriched_field_name: str): """ @@ -311,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): """ @@ -465,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 @@ -479,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], @@ -505,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 @@ -693,6 +744,13 @@ def compute_enriched_shear_modulus(self, shear_modulus_model): self.enr_shear_modulus.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): """ Apply plasticity treatment if criterion is activated : @@ -739,6 +797,11 @@ def compute_enriched_yield_stress(self, yield_stress_model): 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): """ Compute the new local time step. @@ -787,8 +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 + # 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): """ From ce5c06a1e7313c3eb0ebf620c6972b24a090f836 Mon Sep 17 00:00:00 2001 From: Ulrich Houire Date: Tue, 10 May 2022 09:47:18 +0200 Subject: [PATCH 13/16] =?UTF-8?q?Compl=C3=A9ment=20sur=20la=20porosit?= =?UTF-8?q?=C3=A9=20des=20cellules=20enrichies?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- xfv/src/data/data_container.py | 12 +++++++++++- xfv/src/data/porosity_model_props.py | 1 + xfv/src/mesh/mesh1denriched.py | 6 +++++- 3 files changed, 17 insertions(+), 2 deletions(-) diff --git a/xfv/src/data/data_container.py b/xfv/src/data/data_container.py index f74fd69..8d79865 100644 --- a/xfv/src/data/data_container.py +++ b/xfv/src/data/data_container.py @@ -476,10 +476,20 @@ 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) 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)." 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/mesh/mesh1denriched.py b/xfv/src/mesh/mesh1denriched.py index 1a260b6..7b1c6e8 100644 --- 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): """ @@ -598,7 +602,7 @@ def stress_xx_field(self) -> np.array: @property def equivalent_plastic_strain_field(self) -> np.array: """ - Porosity field + equivalent plastic strain field """ return self.cells.equivalent_plastic_strain_field From adaa76cb7195fe8b084295d20d1e57f0a004e7de Mon Sep 17 00:00:00 2001 From: Ulrich Houire Date: Tue, 10 May 2022 09:53:59 +0200 Subject: [PATCH 14/16] =?UTF-8?q?Compl=C3=A9ment=20sur=20la=20porosit?= =?UTF-8?q?=C3=A9=20des=20cellules=20enrichies?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- xfv/src/porosity_model/johnsonmodel.py | 7 +++++-- xfv/src/rupturecriterion/porosity_criterion.py | 2 -- 2 files changed, 5 insertions(+), 4 deletions(-) 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/rupturecriterion/porosity_criterion.py b/xfv/src/rupturecriterion/porosity_criterion.py index 0fdee88..729978f 100644 --- a/xfv/src/rupturecriterion/porosity_criterion.py +++ b/xfv/src/rupturecriterion/porosity_criterion.py @@ -20,6 +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 """ - is_porosity_over_limit_porosity = cells.porosity.new_value >= self.__limit_porosity - cells.porosity.new_value[is_porosity_over_limit_porosity] = self.__limit_porosity return cells.porosity.new_value >= self.__limit_porosity From da225c36de70431db08b804b9b40b2987a4f9e3e Mon Sep 17 00:00:00 2001 From: Ulrich Houire Date: Tue, 10 May 2022 09:55:04 +0200 Subject: [PATCH 15/16] =?UTF-8?q?Incr=C3=A9mentation=20du=20crit=C3=A8re?= =?UTF-8?q?=20double?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- xfv/src/cell/one_dimension_cell.py | 2 +- xfv/src/data/data_container.py | 12 ++++++- xfv/src/data/rupture_criterion_props.py | 14 +++++++-- xfv/src/rupturecriterion/__init__.py | 1 + xfv/src/rupturecriterion/doublecriterion.py | 35 +++++++++++++++++++++ 5 files changed, 60 insertions(+), 4 deletions(-) create mode 100644 xfv/src/rupturecriterion/doublecriterion.py diff --git a/xfv/src/cell/one_dimension_cell.py b/xfv/src/cell/one_dimension_cell.py index a2539c1..47869d1 100644 --- 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) diff --git a/xfv/src/data/data_container.py b/xfv/src/data/data_container.py index 8d79865..b6ddb3a 100644 --- a/xfv/src/data/data_container.py +++ b/xfv/src/data/data_container.py @@ -31,7 +31,7 @@ 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, @@ -486,6 +486,10 @@ def __get_porosity_model_props(matter) -> Optional[Tuple[PorosityModelProps, str 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, @@ -669,6 +673,8 @@ def __get_rheology_props(self, matter) -> Tuple[Optional[ShearModulusProps], coef = coef["InitThermo"] rho_0 = float(coef["initial_density"]) params = matter.get('rheology') + 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) @@ -784,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/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/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) From ce273f87a9bf9d2e66e915b70e82e4ff79e2faa3 Mon Sep 17 00:00:00 2001 From: Ulrich Houire Date: Tue, 10 May 2022 09:55:40 +0200 Subject: [PATCH 16/16] =?UTF-8?q?Compl=C3=A9ments=20sur=20les=20output=20e?= =?UTF-8?q?t=20le=20post-processing?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- xfv/post_processing/profile_figure_in_space.py | 15 +++++++-------- xfv/src/output_manager/outputdatabaseexploit.py | 6 +++--- xfv/src/output_manager/outputmanager.py | 7 ++++++- 3 files changed, 16 insertions(+), 12 deletions(-) diff --git a/xfv/post_processing/profile_figure_in_space.py b/xfv/post_processing/profile_figure_in_space.py index 8557d6d..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='temps '+str(current_time)) - - 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)) + 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/output_manager/outputdatabaseexploit.py b/xfv/src/output_manager/outputdatabaseexploit.py index 61d2b47..b54e3a7 100644 --- a/xfv/src/output_manager/outputdatabaseexploit.py +++ b/xfv/src/output_manager/outputdatabaseexploit.py @@ -30,9 +30,9 @@ 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") diff --git a/xfv/src/output_manager/outputmanager.py b/xfv/src/output_manager/outputmanager.py index 8f0c982..5212c4d 100644 --- a/xfv/src/output_manager/outputmanager.py +++ b/xfv/src/output_manager/outputmanager.py @@ -74,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): """