diff --git a/examples/03-advanced/composite_DCB_test_3D.py b/examples/03-advanced/composite_DCB_test_3D.py new file mode 100644 index 00000000..d61d8b60 --- /dev/null +++ b/examples/03-advanced/composite_DCB_test_3D.py @@ -0,0 +1,274 @@ +""" +Composite double-cantilever beam with a cohesive interface +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The model consists of two composite arms separated by a zero-thickness +cohesive interface. Cohesive elements are inserted only after the prescribed +starter-crack length. The initially unconnected part of the two arms therefore +represents the pre-crack. + +The arms are made from a unidirectional composite whose fibres are aligned +with the beam axis. A mixed-mode bilinear +:class:`fedoo.constitutivelaw.CohesiveLaw` governs the interface, and the +opening displacement is applied incrementally with +:class:`fedoo.problem.NonLinear`. + +The example uses N, mm, and MPa. Its compact mesh demonstrates the cohesive-law +workflow; it is not intended as a mesh-converged fracture benchmark. +""" + +from __future__ import annotations + +import fedoo as fd +import numpy as np + + +############################################################################### +# Geometry and discretization +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# We first define the DCB dimensions, loading, and compact mesh resolution. +# The mesh is deliberately small so that the complete nonlinear example +# remains fast enough for the documentation gallery. + +# Number of nodes along the beam length, width, and through each arm. +nx = 19 +ny = 5 +nz_per_arm = 3 + +length = 60.0 +width = 10.0 +arm_thickness = 1.0 +crack_length = 10.0 +opening = 10 + +fd.ModelingSpace("3D") + + +############################################################################### +# Mesh the two arms +# ~~~~~~~~~~~~~~~~ +# The lower and upper meshes deliberately retain separate nodes at z=0. +# Coincident-but-distinct nodes are necessary because the cohesive element +# interpolates the displacement jump between the two interface faces. +lower = fd.mesh.box_mesh( + nx, + ny, + nz_per_arm, + x_min=0.0, + x_max=length, + y_min=0.0, + y_max=width, + z_min=-arm_thickness, + z_max=0.0, + elm_type="hex8", +) +upper = fd.mesh.box_mesh( + nx, + ny, + nz_per_arm, + x_min=0.0, + x_max=length, + y_min=0.0, + y_max=width, + z_min=0.0, + z_max=arm_thickness, + elm_type="hex8", +) + +# Mesh.stack concatenates the two node lists without merging coincident nodes. +domain = fd.Mesh.stack(lower, upper, name="DCB_domain") +upper_offset = lower.n_nodes + +# Define one volume mesh per arm using the common global node list. Separate +# meshes make it possible to assign distinct materials or ply orientations to +# the arms without altering the interface node numbering. +lower_volume = fd.Mesh( + domain.nodes, + lower.elements, + "hex8", + ndim=3, +) +upper_volume = fd.Mesh( + domain.nodes, + upper.elements + upper_offset, + "hex8", + ndim=3, +) + + +############################################################################### +# Construct the cohesive interface +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# For hex8 elements, nodes 4:8 form the upper face and nodes 0:4 form the +# lower face. Pair the top face of the lower arm with the bottom face of the +# upper arm. A quad4interface element therefore contains eight nodes: +# +# [four nodes on the lower arm, four nodes on the upper arm] +# +n_interface_cells = (nx - 1) * (ny - 1) +lower_interface_faces = lower.elements[-n_interface_cells:, 4:8] +upper_interface_faces = upper.elements[:n_interface_cells, 0:4] + upper_offset +interface_elements = np.hstack((lower_interface_faces, upper_interface_faces)) + +# Do not insert cohesive elements over the starter-crack region. The two arms +# are consequently disconnected for x < crack_length and bonded by cohesive +# elements for x >= crack_length. +interface_centers_x = domain.nodes[interface_elements].mean(axis=1)[:, 0] +interface_elements = interface_elements[interface_centers_x >= crack_length - 1.0e-12] + +interface = fd.Mesh( + domain.nodes, + interface_elements, + "quad4interface", + ndim=3, + name="DCB_cohesive_interface", +) + + +############################################################################### +# Constitutive laws and assemblies +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# Both unidirectional composite arms have fibres along the beam axis (X). +# Separate objects are retained so either arm can easily be changed later. +lower_material = fd.constitutivelaw.CompositeUD( + Vf=0.6, + E_f=250000.0, + E_m=3500.0, + nu_f=0.33, + nu_m=0.3, + angle=0.0, +) +upper_material = fd.constitutivelaw.CompositeUD( + Vf=0.6, + E_f=250000.0, + E_m=3500.0, + nu_f=0.33, + nu_m=0.3, + angle=0.0, +) + +# Bilinear mixed-mode cohesive law. Axis 2 is the local interface-normal +# direction; the quad4interface element supplies the corresponding local frame. +cohesive_law = fd.constitutivelaw.CohesiveLaw( + GIc=0.3, + SImax=60.0, + KI=1.0e4, + GIIc=1.6, + SIImax=None, + KII=1.0e4, + tangent_mode="secant", +) + +lower_assembly = fd.Assembly.create( + fd.weakform.StressEquilibrium(lower_material), + lower_volume, +) +upper_assembly = fd.Assembly.create( + fd.weakform.StressEquilibrium(upper_material), + upper_volume, +) +cohesive_assembly = fd.Assembly.create( + fd.weakform.InterfaceForce(cohesive_law), + interface, +) + +# Keep a separate sum of the two volume assemblies for post-processing. It +# contains only standard hex8 meshes and can therefore be converted directly +# to a PyVista ``MultiMesh``. The cohesive assembly is then added only to the +# global mechanical assembly. +volume_assembly = lower_assembly + upper_assembly +assembly = volume_assembly + cohesive_assembly + + +############################################################################### +# Problem and boundary conditions +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# Material softening makes the problem nonlinear even though large-deformation +# kinematics are disabled. +problem = fd.problem.NonLinear(assembly, nlgeom=False) + +# Load the end faces at x=0. The lower end is fully fixed; the upper end is +# prevented from moving in X and Y and receives the prescribed opening in Z. +lower_left = np.nonzero( + np.isclose(domain.nodes[:, 0], 0.0) + & (domain.nodes[:, 2] < -arm_thickness + 1.0e-12) +)[0] +upper_left = np.nonzero( + np.isclose(domain.nodes[:, 0], 0.0) & (domain.nodes[:, 2] > arm_thickness - 1.0e-12) +)[0] + +problem.bc.add("Dirichlet", lower_left, "Disp", 0.0) +problem.bc.add("Dirichlet", upper_left, "DispX", 0.0) +problem.bc.add("Dirichlet", upper_left, "DispY", 0.0) +problem.bc.add("Dirichlet", upper_left, "DispZ", opening) + + +############################################################################### +# Incremental nonlinear solution +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# The imposed opening is scaled by the problem time from zero to its final +# value. Adaptive increments help the Newton solver cross the damage events. +# Convergence may be slow for this kind of problem. During damage propagation, +# the residual may temporarily increase over several iterations before +# convergence, so the early-divergence check is disabled. +problem.set_nr_criterion(check_early_divergence=False) +problem.nlsolve( + dt=0.05, + tmax=1.0, + update_dt=True, + dt_min=1.0e-5, + tol_nr=5.0e-3, + print_info=0, +) + + +############################################################################### +# Cohesive-zone response +# ~~~~~~~~~~~~~~~~~~~~~~ +# Damage and relative displacement are stored at the interface Gauss points. +# A damage value of one denotes a fully failed cohesive point. +damage = np.asarray(cohesive_assembly.sv["DamageVariable"]) +relative_disp = np.asarray(cohesive_assembly.sv["RelativeDisp"]) + +print("\nDCB cohesive-zone summary") +print("-------------------------") +print(f"n_nodes: {domain.n_nodes}") +print("n_volume_elements: " f"{lower_volume.n_elements + upper_volume.n_elements}") +print(f"n_cohesive_elements: {interface.n_elements}") +print(f"max_damage: {damage.max()}") +print(f"damaged_gauss_points: {np.count_nonzero(damage > 0.0)}") +print("failed_gauss_points: " f"{np.count_nonzero(damage >= 1.0 - 1.0e-10)}") +print(f"max_opening: {relative_disp[2].max()}") + + +############################################################################### +# Plot the deformed DCB +# ~~~~~~~~~~~~~~~~~~~~~ +# Request results from the volume sum explicitly. Calling +# problem.get_results(...) without an assembly would select the summed +# assembly, which also contains the eight-node ``quad4interface`` elements. +# These elements are valid for the mechanical calculation but have no direct +# VTK/PyVista cell equivalent. +results = problem.get_results(volume_assembly, ["Stress", "Strain", "Disp"]) + +# The AssemblySum result is a single MultiMesh dataset containing both arms. +# Element fields are retained independently on each submesh, whereas a +# recovered nodal stress field would be treated as one shared global field and +# would keep only the first arm's recovery. Displacement remains a shared nodal +# field and is used to display the final deformed geometry. +plotter = results.plot( + "Stress", + component="vm", + data_type="Node", + show_edges=True, + show=False, + title="DCB - von Mises stress", + scalar_bar_args={"interactive": False}, +) + +# Use an isometric camera from the negative-X side. This reverses the visual +# direction of the beam axis without changing the mesh, crack position, or +# boundary conditions. +plotter.view_vector((-1, 1, 1), viewup=(0, 0, 1)) +plotter.reset_camera() +plotter.show() diff --git a/examples/FE2/test_fe2.py b/examples/FE2/test_fe2.py new file mode 100644 index 00000000..2a7d004a --- /dev/null +++ b/examples/FE2/test_fe2.py @@ -0,0 +1,51 @@ +# FE2 light test file + +import numpy as np + +import fedoo as fd +from fedoo.util.voigt_tensors import StrainTensorList + + +def test_fe2_uses_problem_level_mean_strain_dofs(): + fd.ModelingSpace("3D") + + micro_mesh = fd.mesh.box_mesh(nx=2, ny=2, nz=2) + micro_material = fd.constitutivelaw.ElasticIsotrop(200_000.0, 0.3) + micro_weakform = fd.weakform.StressEquilibrium(micro_material) + micro_assembly = fd.Assembly.create(micro_weakform, micro_mesh, n_elm_gp=1) + + fe2 = fd.constitutivelaw.FE2(micro_assembly) + macro_mesh = fd.mesh.box_mesh(nx=2, ny=2, nz=2) + macro_weakform = fd.weakform.StressEquilibrium(fe2) + macro_assembly = fd.Assembly.create(macro_weakform, macro_mesh, n_elm_gp=1) + problem = fd.problem.NonLinear(macro_assembly) + + problem.initialize() + + assert "_StrainNodes" not in micro_mesh.node_sets + assert "MeanStrain" in fe2.list_problem[0].global_dof._vector + expected = micro_material.get_elastic_matrix() + np.testing.assert_allclose( + macro_assembly.sv["TangentMatrix"][:, :, 0], + expected, + rtol=1e-10, + atol=1e-8, + ) + + imposed_strain = np.array([1e-4, 0, 0, 0, 0, 0]) + macro_assembly.sv["Strain"] = StrainTensorList(imposed_strain[:, None]) + problem.dtime = 1.0 + fe2._update_pb(0, macro_assembly, problem) + + np.testing.assert_allclose( + macro_assembly.sv["Stress"].asarray()[:, 0], + expected @ imposed_strain, + rtol=1e-8, + atol=1e-8, + ) + + +if __name__ == "__main__": + import pytest + + pytest.main([__file__]) diff --git a/examples/plasticity/plastic_bending_3D.py b/examples/plasticity/plastic_bending_3D.py index 4c01db63..b2aa7369 100644 --- a/examples/plasticity/plastic_bending_3D.py +++ b/examples/plasticity/plastic_bending_3D.py @@ -59,7 +59,7 @@ material = fd.constitutivelaw.ElastoPlasticity( E, nu, Re, name="constitutivelaw" ) - material.SetHardeningFunction("power", H=k, beta=m) + material.set_hardening_function("power", h=k, beta=m) else: material = fd.constitutivelaw.ElasticIsotrop(E, nu, name="constitutivelaw") diff --git a/examples/plasticity/rotation_test.py b/examples/plasticity/rotation_test.py index fff0d114..90f3ba03 100644 --- a/examples/plasticity/rotation_test.py +++ b/examples/plasticity/rotation_test.py @@ -56,7 +56,7 @@ material = fd.constitutivelaw.ElastoPlasticity( E, nu, Re, name="ConstitutiveLaw" ) - material.SetHardeningFunction("power", H=k, beta=m) + material.set_hardening_function("power", h=k, beta=m) else: material = fd.constitutivelaw.ElasticIsotrop(E, nu, name="ConstitutiveLaw") diff --git a/examples/plasticity/shear_test.py b/examples/plasticity/shear_test.py index 3957519b..5df99ccf 100644 --- a/examples/plasticity/shear_test.py +++ b/examples/plasticity/shear_test.py @@ -55,7 +55,7 @@ material = fd.constitutivelaw.ElastoPlasticity( E, nu, Re, name="ConstitutiveLaw" ) - material.SetHardeningFunction("power", H=k, beta=m) + material.set_hardening_function("power", h=k, beta=m) elif mat == 3: props = np.array([3.0, 0.5e2]) material = fd.constitutivelaw.Simcoon("NEOHC", props, name="ConstitutiveLaw") diff --git a/examples/plasticity/torsion_test.py b/examples/plasticity/torsion_test.py index 0c8b17b0..e6443b76 100644 --- a/examples/plasticity/torsion_test.py +++ b/examples/plasticity/torsion_test.py @@ -52,7 +52,7 @@ material = fd.constitutivelaw.ElastoPlasticity( E, nu, Re, name="ConstitutiveLaw" ) - material.SetHardeningFunction("power", H=k, beta=m) + material.set_hardening_function("power", h=k, beta=m) else: material = fd.constitutivelaw.ElasticIsotrop(E, nu, name="ConstitutiveLaw") diff --git a/examples/shell_elements/LaminateExample.py b/examples/shell_elements/LaminateExample.py index 905cc6e1..1bf257ba 100644 --- a/examples/shell_elements/LaminateExample.py +++ b/examples/shell_elements/LaminateExample.py @@ -79,7 +79,7 @@ pb.save_results() # save in vtk # plot the stress distribution -z, StressDistribution = fd.ConstitutiveLaw["PlateSection"].GetStressDistribution( +z, StressDistribution = fd.ConstitutiveLaw["PlateSection"].get_stress_distribution( fd.Assembly["plate"], 200 ) plt.plot(StressDistribution[0], z) diff --git a/examples/shell_elements/simplePlateElementExample.py b/examples/shell_elements/simplePlateElementExample.py index 686413c7..576ad6bf 100644 --- a/examples/shell_elements/simplePlateElementExample.py +++ b/examples/shell_elements/simplePlateElementExample.py @@ -81,7 +81,7 @@ z, StressDistribution = fd.ConstitutiveLaw.get_all()[ "PlateSection" -].GetStressDistribution(assemb, 20) +].get_stress_distribution(assemb, 20) plt.plot(StressDistribution[0], z) # plot the von mises stress on the upper face (position = +1) diff --git a/examples/shell_elements/spherical_shell_compression_plastic.py b/examples/shell_elements/spherical_shell_compression_plastic.py new file mode 100644 index 00000000..342146e3 --- /dev/null +++ b/examples/shell_elements/spherical_shell_compression_plastic.py @@ -0,0 +1,97 @@ +"""Small elastoplastic compression test of a ping-pong ball. + +This reduced example exercises ``ShellHomogeneousNonLinear`` with the +Simcoon ``EPICP`` plane-stress constitutive law. The shell kinematics are +corotational, while strains passed to the material points remain small. +""" + +import numpy as np +import pyvista as pv + +import fedoo as fd + + +def run_simulation( + mesh_resolution=8, + pressure=2.0, + n_steps=2, + print_info=1, +): + """Run a deliberately small pressure-compression simulation.""" + young_modulus = 2_000.0 # MPa + poisson_ratio = 0.37 + yield_stress = 35.0 # MPa + hardening_modulus = 100.0 # MPa + hardening_exponent = 0.3 + radius = 20.0 # mm + thickness = 0.45 # mm + + fd.ModelingSpace("3D") + sphere = pv.Sphere( + radius=radius, + theta_resolution=mesh_resolution, + phi_resolution=mesh_resolution, + ) + mesh = fd.Mesh.from_pyvista(sphere) + + properties = np.array( + [ + young_modulus, + poisson_ratio, + 0.0, + yield_stress, + hardening_modulus, + hardening_exponent, + ] + ) + material = fd.constitutivelaw.Simcoon("EPICP", properties) + material.tangent_mode = 1 + shell = fd.constitutivelaw.ShellHomogeneousNonLinear( + material, + thickness, + n_thickness_points=3, + k=5 / 6, + ) + + weakform = fd.weakform.PlateEquilibrium(shell, nlgeom=True) + assembly = fd.Assembly.create(weakform, mesh) + + loaded_elements = mesh.find_elements( + f"Z>{mesh.bounding_box.zmax-3} or " f"Z<{mesh.bounding_box.zmin+3}" + ) + pressure_load = fd.constraint.Pressure( + mesh.extract_elements(loaded_elements), + pressure, + ) + + problem = fd.problem.NonLinear(assembly, nlgeom=True) + problem.bc.add(pressure_load) + + nodes = mesh.nodes + node_a = int(np.argmin(nodes[:, 0])) + node_b = int(np.argmax(nodes[:, 0])) + node_c = int(np.argmax(nodes[:, 1])) + problem.bc.add("Dirichlet", node_a, "Disp", 0) + problem.bc.add("Dirichlet", node_b, ["DispY", "DispZ"], 0) + problem.bc.add("Dirichlet", node_c, "DispZ", 0) + problem.set_nr_criterion( + "Displacement", + tol=5e-3, + max_subiter=15, + adaptive_stiffness=True, + ) + problem.nlsolve( + dt=1 / n_steps, + tmax=1, + update_dt=False, + print_info=print_info, + ) + return problem, assembly, shell + + +if __name__ == "__main__": + problem, assembly, shell = run_simulation() + displacement = problem.get_disp() + deformed = assembly.mesh.to_pyvista() + deformed.points += displacement.T + deformed.plot(show_edges=True) diff --git a/examples/thermal/thermo_meca_weak_coupling_3D.py b/examples/thermal/thermo_meca_weak_coupling_3D.py index bcd1aab1..9ab5d3a7 100644 --- a/examples/thermal/thermo_meca_weak_coupling_3D.py +++ b/examples/thermal/thermo_meca_weak_coupling_3D.py @@ -74,7 +74,7 @@ mechancial_law = fd.constitutivelaw.ElastoPlasticity( E, nu, Re, name="MechanicalLaw" ) - mechancial_law.SetHardeningFunction("power", H=k, beta=m) + mechancial_law.set_hardening_function("power", h=k, beta=m) else: mechancial_law = fd.constitutivelaw.ElasticIsotrop(E, nu, name="MechanicalLaw") diff --git a/fedoo/constitutivelaw/__init__.py b/fedoo/constitutivelaw/__init__.py index 2949895c..983cff07 100644 --- a/fedoo/constitutivelaw/__init__.py +++ b/fedoo/constitutivelaw/__init__.py @@ -66,7 +66,9 @@ :template: custom-class-template.rst ShellLaminate + ShellLaminateNonLinear ShellHomogeneous + ShellHomogeneousNonLinear Thermal constitutive law ====================================== @@ -83,7 +85,6 @@ from .beam import BeamCircular, BeamPipe, BeamProperties, BeamRectangular from .cohesivelaw import CohesiveLaw -from .cohesivelaw_mod import CohesiveLaw_mod from .composite_ud import CompositeUD from .elastic_anisotropic import ElasticAnisotropic from .elastic_isotrop import ElasticIsotrop @@ -91,13 +92,18 @@ from .elasto_plasticity import ElastoPlasticity from .fe2 import FE2 from .heterogeneous import Heterogeneous -from .shell import ShellBase, ShellHomogeneous, ShellLaminate +from .shell import ( + ShellBase, + ShellHomogeneous, + ShellHomogeneousNonLinear, + ShellLaminate, + ShellLaminateNonLinear, +) from .permeability import HolmesMowPermeability, KozenyCarmanPermeability from .poro_fluid import PoroFluidProperties from .simcoon_umat import Simcoon from .spring import Spring from .thermal_prop import ThermalProperties -from .viso_elastic_orthotropic import ViscoElasticComposites __all__ = [ "BeamCircular", @@ -105,7 +111,6 @@ "BeamProperties", "BeamRectangular", "CohesiveLaw", - "CohesiveLaw_mod", "CompositeUD", "ElasticAnisotropic", "ElasticIsotrop", @@ -115,12 +120,13 @@ "Heterogeneous", "ShellBase", "ShellHomogeneous", + "ShellHomogeneousNonLinear", "ShellLaminate", + "ShellLaminateNonLinear", "HolmesMowPermeability", "KozenyCarmanPermeability", "PoroFluidProperties", "Simcoon", "Spring", "ThermalProperties", - "ViscoElasticComposites", ] diff --git a/fedoo/constitutivelaw/cohesivelaw.py b/fedoo/constitutivelaw/cohesivelaw.py index 52e3ff42..0a0c4346 100644 --- a/fedoo/constitutivelaw/cohesivelaw.py +++ b/fedoo/constitutivelaw/cohesivelaw.py @@ -4,7 +4,6 @@ from fedoo.core.base import ConstitutiveLaw from fedoo.core.base import AssemblyBase import numpy as np -from numpy import linalg class CohesiveLaw(Spring): @@ -34,11 +33,25 @@ class CohesiveLaw(Spring): The axis is defined in local coordinate system. name: str, optional The name of the constitutive law + tangent_mode: {"secant", "consistent"}, default="secant" + Tangent stiffness used by the nonlinear solver. The secant tangent uses + the current damaged stiffness and is the more robust default. The + consistent tangent includes the derivative of damage during active + loading and can be selected explicitly. """ # Use with WeakForm.InterfaceForce def __init__( - self, GIc=0.3, SImax=60, KI=1e4, GIIc=1.6, SIImax=None, KII=5e4, axis=2, name="" + self, + GIc=0.3, + SImax=60, + KI=1e4, + GIIc=1.6, + SIImax=None, + KII=5e4, + axis=2, + name="", + tangent_mode="secant", ): # GIc la ténacité (l'énergie à la rupture = l'aire sous la courbe du modèle en N/mm) # SImax = 60. # la contrainte normale maximale de l'interface (MPa) @@ -55,6 +68,14 @@ def __init__( # axis = axe dans le repère local perpendiculaire au plan (will be deprecated because = 2 by convention) ConstitutiveLaw.__init__(self, name) # heritage + if axis not in (0, 1, 2): + raise ValueError("axis must be 0, 1, or 2.") + + tangent_mode = tangent_mode.lower() + if tangent_mode not in ("consistent", "secant"): + raise ValueError("tangent_mode must be either 'consistent' or 'secant'.") + self.tangent_mode = tangent_mode + self.parameters = { "GIc": GIc, "SImax": SImax, @@ -65,6 +86,26 @@ def __init__( "axis": axis, } + delta_0_I = SImax / KI + delta_m_I = 2 * GIc / SImax + if SIImax is None: + SIImax_check = SImax * np.sqrt(GIIc / GIc) + else: + SIImax_check = SIImax + delta_0_II = SIImax_check / KII + delta_m_II = 2 * GIIc / SIImax_check + if delta_m_I <= delta_0_I or delta_m_II <= delta_0_II: + raise ValueError( + "Cohesive parameters must satisfy delta_m > delta_0 " + "in both fracture modes." + ) + + # Mode-onset / failure separations are constant after construction; + # cache them so the damage update does not recompute them per Gauss point. + self._delta_0_I = delta_0_I + self._delta_0_II = delta_0_II + self._delta_m_II = delta_m_II + def initialize(self, assembly, pb): assembly.sv["InterfaceStress"] = 0 # Interface Stress assembly.sv["DamageVariable"] = 0 # damage variable @@ -76,7 +117,8 @@ def initialize(self, assembly, pb): ) assembly.sv["TangentMatrix"] = self.get_K(assembly) - def get_tangent_matrix(self, assembly): + def get_secant_matrix(self, assembly): + """Return the current damaged secant stiffness in local coordinates.""" Umd = 1 - assembly.sv["DamageVariable"] UmdI = 1 - assembly.sv["DamageVariableOpening"] @@ -87,16 +129,45 @@ def get_tangent_matrix(self, assembly): Kdiag = [Kt if i != axis else Kn for i in range(3)] return [[Kdiag[0], 0, 0], [0, Kdiag[1], 0], [0, 0, Kdiag[2]]] - # return [[Kdiag[0], 0, 0], [0, Kdiag[1], 0], [0,0,Kdiag[2]]] - # if get_Dimension() == "3D": # tester si marche avec contrainte plane ou def plane - # Kdiag = [Umd*self.parameters['KII'] if i != axis else UmdI*self.parameters['KI'] for i in range(3)] - # return [[Kdiag[0], 0, 0], [0, Kdiag[1], 0], [0,0,Kdiag[2]]] - # else: - # Kdiag = [Umd*self.parameters['KII'] if i != axis else UmdI*self.parameters['KI'] for i in range(2)] - # return [[Kdiag[0], 0], [0, Kdiag[1]]] + def get_tangent_matrix(self, assembly, delta=None, damage_gradient=None): + """Return the selected tangent stiffness in local coordinates. + + The consistent correction is active only while damage grows beyond + its committed value. During unloading, reloading below the historical + maximum, and after complete failure, ``damage_gradient`` is zero and + this method returns the damaged secant stiffness. + """ + secant = self.get_secant_matrix(assembly) + if self.tangent_mode == "secant" or delta is None: + return secant + + if damage_gradient is None: + _, _, damage_gradient = self._compute_damage(assembly, delta) - def get_K(self, assembly): - return self.local2global_K(self.get_tangent_matrix(assembly)) + delta = np.asarray(np.broadcast_arrays(*delta), dtype=float) + n_components = delta.shape[0] + opening = delta[self.parameters["axis"]] > 0 + axis = self.parameters["axis"] + stiffness = [ + self.parameters["KII"] if i != axis else self.parameters["KI"] + for i in range(n_components) + ] + + tangent = [[None for _ in range(n_components)] for _ in range(n_components)] + for i in range(n_components): + row_gradient = damage_gradient + if i == axis: + # Normal contact remains elastic in compression. + row_gradient = damage_gradient * opening + for j in range(n_components): + correction = stiffness[i] * delta[i] * row_gradient[j] + tangent[i][j] = secant[i][j] - correction + return tangent + + def get_K(self, assembly, delta=None, damage_gradient=None): + return self.local2global_K( + self.get_tangent_matrix(assembly, delta, damage_gradient) + ) def set_damage(self, assembly, value, irreversible=True): """ @@ -125,11 +196,15 @@ def update_damage(self, assembly, U, irreversible=False, type_data="PG"): if isinstance(assembly, str): assembly = AssemblyBase.get_all()[assembly] - op_delta = assembly.space.op_disp() + # In an updated-Lagrangian analysis, the current assembly carries the + # deformed interface geometry and its updated local coordinate system. + # For a geometrically linear analysis, assembly.current is assembly. + result_assembly = getattr(assembly, "current", assembly) + op_delta = result_assembly.space.op_disp() if type_data == "Node": - delta = [assembly.get_node_results(op, U) for op in op_delta] + delta = [result_assembly.get_node_results(op, U) for op in op_delta] else: - delta = [assembly.get_gp_results(op, U) for op in op_delta] + delta = [result_assembly.get_gp_results(op, U) for op in op_delta] self._update_damage(assembly, delta) @@ -138,168 +213,216 @@ def update_damage(self, assembly, U, irreversible=False, type_data="PG"): "DamageVariable" ].copy() - def _update_damage(self, assembly, delta): - alpha = 2 # for the power low - if ( - np.isscalar(assembly.sv["DamageVariable"]) - and assembly.sv["DamageVariable"] == 0 - ): - assembly.sv["DamageVariable"] = 0 * delta[0] - if ( - np.isscalar(assembly.sv["DamageVariableOpening"]) - and assembly.sv["DamageVariableOpening"] == 0 - ): - assembly.sv["DamageVariableOpening"] = 0 * delta[0] + def _compute_damage(self, assembly, delta): + """Compute damage and its derivative with respect to separation. - delta_n = delta[self.parameters["axis"]] - delta_t = [d for i, d in enumerate(delta) if i != self.parameters["axis"]] - if len(delta_t) == 1: - delta_t = delta_t[0] - else: - delta_t = np.sqrt(delta_t[0] ** 2 + delta_t[1] ** 2) - - # mode I - delta_0_I = ( - self.parameters["SImax"] / self.parameters["KI"] - ) # critical relative displacement (begining of the damage) - delta_m_I = ( - 2 * self.parameters["GIc"] / self.parameters["SImax"] - ) # maximal relative displacement (total failure) - - # mode II - SIImax = self.parameters["SIImax"] - if SIImax == None: - SIImax = self.parameters["SImax"] * np.sqrt( - self.parameters["GIIc"] / self.parameters["GIc"] - ) # value by default used mainly to treat mode I dominant problems - delta_0_II = SIImax / self.parameters["KII"] - delta_m_II = 2 * self.parameters["GIIc"] / SIImax - - t0 = 0.0 * delta_n - tm = 0.0 * delta_n - dta = 0.0 * delta_n - - test = delta_n > 0 # test if traction loading (opening mode) - ind_traction = np.nonzero(test)[ - 0 - ] # indice of value where delta_n > 0 ie traction loading - ind_compr = np.nonzero(test - 1)[0] - - # beta = 0*delta_n - beta = ( - delta_t[ind_traction] / delta_n[ind_traction] - ) # le rapport de mixité de mode - - t0[ind_traction] = (delta_0_II * delta_0_I) * ( - np.sqrt((1 + (beta**2)) / ((delta_0_II**2) + ((beta * delta_0_I) ** 2))) - ) # Critical relative displacement in mixed mode - t0[ind_compr] = ( - delta_0_II # Critical relative displacement in mixed mode (only mode II) - ) + The derivative is the algorithmic derivative: it is nonzero only when + trial damage exceeds the irreversible damage stored at the start of the + increment. Unloading and elastic reloading therefore use the damaged + secant stiffness. + """ + alpha = 2.0 + delta = np.asarray(np.broadcast_arrays(*delta), dtype=float) + n_components = delta.shape[0] + if n_components not in (2, 3): + raise ValueError( + "CohesiveLaw expects two (2D) or three (3D) separation components." + ) + axis = self.parameters["axis"] + if axis >= n_components: + raise ValueError( + f"CohesiveLaw normal 'axis'={axis} is out of range for " + f"{n_components} separation components; set axis to 0 or 1 in 2D." + ) - tm[ind_traction] = (2 * ((1 + beta) ** 2) / t0[ind_traction]) * ( - ( - ((self.parameters["KI"] / self.parameters["GIc"]) ** alpha) - + ( - ((self.parameters["KII"] * beta**2) / self.parameters["GIIc"]) - ** alpha - ) + state_shape = delta.shape[1:] + delta = delta.reshape(n_components, -1) + n_points = delta.shape[1] + tangential_axes = [i for i in range(n_components) if i != axis] + + delta_n = delta[axis] + delta_t_vector = delta[tangential_axes] + delta_t = np.sqrt(np.sum(delta_t_vector**2, axis=0)) + + delta_0_I = self._delta_0_I + delta_0_II = self._delta_0_II + delta_m_II = self._delta_m_II + + opening = delta_n > 0.0 + compression = ~opening + t0 = np.empty(n_points) + tm = np.empty(n_points) + equivalent_delta = np.empty(n_points) + beta = np.zeros(n_points) + t0_beta = np.zeros(n_points) + tm_beta = np.zeros(n_points) + + if np.any(opening): + beta[opening] = delta_t[opening] / delta_n[opening] + beta_opening = beta[opening] + + denominator = delta_0_II**2 + (beta_opening * delta_0_I) ** 2 + t0[opening] = ( + delta_0_II * delta_0_I * np.sqrt((1.0 + beta_opening**2) / denominator) + ) + + mode_I_term = self.parameters["KI"] / self.parameters["GIc"] + mode_II_term = ( + self.parameters["KII"] * beta_opening**2 / self.parameters["GIIc"] + ) + power_sum = mode_I_term**alpha + mode_II_term**alpha + tm[opening] = ( + 2.0 + * (1.0 + beta_opening) ** 2 + / t0[opening] + * power_sum ** (-1.0 / alpha) + ) + equivalent_delta[opening] = np.sqrt( + delta_t[opening] ** 2 + delta_n[opening] ** 2 + ) + + t0_beta[opening] = t0[opening] * ( + beta_opening / (1.0 + beta_opening**2) + - beta_opening * delta_0_I**2 / denominator + ) + power_sum_beta = ( + alpha + * mode_II_term ** (alpha - 1.0) + * (2.0 * self.parameters["KII"] * beta_opening) + / self.parameters["GIIc"] + ) + tm_beta[opening] = tm[opening] * ( + 2.0 / (1.0 + beta_opening) + - t0_beta[opening] / t0[opening] + - power_sum_beta / (alpha * power_sum) ) - ** (-1 / alpha) - ) # Maximal relative displacement in mixed mode (power low criterion) - tm[ind_compr] = ( - delta_m_II # Maximal relative displacement in mixed mode (power low criterion) - ) - dta[ind_traction] = np.sqrt( - delta_t[ind_traction] ** 2 + delta_n[ind_traction] ** 2 - ) # Actual relative displacement in mixed mode - dta[ind_compr] = delta_t[ - ind_compr - ] # Actual relative displacement in mixed mode + t0[compression] = delta_0_II + tm[compression] = delta_m_II + equivalent_delta[compression] = delta_t[compression] - # --------------------------------------------------------------------------------------------------------------- - # La variable d'endommagement "d" - # --------------------------------------------------------------------------------------------------------------- - d = (dta >= tm).astype(float) # initialize d to 1 if dta>tm and else d=0 - test = np.nonzero((dta > t0) * (dta < tm))[ - 0 - ] # indices where dta>t0 and dta= tm).astype(float) + softening = (equivalent_delta > t0) & (equivalent_delta < tm) + softening_factor = np.zeros(n_points) + softening_factor[softening] = tm[softening] / (tm[softening] - t0[softening]) + trial_damage[softening] = softening_factor[softening] * ( + 1.0 - t0[softening] / equivalent_delta[softening] + ) - d[test] = (tm[test] / (tm[test] - t0[test])) * (1 - (t0[test] / dta[test])) + irreversible = np.asarray( + assembly.sv["DamageVariableIrreversible"], dtype=float + ) + irreversible = np.broadcast_to(irreversible, state_shape).reshape(-1) + damage = np.maximum(irreversible, trial_damage) + + # max(d_irreversible, d_trial) makes this derivative vanish unless + # trial damage is actively increasing. + active = softening & (trial_damage > irreversible) + damage_gradient = np.zeros((n_components, n_points)) + + active_compression = active & compression + if np.any(active_compression): + r = equivalent_delta[active_compression] + damage_r = ( + softening_factor[active_compression] * t0[active_compression] / r**2 + ) + for component, tangential_axis in enumerate(tangential_axes): + damage_gradient[tangential_axis, active_compression] = ( + damage_r * delta_t_vector[component, active_compression] / r + ) - if ( - np.isscalar(assembly.sv["DamageVariableIrreversible"]) - and assembly.sv["DamageVariableIrreversible"] == 0 - ): - assembly.sv["DamageVariable"] = ( - d # I don't know why assembly.sv['DamageVariable'] = d end up in bads values + active_opening = active & opening + if np.any(active_opening): + r = equivalent_delta[active_opening] + t0_active = t0[active_opening] + tm_active = tm[active_opening] + factor = softening_factor[active_opening] + + factor_beta = ( + tm_active * t0_beta[active_opening] + - tm_beta[active_opening] * t0_active + ) / (tm_active - t0_active) ** 2 + damage_r = factor * t0_active / r**2 + damage_beta = factor_beta * (1.0 - t0_active / r) - ( + factor * t0_beta[active_opening] / r ) - else: - assembly.sv["DamageVariable"] = np.max( - [assembly.sv["DamageVariableIrreversible"], d], axis=0 + + delta_n_active = delta_n[active_opening] + beta_active = beta[active_opening] + damage_gradient[axis, active_opening] = ( + damage_r * delta_n_active / r + - damage_beta * beta_active / delta_n_active ) - assembly.sv["DamageVariableOpening"] = ( - (delta_n > 0) * assembly.sv["DamageVariable"] - ) # for opening the damage in considered to 0 when the relative displacement is negative (conctact) + delta_t_active = delta_t[active_opening] + nonzero_tangent = delta_t_active > 0.0 + for component, tangential_axis in enumerate(tangential_axes): + gradient = damage_r * (delta_t_vector[component, active_opening] / r) + gradient[nonzero_tangent] += damage_beta[nonzero_tangent] * ( + delta_t_vector[component, active_opening][nonzero_tangent] + / ( + delta_t_active[nonzero_tangent] + * delta_n_active[nonzero_tangent] + ) + ) + damage_gradient[tangential_axis, active_opening] = gradient - # verification : the damage variable should be between 0 and 1 - if ( - assembly.sv["DamageVariable"].min() < 0 - or assembly.sv["DamageVariable"].max() > 1 - ): - print("Warning : the value of damage variable is incorrect") + return ( + damage.reshape(state_shape), + (opening * damage).reshape(state_shape), + damage_gradient.reshape((n_components,) + state_shape), + ) + + def _update_damage(self, assembly, delta): + damage, damage_opening, damage_gradient = self._compute_damage(assembly, delta) + assembly.sv["DamageVariable"] = damage + assembly.sv["DamageVariableOpening"] = damage_opening + return damage_gradient def reset(self): pass def set_start(self, assembly, pb): - # Set Irreversible Damage + # Commit damage at the end of the converged increment. At the start of + # the next increment, the exact algorithmic response is unloading from + # that state, so the predictor matrix must be the damaged secant rather + # than the negative active-softening tangent from the previous step. self.update_irreversible_damage(assembly) + assembly.sv["TangentMatrix"] = self.local2global_K( + self.get_secant_matrix(assembly) + ) # def to_start(self, assembly, pb): # #Damage variable will be recompute. NPOthing to be done here (to be checked) # pass - def update(self, assembly, pb): - displacement = pb.get_dof_solution() - K = self.get_K() - assembly.sv["TangentMatrix"] = K - if np.isscalar(displacement) and displacement == 0: - assembly.sv["InterfaceStress"] = assembly.sv["RelativeDisp"] = 0 - else: - op_delta = ( - assembly.space.op_disp() - ) # relative displacement = disp if used with cohesive element - delta = [assembly.get_gp_results(op, displacement) for op in op_delta] - assembly.sv["RelativeDisp"] = delta - - # Compute interface stress - dim = len(delta) - assembly.sv["InterfaceStress"] = [ - sum([delta[j] * K[i][j] for j in range(dim)]) for i in range(dim) - ] # list of 3 objects - def update(self, assembly, pb): displacement = pb.get_dof_solution() if np.isscalar(displacement) and displacement == 0: assembly.sv["InterfaceStress"] = assembly.sv["RelativeDisp"] = 0 - K = self.get_K() + K = self.get_K(assembly) else: - op_delta = ( - assembly.space.op_disp() - ) # relative displacement = disp if used with cohesive element - delta = [assembly.get_gp_results(op, displacement) for op in op_delta] + # InterfaceForce updates assembly.current before this constitutive + # update. Evaluate the jump with its current local frame so damage, + # traction, and the assembled tangent use the same configuration. + result_assembly = getattr(assembly, "current", assembly) + op_delta = result_assembly.space.op_disp() + delta = [ + result_assembly.get_gp_results(op, displacement) for op in op_delta + ] assembly.sv["RelativeDisp"] = delta - self._update_damage(assembly, delta) + damage_gradient = self._update_damage(assembly, delta) dim = len(delta) - K = self.get_K(assembly) + # Traction follows the secant constitutive relation. The + # consistent matrix is its derivative and is used only for the + # Newton Jacobian. + K_secant = self.local2global_K(self.get_secant_matrix(assembly)) assembly.sv["InterfaceStress"] = [ - sum([delta[j] * K[i][j] for j in range(dim)]) for i in range(dim) + sum([delta[j] * K_secant[i][j] for j in range(dim)]) for i in range(dim) ] # list of 3 objects + K = self.get_K(assembly, delta, damage_gradient) assembly.sv["TangentMatrix"] = K diff --git a/fedoo/constitutivelaw/cohesivelaw_mod.py b/fedoo/constitutivelaw/cohesivelaw_mod.py deleted file mode 100644 index 563f5a54..00000000 --- a/fedoo/constitutivelaw/cohesivelaw_mod.py +++ /dev/null @@ -1,398 +0,0 @@ -# derive de ConstitutiveLaw -# Not working law - -from fedoo.constitutivelaw.spring import Spring -from fedoo.core.base import ConstitutiveLaw -from fedoo.core.base import AssemblyBase -import numpy as np -from numpy import linalg - - -class CohesiveLaw_mod(Spring): - # Use with WeakForm.InterfaceForce - def __init__( - self, GIc=0.3, SImax=60, KI=1e4, GIIc=1.6, SIImax=None, KII=5e4, axis=2, name="" - ): - # GIc la ténacité (l'énergie à la rupture = l'aire sous la courbe du modèle en N/mm) - # SImax = 60. # la contrainte normale maximale de l'interface (MPa) - # KI = 1e4 # la raideur des éléments cohésive (la pente du modèle en N/mm3) - # - # # Mode II (12) - # G_IIc = 1.6 - # KII = 5e4 - - # - ##----------------------- la puissance du critère de propagation (cas de critère de Power Law)--------------------------- - # alpha = 2. - # - - ConstitutiveLaw.__init__(self, name) # heritage - self.__DamageVariable = 0 # damage variable - self.__DamageVariableOpening = 0 # DamageVariableOpening is used for the opening mode (mode I). It is equal to DamageVariable in traction and equal to 0 in compression (soft contact law) - self.__DamageVariableIrreversible = ( - 0 # irreversible damage variable used for time evolution - ) - self.__parameters = { - "GIc": GIc, - "SImax": SImax, - "KI": KI, - "GIIc": GIIc, - "SIImax": SIImax, - "KII": KII, - "axis": axis, - } - self.__currentInterfaceStress = None - - def GetKelas(self): # Get elastic rigidity in local coordinates - Umd = 1 - self.__DamageVariable - UmdI = 1 - self.__DamageVariableOpening - - axis = self.__parameters["axis"] - - Kt = Umd * self.__parameters["KII"] - Kn = UmdI * self.__parameters["KI"] - Kdiag = [Kt if i != axis else Kn for i in range(3)] - return [[Kdiag[0], 0, 0], [0, Kdiag[1], 0], [0, 0, Kdiag[2]]] - - # if get_Dimension() == "3D": # tester si marche avec contrainte plane ou def plane - # Kdiag = [Umd*self.__parameters['KII'] if i != axis else UmdI*self.__parameters['KI'] for i in range(3)] - # return [[Kdiag[0], 0, 0], [0, Kdiag[1], 0], [0,0,Kdiag[2]]] - # else: - # Kdiag = [Umd*self.__parameters['KII'] if i != axis else UmdI*self.__parameters['KI'] for i in range(2)] - # return [[Kdiag[0], 0], [0, Kdiag[1]]] - - def get_tangent_matrix(self): # Get tangent moduli - if self.__currentInterfaceStress is None: - return self.GetKelas() - - Umd = 1 - self.__DamageVariable - UmdI = 1 - self.__DamageVariableOpening - - KIelas = UmdI * self.__parameters["KI"] - KIIelas = Umd * self.__parameters["KII"] - - ### Cohesive Zones Data - # mode I - delta_0_I = ( - self.__parameters["SImax"] / self.__parameters["KI"] - ) # critical relative displacement (begining of the damage) - delta_m_I = ( - 2 * self.__parameters["GIc"] / self.__parameters["SImax"] - ) # maximal relative displacement (total failure) - - # mode II - SIImax = self.__parameters["SIImax"] - if SIImax == None: - SIImax = self.__parameters["SImax"] * np.sqrt( - self.__parameters["GIIc"] / self.__parameters["GIc"] - ) # value by default used mainly to treat mode I dominant problems - delta_0_II = SIImax / self.__parameters["KII"] - delta_m_II = 2 * self.__parameters["GIIc"] / SIImax - - # KItangent = -self.__parameters['SImax']/(delta_m_I-delta_0_I) - # KIItangent = -SIImax/(delta_m_II-delta_0_II) - KItangent = 0 - KIItangent = 0 - - # modeI : - test = (self.__DamageVariableOpening - self.__DamageVariableIrreversible) > 0 - test2 = np.logical_not(test) - test = np.logical_and(test, self.__DamageVariableOpening != 1) - # test2 = (self.__DamageVariable == 1) - KI = KIelas * test2 + KItangent * test - - # modeII : - test = (self.__DamageVariable - self.__DamageVariableIrreversible) > 0 - test2 = np.logical_not(test) - test = np.logical_and(test, self.__DamageVariable != 1) - KII = KIIelas * test2 + KIItangent * test - - axis = self.__parameters["axis"] - - Kdiag = [KII if i != axis else KI for i in range(3)] - return [[Kdiag[0], 0, 0], [0, Kdiag[1], 0], [0, 0, Kdiag[2]]] - - # if get_Dimension() == "3D": # tester si marche avec contrainte plane ou def plane - # Kdiag = [KII if i != axis else KI for i in range(3)] - # return [[Kdiag[0], 0, 0], [0, Kdiag[1], 0], [0,0,Kdiag[2]]] - # else: - # Kdiag = [KII if i != axis else KI for i in range(2)] - # return [[Kdiag[0], 0], [0, Kdiag[1]]] - - def set_DamageVariable(self, value): - self.__DamageVariable = value - - def get_DamageVariable(self): - return self.__DamageVariable - - def updateIrreversibleDamage(self): - if np.isscalar(self.__DamageVariable) and self.__DamageVariable == 0: - self.__DamageVariableIrreversible = 0 - else: - self.__DamageVariableIrreversible = self.__DamageVariable.copy() - - #### Not working, need update - def updateDamageVariable( - self, CohesiveAssembly, U, Irreversible=False, typeData="PG" - ): - # Delta is the relative displacement - # OperatorDelta = assembly.space.op_disp() #relative displacement = disp if used with cohesive element - # OperatorDelta, U_vir = get_DispOperator() - if isinstance(CohesiveAssembly, str): - CohesiveAssembly = AssemblyBase.get_all()[CohesiveAssembly] - if typeData == "Node": - delta = [CohesiveAssembly.get_node_results(op, U) for op in OperatorDelta] - else: - delta = [CohesiveAssembly.get_gp_results(op, U) for op in OperatorDelta] - - self.__UpdateDamageVariable(delta) - - if Irreversible == True: - self.__DamageVariableIrreversible = self.__DamageVariable.copy() - - def __UpdateDamageVariable(self, delta): - alpha = 2 # for the power low - if np.isscalar(self.__DamageVariable) and self.__DamageVariable == 0: - self.__DamageVariable = 0 * delta[0] - if ( - np.isscalar(self.__DamageVariableOpening) - and self.__DamageVariableOpening == 0 - ): - self.__DamageVariableOpening = 0 * delta[0] - - # delta_n = delta.pop(self.__parameters['axis']) - # if get_Dimension() == "3D": - # delta_t = np.sqrt(delta[0]**2 + delta[1]**2) - # else: delta_t = delta[0] - - delta_n = delta[self.__parameters["axis"]] - delta_t = [d for i, d in enumerate(delta) if i != self.__parameters["axis"]] - - if len(delta_t) == 1: - delta_t = delta_t[0] - else: - delta_t = np.sqrt(delta_t[0] ** 2 + delta_t[1] ** 2) - - if get_Dimension() == "3D": - delta_t = np.sqrt(delta_t[0] ** 2 + delta_t[1] ** 2) - else: - delta_t = delta_t[0] - - ### Cohesive Zones Data - # mode I - delta_0_I = ( - self.__parameters["SImax"] / self.__parameters["KI"] - ) # critical relative displacement (begining of the damage) - delta_m_I = ( - 2 * self.__parameters["GIc"] / self.__parameters["SImax"] - ) # maximal relative displacement (total failure) - - # mode II - SIImax = self.__parameters["SIImax"] - if SIImax == None: - SIImax = self.__parameters["SImax"] * np.sqrt( - self.__parameters["GIIc"] / self.__parameters["GIc"] - ) # value by default used mainly to treat mode I dominant problems - delta_0_II = SIImax / self.__parameters["KII"] - delta_m_II = 2 * self.__parameters["GIIc"] / SIImax - - # Compute mixed mode relative displacement (actual values and limit values) - - test = delta_n > 0 # test if traction loading (opening mode) - ind_traction = np.nonzero(test)[ - 0 - ] # indice of value where delta_n > 0 ie traction loading - ind_compr = np.nonzero(test - 1)[0] - - t0 = 0.0 * delta_n - tm = 0.0 * delta_n - dta = 0.0 * delta_n - - # beta = 0*delta_n - beta = ( - delta_t[ind_traction] / delta_n[ind_traction] - ) # le rapport de mixité de mode - - t0[ind_traction] = (delta_0_II * delta_0_I) * ( - np.sqrt((1 + (beta**2)) / ((delta_0_II**2) + ((beta * delta_0_I) ** 2))) - ) # Critical relative displacement in mixed mode - t0[ind_compr] = ( - delta_0_II # Critical relative displacement in mixed mode (only mode II) - ) - - tm[ind_traction] = (2 * ((1 + beta) ** 2) / t0[ind_traction]) * ( - ( - ((self.__parameters["KI"] / self.__parameters["GIc"]) ** alpha) - + ( - ((self.__parameters["KII"] * beta**2) / self.__parameters["GIIc"]) - ** alpha - ) - ) - ** (-1 / alpha) - ) # Maximal relative displacement in mixed mode (power low criterion) - tm[ind_compr] = ( - delta_m_II # Maximal relative displacement in mixed mode (power low criterion) - ) - - dta[ind_traction] = np.sqrt( - delta_t[ind_traction] ** 2 + delta_n[ind_traction] ** 2 - ) # Actual relative displacement in mixed mode - dta[ind_compr] = delta_t[ - ind_compr - ] # Actual relative displacement in mixed mode - - # --------------------------------------------------------------------------------------------------------------- - # La variable d'endommagement "d" - # --------------------------------------------------------------------------------------------------------------- - d = (dta >= tm).astype(float) # initialize d to 1 if dta>tm and else d=0 - test = np.nonzero((dta > t0) * (dta < tm))[ - 0 - ] # indices where dta>t0 and dta 0) * self.__DamageVariable - ) # for opening the damage in considered to 0 when the relative displacement is negative (conctact) - - # verification : the damage variable should be between 0 and 1 - if self.__DamageVariable.min() < 0 or self.__DamageVariable.max() > 1: - print("Warning : the value of damage variable is incorrect") - - def NewTimeIncrement(self): - # Set Irreversible Damage - self.updateIrreversibleDamage() - self.__currentSigma = None - - def to_start(self): - # Damage variable and currentInterfaceStress will be recomputed in the next call of GetInterfaceStress - self.__currentInterfaceStress = None - - def reset(self): - """ - reset the constitutive law (time history) - """ - self.__DamageVariable = 0 # damage variable - self.__DamageVariableOpening = 0 # DamageVariableOpening is used for the opening mode (mode I). It is equal to DamageVariable in traction and equal to 0 in compression (soft contact law) - self.__DamageVariableIrreversible = ( - 0 # irreversible damage variable used for time evolution - ) - - def GetInterfaceStress(self, Delta, time=None): - # Delta is the relative displacement vector - self.__UpdateDamageVariable(Delta) - self.__currentInterfaceStress = Spring.GetInterfaceStress(self, Delta, time) - return self.__currentInterfaceStress - - -# def SetLocalFrame(self, localFrame): -# raise NameError("Not implemented: localFrame are not implemented in the context of cohesive laws") - - -# def __UpdateDamageVariable_old(self, delta): -# #--------------------------------------------------------------------------------------------------------- -# ################# interface 90°/0° (Lower interface) ######################## -# #--------------------------------------------------------------------------------------------------------- - - -# alpha = 2 #for the power low -# if np.isscalar(self.__DamageVariable) and self.__DamageVariable == 0: self.__DamageVariable = 0*delta[0] -# if np.isscalar(self.__DamageVariableOpening) and self.__DamageVariableOpening == 0: self.__DamageVariableOpening = 0*delta[0] - -# # delta_n = delta.pop(self.__parameters['axis']) -# # delta_t = np.sqrt(delta[0]**2 + delta[1]**2) -# delta_n = delta[self.__parameters['axis']] -# delta_t = [d for i,d in enumerate(delta) if i != self.__parameters['axis'] ] -# if get_Dimension() == "3D": -# delta_t = np.sqrt(delta_t[0]**2 + delta_t[1]**2) -# else: delta_t = delta_t[0] - -# # mode I -# delta_0_I = self.__parameters['SImax'] / self.__parameters['KI'] # critical relative displacement (begining of the damage) -# delta_m_I = 2*self.__parameters['GIc'] / self.__parameters['SImax'] # maximal relative displacement (total failure) - -# # mode II -# SIImax = self.__parameters['SIImax'] -# if SIImax == None: SIImax = self.__parameters['SImax'] * np.sqrt(self.__parameters['GIIc'] / self.__parameters['GIc']) #value by default used mainly to treat mode I dominant problems -# delta_0_II = SIImax / self.__parameters['KII'] -# delta_m_II = 2*self.__parameters['GIIc'] / SIImax - -# for i in range (len(delta_n)): -# if delta_n[i] > 0 : -# beta= delta_t[i] / (delta_n[i]) # le rapport de mixité de mode - -# t0= (delta_0_II * delta_0_I) * (np.sqrt((1+ (beta**2)) / ((delta_0_II**2)+((beta*delta_0_I)**2)))) # Critical relative displacement in mixed mode - -# tm= (2*((1+ beta)**2)/t0) * ((((self.__parameters['KI']/self.__parameters['GIc'])**alpha) + \ -# (((self.__parameters['KII']*beta**2)/self.__parameters['GIIc'])**alpha))**(-1/alpha)) #Maximal relative displacement in mixed mode (power low criterion) - -# dta= np.sqrt(delta_t[i]**2 + delta_n[i]**2) # Actual relative displacement in mixed mode - -# else : #only mode II -# t0= delta_0_II # Critical relatie displacement in mixed mode -# print(delta_0_II) -# tm= delta_m_II # Maximal relative displacement in mixed mode (power low criterion) - -# dta= delta_t[i] # Actual relative displacement in mixed mode - -# #--------------------------------------------------------------------------------------------------------------- -# # La variable d'endommagement "d" -# #--------------------------------------------------------------------------------------------------------------- -# if dta <= t0: -# di = 0 -# elif dta > t0 and dta < tm: -# di = (tm / (tm - t0)) * (1 - (t0 / dta)) -# else: -# di = 1 - -# if (self.__DamageVariableIrreversible is 0) or (di > self.__DamageVariableIrreversible[i]): -# self.__DamageVariable[i] = di -# else: self.__DamageVariable[i] = self.__DamageVariableIrreversible[i] - -# if delta_n[i] > 0 : -# self.__DamageVariableOpening [i] = self.__DamageVariable[i] -# else : -# self.__DamageVariableOpening [i] = 0 - -# # verification : the damage variable should be between 0 and 1 -# if self.__DamageVariable.min() < 0 or self.__DamageVariable.max() > 1 : -# print ("Warning : the value of damage variable is incorrect") - - -# if __name__=="__main__": -# ModelingSpace("3D") -# GIc = 0.3 ; SImax = 60 -# # delta_I_max = 2*GIc/SImax -# delta_I_max = 0.04 -# nb_iter = 100 -# sig = [] -# delta_plot = [] -# law = CohesiveLaw(GIc=GIc, SImax = SImax, KI = 1e4, GIIc = 1, SIImax=60, KII=1e4, axis = 0) -# for delta_z in np.arange(0,delta_I_max,delta_I_max/nb_iter): -# delta = [np.array([0]), np.array([0]), np.array([delta_z])] -# sig.append(law.GetInterfaceStress(delta)[2]) -# law.updateIrreversibleDamage() -# delta_plot.append(delta_z) -# # print(law.get_DamageVariable()) - -# # for delta_z in np.arange(delta_I_max,-delta_I_max,-delta_I_max/nb_iter): -# # delta = [np.array([0]), np.array([0]), np.array([delta_z])] -# # sig.append(law.GetInterfaceStress(delta)[2]) -# # law.updateIrreversibleDamage() -# # delta_plot.append(delta_z) - -# import matplotlib.pyplot as plt - -# plt.plot(delta_plot, sig) diff --git a/fedoo/constitutivelaw/elasto_plasticity.py b/fedoo/constitutivelaw/elasto_plasticity.py index 5a80af7a..aa56e53c 100644 --- a/fedoo/constitutivelaw/elasto_plasticity.py +++ b/fedoo/constitutivelaw/elasto_plasticity.py @@ -1,376 +1,456 @@ -# derive de ConstitutiveLaw -# The elastoplastic law should be used with an InternalForce WeakForm - -from fedoo.core.mechanical3d import Mechanical3D -from fedoo.util.voigt_tensors import StressTensorList, StrainTensorList +"""Small-strain J2 plasticity with isotropic hardening.""" import numpy as np +from simcoon import Rotation as SimRotation +from fedoo.core.mechanical3d import Mechanical3D +from fedoo.util.voigt_tensors import StrainTensorList, StressTensorList -class ElastoPlasticity(Mechanical3D): - """ - Elasto-Plastic constitutive law. - - This law is based on the assumption of isotropic hardening with the - Von-Mises plasticity criterion. After creating an ElastoPlasticity object, - the hardening function must be set with the Method 'SetHardeningFunction' - This constitutive Law should be associated with - :mod:`fedoo.weakform.StressEquilibrium` - - Parameters - ---------- - YoungModulus: scalars or arrays of gauss point values - Young modulus - PoissonRatio: scalars or arrays of gauss point values - Poisson's Ratio - YieldStress: scalars or arrays of gauss point values - Yield Stress Value - name: str, optional - The name of the constitutive law - """ - def __init__(self, YoungModulus, PoissonRatio, YieldStress, name=""): - # only scalar values of YoungModulus and PoissonRatio are possible - Mechanical3D.__init__(self, name) # heritage +class ElastoPlasticity(Mechanical3D): + """Elasto-plastic constitutive law with isotropic hardening. - self.__YoungModulus = YoungModulus - self.__PoissonRatio = PoissonRatio - self.__YieldStress = YieldStress + The stress integration uses a vectorized radial-return algorithm for the + von Mises yield criterion. In a finite-strain updated-Lagrangian analysis, + the law operates on the corotated strain increment prepared by + :class:`fedoo.weakform.StressEquilibrium`. - self.__P = None # irrevesrible plasticity - self.__currentP = None # current iteration plasticity (reversible) - self.__PlasticStrainTensor = None - self.__currentPlasticStrainTensor = None - self.__currentSigma = 0 # lissStressTensor object describing the last computed stress (GetStress method) - self.__currentGradDisp = None - self.__TangeantModuli = None + The returned elastoplastic tangent is the continuum tangent evaluated at + the updated stress. It is not the algorithmically consistent tangent of + the discrete radial-return integration. Consequently, global Newton + convergence is not guaranteed to be quadratic for finite plastic + increments or non-proportional loading paths. - self.__tol = 1e-6 # tolerance of Newton Raphson used to get the updated plasticity state (constutive law alogorithm) - self.nlgeom = False # will be updated with the Initialize function + .. warning:: - def GetYoungModulus(self): - return self.__YoungModulus + This Python implementation is intended only for pedagogical use and + as a readable reference implementation. For computational analyses, + prefer :class:`fedoo.constitutivelaw.Simcoon`, which provides optimized + constitutive updates and algorithmically consistent tangent options. - def GetPoissonRatio(self): - return self.__PoissonRatio + Parameters + ---------- + young_modulus : float + Young modulus. + poisson_ratio : float + Poisson ratio. + yield_stress : float + Initial yield stress. + name : str, optional + Name of the constitutive law. + """ - def GetYieldStress(self): - return self.__YieldStress + def __init__( + self, + young_modulus, + poisson_ratio, + yield_stress, + name="", + ): + super().__init__(name) + self.young_modulus = young_modulus + self.poisson_ratio = poisson_ratio + self.yield_stress = yield_stress + self.return_mapping_tolerance = 1e-6 + self.max_return_mapping_iterations = 50 + + self._hardening_function = None + self._hardening_function_derivative = None + self._current_stress = None + self._current_plastic_strain = None + self._current_plasticity = None + self._current_tangent = None + + @property + def shear_modulus(self): + """Shear modulus.""" + return self.young_modulus / (2.0 * (1.0 + self.poisson_ratio)) + + def set_return_mapping_tolerance(self, tolerance): + """Set the absolute tolerance used by the local return mapping.""" + if tolerance <= 0: + raise ValueError("The return-mapping tolerance must be positive.") + self.return_mapping_tolerance = tolerance + + def get_elastic_matrix(self, dimension="3D"): + """Return the isotropic elastic matrix in engineering Voigt form.""" + if dimension == "2Dstress": + raise NotImplementedError( + "ElastoPlasticity does not implement the plane-stress " + "('2Dstress') return mapping. For plane-stress models or " + "through-thickness shell plasticity, use a Simcoon plasticity " + 'law, e.g. fedoo.constitutivelaw.Simcoon("EPICP", props), ' + "which supports the 2Dstress assumption." + ) - def SetNewtonRaphsonTolerance(self, tol): - """ - Set the tolerance of the Newton Raphson algorithm used to get the updated plasticity state (constutive law alogorithm) - """ - self.__tol = tol - - def GetHelas(self): - H = np.zeros((6, 6), dtype="object") - E = self.__YoungModulus - nu = self.__PoissonRatio - - H[0, 0] = H[1, 1] = H[2, 2] = E * ( - 1.0 / (1 + nu) + nu / ((1.0 + nu) * (1 - 2 * nu)) - ) # H1 = 2*mu+lamb - H[0, 1] = H[0, 2] = H[1, 2] = E * (nu / ((1 + nu) * (1 - 2 * nu))) # H2 = lamb - H[3, 3] = H[4, 4] = H[5, 5] = 0.5 * E / (1 + nu) # H3 = mu - H[1, 0] = H[0, 1] - H[2, 0] = H[0, 2] - H[2, 1] = H[1, 2] # symétrie - - return H - - def HardeningFunction(self, p): - raise NameError( - "Hardening function not defined. Use the method SetHardeningFunction" + young_modulus = self.young_modulus + poisson_ratio = self.poisson_ratio + dtype = ( + float + if np.isscalar(young_modulus) and np.isscalar(poisson_ratio) + else object ) - - def HardeningFunctionDerivative(self, p): - raise NameError( - "Hardening function not defined. Use the method SetHardeningFunction" + elastic_matrix = np.zeros((6, 6), dtype=dtype) + lame_parameter = ( + young_modulus + * poisson_ratio + / ((1 + poisson_ratio) * (1 - 2 * poisson_ratio)) ) - - def SetHardeningFunction(self, FunctionType, **kargs): + elastic_matrix[:3, :3] = lame_parameter + elastic_matrix[0, 0] += 2 * self.shear_modulus + elastic_matrix[1, 1] += 2 * self.shear_modulus + elastic_matrix[2, 2] += 2 * self.shear_modulus + elastic_matrix[3, 3] = self.shear_modulus + elastic_matrix[4, 4] = self.shear_modulus + elastic_matrix[5, 5] = self.shear_modulus + return elastic_matrix + + def set_hardening_function(self, function_type, **kwargs): + """Define the isotropic hardening function. + + ``function_type="power"`` defines ``R(p) = h * p**beta``. + ``function_type="user"`` accepts ``hardening_function`` and + ``hardening_function_derivative`` callables. """ - Define the hardening function of the ElastoPlasticity law. - FunctionType is the type of hardening function. - - For now, the only defined hardening function is a power law. - * F = H*p^{beta} were p is the cumuled plasticity - - Other type of hardening function may be added in future versions. - - Parameters - ---------- - FunctionType: str - Type of hardening function. - For now, the only possible value is 'power' for power law. - H (keyword argument): scalar - beta(keyword argument): scalar - name: str, optional - The name of the constitutive law + function_type = function_type.lower() + + if function_type == "power": + if "h" not in kwargs: + raise TypeError("Keyword argument 'h' is required.") + if "beta" not in kwargs: + raise TypeError("Keyword argument 'beta' is required.") + + hardening_modulus = kwargs["h"] + beta = kwargs["beta"] + if hardening_modulus < 0: + raise ValueError("The hardening modulus h must be non-negative.") + if beta <= 0: + raise ValueError("The hardening exponent beta must be positive.") + + def hardening_function(plasticity): + return hardening_modulus * np.asarray(plasticity) ** beta + + def hardening_function_derivative(plasticity): + with np.errstate(divide="ignore", invalid="ignore"): + return ( + beta * hardening_modulus * np.asarray(plasticity) ** (beta - 1) + ) - """ - if FunctionType.lower() == "power": - H = None - beta = None - for item in kargs: - if item.lower() == "h": - H = kargs[item] - if item.lower() == "beta": - beta = kargs[item] - - if H is None: - raise NameError("Keyword arguments 'H' missing") - if beta is None: - raise NameError("Keyword arguments 'beta' missing") - - def HardeningFunction(p): - return H * p**beta - - def HardeningFunctionDerivative(p): - return np.nan_to_num( - beta * H * p ** (beta - 1), posinf=1 - ) # replace inf value by 1 - - elif FunctionType.lower() == "user": - HardeningFunction = None - HardeningFunctionDerivative = None - for item in kargs: - if item.lower() == "hardeningfunction": - HardeningFunction = kargs[item] - if item.lower() == "hardeningfunctionderivative": - HardeningFunctionDerivative = kargs[item] - - if HardeningFunction is None: - raise NameError("Keyword arguments 'HardeningFunction' missing") - if HardeningFunctionDerivative is None: - raise NameError( - "Keyword arguments 'HardeningFunctionDerivative' missing" + elif function_type == "user": + hardening_function = kwargs.get("hardening_function") + hardening_function_derivative = kwargs.get("hardening_function_derivative") + if hardening_function is None: + raise TypeError("Keyword argument 'hardening_function' is required.") + if hardening_function_derivative is None: + raise TypeError( + "Keyword argument 'hardening_function_derivative' " "is required." ) + else: + raise ValueError("function_type must be either 'power' or 'user'.") - self.HardeningFunction = HardeningFunction - self.HardeningFunctionDerivative = HardeningFunctionDerivative - - def YieldFunction(self, Stress, p): - return Stress.vonMises() - self.__YieldStress - self.HardeningFunction(p) + self._hardening_function = hardening_function + self._hardening_function_derivative = hardening_function_derivative - def YieldFunctionDerivativeSigma(self, sigma): - """ - Derivative of the Yield Function with respect to the stress tensor defined in sigma - sigma should be a StressTensorList object - """ - return StressTensorList( - (3 / 2) * np.array(sigma.deviatoric()) / sigma.vonMises() - ).toStrain() + def hardening_function(self, plasticity): + """Evaluate the isotropic hardening stress.""" + if self._hardening_function is None: + raise RuntimeError( + "The hardening function has not been defined. Call " + "set_hardening_function first." + ) + return self._hardening_function(plasticity) + + def hardening_function_derivative(self, plasticity): + """Evaluate the derivative of the isotropic hardening stress.""" + if self._hardening_function_derivative is None: + raise RuntimeError( + "The hardening function has not been defined. Call " + "set_hardening_function first." + ) + return self._hardening_function_derivative(plasticity) - def GetPlasticity(self): - return self.__currentP + def yield_function(self, stress, plasticity): + """Evaluate the von Mises yield function.""" + return ( + stress.von_mises() - self.yield_stress - self.hardening_function(plasticity) + ) - def get_strain(self, **kargs): - return self.__currentPlasticStrainTensor + def yield_function_derivative(self, stress): + """Differentiate the yield function with respect to stress.""" + equivalent_stress = np.asarray(stress.von_mises()) + if np.any(equivalent_stress == 0): + raise ZeroDivisionError( + "The yield direction is undefined at zero deviatoric stress." + ) + return StressTensorList( + 1.5 * np.asarray(stress.deviatoric()) / equivalent_stress + ).to_strain() + + @staticmethod + def _as_point_array(values): + array = np.asarray(values, dtype=float) + if array.ndim == 1: + array = array.reshape(6, 1) + if array.ndim != 2 or array.shape[0] != 6: + raise ValueError("Expected a Voigt array with shape (6, n_points).") + return array + + @staticmethod + def _deviatoric(stress): + deviator = stress.copy() + mean_stress = np.mean(stress[:3], axis=0) + deviator[:3] -= mean_stress + return deviator + + @staticmethod + def _flow_direction(stress): + stress_tensor = StressTensorList(stress) + equivalent_stress = np.asarray(stress_tensor.von_mises()) + direction = 1.5 * np.asarray(stress_tensor.deviatoric()) / equivalent_stress + direction[3:] *= 2.0 + return direction + + def _plastic_increment(self, equivalent_trial_stress, plasticity_old): + """Solve all active radial-return consistency equations together.""" + equivalent_trial_stress = np.asarray(equivalent_trial_stress, dtype=float) + plasticity_old = np.asarray(plasticity_old, dtype=float) + initial_residual = ( + equivalent_trial_stress + - self.yield_stress + - self.hardening_function(plasticity_old) + ) + lower_bound = np.zeros_like(initial_residual) + upper_bound = initial_residual / (3.0 * self.shear_modulus) + + def residual(plastic_increment): + return ( + equivalent_trial_stress + - 3.0 * self.shear_modulus * plastic_increment + - self.yield_stress + - self.hardening_function(plasticity_old + plastic_increment) + ) - def get_pk2(self): - return self.__currentSigma + for _ in range(self.max_return_mapping_iterations): + points_to_expand = residual(upper_bound) > 0 + if not np.any(points_to_expand): + break + upper_bound[points_to_expand] *= 2.0 + else: + raise RuntimeError("Unable to bracket every plastic multiplier.") - def get_cauchy(self, **kargs): # same as GetPKII - # alias of GetPKII mainly use for small strain displacement problems - return self.__currentSigma + plastic_increment = 0.5 * upper_bound + converged = np.zeros_like(plastic_increment, dtype=bool) - def get_stress(self, **kargs): # same as GetPKII - # alias of GetPKII mainly use for small strain displacement problems - return self.__currentSigma + for _ in range(self.max_return_mapping_iterations): + value = residual(plastic_increment) + converged |= np.abs(value) <= self.return_mapping_tolerance + if np.all(converged): + return plastic_increment - def get_disp_grad(self): - return self.__currentGradDisp + active = ~converged + positive_residual = active & (value > 0) + negative_residual = active & ~positive_residual + lower_bound[positive_residual] = plastic_increment[positive_residual] + upper_bound[negative_residual] = plastic_increment[negative_residual] - def get_tangent_matrix(self): - return self.__TangeantModuli + hardening_slope = self.hardening_function_derivative( + plasticity_old + plastic_increment + ) + denominator = 3.0 * self.shear_modulus + hardening_slope + candidate = plastic_increment + value / denominator + invalid_candidate = ( + ~np.isfinite(candidate) + | (denominator <= 0) + | (candidate <= lower_bound) + | (candidate >= upper_bound) + ) + candidate[invalid_candidate] = 0.5 * ( + lower_bound[invalid_candidate] + upper_bound[invalid_candidate] + ) + plastic_increment[active] = candidate[active] - # def GetStressOperator(self, localFrame=None): - # H = self.GetH() + unconverged_residual = np.max(np.abs(residual(plastic_increment)[~converged])) + raise RuntimeError( + "The radial-return algorithm did not converge after " + f"{self.max_return_mapping_iterations} iterations; maximum " + f"residual is {unconverged_residual:.6g}." + ) - # eps, eps_vir = GetStrainOperator(self.__currentGradDisp) - # sigma = [sum([0 if eps[j] is 0 else eps[j]*H[i][j] for j in range(6)]) for i in range(6)] + def _integrate(self, total_strain, plasticity_old, plastic_strain_old): + """Integrate all material points from their committed state.""" + strain = self._as_point_array(total_strain) + plastic_strain_old = self._as_point_array(plastic_strain_old) + plasticity_old = np.asarray(plasticity_old, dtype=float).reshape(-1) + n_points = strain.shape[1] + if plastic_strain_old.shape[1] != n_points or len(plasticity_old) != n_points: + raise ValueError("Inconsistent number of constitutive-law points.") + + elastic_matrix = np.asarray(self.get_elastic_matrix(), dtype=float) + trial_stress = elastic_matrix @ (strain - plastic_strain_old) + equivalent_trial_stress = np.asarray(StressTensorList(trial_stress).von_mises()) + trial_yield = ( + equivalent_trial_stress + - self.yield_stress + - self.hardening_function(plasticity_old) + ) + plastic_points = trial_yield > self.return_mapping_tolerance - # return sigma # list de 6 objets de type DiffOp + stress = trial_stress.copy() + plastic_strain = plastic_strain_old.copy() + plasticity = plasticity_old.copy() + tangent = np.repeat(elastic_matrix[:, :, None], n_points, axis=2) - def NewTimeIncrement(self): - # Set Irreversible Plasticity - if self.__P is not None: - self.__P = self.__currentP.copy() - self.__PlasticStrainTensor = self.__currentPlasticStrainTensor.copy() - self.__TangeantModuli = self.GetHelas() + if np.any(plastic_points): + plastic_increment = self._plastic_increment( + equivalent_trial_stress[plastic_points], + plasticity_old[plastic_points], + ) + plasticity[plastic_points] += plastic_increment + + trial_deviator = self._deviatoric(trial_stress[:, plastic_points]) + radial_factor = ( + 1.0 + - 3.0 + * self.shear_modulus + * plastic_increment + / equivalent_trial_stress[plastic_points] + ) + stress[:, plastic_points] = ( + trial_stress[:, plastic_points] + - trial_deviator + + radial_factor[None, :] * trial_deviator + ) - def to_start(self): - if self.__P is None: - self.__currentP = None - self.__currentPlasticStrainTensor = None - else: - self.__currentP = self.__P.copy() - self.__currentPlasticStrainTensor = self.__PlasticStrainTensor.copy() - self.__TangeantModuli = self.GetHelas() + direction = self._flow_direction(stress[:, plastic_points]) + plastic_strain[:, plastic_points] += direction * plastic_increment[None, :] - def reset(self): - """ - reset the constitutive law (time history) - """ - self.__P = None # irrevesrible plasticity - self.__currentP = None # current iteration plasticity (reversible) - self.__PlasticStrainTensor = None - self.__currentPlasticStrainTensor = None - self.__currentSigma = 0 # lissStressTensor object describing the last computed stress (GetStress method) - self.__TangeantModuli = self.GetHelas() - - def initialize(self, assembly, pb, t0=0.0, nlgeom=False): - if self._dimension is None: - self._dimension = assembly.space.get_dimension() - self.NewTimeIncrement() - self.nlgeom = nlgeom - - def update(self, assembly, pb, time): - displacement = pb.get_dof_solution() - - if np.isscalar(displacement) and displacement == 0: - self.__currentGradDisp = 0 - self.__currentSigma = 0 - else: - self.__currentGradDisp = assembly.get_grad_disp(displacement, "GaussPoint") - GradValues = self.__currentGradDisp - if self.nlgeom == False: - Strain = [GradValues[i][i] for i in range(3)] - Strain += [ - GradValues[0][1] + GradValues[1][0], - GradValues[0][2] + GradValues[2][0], - GradValues[1][2] + GradValues[2][1], - ] - else: - Strain = [ - GradValues[i][i] - + 0.5 * sum([GradValues[k][i] ** 2 for k in range(3)]) - for i in range(3) - ] - Strain += [ - GradValues[0][1] - + GradValues[1][0] - + sum([GradValues[k][0] * GradValues[k][1] for k in range(3)]) - ] - Strain += [ - GradValues[0][2] - + GradValues[2][0] - + sum([GradValues[k][0] * GradValues[k][2] for k in range(3)]) - ] - Strain += [ - GradValues[1][2] - + GradValues[2][1] - + sum([GradValues[k][1] * GradValues[k][2] for k in range(3)]) - ] - - TotalStrain = StrainTensorList(Strain) - self.ComputeStress( - TotalStrain, time - ) # compute the total stress in self.__currentSigma - - # print(self.__currentP) - - dphi_dp = self.HardeningFunctionDerivative(self.__currentP) - dphi_dsigma = self.YieldFunctionDerivativeSigma(self.__currentSigma) - Lambda = dphi_dsigma # for isotropic hardening only - test = self.YieldFunction(self.__currentSigma, self.__P) > self.__tol - - Helas = self.GetHelas() - ##### Compute new tangeant moduli - B = sum( - [ - sum([dphi_dsigma[j] * Helas[i][j] for j in range(6)]) * Lambda[i] - for i in range(6) - ] + hardening_slope = self.hardening_function_derivative( + plasticity[plastic_points] ) - Ap = B - dphi_dp - CL = [ - sum([Lambda[j] * Helas[i][j] for j in range(6)]) for i in range(6) - ] # [C:Lambda] - Peps = [ - sum([dphi_dsigma[i] * Helas[i][j] for i in range(6)]) / Ap - for j in range(6) - ] # Peps - # TangeantModuli = [[Helas[i][j] - CL[i]*Peps[j] for j in range(6)] for i in range(6)] - self.__TangeantModuli = [ - [Helas[i][j] - (CL[i] * Peps[j] * test) for j in range(6)] - for i in range(6) - ] - ##### end Compute new tangeant moduli - - def ComputeStress(self, StrainTensor, time=None): - # time not used here because this law require no time effect - # initilialize values plasticity variables if required - if self.__P is None: - self.__P = np.zeros(len(StrainTensor[0])) - self.__currentP = np.zeros(len(StrainTensor[0])) - if self.__PlasticStrainTensor is None: - self.__PlasticStrainTensor = StrainTensorList( - np.zeros((6, len(StrainTensor[0]))) + elastic_flow = elastic_matrix @ direction + denominator = ( + np.einsum("ip,ip->p", direction, elastic_flow) + hardening_slope ) - self.__currentPlasticStrainTensor = StrainTensorList( - np.zeros((6, len(StrainTensor[0]))) + tangent[:, :, plastic_points] -= ( + np.einsum("ip,jp->ijp", elastic_flow, elastic_flow) + / denominator[None, None, :] ) - H = ( - self.GetHelas() - ) # no change of basis because only isotropic behavior are considered - sigma = StressTensorList( - [ - sum( - [ - (StrainTensor[j] - self.__PlasticStrainTensor[j]) * H[i][j] - for j in range(6) - ] - ) - for i in range(6) - ] + return ( + StressTensorList(stress), + StrainTensorList(plastic_strain), + plasticity, + tangent, ) - test = self.YieldFunction(sigma, self.__P) > self.__tol - # print(sum(test)/len(test)*100) - - sigmaFull = np.array(sigma).T - Ep = np.array(self.__PlasticStrainTensor).T - # Ep = np.array(self.__currentPlasticStrainTensor).T - - for pg in range(len(sigmaFull)): - if test[pg] > 0: - sigma = StressTensorList(sigmaFull[pg]) - p = self.__P[pg] - iter = 0 - while abs(self.YieldFunction(sigma, p)) > self.__tol: - dphi_dp = self.HardeningFunctionDerivative(p) - dphi_dsigma = np.array(self.YieldFunctionDerivativeSigma(sigma)) - - Lambda = dphi_dsigma # for associated plasticity - B = sum( - [ - sum([dphi_dsigma[j] * H[i][j] for j in range(6)]) - * Lambda[i] - for i in range(6) - ] - ) - dp = self.YieldFunction(sigma, p) / (B - dphi_dp) - p += dp - Ep[pg] += Lambda * dp - sigma = StressTensorList( - [ - sum( - [ - (StrainTensor[j][pg] - Ep[pg][j]) * H[i][j] - for j in range(6) - ] - ) - for i in range(6) - ] - ) - self.__currentP[pg] = p - sigmaFull[pg] = sigma + def compute_stress( + self, + total_strain, + plasticity_old=None, + plastic_strain_old=None, + ): + """Integrate a material-point state without an assembly. + + This helper is stateless with respect to history: callers advancing + several increments must pass the previously returned plasticity and + plastic strain explicitly. + """ + strain = self._as_point_array(total_strain) + n_points = strain.shape[1] + if plasticity_old is None: + plasticity_old = np.zeros(n_points) + if plastic_strain_old is None: + plastic_strain_old = np.zeros((6, n_points)) + + ( + self._current_stress, + self._current_plastic_strain, + self._current_plasticity, + self._current_tangent, + ) = self._integrate(strain, plasticity_old, plastic_strain_old) + return self._current_stress + + def get_plasticity(self): + """Return plasticity from the most recent integration.""" + return self._current_plasticity + + def get_plastic_strain(self): + """Return plastic strain from the most recent integration.""" + return self._current_plastic_strain + + def get_stress(self): + """Return stress from the most recent integration.""" + return self._current_stress + + def get_tangent_matrix(self, assembly=None, dimension=None): + """Return the current tangent, or the elastic tangent initially.""" + if assembly is not None: + if dimension is None: + dimension = assembly.space.get_dimension() + if "TangentMatrix" in assembly.sv: + return assembly.sv["TangentMatrix"] + if self._current_tangent is not None: + return self._current_tangent + return self.get_elastic_matrix(dimension or "3D") + + def initialize(self, assembly, pb): + """Initialize finite-element state variables.""" + self._dimension = assembly.space.get_dimension() + elastic_matrix = np.asarray( + self.get_elastic_matrix(self._dimension), dtype=float + ) + n_points = assembly.n_gauss_points + assembly.sv["P"] = np.zeros(n_points) + assembly.sv["EP"] = StrainTensorList(np.zeros((6, n_points), order="F")) + assembly.sv["TangentMatrix"] = np.repeat( + elastic_matrix[:, :, None], n_points, axis=2 + ) + self.is_initialized = True - self.__currentPlasticStrainTensor = StrainTensorList(Ep.T) - self.__currentSigma = StressTensorList(sigmaFull.T) # list of 6 objets + def update(self, assembly, pb): + """Update stress, plastic state, and tangent for the current iterate.""" + if "DStrain" in assembly.sv: + total_strain = assembly.sv["Strain"] + assembly.sv["DStrain"] + else: + total_strain = assembly.sv["Strain"] + + start_state = getattr(assembly, "sv_start", assembly.sv) + ( + self._current_stress, + self._current_plastic_strain, + self._current_plasticity, + self._current_tangent, + ) = self._integrate( + total_strain, + start_state["P"], + start_state["EP"], + ) + assembly.sv["Stress"] = self._current_stress + assembly.sv["EP"] = self._current_plastic_strain + assembly.sv["P"] = self._current_plasticity + assembly.sv["TangentMatrix"] = self._current_tangent + + def set_start(self, assembly, pb): + """Commit the converged state and prepare the elastic predictor.""" + if assembly._nlgeom and "DR" in assembly.sv: + rotation = SimRotation.from_matrix(assembly.sv["DR"].transpose(2, 0, 1)) + assembly.sv["EP"] = StrainTensorList( + rotation.apply_strain(assembly.sv["EP"].asarray()) + ) - return self.__currentSigma + elastic_matrix = np.asarray( + self.get_elastic_matrix(self._dimension), dtype=float + ) + assembly.sv["TangentMatrix"] = np.repeat( + elastic_matrix[:, :, None], + assembly.n_gauss_points, + axis=2, + ) + + def reset(self): + """Reset cached results; assembly history is managed by Fedoo.""" + super().reset() + self._current_stress = None + self._current_plastic_strain = None + self._current_plasticity = None + self._current_tangent = None diff --git a/fedoo/constitutivelaw/fe2.py b/fedoo/constitutivelaw/fe2.py index 075e2e15..8864bda2 100644 --- a/fedoo/constitutivelaw/fe2.py +++ b/fedoo/constitutivelaw/fe2.py @@ -1,272 +1,129 @@ -# derive de ConstitutiveLaw -# This law should be used with an StressEquilibrium WeakForm +"""FE² constitutive law.""" -from fedoo.core.mechanical3d import Mechanical3D -from fedoo.weakform.stress_equilibrium import StressEquilibrium +import numpy as np + +from fedoo.constraint.periodic_bc import PeriodicBC from fedoo.core.assembly import Assembly -from fedoo.problem.non_linear import NonLinear -from fedoo.util.voigt_tensors import StressTensorList, StrainTensorList -from fedoo.constraint.periodic_bc import ( - PeriodicBC, -) # , DefinePeriodicBoundaryConditionNonPerioMesh +from fedoo.core.mechanical3d import Mechanical3D from fedoo.homogen.tangent_stiffness import ( - get_tangent_stiffness, get_homogenized_stiffness, + get_tangent_stiffness, ) -import numpy as np -import multiprocessing +from fedoo.problem.non_linear import NonLinear +from fedoo.util.voigt_tensors import StressTensorList, StrainTensorList class FE2(Mechanical3D): - """ - ConstitutiveLaw that solve a Finite Element Problem at each point of gauss - in the contexte of the so called "FE²" method. + """FE² constitutive law based on periodic microscopic problems. + + A microscopic finite-element problem is solved at each macroscopic + integration point. The macroscopic strain is prescribed through the + problem-level ``MeanStrain`` degrees of freedom created by + :class:`PeriodicBC`. Parameters ---------- - assemb: Assembly or Assembly name (str), or list of Assembly - (with len(list) = number of integration points). - Assembly that correspond to the microscopic problem - name: str, optional - The name of the constitutive law + assemb : Assembly, str, or list of Assembly + Microscopic assembly, its registered name, or one microscopic + assembly per macroscopic integration point. + name : str, optional + Name of the constitutive law. """ def __init__(self, assemb, name=""): - # props is a nparray containing all the material variables - # nstatev is a nparray containing all the material variables if isinstance(assemb, str): assemb = Assembly.get_all()[assemb] - Mechanical3D.__init__(self, name) # heritage + super().__init__(name) if isinstance(assemb, list): self.__assembly = [ - Assembly.get_all()[a] if isinstance(a, str) else a for a in assemb + Assembly.get_all()[item] if isinstance(item, str) else item + for item in assemb ] - self.__mesh = [a.mesh for a in self.__assembly] + self.__mesh = [item.mesh for item in self.__assembly] else: - self.__mesh = assemb.mesh self.__assembly = assemb + self.__mesh = assemb.mesh self.list_problem = None - - self.use_elastic_lt = True # option to use the elastic tangeant matrix (in principle = initial tangent matrix) at the begining of each time step - - # self.__currentGradDisp = self.__initialGradDisp = 0 - - # def get_pk2(self): - # return StressTensorList(self.__stress) - - # # def get_kirchhoff(self): - # # return StressTensorList(self.Kirchhoff.T) - - # # def get_cauchy(self): - # # return StressTensorList(self.Cauchy.T) - - # def get_strain(self, **kargs): - # return StrainTensorList(self.__strain) - - # # def get_statev(self): - # # return self.statev.T - - # def get_stress(self, **kargs): #same as GetPKII (used for small def) - # return StressTensorList(self.__stress) - - # # def GetHelas (self): - # # # if self.__L is None: - # # # self.RunUmat(np.eye(3).T.reshape(1,3,3), np.eye(3).T.reshape(1,3,3), time=0., dtime=1.) - - # # return np.squeeze(self.L.transpose(1,2,0)) - - # def get_wm(self): - # return self.__Wm - - # def get_disp_grad(self): - # if np.isscalar(self.__currentGradDisp) and self.__currentGradDisp == 0: return 0 - # else: return self.__currentGradDisp - - # def get_tangent_matrix(self): - - # H = np.squeeze(self.Lt.transpose(1,2,0)) - # return H - - # def NewTimeIncrement(self): - # # self.set_start() #in set_start -> set tangeant matrix to elastic - - # #save variable at the begining of the Time increment - # self.__initialGradDisp = self.__currentGradDisp - # self.Lt = self.L.copy() - - # def to_start(self): - # # self.to_start() - # self.__currentGradDisp = self.__initialGradDisp - # self.Lt = self.L.copy() - - # def reset(self): - # """ - # reset the constitutive law (time history) - # """ - # #a modifier - # self.__currentGradDisp = self.__initialGradDisp = 0 - # # self.__Statev = None - # self.__currentStress = None #lissStressTensor object describing the last computed stress (GetStress method) - # # self.__currentGradDisp = 0 - # # self.__F0 = None + self.use_elastic_lt = True def initialize(self, assembly, pb): - if self.list_problem is None: # only initialize once - nb_points = assembly.n_gauss_points - - # Definition of the set of nodes for boundary conditions - if not (isinstance(self.__mesh, list)): - self.list_mesh = [self.__mesh for i in range(nb_points)] - self.list_assembly = [self.__assembly.copy() for i in range(nb_points)] - else: - self.list_mesh = self.__mesh - self.list_assembly = self.__assembly - - self.list_problem = [] - self._list_volume = np.empty(nb_points) - self._list_center = np.empty(nb_points, dtype=int) - # self.L = np.empty((nb_points,6,6)) - assembly.sv["TangentMatrix"] = np.empty((6, 6, nb_points)) - - print("-- Initialize micro problems --") - for i in range(nb_points): - print("\r", str(i + 1), "/", str(nb_points), end="") - crd = self.list_mesh[i].nodes - type_el = self.list_mesh[i].elm_type - xmax = np.max(crd[:, 0]) - xmin = np.min(crd[:, 0]) - ymax = np.max(crd[:, 1]) - ymin = np.min(crd[:, 1]) - zmax = np.max(crd[:, 2]) - zmin = np.min(crd[:, 2]) - - crd_center = ( - np.array([xmin, ymin, zmin]) + np.array([xmax, ymax, zmax]) - ) / 2 - self._list_volume[i] = ( - (xmax - xmin) * (ymax - ymin) * (zmax - zmin) - ) # total volume of the domain - - if "_StrainNodes" in self.list_mesh[i].node_sets: - strain_nodes = self.list_mesh[i].node_sets["_StrainNodes"] - else: - strain_nodes = self.list_mesh[i].add_nodes( - crd_center, 2 - ) # add virtual nodes for macro strain - self.list_mesh[i].add_node_set(strain_nodes, "_StrainNodes") - - self._list_center[i] = np.linalg.norm( - crd[:-2] - crd_center, axis=1 - ).argmin() - # list_material.append(self.__constitutivelaw.copy()) - - # Type of problem - self.list_problem.append( - NonLinear(self.list_assembly[i], name="_fe2_cell_" + str(i)) - ) - pb_micro = self.list_problem[-1] - meshperio = True - - # Shall add other conditions later on - pb_micro.bc.add( - PeriodicBC( - [ - strain_nodes[0], - strain_nodes[0], - strain_nodes[0], - strain_nodes[1], - strain_nodes[1], - strain_nodes[1], - ], - ["DispX", "DispY", "DispZ", "DispX", "DispY", "DispZ"], - dim=3, - meshperio=meshperio, - name="_fe2_cell_" + str(i), - ) + if self.list_problem is not None: + return + + nb_points = assembly.n_gauss_points + if isinstance(self.__mesh, list): + if len(self.__assembly) != nb_points: + raise ValueError( + "A list of microscopic assemblies must contain one " + "assembly per macroscopic integration point." ) - - pb_micro.bc.add("Dirichlet", [self._list_center[i]], "Disp", 0) - # self.list_assembly[i].initialize() - assembly.sv["TangentMatrix"][:, :, i] = get_homogenized_stiffness( - self.list_assembly[i] + self.list_mesh = self.__mesh + self.list_assembly = self.__assembly + else: + self.list_mesh = [self.__mesh for _ in range(nb_points)] + self.list_assembly = [self.__assembly.copy() for _ in range(nb_points)] + + self.list_problem = [] + self._list_volume = np.empty(nb_points) + assembly.sv["TangentMatrix"] = np.empty((6, 6, nb_points)) + + print("-- Initialize micro problems --") + for index in range(nb_points): + print("\r", index + 1, "/", nb_points, end="") + coordinates = self.list_mesh[index].nodes + lower = np.min(coordinates, axis=0) + upper = np.max(coordinates, axis=0) + box_center = (lower + upper) / 2 + self._list_volume[index] = np.prod(upper - lower) + center_node = np.linalg.norm(coordinates - box_center, axis=1).argmin() + + micro_problem = NonLinear( + self.list_assembly[index], + name=f"_fe2_cell_{index}", + ) + self.list_problem.append(micro_problem) + micro_problem.bc.add( + PeriodicBC( + "small_strain", + dim=3, + name=f"_fe2_cell_{index}", ) + ) + micro_problem.bc.add("Dirichlet", [center_node], "Disp", 0) + assembly.sv["TangentMatrix"][:, :, index] = get_homogenized_stiffness( + self.list_assembly[index] + ) - pb.make_active() - if self.use_elastic_lt: - assembly.sv["ElasticMatrix"] = assembly.sv["TangentMatrix"].copy() - - assembly.sv["Strain"] = StrainTensorList(np.zeros((6, nb_points))) - assembly.sv["Stress"] = StressTensorList(np.zeros((6, nb_points))) - assembly.sv["Wm"] = np.zeros((4, nb_points)) + pb.make_active() + if self.use_elastic_lt: + assembly.sv["ElasticMatrix"] = assembly.sv["TangentMatrix"].copy() - print("") + assembly.sv["Strain"] = StrainTensorList(np.zeros((6, nb_points))) + assembly.sv["Stress"] = StressTensorList(np.zeros((6, nb_points))) + assembly.sv["Wm"] = np.zeros((4, nb_points)) + print("") def set_start(self, assembly, pb): if self.use_elastic_lt: - assembly.sv["TangentMatrix"] = assembly.sv["ElasticMatrix"] + assembly.sv["TangentMatrix"] = assembly.sv["ElasticMatrix"].copy() - def _update_pb(self, id_pb, assembly_macro, pb_macro): + def _update_pb(self, index, assembly_macro, pb_macro): strain = assembly_macro.sv["Strain"] strain_start = assembly_macro.sv_start["Strain"] - nb_points = len(self.list_problem) - pb = self.list_problem[id_pb] - - print("\r", str(id_pb + 1), "/", str(nb_points), end="") - strain_nodes = self.list_mesh[id_pb].node_sets["_StrainNodes"] + micro_problem = self.list_problem[index] - pb.bc.remove("Strain") - pb.bc.add( - "Dirichlet", - [strain_nodes[0]], - "DispX", - strain[0][id_pb], - start_value=strain_start[0][id_pb], - name="Strain", - ) # EpsXX - pb.bc.add( + print("\r", index + 1, "/", len(self.list_problem), end="") + micro_problem.bc.remove("Strain") + micro_problem.bc.add( "Dirichlet", - [strain_nodes[0]], - "DispY", - strain[1][id_pb], - start_value=strain_start[1][id_pb], + "MeanStrain", + strain.asarray()[:, index], + start_value=strain_start.asarray()[:, index], name="Strain", - ) # EpsYY - pb.bc.add( - "Dirichlet", - [strain_nodes[0]], - "DispZ", - strain[2][id_pb], - start_value=strain_start[2][id_pb], - name="Strain", - ) # EpsZZ - pb.bc.add( - "Dirichlet", - [strain_nodes[1]], - "DispX", - strain[3][id_pb], - start_value=strain_start[3][id_pb], - name="Strain", - ) # EpsXY - pb.bc.add( - "Dirichlet", - [strain_nodes[1]], - "DispY", - strain[4][id_pb], - start_value=strain_start[4][id_pb], - name="Strain", - ) # EpsXZ - pb.bc.add( - "Dirichlet", - [strain_nodes[1]], - "DispZ", - strain[5][id_pb], - start_value=strain_start[5][id_pb], - name="Strain", - ) # EpsYZ - - pb.nlsolve( + ) + micro_problem.nlsolve( dt=pb_macro.dtime, tmax=pb_macro.dtime, update_dt=True, @@ -274,36 +131,27 @@ def _update_pb(self, id_pb, assembly_macro, pb_macro): print_info=0, ) - assembly_macro.sv["TangentMatrix"][:, :, id_pb] = get_tangent_stiffness(pb.name) + assembly_macro.sv["TangentMatrix"][:, :, index] = get_tangent_stiffness( + micro_problem.name + ) - stress_field = self.list_assembly[id_pb].sv["Stress"] # computed micro stress - # integrate micro stress to get the macro one - assembly_macro.sv["Stress"].asarray()[:, id_pb] = np.array( + micro_assembly = self.list_assembly[index] + stress_field = micro_assembly.sv["Stress"] + assembly_macro.sv["Stress"].asarray()[:, index] = np.array( [ - 1 - / self._list_volume[id_pb] - * self.list_assembly[id_pb].integrate_field(stress_field[i]) - for i in range(6) + micro_assembly.integrate_field(component) / self._list_volume[index] + for component in stress_field ] ) - Wm_field = self.list_assembly[id_pb].sv["Wm"] - assembly_macro.sv["Wm"][:, id_pb] = ( - 1 / self._list_volume[id_pb] - ) * self.list_assembly[id_pb].integrate_field(Wm_field) + energy_field = micro_assembly.sv.get("Wm") + if energy_field is not None: + assembly_macro.sv["Wm"][:, index] = ( + micro_assembly.integrate_field(energy_field) / self._list_volume[index] + ) def update(self, assembly, pb): - displacement = pb.get_dof_solution() - - # resolution of the micro problem at each gauss points - nb_points = len(self.list_problem) - print("-- Update micro cells --") - - # with multiprocessing.Pool(4) as pool: - # pool.map(self._update_pb, range(nb_points)) - - for id_pb in range(nb_points): - self._update_pb(id_pb, assembly, pb) - + for index in range(len(self.list_problem)): + self._update_pb(index, assembly, pb) print("") diff --git a/fedoo/constitutivelaw/shell.py b/fedoo/constitutivelaw/shell.py index 04c35b34..3d59f13f 100644 --- a/fedoo/constitutivelaw/shell.py +++ b/fedoo/constitutivelaw/shell.py @@ -4,6 +4,7 @@ from fedoo.core.base import ConstitutiveLaw from fedoo.util.voigt_tensors import StressTensorList, StrainTensorList +import copy import numpy as np @@ -110,9 +111,9 @@ def get_strain(self, assembly, **kargs): if np.isscalar(ShellStrain) and ShellStrain == 0: zeros = np.zeros(assembly.n_gauss_points) return StrainTensorList([zeros.copy() for _ in range(6)]) - Strain[0] = ShellStrain[0] + z * ShellStrain[4] # epsXX -> membrane and bending - Strain[1] = ShellStrain[1] - z * ShellStrain[3] # epsYY -> membrane and bending - Strain[3] = ShellStrain[2] # 2epsXY + Strain[0] = ShellStrain[0] + z * ShellStrain[3] # epsXX -> membrane and bending + Strain[1] = ShellStrain[1] + z * ShellStrain[4] # epsYY -> membrane and bending + Strain[3] = ShellStrain[2] + z * ShellStrain[5] # 2epsXY -> membrane and twist Strain[4:6] = ShellStrain[6:8] # 2epsXZ and 2epsYZ -> shear return Strain @@ -194,21 +195,23 @@ def get_stress(self, assembly, **kargs): return StressTensorList(Stress) - def GetStressDistribution(self, assembly, pg, resolution=100): + def get_stress_distribution(self, assembly, pg, resolution=100): h = self.thickness z = np.arange(-h / 2, h / 2, h / resolution) Strain = StrainTensorList([0 for i in range(6)]) ShellStrain = assembly.sv["ShellStrain"] Strain[0] = ( - ShellStrain[0][pg] + z * ShellStrain[4][pg] + ShellStrain[0][pg] + z * ShellStrain[3][pg] ) # epsXX -> membrane and bending Strain[1] = ( - ShellStrain[1][pg] - z * ShellStrain[3][pg] + ShellStrain[1][pg] + z * ShellStrain[4][pg] ) # epsYY -> membrane and bending - Strain[3] = ShellStrain[2][pg] * np.ones_like(z) # 2epsXY + Strain[3] = ( + ShellStrain[2][pg] + z * ShellStrain[5][pg] + ) # 2epsXY -> membrane and twist Strain[4] = ShellStrain[6][pg] * np.ones_like(z) # 2epsXZ -> shear - Strain[5] = ShellStrain[6][pg] * np.ones_like(z) # 2epsYZ -> shear + Strain[5] = ShellStrain[7][pg] * np.ones_like(z) # 2epsYZ -> shear Hplane = self.material.get_elastic_matrix( "2Dstress" @@ -244,6 +247,628 @@ def GetStressDistribution(self, assembly, pg, resolution=100): return z, Stress +class _ShellMaterialPointAssembly: + """Minimal assembly interface used by through-thickness material points.""" + + class _MaterialSpace: + def __init__(self, dimension): + self._dimension = dimension + + def get_dimension(self): + return self._dimension + + @property + def ndim(self): + return 2 if self._dimension == "2Dstress" else 3 + + @staticmethod + def list_variables(): + return () + + def __init__(self, n_gauss_points, dimension="2Dstress"): + self.n_gauss_points = n_gauss_points + self.space = self._MaterialSpace(dimension) + self.sv = { + "Strain": StrainTensorList(np.zeros((6, n_gauss_points))), + "Stress": StressTensorList(np.zeros((6, n_gauss_points))), + } + self.sv_start = dict(self.sv) + self.sv_component = {} + self._nlgeom = False + + def convert_data(self, data, *args, **kwargs): + """Return pointwise material data without mesh interpolation.""" + return data + + +def _copy_state(state): + """Deep-ish copy of a state-variable dict (shared by the nonlinear shells).""" + copied = {} + for key, value in state.items(): + if hasattr(value, "copy"): + copied[key] = value.copy() + else: + copied[key] = copy.deepcopy(value) + return copied + + +def _tangent_array(tangent, n_points): + """Normalize scalar/pointwise tangent entries to a dense (6, 6, n) array.""" + dense = np.empty((6, 6, n_points), dtype=float) + for row in range(6): + for column in range(6): + dense[row, column] = tangent[row][column] + return dense + + +def _strain_matrix(z): + """Membrane+bending strain-interpolation matrix at through-thickness ``z``.""" + matrix = np.zeros((6, 8)) + matrix[0, 0] = 1 + matrix[0, 3] = z + matrix[1, 1] = 1 + matrix[1, 4] = z + matrix[3, 2] = 1 + matrix[3, 5] = z + return matrix + + +class ShellHomogeneousNonLinear(ShellBase): + """Homogeneous shell integrated from nonlinear plane-stress material points. + + The membrane and bending response is obtained by calling a private copy + of ``material`` at Gauss--Legendre points through the thickness. The + material must support the standard Fedoo ``"2Dstress"`` constitutive-law + interface. Transverse shear retains the elastic Reissner--Mindlin + treatment used by :class:`ShellHomogeneous`. + + This implementation is restricted to small-strain material laws. For + plasticity (or any law needing a plane-stress return mapping) use a Simcoon + law, e.g. ``fedoo.constitutivelaw.Simcoon("EPICP", props)``; the pedagogical + :class:`ElastoPlasticity` is 3D-only and does not support ``"2Dstress"``. + + Parameters + ---------- + material : ConstitutiveLaw or str + Plane-stress-compatible material law or its registered name. + thickness : float + Shell thickness. + n_thickness_points : int, default=5 + Number of Gauss--Legendre points through the thickness. + k : float, default=1 + Transverse shear correction factor. + name : str, optional + Name of the shell constitutive law. + """ + + def __init__( + self, + material, + thickness, + n_thickness_points=5, + k=1, + name="", + ): + if isinstance(material, str): + material = ConstitutiveLaw.get_all()[material] + if n_thickness_points < 1: + raise ValueError("n_thickness_points must be at least one.") + + super().__init__(thickness, k, name) + self.material = copy.deepcopy(material) + self._shear_material = copy.deepcopy(material) + self.n_thickness_points = n_thickness_points + points, weights = np.polynomial.legendre.leggauss(n_thickness_points) + self._z = points * thickness / 2 + self._weights = weights * thickness / 2 + self._material_assembly = None + self._n_shell_points = None + self._shear_matrix = None + + def compute_area_density(self): + density = self._material_density(self.material, self.name) + return density * self.thickness + + def compute_rotary_density(self): + density = self._material_density(self.material, self.name) + return density * self.thickness**3 / 12.0 + + def _elastic_shear_matrix(self): + if self._shear_matrix is not None: + return self._shear_matrix + if not hasattr(self.material, "get_elastic_matrix"): + raise RuntimeError( + "The shell law must be initialized before its transverse " + "shear stiffness can be requested." + ) + elastic = np.asarray(self.material.get_elastic_matrix("3D"), dtype=float) + return elastic[np.ix_([4, 5], [4, 5])] + + def get_shell_stiffness_matrix(self): + if self._material_assembly is not None: + return self._shell_tangent + + plane = np.asarray(self.material.get_elastic_matrix("2Dstress"), dtype=float) + tangent = np.zeros((8, 8)) + for z, weight in zip(self._z, self._weights): + strain_matrix = _strain_matrix(z) + tangent += weight * strain_matrix.T @ plane @ strain_matrix + tangent[6:8, 6:8] = self.k * self.thickness * self._elastic_shear_matrix() + return tangent + + def initialize(self, assembly, pb): + self._n_shell_points = assembly.n_gauss_points + n_material_points = self.n_thickness_points * self._n_shell_points + self._material_assembly = _ShellMaterialPointAssembly(n_material_points) + shear_assembly = _ShellMaterialPointAssembly(1, dimension="3D") + self.material.reset() + self._shear_material.reset() + self.material.initialize(self._material_assembly, pb) + self._shear_material.initialize(shear_assembly, pb) + shear_tangent = _tangent_array(shear_assembly.sv["TangentMatrix"], 1) + shear_tangent = shear_tangent[:, :, 0] + self._shear_matrix = shear_tangent[np.ix_([4, 5], [4, 5])] + self._material_assembly.sv_start = _copy_state(self._material_assembly.sv) + + assembly.sv["_ShellStiffnessMatrix"] = self._integrate_tangent() + assembly.sv["ShellStress"] = 0 + + def _material_strain(self, shell_strain): + shell_array = np.asarray( + [ + ( + np.zeros(self._n_shell_points) + if np.isscalar(component) + else component + ) + for component in shell_strain + ] + ) + strains = [_strain_matrix(z) @ shell_array for z in self._z] + return StrainTensorList(np.concatenate(strains, axis=1)) + + def _integrate_tangent(self): + material_tangent = _tangent_array( + self._material_assembly.sv["TangentMatrix"], + self.n_thickness_points * self._n_shell_points, + ) + + shell_tangent = np.zeros((8, 8, self._n_shell_points)) + for index, (z, weight) in enumerate(zip(self._z, self._weights)): + point_slice = slice( + index * self._n_shell_points, + (index + 1) * self._n_shell_points, + ) + strain_matrix = _strain_matrix(z) + shell_tangent += weight * np.einsum( + "ia,ijp,jb->abp", + strain_matrix, + material_tangent[:, :, point_slice], + strain_matrix, + ) + + shell_tangent[6:8, 6:8] = ( + self.k * self.thickness * self._elastic_shear_matrix()[:, :, None] + ) + self._shell_tangent = shell_tangent + return shell_tangent + + def _integrate_stress(self, shell_strain): + material_stress = self._material_assembly.sv["Stress"].asarray() + resultants = np.zeros((8, self._n_shell_points)) + for index, (z, weight) in enumerate(zip(self._z, self._weights)): + point_slice = slice( + index * self._n_shell_points, + (index + 1) * self._n_shell_points, + ) + resultants += weight * ( + _strain_matrix(z).T @ material_stress[:, point_slice] + ) + + shear_strain = np.asarray( + [ + ( + np.zeros(self._n_shell_points) + if np.isscalar(shell_strain[index]) + else shell_strain[index] + ) + for index in [6, 7] + ] + ) + resultants[6:8] = ( + self.k * self.thickness * self._elastic_shear_matrix() @ shear_strain + ) + return type(shell_strain)(resultants) + + def update(self, assembly, pb): + shell_strain = assembly.sv["ShellStrain"] + if np.isscalar(shell_strain) and shell_strain == 0: + shell_strain = [np.zeros(self._n_shell_points) for _ in range(8)] + + self._material_assembly.sv["Strain"] = self._material_strain(shell_strain) + self.material.update(self._material_assembly, pb) + assembly.sv["_ShellStiffnessMatrix"] = self._integrate_tangent() + assembly.sv["ShellStress"] = self._integrate_stress(shell_strain) + + def set_start(self, assembly, pb): + self.material.set_start(self._material_assembly, pb) + self._material_assembly.sv_start = _copy_state(self._material_assembly.sv) + + def to_start(self, assembly, pb): + self.material.to_start(self._material_assembly, pb) + self._material_assembly.sv = _copy_state(self._material_assembly.sv_start) + + def get_stress_distribution(self, assembly, pg, resolution=None): + """Return stresses through the thickness at one shell Gauss point. + + By default, stresses are returned at the constitutive integration + points. If ``resolution`` is given, every stress component is + interpolated onto that many equally spaced positions for + visualization. + """ + if self._material_assembly is None: + raise RuntimeError("The shell law has not been initialized.") + if not 0 <= pg < self._n_shell_points: + raise IndexError("pg is outside the shell Gauss-point range.") + + stress = ( + self._material_assembly.sv["Stress"] + .asarray()[ + :, + pg :: self._n_shell_points, + ] + .copy() + ) + shell_strain = assembly.sv["ShellStrain"] + shear_strain = np.array( + [ + (0.0 if np.isscalar(shell_strain[index]) else shell_strain[index][pg]) + for index in [6, 7] + ] + ) + stress[4:6] = self._elastic_shear_matrix() @ shear_strain[:, None] + + z = self._z.copy() + if resolution is not None: + if resolution < 2: + raise ValueError("resolution must be at least two.") + output_z = np.linspace( + -self.thickness / 2, + self.thickness / 2, + resolution, + ) + stress = np.array( + [np.interp(output_z, z, component) for component in stress] + ) + z = output_z + return z, StressTensorList(stress) + + +class ShellLaminateNonLinear(ShellBase): + """Layered shell integrated from nonlinear plane-stress material points. + + Each layer owns an independent copy of its material law and its internal + variables. Membrane and bending stresses and tangents are integrated at + Gauss--Legendre points in every layer. Transverse shear is treated + elastically using the three-dimensional initial tangent of each material + and the shell correction factor ``k``. + + This implementation is restricted to small-strain material laws. For + plasticity (or any law needing a plane-stress return mapping) use a Simcoon + law, e.g. ``fedoo.constitutivelaw.Simcoon("EPICP", props)``; the pedagogical + :class:`ElastoPlasticity` is 3D-only and does not support ``"2Dstress"``. + + Parameters + ---------- + list_mat : sequence of ConstitutiveLaw or str + Material law, or registered material name, for every layer, ordered + from the bottom to the top surface. + list_thickness : sequence of float + Thickness of every layer. + n_thickness_points : int or sequence of int, default=3 + Number of Gauss--Legendre points in each layer. + k : float, default=1 + Transverse shear correction factor. + name : str, optional + Name of the shell constitutive law. + """ + + def __init__( + self, + list_mat, + list_thickness, + n_thickness_points=3, + k=1, + name="", + ): + if len(list_mat) != len(list_thickness): + raise ValueError("list_mat and list_thickness must have the same length.") + if len(list_mat) == 0: + raise ValueError("At least one laminate layer is required.") + if np.any(np.asarray(list_thickness) <= 0): + raise ValueError("Every layer thickness must be positive.") + + materials = [ + ( + ConstitutiveLaw.get_all()[material] + if isinstance(material, str) + else material + ) + for material in list_mat + ] + if np.isscalar(n_thickness_points): + point_counts = [int(n_thickness_points)] * len(materials) + else: + point_counts = [int(value) for value in n_thickness_points] + if len(point_counts) != len(materials): + raise ValueError( + "n_thickness_points must be an integer or contain one " + "value per layer." + ) + if any(value < 1 for value in point_counts): + raise ValueError("Every layer must have at least one thickness point.") + + thickness = float(np.sum(list_thickness)) + super().__init__(thickness, k, name) + self.materials = [copy.deepcopy(material) for material in materials] + self._shear_materials = [copy.deepcopy(material) for material in materials] + self.list_thickness = np.asarray(list_thickness, dtype=float) + self.n_thickness_points = point_counts + self._interfaces = ( + np.concatenate(([0.0], np.cumsum(self.list_thickness))) - thickness / 2 + ) + self._layer_z = [] + self._layer_weights = [] + for index, point_count in enumerate(point_counts): + points, weights = np.polynomial.legendre.leggauss(point_count) + lower = self._interfaces[index] + upper = self._interfaces[index + 1] + self._layer_z.append((lower + upper) / 2 + points * (upper - lower) / 2) + self._layer_weights.append(weights * (upper - lower) / 2) + + self._material_assemblies = None + self._shear_matrices = None + self._n_shell_points = None + + def compute_area_density(self): + return sum( + self._material_density(material, self.name) * thickness + for material, thickness in zip(self.materials, self.list_thickness) + ) + + def compute_rotary_density(self): + return sum( + self._material_density(material, self.name) + * (self._interfaces[index + 1] ** 3 - self._interfaces[index] ** 3) + / 3 + for index, material in enumerate(self.materials) + ) + + def _initial_shear_matrix(self, index): + if self._shear_matrices is not None: + return self._shear_matrices[index] + material = self.materials[index] + if not hasattr(material, "get_elastic_matrix"): + raise RuntimeError( + "The laminate must be initialized before its transverse " + "shear stiffness can be requested." + ) + elastic = np.asarray(material.get_elastic_matrix("3D"), dtype=float) + return elastic[np.ix_([4, 5], [4, 5])] + + def get_shell_stiffness_matrix(self): + if self._material_assemblies is not None: + return self._shell_tangent + + tangent = np.zeros((8, 8)) + for index, material in enumerate(self.materials): + plane = np.asarray(material.get_elastic_matrix("2Dstress"), dtype=float) + for z, weight in zip(self._layer_z[index], self._layer_weights[index]): + strain_matrix = _strain_matrix(z) + tangent += weight * strain_matrix.T @ plane @ strain_matrix + tangent[6:8, 6:8] += ( + self.k * self.list_thickness[index] * self._initial_shear_matrix(index) + ) + return tangent + + def initialize(self, assembly, pb): + self._n_shell_points = assembly.n_gauss_points + self._material_assemblies = [] + self._shear_matrices = [] + + for material, shear_material, point_count in zip( + self.materials, + self._shear_materials, + self.n_thickness_points, + ): + material_assembly = _ShellMaterialPointAssembly( + point_count * self._n_shell_points + ) + shear_assembly = _ShellMaterialPointAssembly(1, dimension="3D") + material.reset() + shear_material.reset() + material.initialize(material_assembly, pb) + shear_material.initialize(shear_assembly, pb) + + shear_tangent = _tangent_array(shear_assembly.sv["TangentMatrix"], 1) + shear_tangent = shear_tangent[:, :, 0] + self._shear_matrices.append(shear_tangent[np.ix_([4, 5], [4, 5])]) + material_assembly.sv_start = _copy_state(material_assembly.sv) + self._material_assemblies.append(material_assembly) + + assembly.sv["_ShellStiffnessMatrix"] = self._integrate_tangent() + assembly.sv["ShellStress"] = 0 + + def _material_strain(self, shell_strain, layer): + shell_array = np.asarray( + [ + ( + np.zeros(self._n_shell_points) + if np.isscalar(component) + else component + ) + for component in shell_strain + ] + ) + strains = [_strain_matrix(z) @ shell_array for z in self._layer_z[layer]] + return StrainTensorList(np.concatenate(strains, axis=1)) + + def _integrate_tangent(self): + shell_tangent = np.zeros((8, 8, self._n_shell_points)) + for layer, material_assembly in enumerate(self._material_assemblies): + material_tangent = _tangent_array( + material_assembly.sv["TangentMatrix"], + self.n_thickness_points[layer] * self._n_shell_points, + ) + + for point, (z, weight) in enumerate( + zip( + self._layer_z[layer], + self._layer_weights[layer], + ) + ): + point_slice = slice( + point * self._n_shell_points, + (point + 1) * self._n_shell_points, + ) + strain_matrix = _strain_matrix(z) + shell_tangent += weight * np.einsum( + "ia,ijp,jb->abp", + strain_matrix, + material_tangent[:, :, point_slice], + strain_matrix, + ) + + shell_tangent[6:8, 6:8] += ( + self.k + * self.list_thickness[layer] + * self._shear_matrices[layer][:, :, None] + ) + + self._shell_tangent = shell_tangent + return shell_tangent + + def _integrate_stress(self, shell_strain): + resultants = np.zeros((8, self._n_shell_points)) + for layer, material_assembly in enumerate(self._material_assemblies): + material_stress = material_assembly.sv["Stress"].asarray() + for point, (z, weight) in enumerate( + zip( + self._layer_z[layer], + self._layer_weights[layer], + ) + ): + point_slice = slice( + point * self._n_shell_points, + (point + 1) * self._n_shell_points, + ) + resultants += weight * ( + _strain_matrix(z).T @ material_stress[:, point_slice] + ) + + shear_strain = np.asarray( + [ + ( + np.zeros(self._n_shell_points) + if np.isscalar(shell_strain[index]) + else shell_strain[index] + ) + for index in [6, 7] + ] + ) + shear_stiffness = sum( + thickness * matrix + for thickness, matrix in zip(self.list_thickness, self._shear_matrices) + ) + resultants[6:8] = self.k * shear_stiffness @ shear_strain + return type(shell_strain)(resultants) + + def update(self, assembly, pb): + shell_strain = assembly.sv["ShellStrain"] + if np.isscalar(shell_strain) and shell_strain == 0: + shell_strain = [np.zeros(self._n_shell_points) for _ in range(8)] + + for layer, (material, material_assembly) in enumerate( + zip(self.materials, self._material_assemblies) + ): + material_assembly.sv["Strain"] = self._material_strain(shell_strain, layer) + material.update(material_assembly, pb) + + assembly.sv["_ShellStiffnessMatrix"] = self._integrate_tangent() + assembly.sv["ShellStress"] = self._integrate_stress(shell_strain) + + def set_start(self, assembly, pb): + for material, material_assembly in zip( + self.materials, self._material_assemblies + ): + material.set_start(material_assembly, pb) + material_assembly.sv_start = _copy_state(material_assembly.sv) + + def to_start(self, assembly, pb): + for material, material_assembly in zip( + self.materials, self._material_assemblies + ): + material.to_start(material_assembly, pb) + material_assembly.sv = _copy_state(material_assembly.sv_start) + + def get_stress_distribution(self, assembly, pg, resolution=None): + """Return layerwise stresses at one shell Gauss point. + + Without ``resolution``, values are returned at the actual material + integration points. With ``resolution``, each layer is sampled + independently so stress jumps at material interfaces are preserved. + """ + if self._material_assemblies is None: + raise RuntimeError("The shell law has not been initialized.") + if not 0 <= pg < self._n_shell_points: + raise IndexError("pg is outside the shell Gauss-point range.") + if resolution is not None and resolution < 2: + raise ValueError("resolution must be at least two.") + + shell_strain = assembly.sv["ShellStrain"] + shear_strain = np.array( + [ + (0.0 if np.isscalar(shell_strain[index]) else shell_strain[index][pg]) + for index in [6, 7] + ] + ) + z_values = [] + stress_values = [] + for layer, material_assembly in enumerate(self._material_assemblies): + layer_stress = ( + material_assembly.sv["Stress"] + .asarray()[ + :, + pg :: self._n_shell_points, + ] + .copy() + ) + layer_stress[4:6] = self._shear_matrices[layer] @ shear_strain[:, None] + layer_z = self._layer_z[layer] + + if resolution is not None: + output_z = np.linspace( + self._interfaces[layer], + self._interfaces[layer + 1], + resolution, + ) + layer_stress = np.array( + [ + np.interp(output_z, layer_z, component) + for component in layer_stress + ] + ) + layer_z = output_z + + z_values.append(layer_z) + stress_values.append(layer_stress) + + return ( + np.concatenate(z_values), + StressTensorList(np.concatenate(stress_values, axis=1)), + ) + + class ShellLaminate(ShellBase): def __init__(self, listMat, list_thickness, k=1, name=""): # assert get_Dimension() == '3D', "No 2D model for a shell kinematic. Choose '3D' problem dimension." @@ -379,7 +1004,7 @@ def get_stress(self, assembly, **kargs): return StressTensorList(Stress) - def GetStressDistribution(self, assembly, pg, resolution=100): + def get_stress_distribution(self, assembly, pg, resolution=100): h = self.thickness z = np.linspace(-h / 2, h / 2, resolution) @@ -387,14 +1012,16 @@ def GetStressDistribution(self, assembly, pg, resolution=100): ShellStrain = assembly.sv["ShellStrain"] Strain[0] = ( - ShellStrain[0][pg] + z * ShellStrain[4][pg] + ShellStrain[0][pg] + z * ShellStrain[3][pg] ) # epsXX -> membrane and bending Strain[1] = ( - ShellStrain[1][pg] - z * ShellStrain[3][pg] + ShellStrain[1][pg] + z * ShellStrain[4][pg] ) # epsYY -> membrane and bending - Strain[3] = ShellStrain[2][pg] * np.ones_like(z) # 2epsXY + Strain[3] = ( + ShellStrain[2][pg] + z * ShellStrain[5][pg] + ) # 2epsXY -> membrane and twist Strain[4] = ShellStrain[6][pg] * np.ones_like(z) # 2epsXZ -> shear - Strain[5] = ShellStrain[6][pg] * np.ones_like(z) # 2epsYZ -> shear + Strain[5] = ShellStrain[7][pg] * np.ones_like(z) # 2epsYZ -> shear layer_z = [ list((pos - self.__layer) <= 0).index(True) - 1 for pos in z diff --git a/fedoo/constitutivelaw/viso_elastic_orthotropic.py b/fedoo/constitutivelaw/viso_elastic_orthotropic.py deleted file mode 100644 index 7a136083..00000000 --- a/fedoo/constitutivelaw/viso_elastic_orthotropic.py +++ /dev/null @@ -1,112 +0,0 @@ -# derive de ConstitutiveLaw -####WARNING: not working constitutive law - -from fedoo.core.mechanical3d import Mechanical3D -import scipy as sp - - -class ViscoElasticComposites(Mechanical3D): - def __init__( - self, - EL, - ET, - GLT, - GTT, - nuLT, - nuTT, - CL=0, - CT=0, - CLT=0, - RefStrainRate=1, - SLc_T=None, - SLc_C=None, - SYc_T=None, - SYc_C=None, - SZc_T=None, - SZc_C=None, - SLYc=None, - SLZc=None, - name="", - ): - Mechanical3D.__init__(self, name) # heritage - - self.__parameters = { - "EL": EL, - "ET": ET, - "GLT": GLT, - "GTT": GTT, - "nuLT": nuLT, - "nuTT": nuTT, - "CL": CL, - "CT": CT, - "CLT": CLT, - "RefStrainRate": RefStrainRate, - "SLc_T": SLc_T, - "SLc_C": SLc_C, - "SYc_T": SYc_T, - "SYc_C": SYc_C, - "SZc_T": SZc_T, - "SZc_C": SZc_C, - "SLYc": SLYc, - "SLZc": SLZc, - } - - self.__DamageVariable = [0, 0, 0, 0, 0, 0] - - def SetStrainRate(self, StrainRate): - self.__StrainRate = StrainRate - - def get_stress(self, localFrame=None): # methode virtuel - # tester si contrainte plane ou def plane - # if get_Dimension() == "2Dstress": - # print('ViscoElasticComposites law for 2Dstress is not implemented') - # return NotImplemented - - for key in self.__parameters: - exec(key + '= self.__parameters["' + key + '"]') - StrainRate = self.__StrainRate - - if isinstance(EL, (float, int, np.number)): - H = sp.empty((6, 6)) - elif isinstance(EL, (sp.ndarray, list)): - H = sp.zeros((6, 6, len(EL))) - else: - H = sp.zeros((6, 6), dtype="object") - - d1, d2, d3, d4, d5, d6 = self.__DamageVariables - - StrainRateEffect = [ - (StrainRate[i] >= RefStrainRate) * np.log(StrainRate[i] / RefStrainRate) - for i in range(6) - ] - - nuTL = nuLT * ET / EL - k = 1 - nuTT**2 - 2 * nuLT * nuTL - 2 * nuLT * nuTT * nuTL - H[0, 0] = ((1 - d1) * EL * (1 + CL * StrainRateEffect[0])) * (1 - nuTT**2) / k - H[1, 1] = ( - ((1 - d2) * ET * (1 + CT * StrainRateEffect[1])) * (1 - nuLT * nuTL) / k - ) - H[2, 2] = ( - ((1 - d3) * ET * (1 + CT * StrainRateEffect[2])) * (1 - nuLT * nuTL) / k - ) - H[0, 1] = H[1, 0] = H[0, 2] = H[2, 0] = ( - ((1 - d1) * EL * (1 + CL * StrainRateEffect[0])) * (nuTT * nuTL + nuTL) / k - ) - H[1, 2] = H[2, 1] = ( - ((1 - d2) * ET * (1 + CT * StrainRateEffect[1])) * (nuLT * nuTL + nuTT) / k - ) - H[3, 3] = (1 - d4) * GLT * (1 + CLT * StrainRateEffect[3]) - H[4, 4] = (1 - d5) * GLT * (1 + CLT * StrainRateEffect[4]) - H[5, 5] = (1 - d6) * GTT * (1 + CLT * StrainRateEffect[5]) - - H = self._ConsitutiveLaw__ChangeBasisH(H) - - # eps, eps_vir = GetStrainOperator() - sigma = [sum([eps[j] * H[i][j] for j in range(6)]) for i in range(6)] - - return sigma # list de 6 objets de type DiffOp - - def updateDamage(self): - for key in self.__parameters: - exec(key + '= self.__parameters["' + key + '"]') - return NotImplemented diff --git a/fedoo/core/boundary_conditions.py b/fedoo/core/boundary_conditions.py index c081b52b..6972ed23 100644 --- a/fedoo/core/boundary_conditions.py +++ b/fedoo/core/boundary_conditions.py @@ -440,6 +440,10 @@ def create( if isinstance(variable, list): if np.isscalar(value): value = [value for var in variable] + if start_value is None or np.isscalar(start_value): + start_values = [start_value for var in variable] + else: + start_values = start_value return ListBC( [ @@ -449,7 +453,7 @@ def create( var, value[i], time_func=time_func, - start_value=start_value, + start_value=start_values[i], name=name, ) for i, var in enumerate(variable) diff --git a/fedoo/lib_elements/element_list.py b/fedoo/lib_elements/element_list.py index 7e8f1150..7ccbbb61 100644 --- a/fedoo/lib_elements/element_list.py +++ b/fedoo/lib_elements/element_list.py @@ -102,6 +102,9 @@ def get_node_elm_coordinates(element, nNd_elm=None): return np.c_[[0.0, 1.0]] elif nNd_elm == 3: return np.c_[[0.0, 1.0, 0.5]] + elif element == "lin2interface": + if nNd_elm == 4: + return np.c_[[0.0, 1.0, 0.0, 1.0]] elif element in ["tri3", "tri6", "tri3bubble"]: if nNd_elm == 3: return np.c_[[0.0, 1.0, 0.0], [0.0, 0.0, 1.0]] diff --git a/fedoo/weakform/interface_force.py b/fedoo/weakform/interface_force.py index 8cec444f..33ddbd78 100644 --- a/fedoo/weakform/interface_force.py +++ b/fedoo/weakform/interface_force.py @@ -9,7 +9,7 @@ class InterfaceForce(WeakFormBase): Weak formulation of the interface equilibrium equation. * Require an interface constitutive law such as :mod:`fedoo.constitutivelaw.CohesiveLaw` or :mod:`fedoo.constitutivelaw.Spring` - * Geometrical non linearities not implemented + * Updated-Lagrangian geometrical nonlinearities are supported. Parameters ---------- @@ -18,12 +18,14 @@ class InterfaceForce(WeakFormBase): (:mod:`fedoo.constitutivelaw`) name: str, optional name of the WeakForm - nlgeom: bool (default = False) - For future development - If True, return a NotImplemented Error + nlgeom: {None, bool, "UL"}, optional + Geometrical-nonlinearity formulation. If ``None`` (default), use the + value defined by the associated problem. ``True`` and ``"UL"`` select + the updated-Lagrangian formulation. Total Lagrangian (``"TL"``) is not + implemented for interface forces. """ - def __init__(self, constitutivelaw, name="", nlgeom=False, space=None): + def __init__(self, constitutivelaw, name="", nlgeom=None, space=None): if isinstance(constitutivelaw, str): constitutivelaw = ConstitutiveLaw[constitutivelaw] @@ -39,27 +41,17 @@ def __init__(self, constitutivelaw, name="", nlgeom=False, space=None): self.constitutivelaw = constitutivelaw - self.nlgeom = nlgeom # geometric non linearities -> False, True, 'UL' or 'TL' (True or 'UL': updated lagrangian - 'TL': total lagrangian) - """Method used to treat the geometric non linearities. - * Set to False if geometric non linarities are ignored (default). - * Set to True or 'UL' to use the updated lagrangian method (update the mesh) - * Set to 'TL' to use the total lagrangian method (base on the initial mesh with initial displacement effet) - """ + self.nlgeom = nlgeom self.assembly_options["assume_sym"] = False # symetric ? def initialize(self, assembly, pb): - if self.nlgeom: - if self.nlgeom is True: - self.nlgeom = "UL" - elif isinstance(self.nlgeom, str): - self.nlgeom = self.nlgeom.upper() - if self.nlgeom != "UL": - raise NotImplementedError( - f"{self.nlgeom} nlgeom not implemented for Interface force." - ) - else: - raise TypeError("nlgeom should be in {'TL', 'UL', True, False}") + self._initialize_nlgeom(assembly, pb) + self.nlgeom = assembly._nlgeom + if self.nlgeom == "TL": + raise NotImplementedError( + "TL nlgeom is not implemented for InterfaceForce. Use 'UL'." + ) def update(self, assembly, pb): # function called when the problem is updated (NR loop or time increment) diff --git a/tests/test_Laminate.py b/tests/test_Laminate.py index c3c0a6b8..96e4dd24 100644 --- a/tests/test_Laminate.py +++ b/tests/test_Laminate.py @@ -79,7 +79,9 @@ def test_laminate(): assert np.abs(pb.get_disp("DispZ")[node_right_center] + 25.698414360666018) < 1e-7 # plot the stress distribution - # z, StressDistribution = fd.ConstitutiveLaw['PlateSection'].GetStressDistribution(fd.Assembly['plate'],200) + # z, stress = fd.ConstitutiveLaw["PlateSection"].get_stress_distribution( + # fd.Assembly["plate"], 200 + # ) # plt.plot(StressDistribution[0], z) diff --git a/tests/test_cohesive_interface_2d.py b/tests/test_cohesive_interface_2d.py new file mode 100644 index 00000000..095dd10a --- /dev/null +++ b/tests/test_cohesive_interface_2d.py @@ -0,0 +1,46 @@ +import numpy as np + +import fedoo as fd + + +def test_lin2interface_with_cohesive_law(): + fd.ModelingSpace("2D") + + nodes = np.array( + [ + [0.0, 0.0], + [1.0, 0.0], + [0.0, 0.0], + [1.0, 0.0], + ] + ) + mesh = fd.Mesh( + nodes, + np.array([[0, 1, 2, 3]]), + "lin2interface", + ) + law = fd.constitutivelaw.CohesiveLaw( + axis=1, + tangent_mode="consistent", + ) + assembly = fd.Assembly.create( + fd.weakform.InterfaceForce(law), + mesh, + ) + problem = fd.problem.NonLinear(assembly) + problem.bc.add("Dirichlet", [0, 1], "Disp", 0.0) + problem.bc.add("Dirichlet", [2, 3], "DispX", 0.0) + problem.bc.add("Dirichlet", [2, 3], "DispY", 0.008) + + problem.nlsolve( + dt=1.0, + tmax=1.0, + tol_nr=1.0e-8, + print_info=0, + ) + + assert np.allclose(problem.get_disp()[1, 2:], 0.008) + assert np.allclose(assembly.sv["DamageVariable"], 0.625) + tangent = assembly.sv["TangentMatrix"] + assert np.asarray(tangent[0][0]).shape == (2,) + assert np.asarray(tangent[1][1]).shape == (2,) diff --git a/tests/test_cohesivelaw.py b/tests/test_cohesivelaw.py new file mode 100644 index 00000000..afe497e3 --- /dev/null +++ b/tests/test_cohesivelaw.py @@ -0,0 +1,236 @@ +from types import SimpleNamespace + +import numpy as np + +from fedoo.constitutivelaw import CohesiveLaw + + +def _matrix_at_material_point(matrix): + return np.array( + [[np.asarray(value).reshape(-1)[0] for value in row] for row in matrix], + dtype=float, + ) + + +def _evaluate_traction_and_tangent(law, delta, irreversible=0.0): + """Evaluate one material point with a fixed committed damage state.""" + assembly = SimpleNamespace( + sv={ + "DamageVariable": irreversible, + "DamageVariableOpening": irreversible, + "DamageVariableIrreversible": irreversible, + } + ) + delta_arrays = [np.array([component], dtype=float) for component in delta] + damage_gradient = law._update_damage(assembly, delta_arrays) + + secant = _matrix_at_material_point(law.get_secant_matrix(assembly)) + tangent = _matrix_at_material_point( + law.get_tangent_matrix(assembly, delta_arrays, damage_gradient) + ) + traction = secant @ np.asarray(delta, dtype=float) + return traction, tangent + + +def _finite_difference_tangent(law, delta, irreversible=0.0, step=1.0e-8): + delta = np.asarray(delta, dtype=float) + tangent = np.empty((3, 3)) + for component in range(3): + perturbation = np.zeros(3) + perturbation[component] = step + traction_plus, _ = _evaluate_traction_and_tangent( + law, delta + perturbation, irreversible + ) + traction_minus, _ = _evaluate_traction_and_tangent( + law, delta - perturbation, irreversible + ) + tangent[:, component] = (traction_plus - traction_minus) / (2.0 * step) + return tangent + + +def test_cohesive_law_zero_displacement_update(): + law = CohesiveLaw(KI=1.0e4, KII=2.0e4) + assembly = SimpleNamespace(sv={}) + problem = SimpleNamespace(get_dof_solution=lambda: 0) + + law.initialize(assembly, problem) + law.update(assembly, problem) + + assert assembly.sv["InterfaceStress"] == 0 + assert assembly.sv["RelativeDisp"] == 0 + assert np.diag(assembly.sv["TangentMatrix"]).tolist() == [ + 2.0e4, + 2.0e4, + 1.0e4, + ] + + +def test_cohesive_law_update_uses_current_assembly_local_frame(): + law = CohesiveLaw(KI=1.0e4, KII=2.0e4) + delta = [ + np.array([0.001]), + np.array([0.002]), + np.array([0.003]), + ] + operators = [object(), object(), object()] + current = SimpleNamespace( + space=SimpleNamespace(op_disp=lambda: operators), + get_gp_results=lambda operator, displacement: delta[operators.index(operator)], + ) + assembly = SimpleNamespace( + sv={ + "DamageVariable": 0, + "DamageVariableOpening": 0, + "DamageVariableIrreversible": 0, + }, + current=current, + ) + problem = SimpleNamespace(get_dof_solution=lambda: np.ones(1)) + + law.update(assembly, problem) + + assert all( + computed is expected + for computed, expected in zip(assembly.sv["RelativeDisp"], delta) + ) + assert np.allclose( + np.asarray(assembly.sv["InterfaceStress"]).reshape(3), + [20.0, 40.0, 30.0], + ) + + +def test_cohesive_law_mode_i_damage_and_compression_contact(): + law = CohesiveLaw(GIc=0.3, SImax=60.0, KI=1.0e4, axis=2) + assembly = SimpleNamespace( + sv={ + "DamageVariable": 0, + "DamageVariableOpening": 0, + "DamageVariableIrreversible": 0, + } + ) + + # Damage starts at delta_0 = SImax / KI = 0.006 and reaches one at + # delta_m = 2 GIc / SImax = 0.01. + law._update_damage( + assembly, + [np.zeros(3), np.zeros(3), np.array([0.0, 0.008, 0.011])], + ) + + assert np.allclose(assembly.sv["DamageVariable"], [0.0, 0.625, 1.0]) + assert np.allclose( + assembly.sv["DamageVariableOpening"], assembly.sv["DamageVariable"] + ) + fully_open_stiffness = law.get_tangent_matrix(assembly) + assert np.allclose( + [fully_open_stiffness[i][i][2] for i in range(3)], + 0.0, + ) + + law.update_irreversible_damage(assembly) + law._update_damage( + assembly, + [np.zeros(3), np.zeros(3), np.array([-0.001, -0.001, -0.001])], + ) + + # Damage is irreversible, while normal stiffness is restored in + # compression to provide the cohesive law's soft-contact response. + assert np.allclose(assembly.sv["DamageVariable"], [0.0, 0.625, 1.0]) + assert np.allclose(assembly.sv["DamageVariableOpening"], 0.0) + assert np.allclose( + law.get_tangent_matrix(assembly)[2][2], + np.full(3, 1.0e4), + ) + closed_stiffness = law.get_tangent_matrix(assembly) + assert np.allclose( + [closed_stiffness[i][i][2] for i in range(3)], + [0.0, 0.0, 1.0e4], + ) + + +def test_cohesive_law_consistent_tangent_in_mixed_mode(): + law = CohesiveLaw(tangent_mode="consistent") + delta = np.array([0.003, 0.002, 0.007]) + + _, tangent = _evaluate_traction_and_tangent(law, delta) + finite_difference = _finite_difference_tangent(law, delta) + + assert np.allclose(tangent, finite_difference, rtol=2.0e-6, atol=2.0e-4) + + +def test_cohesive_law_consistent_tangent_in_pure_mode_i(): + law = CohesiveLaw(tangent_mode="consistent") + delta = np.array([0.0, 0.0, 0.008]) + + _, tangent = _evaluate_traction_and_tangent(law, delta) + finite_difference = _finite_difference_tangent(law, delta) + + # The tangential norm has a cusp at exactly zero; its centred finite + # difference converges one-sidedly and therefore needs a slightly looser + # tolerance than the smooth mixed-mode point. + assert np.allclose(tangent, finite_difference, rtol=1.0e-5, atol=2.0e-4) + assert np.isclose(tangent[2, 2], -1.5e4) + + +def test_cohesive_law_consistent_tangent_in_mode_ii_compression(): + law = CohesiveLaw(tangent_mode="consistent") + delta = np.array([0.008, 0.003, -0.001]) + + _, tangent = _evaluate_traction_and_tangent(law, delta) + finite_difference = _finite_difference_tangent(law, delta) + + assert np.allclose(tangent, finite_difference, rtol=2.0e-6, atol=2.0e-4) + assert tangent[2, 2] == law.parameters["KI"] + + +def test_cohesive_law_consistent_tangent_uses_secant_during_unloading(): + law = CohesiveLaw(tangent_mode="consistent") + delta = np.array([0.0, 0.0, 0.007]) + + _, tangent = _evaluate_traction_and_tangent(law, delta, irreversible=0.6) + finite_difference = _finite_difference_tangent(law, delta, irreversible=0.6) + + assert np.allclose(tangent, finite_difference) + assert np.allclose(np.diag(tangent), [2.0e4, 2.0e4, 4.0e3]) + + +def test_cohesive_law_commits_damage_with_secant_predictor(): + law = CohesiveLaw(tangent_mode="consistent") + assembly = SimpleNamespace( + sv={ + "DamageVariable": np.array([0.625]), + "DamageVariableOpening": np.array([0.625]), + "DamageVariableIrreversible": 0.0, + "TangentMatrix": [[0.0, 0.0, 0.0]] * 3, + } + ) + + law.set_start(assembly, None) + + assert np.allclose( + np.diag(_matrix_at_material_point(assembly.sv["TangentMatrix"])), + [1.875e4, 1.875e4, 3.75e3], + ) + assert np.allclose(assembly.sv["DamageVariableIrreversible"], np.array([0.625])) + + +def test_cohesive_law_secant_mode_and_validation(): + assert CohesiveLaw().tangent_mode == "secant" + + law = CohesiveLaw(tangent_mode="secant") + delta = np.array([0.003, 0.002, 0.007]) + _, tangent = _evaluate_traction_and_tangent(law, delta) + + _, consistent_tangent = _evaluate_traction_and_tangent( + CohesiveLaw(tangent_mode="consistent"), delta + ) + assert np.allclose(tangent, np.diag(np.diag(tangent))) + assert not np.allclose(tangent, consistent_tangent) + + with np.testing.assert_raises_regex(ValueError, "tangent_mode"): + CohesiveLaw(tangent_mode="invalid") + + +if __name__ == "__main__": + import pytest + + pytest.main([__file__]) diff --git a/tests/test_elasto_plasticity.py b/tests/test_elasto_plasticity.py new file mode 100644 index 00000000..588db2e8 --- /dev/null +++ b/tests/test_elasto_plasticity.py @@ -0,0 +1,98 @@ +import numpy as np + +from fedoo.constitutivelaw import ElastoPlasticity +from fedoo.util.voigt_tensors import StrainTensorList + + +E = 200e3 +NU = 0.3 +SIGMA_Y = 300.0 +H = 1000.0 +BETA = 0.3 + + +def make_material(): + material = ElastoPlasticity(E, NU, SIGMA_Y) + material.set_hardening_function("power", h=H, beta=BETA) + return material + + +def test_radial_return_matches_simcoon_simple_shear_reference(): + material = make_material() + plasticity = np.zeros(1) + plastic_strain = np.zeros((6, 1)) + + for gamma in np.linspace(0.0, 0.02, 21)[1:]: + strain = np.zeros((6, 1)) + strain[3, 0] = gamma + stress = material.compute_stress( + StrainTensorList(strain), + plasticity, + plastic_strain, + ) + plasticity = material.get_plasticity().copy() + plastic_strain = material.get_plastic_strain().asarray(copy=True) + + # Reference values produced by Simcoon's EPICP law for the same 20 + # proportional shear increments. + np.testing.assert_allclose(stress[3], [314.581143696], rtol=1e-9) + np.testing.assert_allclose(material.get_plasticity(), [0.00918589977986], rtol=1e-9) + np.testing.assert_allclose( + material.yield_function(stress, material.get_plasticity()), + [0.0], + atol=1e-6, + ) + + +def test_plastic_tangent_uses_positive_hardening_denominator(): + material = make_material() + strain = np.zeros((6, 1)) + strain[3, 0] = 0.02 + material.compute_stress(StrainTensorList(strain)) + + plasticity = material.get_plasticity()[0] + hardening_slope = BETA * H * plasticity ** (BETA - 1) + shear_modulus = E / (2 * (1 + NU)) + expected_shear_tangent = ( + shear_modulus * hardening_slope / (3 * shear_modulus + hardening_slope) + ) + + tangent = material.get_tangent_matrix() + np.testing.assert_allclose(tangent[3, 3, 0], expected_shear_tangent, rtol=1e-12) + + +def test_vectorized_integration_matches_independent_material_points(): + gammas = np.array([0.0, 0.001, 0.003, 0.01, 0.02]) + strains = np.zeros((6, len(gammas))) + strains[3] = gammas + + vectorized_material = make_material() + vectorized_stress = vectorized_material.compute_stress(StrainTensorList(strains)) + + independent_stress = np.empty_like(strains) + for point, gamma in enumerate(gammas): + material = make_material() + strain = np.zeros((6, 1)) + strain[3, 0] = gamma + independent_stress[:, point] = material.compute_stress( + StrainTensorList(strain) + ).asarray()[:, 0] + + np.testing.assert_allclose( + vectorized_stress.asarray(), + independent_stress, + rtol=1e-12, + atol=1e-12, + ) + + +def test_legacy_commit_api_was_removed(): + material = make_material() + assert not hasattr(material, "NewTimeIncrement") + assert not hasattr(material, "ComputeStress") + + +if __name__ == "__main__": + import pytest + + pytest.main([__file__]) diff --git a/tests/test_interface_force_nlgeom.py b/tests/test_interface_force_nlgeom.py new file mode 100644 index 00000000..9bbf7b2d --- /dev/null +++ b/tests/test_interface_force_nlgeom.py @@ -0,0 +1,66 @@ +from types import SimpleNamespace + +import pytest + +import fedoo as fd + + +class _AssemblyStub: + def __init__(self): + self._nlgeom = False + self.updated_displacement = None + + def set_disp(self, displacement): + self.updated_displacement = displacement + + +def _interface_force(nlgeom=None): + fd.ModelingSpace("3D") + law = fd.constitutivelaw.Spring(Kx=1.0, Ky=1.0, Kz=1.0) + return fd.weakform.InterfaceForce(law, nlgeom=nlgeom) + + +def test_interface_force_inherits_problem_nlgeom(): + weakform = _interface_force() + assembly = _AssemblyStub() + displacement = object() + problem = SimpleNamespace( + nlgeom=True, + get_disp=lambda: displacement, + ) + + weakform.initialize(assembly, problem) + weakform.update(assembly, problem) + + assert assembly._nlgeom == "UL" + assert weakform.nlgeom == "UL" + assert assembly.updated_displacement is displacement + + +def test_interface_force_explicit_false_overrides_problem_nlgeom(): + weakform = _interface_force(nlgeom=False) + assembly = _AssemblyStub() + problem = SimpleNamespace( + nlgeom=True, + get_disp=lambda: object(), + ) + + weakform.initialize(assembly, problem) + weakform.update(assembly, problem) + + assert assembly._nlgeom is False + assert weakform.nlgeom is False + assert assembly.updated_displacement is None + + +def test_interface_force_rejects_total_lagrangian_formulation(): + weakform = _interface_force(nlgeom="TL") + assembly = _AssemblyStub() + problem = SimpleNamespace(nlgeom=False) + + with pytest.raises(NotImplementedError, match="InterfaceForce"): + weakform.initialize(assembly, problem) + + +if __name__ == "__main__": + pytest.main([__file__])