Axisymmetric (2Daxi) fixes and improvements

examples/03-advanced/thick_cylinder_inflation_axi.py

@@ -0,0 +1,123 @@
+"""
+Thick cylinder inflation: F_theta-theta verification (2Daxi)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Pinning benchmark for the canonical hoop deformation gradient
+:math:`F_{\\theta\\theta} = r_{\\text{current}} / R_{\\text{reference}}`
+in a finite-strain updated-lagrangian axisymmetric problem.
+
+A thick-walled annulus is inflated by prescribing a radial Dirichlet
+displacement on the outer boundary (the inner boundary is clamped
+radially). The simulation runs in ``2Daxi`` with ``nlgeom='UL'``, then
+the example reads the deformation gradient stored on the assembly's
+state vector and asserts that the hoop component matches the textbook
+relation :math:`F_{\\theta\\theta} = r/R` (Bonet & Wood, Box 8.3;
+Holzapfel Sec. 2.5).
+
+The assertion is what protects against regressions in
+``_comp_grad_disp``: an earlier fedoo implementation divided ``u_r``
+by ``r_current`` instead of the reference radius ``R``, giving the
+small-strain limit ``1 + u_r/r`` rather than the correct ``r/R``. At
+the inflation level used here (10-30% hoop stretch) the two formulas
+disagree by several percent, so this example fails loudly if the bug
+returns.
+
+This example also demonstrates how to extract the per-Gauss-point
+deformation gradient from ``assembly.sv``.
+"""
+
+import fedoo as fd
+import numpy as np
+
+# ---------------------------------------------------------------------------
+# 1. Mesh: a thick annulus in the (r, z) plane.
+# ---------------------------------------------------------------------------
+# r in [A, B], z in [0, H]. The symmetry axis is Y in 2D (i.e. r = X column).
+A, B = 1.0, 2.0
+H = 0.5
+
+fd.ModelingSpace("2Daxi")
+mesh = fd.mesh.rectangle_mesh(nx=21, ny=11, x_min=A, x_max=B, y_min=0.0, y_max=H)
+
+# ---------------------------------------------------------------------------
+# 2. Linear elastic isotropic constitutive law.
+# ---------------------------------------------------------------------------
+# Linear elasticity is sufficient: the F_theta-theta = r/R relation is a pure
+# kinematic identity, independent of the constitutive law. We only need the
+# solver to converge to the prescribed boundary motion.
+E = 1.0e3
+nu = 0.3
+material = fd.constitutivelaw.ElasticIsotrop(E, nu)
+
+wf = fd.weakform.StressEquilibrium(material)
+assembly = fd.Assembly.create(wf, mesh)
+
+# ---------------------------------------------------------------------------
+# 3. Updated-Lagrangian non-linear problem.
+# ---------------------------------------------------------------------------
+pb = fd.problem.NonLinear(assembly, nlgeom="UL")
+pb.set_nr_criterion("Displacement", tol=1e-6, max_subiter=15)
+
+# Inner radius: clamped radially. Top and bottom: clamped axially to keep
+# the problem axisymmetric and 2D (no plane-strain ambiguity).
+left = mesh.node_sets["left"]  # r = A
+right = mesh.node_sets["right"]  # r = B
+bottom = mesh.node_sets["bottom"]
+top = mesh.node_sets["top"]
+
+# Inner radius radially fixed -> drives a non-uniform F_theta-theta(R).
+pb.bc.add("Dirichlet", left, "DispX", 0.0)
+# Symmetry / axial confinement to keep things 2D.
+pb.bc.add("Dirichlet", bottom, "DispY", 0.0)
+pb.bc.add("Dirichlet", top, "DispY", 0.0)
+# Outer radius pulled outward by delta -> a rich F_theta-theta field.
+delta = 0.6  # 30% radial expansion at r = B
+pb.bc.add("Dirichlet", right, "DispX", delta)
+
+pb.nlsolve(dt=0.1, tmax=1.0, update_dt=True, print_info=0)
+
+# ---------------------------------------------------------------------------
+# 4. Verify F_theta-theta = r_current / R_reference at every gauss point.
+# ---------------------------------------------------------------------------
+# F is stored as (3, 3, n_gauss_points), Fortran-ordered.
+F = assembly.sv["F"]
+F_tt = F[2, 2, :]
+
+# R0 is the reference radius captured at problem initialize.
+R0 = assembly.sv["_R0_gausspoints"]
+
+# Current r at gauss points: same interpolation, but on the deformed mesh.
+r_current = assembly.current.mesh.convert_data(
+    assembly.current.mesh.nodes[:, 0],
+    "Node",
+    "GaussPoint",
+    n_elm_gp=assembly.n_elm_gp,
+)
+
+F_tt_expected = r_current / R0
+
+err_max = np.max(np.abs(F_tt - F_tt_expected))
+print(f"Inflation: delta = {delta}, max |F_tt - r/R| = {err_max:.3e}")
+print(f"Hoop stretch range: [{F_tt.min():.4f}, {F_tt.max():.4f}]")
+print(f"Reference radius range: [{R0.min():.4f}, {R0.max():.4f}]")
+
+# Tight tolerance: F_theta-theta is computed *exactly* from u_r and R0 by
+# `_comp_grad_disp`; any discrepancy beyond fp noise indicates a regression.
+assert err_max < 1e-10, (
+    f"F_theta-theta drifted from r/R by {err_max:.3e}; "
+    "this is a regression in fedoo/weakform/stress_equilibrium._comp_grad_disp."
+)
+
+# ---------------------------------------------------------------------------
+# 5. Negative control: the *buggy* formula 1 + u_r / r_current would give a
+#    different (incorrect) F_theta-theta. Show how badly it would have
+#    differed at the stretches reached here.
+# ---------------------------------------------------------------------------
+u_r_gp = r_current - R0
+F_tt_buggy = 1.0 + u_r_gp / r_current  # the pre-fix formula
+disagreement = np.max(np.abs(F_tt_buggy - F_tt_expected))
+print(
+    "Disagreement between the (correct) r/R and the (buggy) 1 + u_r/r "
+    f"formulas at this stretch: {disagreement:.3e} "
+    f"({100 * disagreement / F_tt_expected.max():.2f}% of peak hoop stretch)."
+)

