Refactor and improve several constitutive laws

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

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

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

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

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

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

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)

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)