examples/03-advanced/tube_compression_ipc.py

@@ -0,0 +1,97 @@
+"""
+Compression of a tube using 2D axisymmetric model — IPC self-contact
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Same problem as ``tube_compression.py`` but with ``IPCSelfContact``
+instead of the penalty ``SelfContact``. Used to benchmark the
+2π·r-weighted IPC formulation against the penalty reference.
+
+Run after ``tube_compression.py`` so both result folders exist; this
+script prints reduced metrics from each run for side-by-side
+comparison.
+"""
+
+import os
+
+import fedoo as fd
+import numpy as np
+
+fd.ModelingSpace("2Daxi")
+mesh = fd.mesh.rectangle_mesh(5, 240, 23, 25, 0, 180)
+
+sigma_y = 300
+k = 1000
+m = 0.3
+E = 200e3
+nu = 0.3
+props = np.array([E, nu, 1e-5, sigma_y, k, m])
+material = fd.constitutivelaw.Simcoon("EPICP", props)
+
+NLGEOM = "UL"
+wf = fd.weakform.StressEquilibrium(material, nlgeom=NLGEOM)
+solid_assembly = fd.Assembly.create(wf, mesh)
+
+# --- IPC self-contact (replaces penalty SelfContact) ---
+contact = fd.constraint.IPCSelfContact(
+    mesh,
+    dhat=1e-3,
+    dhat_is_relative=True,
+    use_ccd=True,
+)
+
+assembly = fd.Assembly.sum(solid_assembly, contact)
+
+pb = fd.problem.NonLinear(assembly, nlgeom=NLGEOM)
+pb.set_nr_criterion(
+    "Displacement",
+    tol=1e-2,
+    max_subiter=20,
+    adaptive_stiffness=True,
+)
+
+if not os.path.isdir("results"):
+    os.mkdir("results")
+res_ipc = pb.add_output(
+    "results/tube_compression_ipc",
+    solid_assembly,
+    ["Disp", "Stress", "Strain", "P"],
+)
+
+bottom = mesh.node_sets["bottom"]
+top = mesh.node_sets["top"]
+
+pb.bc.add("Dirichlet", bottom, "Disp", 0)
+pb.bc.add("Dirichlet", top, "Disp", [0, -150])
+pb.add_line_search()
+pb.nlsolve(dt=0.01, tmax=1, update_dt=True, print_info=0, dt_min=1e-8)
+
+###############################################################################
+# Reduced numeric metrics — peak axial stress, peak equivalent plastic strain,
+# and final top displacement at the last saved iteration.
+
+
+def _summarise(res, label):
+    res.load(-1)
+    stress_yy = np.asarray(res.get_data("Stress", component="YY", data_type="Node"))
+    p = np.asarray(res.get_data("P", data_type="Node"))
+    disp_y = np.asarray(res.get_data("Disp", component="Y", data_type="Node"))
+    print(f"--- {label} ---")
+    print(f"  peak |Stress YY|     = {np.abs(stress_yy).max():.4e}")
+    print(f"  peak P (eq. plastic) = {p.max():.4e}")
+    print(f"  min Disp Y           = {disp_y.min():.4e}")
+    print(f"  max Disp Y           = {disp_y.max():.4e}")
+
+
+print()
+_summarise(res_ipc, "IPC (this run)")
+
+# Compare against the penalty reference if it has been run.
+penalty_path = "results/tube_compressoin.fdz"  # original spelling kept as in ref script
+if os.path.isfile(penalty_path):
+    res_penalty = fd.DataSet.read(penalty_path)
+    _summarise(res_penalty, "Penalty (reference)")
+else:
+    print(
+        f"\n(Penalty reference not found at '{penalty_path}'. Run "
+        "tube_compression.py first to populate it for comparison.)"
+    )

fedoo/constitutivelaw/composite_ud.py

@@ -34,6 +34,14 @@ class CompositeUD(ElasticAnisotropic):
         Z direction (if defined, the local material coordinates are used).
     name: str, optional
         The name of the constitutive law
+
+    Notes
+    -----
+    In ``2Daxi``, slot 2 is the hoop direction (not Cartesian z). The
+    fibers default to the X (= radial) direction; rotating with ``angle``
+    around the implicit "Z" axis effectively rotates between r and the
+    hoop direction. See :class:`fedoo.core.mechanical3d.Mechanical3D` for
+    the slot-ordering convention.
     """
 
     def __init__(

fedoo/constitutivelaw/elastic_anisotropic.py

@@ -21,6 +21,14 @@ class ElasticAnisotropic(Mechanical3D):
         If H is a list of gauss point values, the shape shoud be H.shape = (6,6,NumberOfGaussPoints)
     name : str, optional
         The name of the constitutive law
+
+    Notes
+    -----
+    H is interpreted in the slot ordering documented on
+    :class:`fedoo.core.mechanical3d.Mechanical3D`. In ``2Daxi``, slot 2 is
+    the **hoop** direction (not Cartesian z): rows/columns labelled "3" of
+    H apply to the hoop response. Define the rigidity matrix accordingly,
+    or use a local frame.
     """
 
     def __init__(self, H, name=""):

fedoo/constitutivelaw/elastic_orthotropic.py

@@ -27,6 +27,13 @@ class ElasticOrthotropic(ElasticAnisotropic):
         Shear modulus
     nuYZ, nuXZ, nuXY: scalars or arrays of gauss point values
         Poisson's ratio
+
+    Notes
+    -----
+    Material directions X / Y / Z map to fedoo's slot ordering. In ``2Daxi``,
+    slot 2 carries the hoop response, so **EZ, GYZ, GXZ, nuYZ, nuXZ refer
+    to the hoop direction**, not Cartesian z. See
+    :class:`fedoo.core.mechanical3d.Mechanical3D` for details.
     """
 
     def __init__(self, Ex, Ey, Ez, Gyz, Gxz, Gxy, nuyz, nuxz, nuxy, name=""):

fedoo/constitutivelaw/simcoon_umat.py

@@ -21,6 +21,27 @@ class Simcoon(Mechanical3D):
         The constitive laws properties
     name : str
         The name of the constitutive law
+
+    Notes
+    -----
+    UMAT compatibility with the ``2Daxi`` ModelingSpace:
+
+    * **Isotropic UMATs** (e.g. ``ELI``, ``NEOH``, ``MOON``, ``EPICP`` J2
+      plasticity): supported. Hooke's response and J2 invariants are
+      invariant under the slot remapping fedoo applies in 2Daxi
+      (cf. :class:`fedoo.core.mechanical3d.Mechanical3D`).
+
+    * **Orthotropic / anisotropic UMATs** (e.g. ``ELIO``, composite
+      Mori-Tanaka, anisotropic damage / plasticity): user-supplied
+      material direction "3" is silently the **hoop** direction in
+      2Daxi (because slot 2 of the 6-vector carries ε_θθ). Define the
+      stiffness / hardening parameters with this convention or the
+      response will not match the intended material orientation.
+
+    * Hyperelastic laws are gated on plane stress (see
+      ``_Lt_from_F`` branch); they remain compatible with 2Daxi at
+      finite strain provided the F[θθ] = r/R fix is in effect (see
+      :func:`fedoo.weakform.stress_equilibrium._comp_grad_disp`).
     """
 
     def __init__(self, umat_name, props, name=""):

fedoo/constraint/contact.py

@@ -297,7 +297,7 @@ def contact_search(self, contact_list={}, update_contact=True):
         contact_elements = []
         contact_g = []
 
-        is_axi = self.space._dimension == "2Daxi"
+        is_axi = self.space.is_axisymmetric
         if is_axi:
             contact_r = []  # radial coordinate for 2πr weighting