diff --git a/examples/03-advanced/thick_cylinder_inflation_axi.py b/examples/03-advanced/thick_cylinder_inflation_axi.py new file mode 100644 index 00000000..91a96d9e --- /dev/null +++ b/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)." +) diff --git a/examples/03-advanced/tube_compression_ipc.py b/examples/03-advanced/tube_compression_ipc.py new file mode 100644 index 00000000..1a2b3d0f --- /dev/null +++ b/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.)" + ) diff --git a/fedoo/constitutivelaw/composite_ud.py b/fedoo/constitutivelaw/composite_ud.py index 6283ece8..6146b58e 100644 --- a/fedoo/constitutivelaw/composite_ud.py +++ b/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__( diff --git a/fedoo/constitutivelaw/elastic_anisotropic.py b/fedoo/constitutivelaw/elastic_anisotropic.py index 1c8a4406..496afdfb 100644 --- a/fedoo/constitutivelaw/elastic_anisotropic.py +++ b/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=""): diff --git a/fedoo/constitutivelaw/elastic_orthotropic.py b/fedoo/constitutivelaw/elastic_orthotropic.py index cf7642f2..644d205e 100644 --- a/fedoo/constitutivelaw/elastic_orthotropic.py +++ b/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=""): diff --git a/fedoo/constitutivelaw/simcoon_umat.py b/fedoo/constitutivelaw/simcoon_umat.py index a1fbf626..445e687d 100644 --- a/fedoo/constitutivelaw/simcoon_umat.py +++ b/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=""): diff --git a/fedoo/constraint/contact.py b/fedoo/constraint/contact.py index d17c9f1f..430c6ed8 100644 --- a/fedoo/constraint/contact.py +++ b/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 diff --git a/fedoo/constraint/ipc_contact.py b/fedoo/constraint/ipc_contact.py index 762ef186..c4d6c8c1 100644 --- a/fedoo/constraint/ipc_contact.py +++ b/fedoo/constraint/ipc_contact.py @@ -77,6 +77,20 @@ class IPCContact(AssemblyBase): nodes must use ``use_ccd=True`` instead. An error is raised at initialization if ``use_ogc=True`` is used with a solid mesh. + **Axisymmetric (2Daxi)** — Under a ``"2Daxi"`` ModelingSpace the + barrier gradient and hessian returned by ipctk are scaled by the + circumferential integration weight ``2*pi*r`` per surface vertex + before being scattered into the global system. This is the IPC + analogue of fedoo's per-Gauss-point ``axi_volume_weight`` for weak + forms. The scaling is applied per vertex (using the reference + radius ``R₀``); it is exact when both endpoints of a contact pair + sit at similar radii, and slightly asymmetric otherwise. Vertices + on the symmetry axis (r = 0) get weight 0, which is geometrically + correct (a "ring" of zero radius collapses to a point). In 2Daxi + the planar (r, z) distance equals the true 3D minimum distance + between the corresponding circles (closest points share the + azimuth), so the barrier distance and PSD projection are unchanged. + Parameters ---------- mesh : fedoo.Mesh @@ -234,6 +248,10 @@ def __init__( self._ogc_trust_region = None self._ogc_vertices_at_start = None + # Axisymmetric (2Daxi) integration weights — built in initialize() + self._axi_weight = None # 1D array (n_surf*ndim,) interleaved 2*pi*r + self._axi_diag = None # sparse diag(_axi_weight) for hessian scaling + self.sv = {} self.sv_start = {} @@ -410,6 +428,12 @@ def _compute_ipc_contributions(self, vertices, compute="all"): global_vector = P @ ipctk_gradient global_matrix = P @ ipctk_hessian @ P.T + Under a 2Daxi ModelingSpace, the per-surface-DOF gradient/hessian + are first scaled by ``2*pi*r`` (per surface vertex, broadcast over + coordinates) so that the planar barrier energy is integrated + around the symmetry axis — the IPC analogue of fedoo's + ``axi_volume_weight`` for weak forms. + Sets self.global_matrix and self.global_vector. """ n_mesh_dof = self.space.nvar * self.mesh.n_nodes @@ -421,12 +445,16 @@ def _compute_ipc_contributions(self, vertices, compute="all"): ipctk = _import_ipctk() P = self._scatter_matrix + axi_w = self._axi_weight # None outside 2Daxi + axi_D = self._axi_diag # None outside 2Daxi # Barrier contributions if compute != "matrix": grad_surf = self._barrier_potential.gradient( self._collisions, self._collision_mesh, vertices ) + if axi_w is not None: + grad_surf = axi_w * grad_surf # ipctk gradient points toward increasing barrier (toward contact). # The repulsive force is -gradient. In fedoo, global_vector is # added to RHS: K*dX = B + D, so D = -kappa * gradient. @@ -449,6 +477,8 @@ def _compute_ipc_contributions(self, vertices, compute="all"): vertices, project_hessian_to_psd=ipctk.PSDProjectionMethod.NONE, ) + if axi_D is not None: + hess_surf = axi_D @ hess_surf @ axi_D self.global_matrix = self._kappa * (P @ hess_surf @ P.T) # Friction contributions @@ -463,6 +493,8 @@ def _compute_ipc_contributions(self, vertices, compute="all"): self._collision_mesh, vertices, ) + if axi_w is not None: + fric_grad_surf = axi_w * fric_grad_surf self.global_vector += -(P @ fric_grad_surf) if compute != "vector": @@ -472,6 +504,8 @@ def _compute_ipc_contributions(self, vertices, compute="all"): vertices, project_hessian_to_psd=ipctk.PSDProjectionMethod.CLAMP, ) + if axi_D is not None: + fric_hess_surf = axi_D @ fric_hess_surf @ axi_D self.global_matrix += P @ fric_hess_surf @ P.T # Pad with zeros to account for extra global DOFs (e.g. PeriodicBC) @@ -803,6 +837,15 @@ def initialize(self, pb): self._rest_positions.max(axis=0) - self._rest_positions.min(axis=0) ) + # Build axisymmetric (2π·r) integration weight on surface DOFs. + # Applied per surface vertex, broadcast over the ipctk interleaved + # layout [x0, y0, x1, y1, ...]. Reference radius R₀ is used, + # consistent with the rest of the 2Daxi pipeline. + if self.space.is_axisymmetric: + r_surf = self._rest_positions[:, 0] + self._axi_weight = np.repeat(2.0 * np.pi * r_surf, ndim) + self._axi_diag = sparse.diags(self._axi_weight) + # Compute actual dhat if self._dhat_is_relative: self._actual_dhat = self._dhat * self._bbox_diag diff --git a/fedoo/constraint/periodic_bc.py b/fedoo/constraint/periodic_bc.py index 318bb97a..e73460c7 100644 --- a/fedoo/constraint/periodic_bc.py +++ b/fedoo/constraint/periodic_bc.py @@ -1275,6 +1275,12 @@ def _add_additional_rot_dof(self, problem, res, get_boundary): ) def initialize(self, problem, dic_closest_points_on_boundaries=None): + if problem.space.is_axisymmetric: + raise NotImplementedError( + "PeriodicBC does not support the '2Daxi' ModelingSpace: " + "translational periodicity in the (r, z) plane is not " + "well-defined for an axisymmetric body of revolution." + ) # Use stored dictionary if none provided explicitly if dic_closest_points_on_boundaries is None: dic_closest_points_on_boundaries = self._dic_closest_points_on_boundaries diff --git a/fedoo/constraint/rigid_tie.py b/fedoo/constraint/rigid_tie.py index 8cccacb8..9e42dd3d 100644 --- a/fedoo/constraint/rigid_tie.py +++ b/fedoo/constraint/rigid_tie.py @@ -1,5 +1,7 @@ """Rigid Tie constraint.""" +import warnings + import numpy as np from fedoo.core.boundary_conditions import BCBase, MPC, ListBC @@ -98,6 +100,16 @@ def __repr__(self): return "\n".join(list_str) def initialize(self, problem): + if problem.space.is_axisymmetric: + warnings.warn( + "RigidTie under the '2Daxi' ModelingSpace: the constraint " + "is applied as pure-kinematics MPCs and will function, but " + "only rigid motions consistent with axisymmetry (axial " + "translation along the symmetry axis) preserve the " + "axisymmetric assumption. Radial translation and any " + "rotation around the in-plane axes break it.", + stacklevel=2, + ) if self.center is None: # initialize the rotation center at center of rigid nodes bounding box nodes_crd = problem.mesh.nodes[self.list_nodes] @@ -340,6 +352,14 @@ def __repr__(self): return "\n".join(list_str) def initialize(self, problem): + if problem.space.is_axisymmetric: + warnings.warn( + "RigidTie2D under the '2Daxi' ModelingSpace: the constraint " + "is applied as pure-kinematics MPCs and will function, but " + "only axial (z) translation preserves axisymmetry. Radial " + "translation and the in-plane RigidRotZ rotation break it.", + stacklevel=2, + ) if self.center is None: # initialize the rotation center at center of rigid nodes bounding box nodes_crd = problem.mesh.nodes[self.list_nodes] diff --git a/fedoo/core/mechanical3d.py b/fedoo/core/mechanical3d.py index 871958d8..711be31e 100644 --- a/fedoo/core/mechanical3d.py +++ b/fedoo/core/mechanical3d.py @@ -7,7 +7,101 @@ class Mechanical3D(ConstitutiveLaw): - """Base class for mechanical constitutive laws.""" + """Base class for mechanical constitutive laws. + + Strain / stress 6-vector slot ordering + -------------------------------------- + Strains and stresses are stored as 6-vectors (Voigt-like) at every Gauss + point. The slot interpretation depends on the active ``ModelingSpace``: + + ========= =================== =========================== =================== + slot 3D 2Dplane (plane strain) 2Daxi + ========= =================== =========================== =================== + 0 ε_xx ε_xx ε_rr + 1 ε_yy ε_yy ε_zz (z = Y axis) + 2 ε_zz 0 (unused) ε_θθ (= u_r / R) + 3 γ_xy γ_xy γ_rz + 4 γ_xz 0 0 (γ_rθ = 0 by sym.) + 5 γ_yz 0 0 (γ_θz = 0 by sym.) + ========= =================== =========================== =================== + + For ``2Dstress`` (plane stress), the slot layout is the same as + ``2Dplane``, but the constitutive law internally computes a nonzero + out-of-plane strain ε_zz so that σ_zz = 0 — the symbolic strain + operator places 0 in slot 2, but post-processed strain output from + a constitutive law that relaxes σ_zz may have a nonzero ε_zz that + is not reflected in slot 2. + + Rationale for the 2Daxi mapping + ------------------------------- + In axisymmetric kinematics, γ_rθ = γ_θz = 0 by symmetry, so γ_rz is the + only nonzero shear — exactly as γ_xy is the only nonzero shear in 2D + plane strain / plane stress. Placing γ_rz in slot 3 keeps slot 3 = "the + in-plane shear" across all 2D regimes, so element B-matrices, the + assembly pipeline, and most constitutive-law code are dimension-agnostic + in their handling of slot 3. The 3D-Voigt slot 2 (normally ε_zz in + Cartesian) is repurposed to carry ε_θθ; the physical ε_zz of 2Daxi + lives in slot 1 because z is along the Y axis (z ≡ Y, r ≡ X). + + Convention + ---------- + For a 2Daxi mesh, ``mesh.nodes[:, 0]`` is the radial coordinate r and + ``mesh.nodes[:, 1]`` is the axial coordinate z. The symmetry axis is the + Y axis in 2D; after ``axisymmetric_extrusion`` the symmetry axis becomes + Z in the resulting 3D mesh. + + Implication for orthotropic / anisotropic constitutive laws + ----------------------------------------------------------- + The 6×6 tangent matrix produced by a constitutive law is consumed by + fedoo with the slot mapping above. For a law whose material directions + are labelled (1, 2, 3): in 2Daxi, axis 1 is r, axis 2 is z, **axis 3 is + the hoop direction**. Authors of orthotropic / anisotropic laws should + declare moduli with this convention in mind; otherwise the user-supplied + "direction 3" stiffness will be silently applied to the hoop response. + + Theory of axisymmetric kinematics + --------------------------------- + Label the reference configuration in cylindrical coordinates + ``(R, Θ, Z)`` and the current configuration ``(r, θ, z)``. An + axisymmetric motion *without twist* satisfies:: + + r = r(R, Z), z = z(R, Z), θ = Θ. + + In the orthonormal cylindrical bases ``(e_R, e_Θ, e_Z)`` and + ``(e_r, e_θ, e_z)`` the deformation gradient takes the block form + + .. math:: + + F = \\frac{\\partial r}{\\partial R}\\, e_r \\!\\otimes\\! e_R + + \\frac{\\partial r}{\\partial Z}\\, e_r \\!\\otimes\\! e_Z + + \\frac{\\partial z}{\\partial R}\\, e_z \\!\\otimes\\! e_R + + \\frac{\\partial z}{\\partial Z}\\, e_z \\!\\otimes\\! e_Z + + \\frac{r}{R}\\, e_θ \\!\\otimes\\! e_Θ. + + Two consequences are load-bearing for the implementation: + + * **Block structure.** All hoop-coupling components vanish: + ``F_rΘ = F_zΘ = F_θR = F_θZ = 0``. Therefore γ_rθ = γ_θz = 0 in + every kinematic measure, which is why slot 3 (γ_rz) is the only + shear component in 2Daxi (slots 4 and 5 are identically zero). + + * **Hoop deformation gradient.** The hoop stretch is + ``λ_θ = F_θΘ = r/R = 1 + u_r/R`` (where ``u_r = r − R`` is the + *total* radial displacement and ``R`` is the *reference* radial + coordinate). At small strain ``r ≈ R`` and this reduces to the + conventional ``ε_θθ = u_r/r``, which is what fedoo's + ``_comp_linear_strain`` and the small-strain weakforms compute via + the symbolic operator ``space.variable("DispX") * (1/r_gp)``. At + finite strain the two definitions differ; fedoo's UL + 2Daxi + pipeline therefore divides ``u_r`` by the *reference* radius + ``R = mesh.nodes[:, 0]`` (captured at problem ``initialize`` as + ``assembly.sv["_R0_gausspoints"]``), not by ``r_current``. + + See Bonet & Wood, *Nonlinear Continuum Mechanics for Finite Element + Analysis* (2008), Box 8.3; Holzapfel, *Nonlinear Solid Mechanics* + (2000), §2.5; Belytschko, Liu & Moran, *Nonlinear Finite Elements + for Continua and Structures* (2014), §4.5. + """ # model of constitutive law for InternalForce Weakform diff --git a/fedoo/core/modelingspace.py b/fedoo/core/modelingspace.py index 9f040375..ab9a33d7 100644 --- a/fedoo/core/modelingspace.py +++ b/fedoo/core/modelingspace.py @@ -84,6 +84,11 @@ def get_dimension(self): """ return self._dimension + @property + def is_axisymmetric(self): + """True if this ModelingSpace uses the 2Daxi axisymmetric assumption.""" + return self._dimension == "2Daxi" + def make_active(self): """Define the modeling space as the active ModelingSpace.""" ModelingSpace._active = self @@ -287,6 +292,22 @@ def op_div_u(self): return self.derivative("DispX", "X") + self.derivative("DispY", "Y") def op_strain(self, InitialGradDisp=None): + """Symbolic strain operator as a list of six DiffOps. + + The returned 6-vector follows the slot ordering documented on + :class:`fedoo.core.mechanical3d.Mechanical3D`. In 2D regimes + slot 2 is set to 0 here; for ``2Daxi`` problems the hoop term + ``eps[2] = u_r / r`` is patched in by each mechanical weakform's + ``get_weak_equation`` (see e.g. ``StressEquilibrium``) using + ``assembly.sv["_R_gausspoints"]``. In 2D plane strain / plane + stress slot 2 stays 0. + + Note for plane stress (``2Dstress``): the symbolic slot 2 is 0 + (kinematic zero), but the constitutive law may compute a + nonzero out-of-plane ε_zz when enforcing the plane-stress + condition σ_zz = 0. Post-processed strain output from such a + law is therefore not guaranteed to have slot 2 = 0. + """ # InitialGradDisp = StrainOperator.__InitialGradDisp if (InitialGradDisp is None) or ( diff --git a/fedoo/post_processing/axi_to_3d.py b/fedoo/post_processing/axi_to_3d.py index e612a345..9dd6a844 100644 --- a/fedoo/post_processing/axi_to_3d.py +++ b/fedoo/post_processing/axi_to_3d.py @@ -5,6 +5,108 @@ from zipfile import Path +def _revolve_axi_field(field, theta, n_nodes_2d, kind=None, field_name=None): + """Revolve a 2D axisymmetric field around theta to a 3D field. + + The 2D (r, z) plane is interpreted using fedoo's 2Daxi conventions: + r = column 0, z = column 1, hoop component in Voigt slot 2 of any + 6-vector. The output 3D Cartesian frame uses (X, Y, Z) with the + symmetry axis along Z, matching ``axisymmetric_extrusion``. + + Parameters + ---------- + field : ndarray, shape (n_components, n_nodes_2d) or (n_nodes_2d,) + Field on the 2D mesh. + theta : ndarray, shape (n_theta,) + Angular sample positions. + n_nodes_2d : int + Number of nodes in the 2D mesh. + kind : {"stress", "strain"}, optional + For 6-tensor fields, declares whether the off-diagonal Voigt slot + stores a tensor component (``"stress"`` -> ``sigma_ij``) or an + engineering shear (``"strain"`` -> ``gamma_ij = 2 * eps_ij``). A + strain field needs an extra factor of 2 in the rotated slot 3 + (``gamma_xy``) because the hoop / radial diagonals are *not* + engineering quantities. ``None`` falls back to the + ``field_name`` heuristic below. + field_name : str, optional + Used only when ``kind`` is None: any name containing ``"Strain"`` + is treated as ``kind="strain"``, otherwise ``kind="stress"``. + + Returns + ------- + ndarray, shape (n_components_out, n_nodes_2d * n_theta) + Revolved field. ``n_components_out`` is 3 for a 2-vector input, + 6 for a 6-tensor input, and equal to ``n_components`` otherwise + (scalars and unrecognised shapes are tiled along theta). + """ + field = np.asarray(field) + n_theta = theta.shape[0] + cos_t = np.cos(theta) + sin_t = np.sin(theta) + + # Promote a 1-D scalar field to 2-D for uniform handling. + if field.ndim == 1: + field = field.reshape(1, -1) + + n_comp = field.shape[0] + + # Convention: the extruded 3D mesh has nodes laid out node-major, + # i.e. the 2D node index varies slowest and theta varies fastest. + # All per-node arrays we build below have shape (n_nodes_2d, n_theta) + # so that ``.ravel()`` produces [n0t0, n0t1, ..., n1t0, n1t1, ...]. + cos_row = cos_t[None, :] # (1, n_theta) + sin_row = sin_t[None, :] + + if n_comp == 2: + # Cylindrical (dr, dz) vector -> Cartesian (X, Y, Z) + # at each angle theta_k: (dr*cos, dr*sin, dz). + dr = field[0][:, None] # (n_nodes_2d, 1) + dz = field[1][:, None] + out = np.empty((3, n_nodes_2d * n_theta), dtype=field.dtype) + out[0] = (dr * cos_row).ravel() + out[1] = (dr * sin_row).ravel() + out[2] = np.broadcast_to(dz, (n_nodes_2d, n_theta)).ravel() + return out + + if n_comp == 6: + # Symmetric second-order tensor in fedoo's 2Daxi slot ordering + # (s_rr, s_zz, s_tt, s_rz, 0, 0) -> Cartesian Voigt (xx, yy, zz, xy, xz, yz). + # Rotation about the symmetry axis Z by theta gives, for the underlying + # tensor: s_xy_tensor = sin*cos (s_rr - s_tt). For strain fields whose + # slot 3 stores an engineering shear gamma = 2*eps, the output slot 3 + # gamma_xy equals 2 * sin*cos * (s_rr - s_tt) = sin(2 theta)*(...). + # Slots 4 and 5 (rz -> xz/yz) carry the same convention as slot 3 of the + # input on both sides, so no extra factor is needed there. + if kind is None: + kind = "strain" if (field_name and "Strain" in field_name) else "stress" + if kind not in ("stress", "strain"): + raise ValueError(f"kind must be 'stress' or 'strain', got {kind!r}.") + diag_to_xy = 2.0 if kind == "strain" else 1.0 + + s_rr = field[0][:, None] # (n_nodes_2d, 1) + s_zz = field[1][:, None] + s_tt = field[2][:, None] + s_rz = field[3][:, None] + c2 = (cos_t**2)[None, :] + s2 = (sin_t**2)[None, :] + sc = (sin_t * cos_t)[None, :] + out = np.empty((6, n_nodes_2d * n_theta), dtype=field.dtype) + out[0] = (c2 * s_rr + s2 * s_tt).ravel() + out[1] = (s2 * s_rr + c2 * s_tt).ravel() + out[2] = np.broadcast_to(s_zz, (n_nodes_2d, n_theta)).ravel() + out[3] = (diag_to_xy * sc * (s_rr - s_tt)).ravel() + out[4] = (cos_row * s_rz).ravel() + out[5] = (sin_row * s_rz).ravel() + return out + + # Scalar fields and any unrecognised shape: tile along theta unchanged + # (node-major ordering: each 2D node value repeated n_theta times). + return (field.reshape(n_comp, n_nodes_2d, 1) * np.ones((1, 1, n_theta))).reshape( + n_comp, n_nodes_2d * n_theta + ) + + def axi_to_3d(axi_data: DataSet | MultiFrameDataSet, n_theta: int = 41, filename=None): """Convert axisymmetric data into a full 3D representation. @@ -77,7 +179,6 @@ def axi_to_3d_dataset(axi_data: DataSet, n_theta: int = 41): mesh = axi_data.mesh full_mesh = axisymmetric_extrusion(mesh, n_theta, merge_nodes=False) theta = np.linspace(0, 2 * np.pi, n_theta) - e_r = np.tile(np.column_stack((np.cos(theta), np.sin(theta))), (mesh.n_nodes, 1)) res3d = DataSet(full_mesh) for field in [ @@ -86,17 +187,12 @@ def axi_to_3d_dataset(axi_data: DataSet, n_theta: int = 41): *axi_data.element_data, ]: if field == "Disp": - # dr,dz = np.tile(res.node_data[field],n_theta) - dr, dz = ( - axi_data.node_data[field].reshape(2, mesh.n_nodes, 1) - * np.ones((1, 1, n_theta)) - ).reshape(2, -1) - res3d.node_data[field] = np.column_stack((dr[:, np.newaxis] * e_r, dz)).T + data2d = axi_data.node_data[field] else: - res3d.node_data[field] = ( - axi_data.get_data(field, data_type="Node").reshape(-1, mesh.n_nodes, 1) - * np.ones((1, 1, n_theta)) - ).reshape(-1, mesh.n_nodes * n_theta) + data2d = axi_data.get_data(field, data_type="Node") + res3d.node_data[field] = _revolve_axi_field( + data2d, theta, mesh.n_nodes, field_name=field + ) res3d.scalar_data = axi_data.scalar_data return res3d @@ -164,10 +260,9 @@ def __init__(self, axi_data: MultiFrameDataSet, n_theta: int = 41): self.n_theta = n_theta mesh = axisymmetric_extrusion(self.mesh2d, n_theta, merge_nodes=False) - theta = np.linspace(0, 2 * np.pi, n_theta) - self._er = np.tile( - np.column_stack((np.cos(theta), np.sin(theta))), (self.mesh2d.n_nodes, 1) - ) + # Precompute the angular sample positions once; reused on every + # frame load by _revolve_axi_field. + self._theta = np.linspace(0, 2 * np.pi, n_theta) DataSet.__init__(self, mesh) if self.loaded_iter is None: @@ -200,22 +295,17 @@ def field_names(self): def _import_disp_to_3d(self): if "Disp" not in self.node_data and "Disp" in self.axi_data.node_data: - dr, dz = ( - self.axi_data.node_data["Disp"].reshape(2, self.mesh2d.n_nodes, 1) - * np.ones((1, 1, self.n_theta)) - ).reshape(2, -1) - self.node_data["Disp"] = np.column_stack( - (dr[:, np.newaxis] * self._er, dz) - ).T + self.node_data["Disp"] = _revolve_axi_field( + self.axi_data.node_data["Disp"], self._theta, self.mesh2d.n_nodes + ) def get_data(self, field, component=None, data_type=None, return_data_type=False): if field not in self.node_data: # import field as 3d node data if not already present data2d = self.axi_data.get_data(field, data_type="Node") - self.node_data[field] = ( - data2d.reshape(-1, self.mesh2d.n_nodes, 1) - * np.ones((1, 1, self.n_theta)) - ).reshape(-1, self.mesh2d.n_nodes * self.n_theta) + self.node_data[field] = _revolve_axi_field( + data2d, self._theta, self.mesh2d.n_nodes, field_name=field + ) return DataSet.get_data(self, field, component, data_type, return_data_type) diff --git a/fedoo/weakform/_axi_utils.py b/fedoo/weakform/_axi_utils.py new file mode 100644 index 00000000..1483f5b6 --- /dev/null +++ b/fedoo/weakform/_axi_utils.py @@ -0,0 +1,41 @@ +"""Internal helpers for 2Daxi weakforms. + +Centralises the ``2*pi*r`` axisymmetric volume / surface integration +weight so that every mechanical, inertial, damping, contact and +interface weakform applies the factor consistently. +""" + +import numpy as np + + +def axi_volume_weight(assembly): + """Return ``2*pi*r`` evaluated at each Gauss point. + + Reuses ``assembly.sv["_R_gausspoints"]`` when it has already been + populated by :class:`fedoo.weakform.StressEquilibrium` (which keeps + it refreshed at the current configuration in UL). Otherwise builds + the radial coordinate from the current mesh on the fly so that + standalone uses (e.g. an ``Inertia`` weakform without a companion + ``StressEquilibrium``) still get the correct weight. + + Parameters + ---------- + assembly : fedoo.Assembly + The assembly the weakform integrates over. Must be associated + with a 2Daxi ``ModelingSpace``. + + Returns + ------- + ndarray, shape (n_gauss_points,) + ``2*pi*r`` at every Gauss point of the current mesh. + """ + rr = assembly.sv.get("_R_gausspoints") + if rr is None: + mesh = assembly.current.mesh + rr = mesh.convert_data( + mesh.nodes[:, 0], + "Node", + "GaussPoint", + n_elm_gp=assembly.n_elm_gp, + ) + return (2 * np.pi) * rr diff --git a/fedoo/weakform/damping_stabilization.py b/fedoo/weakform/damping_stabilization.py index c1889f3b..2cbc1ece 100644 --- a/fedoo/weakform/damping_stabilization.py +++ b/fedoo/weakform/damping_stabilization.py @@ -2,6 +2,7 @@ import numpy as np from fedoo.core.weakform import WeakFormBase +from fedoo.weakform._axi_utils import axi_volume_weight class ArtificialDamping(WeakFormBase): @@ -181,15 +182,15 @@ def get_weak_equation(self, assembly, pb): [a * b * (self._c_stab / dt) for (a, b) in zip(op_var_vir, op_var)] ) + # 2*pi*r weight applied on the tangent so both the tangent + # contribution and the residual derived from it via operator_apply + # inherit it consistently. + if self.space.is_axisymmetric: + tangent_matrix = tangent_matrix * axi_volume_weight(assembly) + # 4. Residual Contribution: c_stab * M* * v_pseudo if not np.array_equal(delta_u, 0): damping_force = assembly.operator_apply(tangent_matrix, delta_u) - - # Axisymmetric correction - # if self.space is not None and getattr(self.space, "_dimension", "") == "2Daxi": - # rr = assembly.sv["_R_gausspoints"] - # damping_force = damping_force * ((2 * np.pi) * rr) - return tangent_matrix + damping_force return tangent_matrix diff --git a/fedoo/weakform/distributed_load.py b/fedoo/weakform/distributed_load.py index 497eb82e..7940cd51 100644 --- a/fedoo/weakform/distributed_load.py +++ b/fedoo/weakform/distributed_load.py @@ -58,7 +58,7 @@ def get_weak_equation(self, assembly: AssemblyBase, pb: ProblemBase) -> DiffOp: vec_u = self.space.op_disp() normals = assembly.sv["Normals"] - if self.space._dimension == "2Daxi": + if self.space.is_axisymmetric: mesh = assembly.current.mesh rr = mesh.convert_data( mesh.nodes[:, 0], @@ -137,7 +137,7 @@ def get_weak_equation(self, assembly: AssemblyBase, pb: ProblemBase) -> DiffOp: """Return the weak equation related to the current state.""" vec_u = self.space.op_disp() - if self.space._dimension == "2Daxi": + if self.space.is_axisymmetric: mesh = assembly.current.mesh rr = mesh.convert_data( mesh.nodes[:, 0], diff --git a/fedoo/weakform/heat_equation.py b/fedoo/weakform/heat_equation.py index 8a349b27..40caf08b 100644 --- a/fedoo/weakform/heat_equation.py +++ b/fedoo/weakform/heat_equation.py @@ -24,7 +24,7 @@ def __init__(self, thermal_constitutivelaw, name=None, nlgeom=False, space=None) name = thermal_constitutivelaw.name WeakFormBase.__init__(self, name, space) - if self.space._dimension == "2Daxi": + if self.space.is_axisymmetric: raise NotImplementedError self.space.new_variable("Temp") # temperature diff --git a/fedoo/weakform/implicit_dynamic.py b/fedoo/weakform/implicit_dynamic.py index 4049daef..f3c7d2a3 100644 --- a/fedoo/weakform/implicit_dynamic.py +++ b/fedoo/weakform/implicit_dynamic.py @@ -81,7 +81,7 @@ def get_weak_equation(self, assembly, pb): # Apply the operator to the nodal v_curr to get the damping force residual if not np.array_equal(residual_val, 0): diff_op += assembly.operator_apply(wf, residual_val.ravel()) - if self.space._dimension == "2Daxi": + if self.space.is_axisymmetric: rr = assembly.sv["_R_gausspoints"] return diff_op * ((2 * np.pi) * rr) else: @@ -136,7 +136,7 @@ def get_weak_equation(self, assembly, pb): damping_force_wf = self.damping_coef * assembly.operator_apply( mat, v_curr.ravel() ) - if self.space._dimension == "2Daxi": + if self.space.is_axisymmetric: rr = assembly.sv["_R_gausspoints"] damping_force_wf = damping_force_wf * ((2 * np.pi) * rr) return scaled_mat + vec + damping_force_wf @@ -268,6 +268,14 @@ def __init__( ): super().__init__(name, space) + if self.space.is_axisymmetric: + raise NotImplementedError( + "ImplicitDynamic2 is not implemented in '2Daxi'. " + "Use the ImplicitDynamic factory (which builds an " + "ImplicitDynamicSum from _NewmarkStiffness + _NewmarkInertia) " + "instead — that path correctly applies the 2*pi*r weight." + ) + if name != "": stiffness_name = name + "_stiffness" else: @@ -400,7 +408,7 @@ def get_weak_equation(self, assembly, pb): initial_stress = assembly.sv["PK2"] else: eps = self.space.op_strain() - if self.space._dimension == "2Daxi": + if self.space.is_axisymmetric: raise NotImplementedError("2Daxi not implemented.") # rr = assembly.sv["_R_gausspoints"] # eps[2] = self.space.variable("DispX") * np.divide( @@ -502,7 +510,7 @@ def get_weak_equation(self, assembly, pb): for i in range(6) ] ) - # if self.space._dimension == "2Daxi": + # if self.space.is_axisymmetric: # diff_op = diff_op * ((2 * np.pi) * rr) return diff_op diff --git a/fedoo/weakform/inertia.py b/fedoo/weakform/inertia.py index e1e7ddfb..2e76415a 100644 --- a/fedoo/weakform/inertia.py +++ b/fedoo/weakform/inertia.py @@ -1,4 +1,5 @@ from fedoo.core.weakform import WeakFormBase +from fedoo.weakform._axi_utils import axi_volume_weight class Inertia(WeakFormBase): @@ -36,7 +37,12 @@ def get_weak_equation(self, assembly, pb): op_dU = self.space.op_disp() # displacement increment (incremental formulation) op_dU_vir = [du.virtual if du != 0 else 0 for du in op_dU] - return sum([a * b * self.density for (a, b) in zip(op_dU_vir, op_dU)]) + diff_op = sum([a * b * self.density for (a, b) in zip(op_dU_vir, op_dU)]) + + if self.space.is_axisymmetric: + diff_op = diff_op * axi_volume_weight(assembly) + + return diff_op class RotaryInertia(WeakFormBase): @@ -62,6 +68,12 @@ class RotaryInertia(WeakFormBase): def __init__(self, rotary_inertia, name="", space=None): WeakFormBase.__init__(self, name, space) + if self.space.is_axisymmetric: + raise NotImplementedError( + "RotaryInertia is not defined in '2Daxi' (no rotational " + "degree of freedom in axisymmetric kinematics)." + ) + if self.space.ndim == 3: self.space.new_variable("RotX") # torsion rotation self.space.new_variable("RotY") diff --git a/fedoo/weakform/interface_force.py b/fedoo/weakform/interface_force.py index 80f90283..8cec444f 100644 --- a/fedoo/weakform/interface_force.py +++ b/fedoo/weakform/interface_force.py @@ -1,5 +1,6 @@ from fedoo.core.weakform import WeakFormBase from fedoo.core.base import ConstitutiveLaw +from fedoo.weakform._axi_utils import axi_volume_weight import numpy as np @@ -122,4 +123,9 @@ def get_weak_equation(self, assembly, pb): ] ) + # An interface in the (r, z) plane is a cylindrical surface with + # area dA = 2*pi*r * ds. + if self.space.is_axisymmetric: + diff_op = diff_op * axi_volume_weight(assembly) + return diff_op diff --git a/fedoo/weakform/stress_equilibrium.py b/fedoo/weakform/stress_equilibrium.py index 0a96f687..e69d1219 100644 --- a/fedoo/weakform/stress_equilibrium.py +++ b/fedoo/weakform/stress_equilibrium.py @@ -94,7 +94,7 @@ def get_weak_equation(self, assembly, pb): "Stress" ] # Stress = Cauchy for updated lagrangian method - if self.space._dimension == "2Daxi": + if self.space.is_axisymmetric: rr = assembly.sv["_R_gausspoints"] # nlgeom = False @@ -133,7 +133,7 @@ def get_weak_equation(self, assembly, pb): ] ) - if self.space._dimension == "2Daxi": + if self.space.is_axisymmetric: DiffOp = DiffOp * ((2 * np.pi) * rr) return DiffOp @@ -159,18 +159,44 @@ def initialize(self, assembly, pb): ) assembly.sv["DispGradient"] = 0 - if self.space._dimension == "2Daxi": - assembly.sv["_R_gausspoints"] = assembly.mesh.convert_data( - assembly.mesh.nodes[:, 0], + if self.space.is_axisymmetric: + # ``assembly.mesh`` is treated as the *reference* configuration: + # any subsequent ``set_disp`` only mutates ``assembly.current.mesh`` + # (see fedoo/core/assembly.py: set_disp). Capture R0 from the + # reference mesh here, once. If the mesh has already been deformed + # before initialize is reached (unusual, e.g. chained problems + # sharing an assembly), the captured R0 will be the deformed + # radius, breaking F_theta-theta = r/R at finite strain. Callers + # should rebuild the assembly before initializing a 2Daxi problem + # in that case. + r_nodes = assembly.mesh.nodes[:, 0] + if r_nodes.min() < 0: + raise ValueError( + "2Daxi requires non-negative radial coordinates " + "(mesh.nodes[:, 0] >= 0). Found " + f"min(r) = {r_nodes.min():.6g}. " + "The convention is r = X (column 0), z = Y (column 1)." + ) + # Reference radial coordinate at gauss points (initial mesh). + # Captured once and never overwritten: used to form the canonical + # hoop deformation gradient F_theta-theta = r_current / R_reference + # in finite-strain UL+axi (Bonet & Wood, Box 8.3). + assembly.sv["_R0_gausspoints"] = assembly.mesh.convert_data( + r_nodes, "Node", "GaussPoint", n_elm_gp=assembly.n_elm_gp, ) + # Current radial coordinate at gauss points. Equal to the reference + # at initialize; refreshed each iteration in _comp_grad_disp to the + # deformed mesh (used for the 2*pi*r weak-form weight and the + # symbolic operator eps[2] = DispX / r at the current config). + assembly.sv["_R_gausspoints"] = assembly.sv["_R0_gausspoints"].copy() if assembly._nlgeom: if assembly._nlgeom == "TL": assembly.sv["PK2"] = 0 - if self.space._dimension == "2Daxi": + if self.space.is_axisymmetric: raise NotImplementedError( "'2Daxi' ModelingSpace is not implemented with \ total lagrangian formulation. Use update \ @@ -470,29 +496,23 @@ def corate(self, value): # function to compute the displacement gradient def _comp_grad_disp(assembly, displacement): grad_values = assembly.get_grad_disp(displacement, "GaussPoint") - if assembly.space._dimension == "2Daxi": - # mesh = assembly.current.mesh - # eps_tt = np.divide( - # pb.get_disp()[0], - # mesh.nodes[:, 0], - # out=np.zeros(mesh.n_nodes), - # where=mesh.nodes[:, 0] != 0, - # ) # put zero if X==0 (division by 0) - # grad_values[2][2] = mesh.convert_data( - # eps_tt, "Node", "GaussPoint" - # ) + if assembly.space.is_axisymmetric: mesh = assembly.current.mesh - rr = mesh.convert_data( + # Refresh r_current at gauss points: used by the symbolic + # operator eps[2] = DispX / r in get_weak_equation and by the + # 2*pi*r weak-form integration weight, both at the current config. + assembly.sv["_R_gausspoints"] = mesh.convert_data( mesh.nodes[:, 0], "Node", "GaussPoint", n_elm_gp=assembly.n_elm_gp, ) - assembly.sv["_R_gausspoints"] = rr - # grad_values[2][2] = mesh.convert_data( - # pb.get_disp()[0]/mesh.nodes[:, 0], 'Node') - + # F_theta-theta = r_current / R_reference, hence + # grad_values[2][2] = u_r / R_reference. See the + # "Theory of axisymmetric kinematics" section of + # :class:`fedoo.core.mechanical3d.Mechanical3D` for the derivation. + R0 = assembly.sv["_R0_gausspoints"] rank_dispx = assembly.space.variable_rank("DispX") n = mesh.n_nodes grad_values[2][2] = np.divide( @@ -502,10 +522,10 @@ def _comp_grad_disp(assembly, displacement): "GaussPoint", n_elm_gp=assembly.n_elm_gp, ), - rr, - out=np.zeros_like(rr), - where=rr != 0, - ) # put zero if X==0 (division by 0) + R0, + out=np.zeros_like(R0), + where=R0 != 0, + ) # zero at the symmetry axis (R0 == 0) assembly.sv["DispGradient"] = grad_values return grad_values diff --git a/fedoo/weakform/stress_equilibrium_bbar.py b/fedoo/weakform/stress_equilibrium_bbar.py index cb051c81..8233054a 100644 --- a/fedoo/weakform/stress_equilibrium_bbar.py +++ b/fedoo/weakform/stress_equilibrium_bbar.py @@ -49,6 +49,12 @@ class StressEquilibriumBbar(StressEquilibrium): def get_weak_equation(self, assembly, pb): """Get the weak equation related to the current problem state.""" if assembly._nlgeom == "TL": # add initial displacement effect + if self.space.is_axisymmetric: + raise NotImplementedError( + "'2Daxi' ModelingSpace is not implemented with the total " + "lagrangian formulation for the b-bar / sub-integration " + "weak form. Use updated lagrangian instead." + ) eps = self.space.op_strain(assembly.sv["DispGradient"]) initial_stress = assembly.sv["PK2"] else: @@ -75,7 +81,7 @@ def get_weak_equation(self, assembly, pb): "Stress" ] # Stress = Cauchy for updated lagrangian method - if self.space._dimension == "2Daxi": + if self.space.is_axisymmetric: rr = assembly.sv["_R_gausspoints"] # nlgeom = False @@ -110,7 +116,7 @@ def get_weak_equation(self, assembly, pb): ] ) - if self.space._dimension == "2Daxi": + if self.space.is_axisymmetric: DiffOp = DiffOp * ((2 * np.pi) * rr) return DiffOp @@ -174,6 +180,14 @@ class StressEquilibriumFbar(StressEquilibrium): def __init__(self, constitutivelaw, name="", nlgeom=None, space=None): super().__init__(constitutivelaw, name, nlgeom, space) + if self.space.is_axisymmetric: + raise NotImplementedError( + "StressEquilibriumFbar is not implemented in '2Daxi': the " + "F-bar volumetric correction requires the hoop component " + "of the strain at the element centroid, which is not yet " + "wired through _op_strain_center. Use StressEquilibrium " + "(standard) or StressEquilibriumMixed in 2Daxi instead." + ) self.fbar = True self.assembly_options["elm_type", "quad4"] = "quad4sri" self.assembly_options["elm_type", "hex8"] = "hex8sri" @@ -181,6 +195,12 @@ def __init__(self, constitutivelaw, name="", nlgeom=None, space=None): def get_weak_equation(self, assembly, pb): """Get the weak equation related to the current problem state.""" if assembly._nlgeom == "TL": # add initial displacement effect + if self.space.is_axisymmetric: + raise NotImplementedError( + "'2Daxi' ModelingSpace is not implemented with the total " + "lagrangian formulation for the f-bar weak form. " + "Use updated lagrangian instead." + ) eps = self.space.op_strain(assembly.sv["DispGradient"]) initial_stress = assembly.sv["PK2"] else: @@ -189,7 +209,7 @@ def get_weak_equation(self, assembly, pb): "Stress" ] # Stress = Cauchy for updated lagrangian method - if self.space._dimension == "2Daxi": + if self.space.is_axisymmetric: rr = assembly.sv["_R_gausspoints"] # nlgeom = False @@ -259,7 +279,7 @@ def get_weak_equation(self, assembly, pb): ] ) - if self.space._dimension == "2Daxi": + if self.space.is_axisymmetric: DiffOp = DiffOp * ((2 * np.pi) * rr) return DiffOp @@ -364,6 +384,13 @@ class HourglassStiffness(WeakFormBase): def __init__(self, stiffness_coef=0.01, name="", nlgeom=False, space=None): WeakFormBase.__init__(self, name, space) + if self.space.is_axisymmetric: + raise NotImplementedError( + "HourglassStiffness is not implemented in '2Daxi': the " + "2*pi*r integration weight on the stabilization term has " + "not been validated. Use full integration (no hourglass " + "stabilization needed) in 2Daxi." + ) self.assembly_options["n_elm_gp"] = 1 self.assembly_options["elm_type", "quad4"] = "quad4hourglass" self.assembly_options["elm_type", "hex8"] = "hex8hourglass" @@ -458,7 +485,7 @@ def get_weak_equation(self, assembly, pb): DiffOp = DiffOp * hourglass_stiffness - # if self.space._dimension == "2Daxi": + # if self.space.is_axisymmetric: # # not sur it works # DiffOp = DiffOp * ((2 * np.pi) * rr) diff --git a/fedoo/weakform/stress_equilibrium_mixed.py b/fedoo/weakform/stress_equilibrium_mixed.py index 8b0e860c..f878dbc3 100644 --- a/fedoo/weakform/stress_equilibrium_mixed.py +++ b/fedoo/weakform/stress_equilibrium_mixed.py @@ -59,6 +59,12 @@ def get_weak_equation(self, assembly, pb): # 1. Kinematics (Strain / Deformation Gradient) # ----------------------------------------------------------- if assembly._nlgeom == "TL": # add initial displacement effect + if self.space.is_axisymmetric: + raise NotImplementedError( + "'2Daxi' ModelingSpace is not implemented with the total " + "lagrangian formulation for the mixed (u, p) weak form. " + "Use updated lagrangian instead." + ) eps = self.space.op_strain(assembly.sv["DispGradient"]) initial_stress = assembly.sv["PK2"] else: @@ -66,7 +72,7 @@ def get_weak_equation(self, assembly, pb): # Stress = Cauchy for updated lagrangian method initial_stress = assembly.sv["Stress"] - if self.space._dimension == "2Daxi": + if self.space.is_axisymmetric: rr = assembly.sv["_R_gausspoints"] eps[2] = self.space.variable("DispX") * np.divide( 1, rr, out=np.zeros_like(rr), where=rr != 0 @@ -222,7 +228,7 @@ def get_weak_equation(self, assembly, pb): DiffOp += P_inc.virtual * (pressure_residual + pressure_stiffness_term) # Axisymmetric volume integration factor - if self.space._dimension == "2Daxi": + if self.space.is_axisymmetric: DiffOp = DiffOp * ((2 * np.pi) * rr) return DiffOp diff --git a/tests/test_axi_phase3.py b/tests/test_axi_phase3.py new file mode 100644 index 00000000..4da2070a --- /dev/null +++ b/tests/test_axi_phase3.py @@ -0,0 +1,78 @@ +"""Regression tests for Phase 3 fixes (Inertia / damping_stabilization / +InterfaceForce 2*pi*r weighting in 2Daxi). + +These tests don't run a full FEM solve. They build a tiny 2Daxi assembly +with a non-uniform radial coordinate, evaluate the symbolic weak form +on a known displacement, and check that the assembled mass / damping / +interface-force scales with the radial coordinate as the canonical +2*pi*r dV factor would predict. + +If a weakform skips the 2*pi*r weight, its global matrix scales like +the uniform mesh integral and the test fails. +""" + +import numpy as np +import pytest + +import fedoo as fd + + +def _build_2daxi_unit_square_mesh(n=2): + """Build a 1x1 quad mesh in (r, z) with r in [1, 2].""" + fd.ModelingSpace("2Daxi") + mesh = fd.mesh.rectangle_mesh(n + 1, n + 1, 1.0, 2.0, 0.0, 1.0) + return mesh + + +def test_inertia_axi_mass_matrix_scales_with_2pi_r(): + """Total mass of a unit-density 2Daxi annulus matches 2*pi * integral(r dA).""" + mesh = _build_2daxi_unit_square_mesh(n=4) + density = 1.0 + + inertia = fd.weakform.Inertia(density) + asm = fd.Assembly.create(inertia, mesh) + # Trigger initialize via a tiny fake problem-style call by assembling + # the global matrix directly. + asm.assemble_global_mat() + M = asm.get_global_matrix() + + # Sum of all entries of the consistent mass matrix gives total mass + # (because mass matrix for a unit displacement field is total mass per DOF + # block; sum of all entries of the diagonal blocks is total volume*density). + # Equivalent test: sum(M) per displacement component. + n_dof = M.shape[0] + n_nodes = mesh.n_nodes + # extract block for DispX (rank 0): rows/cols [0:n_nodes] + M_xx = M[:n_nodes, :n_nodes] + total_mass = float(M_xx.sum()) + + # analytic total volume of an annulus r in [1, 2], z in [0, 1] revolved 360 deg + # V = 2*pi * integral(r dr dz) = 2*pi * (r^2/2 from 1 to 2) * (1) = 2*pi * 1.5 = 3*pi + expected_mass = 2 * np.pi * 0.5 * (2.0**2 - 1.0**2) * 1.0 * density + assert total_mass == pytest.approx(expected_mass, rel=5e-3) + + +def test_inertia_3d_unaffected_by_axi_path(): + """Sanity: 3D Inertia still produces volume-integrated mass without 2*pi*r.""" + fd.ModelingSpace("3D") + mesh = fd.mesh.box_mesh(3, 3, 3, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0) + + density = 2.0 + inertia = fd.weakform.Inertia(density) + asm = fd.Assembly.create(inertia, mesh) + asm.assemble_global_mat() + M = asm.get_global_matrix() + + n_nodes = mesh.n_nodes + M_xx = M[:n_nodes, :n_nodes] + total_mass = float(M_xx.sum()) + # unit cube volume * density + expected = 1.0 * density + assert total_mass == pytest.approx(expected, rel=1e-3) + + +def test_rotary_inertia_axi_raises(): + """RotaryInertia in 2Daxi is not physically defined: must raise.""" + fd.ModelingSpace("2Daxi") + with pytest.raises(NotImplementedError, match="2Daxi"): + fd.weakform.RotaryInertia(1.0) diff --git a/tests/test_axi_simcoon_structural.py b/tests/test_axi_simcoon_structural.py new file mode 100644 index 00000000..133a8933 --- /dev/null +++ b/tests/test_axi_simcoon_structural.py @@ -0,0 +1,247 @@ +"""Structural tests for the simcoon API on axi-block tensors. + +These tests are independent of fedoo's problem assembly: they feed +a deformation gradient with the block structure expected for an +axisymmetric kinematics into simcoon's lower-level routines and +check that the *structural* zeros (slots 4 and 5 of the 6-vector, +columns/rows 2 of any 3x3 returned tensor) are preserved. + +The 6-vector slot ordering used here is the one documented on +``fedoo.core.mechanical3d.Mechanical3D``:: + + (rr, zz, theta-theta, gamma_rz, 0, 0) + +In the 3x3 layout, the hoop direction lives on row/col 2; the +meridional (r, z) plane occupies rows/cols (0, 1). + +If any of these tests fail, simcoon's polar / objective-rate / +rotation routines do *not* preserve the axi block structure +exactly, and fedoo's UL+axi pipeline cannot be relied upon at +finite strain even after the F[theta-theta] hoop-formula fix. +""" + +import numpy as np +import pytest + +simcoon = pytest.importorskip("simcoon") +sim = simcoon +from simcoon import Rotation + + +def _axi_block_F(F_rr=1.2, F_zz=0.95, F_rz=0.1, F_zr=0.05, F_tt=1.5): + """Build a single 3x3 deformation gradient with axi block structure.""" + F = np.eye(3, dtype=float) + F[0, 0] = F_rr + F[1, 1] = F_zz + F[0, 1] = F_rz + F[1, 0] = F_zr + F[2, 2] = F_tt + # off-block-diagonal entries (hoop coupling) must stay exactly zero + return F + + +def test_objective_rate_log_preserves_axi_block_DR(): + """sim.objective_rate('log', I, F_axi, dt) returns DR with axi block.""" + F1 = np.asfortranarray(_axi_block_F().reshape(3, 3, 1)) + F0 = np.asfortranarray(np.eye(3).reshape(3, 3, 1)) + + DStrain, D, DR, Omega = sim.objective_rate("log", F0, F1, 1.0, True) + + DR0 = DR[:, :, 0] + # off-block-diagonal hoop entries must be exactly 0 + assert DR0[0, 2] == 0.0 + assert DR0[1, 2] == 0.0 + assert DR0[2, 0] == 0.0 + assert DR0[2, 1] == 0.0 + # hoop direction does not rotate in axisymmetric kinematics + assert DR0[2, 2] == pytest.approx(1.0, abs=1e-14) + + +def test_objective_rate_log_preserves_axi_block_DStrain(): + """The 6-vector strain rate has zero shear in slots 4 and 5.""" + F1 = np.asfortranarray(_axi_block_F().reshape(3, 3, 1)) + F0 = np.asfortranarray(np.eye(3).reshape(3, 3, 1)) + + DStrain, D, DR, Omega = sim.objective_rate("log", F0, F1, 1.0, True) + + # DStrain is the 6-vector strain increment: slots 4, 5 = gamma_xz, + # gamma_yz must remain identically zero + assert DStrain[4, 0] == 0.0 + assert DStrain[5, 0] == 0.0 + # D is the 3x3 rate-of-deformation tensor: hoop row/col must be zero + # off-block-diagonal (only the diagonal D[2,2] is allowed to be nonzero) + D0 = D[:, :, 0] + assert D0[0, 2] == 0.0 and D0[1, 2] == 0.0 + assert D0[2, 0] == 0.0 and D0[2, 1] == 0.0 + + +def test_objective_rate_log_preserves_axi_block_Omega(): + """The spin tensor has zero entries in the hoop row / column.""" + F1 = np.asfortranarray(_axi_block_F().reshape(3, 3, 1)) + F0 = np.asfortranarray(np.eye(3).reshape(3, 3, 1)) + + _, _, _, Omega = sim.objective_rate("log", F0, F1, 1.0, True) + + Omega0 = Omega[:, :, 0] + assert Omega0[0, 2] == 0.0 + assert Omega0[1, 2] == 0.0 + assert Omega0[2, 0] == 0.0 + assert Omega0[2, 1] == 0.0 + # hoop axis does not spin in axisymmetric motion + assert Omega0[2, 2] == pytest.approx(0.0, abs=1e-14) + + +@pytest.mark.parametrize("corate", ["log", "log_R", "jaumann", "green_naghdi"]) +def test_all_objective_rates_preserve_axi_block(corate): + """Every objective rate fedoo exposes should preserve the axi block.""" + F1 = np.asfortranarray(_axi_block_F().reshape(3, 3, 1)) + F0 = np.asfortranarray(np.eye(3).reshape(3, 3, 1)) + + DStrain, D, DR, Omega = sim.objective_rate(corate, F0, F1, 1.0, True) + + DR0 = DR[:, :, 0] + assert DR0[0, 2] == 0.0 and DR0[1, 2] == 0.0 + assert DR0[2, 0] == 0.0 and DR0[2, 1] == 0.0 + assert DR0[2, 2] == pytest.approx(1.0, abs=1e-14) + assert DStrain[4, 0] == 0.0 and DStrain[5, 0] == 0.0 + + +def test_log_strain_voigt_preserves_axi_zeros(): + """sim.Log_strain in Voigt form has zero gamma_xz / gamma_yz.""" + F = _axi_block_F() + eps = sim.Log_strain(F, True, False) # voigt_form=True + # eps has shape (6, 1) — flatten for easier indexing + eps = np.asarray(eps).ravel() + assert eps[4] == 0.0 + assert eps[5] == 0.0 + + +def test_log_strain_hoop_component_equals_ln_lambda_theta(): + """Slot 2 of the Voigt log strain is ln(F_theta-theta).""" + lam_theta = 1.5 + F = _axi_block_F(F_tt=lam_theta) + eps = np.asarray(sim.Log_strain(F, True, False)).ravel() + assert eps[2] == pytest.approx(np.log(lam_theta), rel=1e-12) + + +def test_log_strain_hoop_independent_of_meridional_block(): + """Hoop log strain does not couple to r-z deformation.""" + lam_theta = 1.3 + eps_a = np.asarray( + sim.Log_strain(_axi_block_F(F_tt=lam_theta), True, False) + ).ravel() + eps_b = np.asarray( + sim.Log_strain( + _axi_block_F(F_rr=2.0, F_zz=0.5, F_rz=0.3, F_zr=0.2, F_tt=lam_theta), + True, + False, + ) + ).ravel() + # slot 2 = ln(lam_theta) regardless of meridional block + assert eps_a[2] == pytest.approx(eps_b[2], rel=1e-12) + + +def test_rotation_apply_strain_preserves_axi_zeros(): + """A meridional rotation applied to an axi strain keeps slots 4, 5 = 0.""" + # build an axi-block DR (rotation by 30 deg in r-z plane) + angle = np.pi / 6 + DR_full = np.eye(3) + DR_full[0, 0] = np.cos(angle) + DR_full[0, 1] = -np.sin(angle) + DR_full[1, 0] = np.sin(angle) + DR_full[1, 1] = np.cos(angle) + + rot = Rotation.from_matrix(DR_full[np.newaxis, :, :]) + + # axi strain: (rr, zz, theta-theta, gamma_rz, 0, 0) + strain_axi = np.array([[0.05], [-0.02], [0.10], [0.03], [0.0], [0.0]], order="F") + rotated = rot.apply_strain(strain_axi) + + rotated = np.asarray(rotated) + # off-meridional shears must remain identically zero + assert rotated[4, 0] == pytest.approx(0.0, abs=1e-14) + assert rotated[5, 0] == pytest.approx(0.0, abs=1e-14) + + +def test_rotation_apply_strain_preserves_hoop_invariance(): + """A meridional rotation does not change the hoop strain component.""" + angle = np.pi / 4 + DR_full = np.eye(3) + DR_full[0, 0] = np.cos(angle) + DR_full[0, 1] = -np.sin(angle) + DR_full[1, 0] = np.sin(angle) + DR_full[1, 1] = np.cos(angle) + + rot = Rotation.from_matrix(DR_full[np.newaxis, :, :]) + + hoop = 0.07 + strain_axi = np.array([[0.10], [-0.03], [hoop], [0.04], [0.0], [0.0]], order="F") + rotated = np.asarray(rot.apply_strain(strain_axi)) + # slot 2 (hoop) is invariant under rotation about the hoop axis + assert rotated[2, 0] == pytest.approx(hoop, rel=1e-12) + + +def test_rotation_apply_stress_preserves_axi_zeros(): + """Same structural property for stress rotation.""" + angle = -np.pi / 3 + DR_full = np.eye(3) + DR_full[0, 0] = np.cos(angle) + DR_full[0, 1] = -np.sin(angle) + DR_full[1, 0] = np.sin(angle) + DR_full[1, 1] = np.cos(angle) + + rot = Rotation.from_matrix(DR_full[np.newaxis, :, :]) + + # axi stress: (sigma_rr, sigma_zz, sigma_theta, sigma_rz, 0, 0) + stress_axi = np.array([[100.0], [50.0], [80.0], [20.0], [0.0], [0.0]], order="F") + rotated = np.asarray(rot.apply_stress(stress_axi)) + assert rotated[4, 0] == pytest.approx(0.0, abs=1e-12) + assert rotated[5, 0] == pytest.approx(0.0, abs=1e-12) + # hoop stress is invariant under in-plane rotation + assert rotated[2, 0] == pytest.approx(80.0, rel=1e-12) + + +# --------------------------------------------------------------------------- +# Phase 6d: internal-variable rotation through DR derived from objective_rate. +# This mirrors fedoo's set_start codepath +# (see stress_equilibrium.StressEquilibrium.set_start at "rot.apply_strain"): +# the DR returned by sim.objective_rate is fed to Rotation.from_matrix and +# applied to the saved strain history. We verify that, given an axi-block F1, +# this round-trip preserves the slot zeros for the strain history. +# --------------------------------------------------------------------------- + + +def test_set_start_rotation_axi_block_preserves_slot_zeros(): + """DR from objective_rate(I, F_axi) rotates an axi strain to an axi strain.""" + F1 = np.asfortranarray(_axi_block_F().reshape(3, 3, 1)) + F0 = np.asfortranarray(np.eye(3).reshape(3, 3, 1)) + _, _, DR, _ = sim.objective_rate("log", F0, F1, 1.0, True) + + rot = Rotation.from_matrix(DR.transpose(2, 0, 1)) + + strain_axi = np.array([[0.05], [-0.02], [0.10], [0.03], [0.0], [0.0]], order="F") + rotated = np.asarray(rot.apply_strain(strain_axi)) + assert rotated[4, 0] == pytest.approx(0.0, abs=1e-12) + assert rotated[5, 0] == pytest.approx(0.0, abs=1e-12) + # hoop strain unchanged because DR has DR[2,2] = 1 and zero hoop coupling + assert rotated[2, 0] == pytest.approx(0.10, rel=1e-12) + + +def test_set_start_rotation_preserves_axi_zeros_for_pure_meridional_F(): + """A purely meridional rotation of F (no hoop stretch) -> hoop component + of any axi strain remains untouched.""" + F1 = np.eye(3) + angle = np.pi / 8 + F1[0, 0] = np.cos(angle) + F1[0, 1] = -np.sin(angle) + F1[1, 0] = np.sin(angle) + F1[1, 1] = np.cos(angle) + F1 = np.asfortranarray(F1.reshape(3, 3, 1)) + F0 = np.asfortranarray(np.eye(3).reshape(3, 3, 1)) + _, _, DR, _ = sim.objective_rate("log", F0, F1, 1.0, True) + + DR0 = DR[:, :, 0] + # DR must remain axi-block + assert DR0[0, 2] == 0.0 and DR0[1, 2] == 0.0 + assert DR0[2, 0] == 0.0 and DR0[2, 1] == 0.0 + assert DR0[2, 2] == pytest.approx(1.0, abs=1e-14) diff --git a/tests/test_axi_to_3d_rotation.py b/tests/test_axi_to_3d_rotation.py new file mode 100644 index 00000000..e0207900 --- /dev/null +++ b/tests/test_axi_to_3d_rotation.py @@ -0,0 +1,207 @@ +"""Tests for the cylindrical->Cartesian rotation of vector / tensor fields +when revolving 2Daxi data into a full 3D representation. + +These tests pin down the analytical rotation rules for the helper +``_revolve_axi_field``: + +* Scalar fields are tiled along theta unchanged. +* 2-vector fields ``(dr, dz)`` map to 3D Cartesian + ``(dr*cos, dr*sin, dz)`` at each theta. +* 6-tensor fields in fedoo's slot ordering ``(rr, zz, tt, rz, 0, 0)`` + rotate about the symmetry axis Z by theta to produce the Cartesian + Voigt 6-vector. + +If a future change accidentally drops the rotation, these tests fail. +""" + +import numpy as np +import pytest + +from fedoo.post_processing.axi_to_3d import _revolve_axi_field + + +def _theta(n=4): + return np.linspace(0, 2 * np.pi, n, endpoint=False) + + +def test_scalar_field_is_tiled(): + n_nodes = 3 + n_theta = 5 + field = np.array([1.0, 2.0, 3.0]) + theta = _theta(n_theta) + out = _revolve_axi_field(field, theta, n_nodes) + assert out.shape == (1, n_nodes * n_theta) + # node-major layout: each 2D node value repeated for every theta. + # Reshape to (n_nodes_2d, n_theta); each row should be constant. + grid = out.reshape(1, n_nodes, n_theta) + for i in range(n_nodes): + assert np.allclose(grid[0, i], field[i]) + + +def test_vector_field_disp(): + n_nodes = 2 + # node 0: (dr, dz) = (1, 0); node 1: (dr, dz) = (0, 5) + field = np.array([[1.0, 0.0], [0.0, 5.0]]) + theta = _theta(4) # 0, pi/2, pi, 3pi/2 + out = _revolve_axi_field(field, theta, n_nodes) + + assert out.shape == (3, n_nodes * 4) + + # node-major layout: index [node_i, theta_k] after reshape (n_nodes, n_theta). + out_X = out[0].reshape(n_nodes, 4) + out_Y = out[1].reshape(n_nodes, 4) + out_Z = out[2].reshape(n_nodes, 4) + + # node 0 at theta = 0: (1*cos0, 1*sin0, 0) = (1, 0, 0) + assert out_X[0, 0] == pytest.approx(1.0) + assert out_Y[0, 0] == pytest.approx(0.0) + assert out_Z[0, 0] == pytest.approx(0.0) + + # node 0 at theta = pi/2: (cos, sin, 0) = (0, 1, 0) + assert out_X[0, 1] == pytest.approx(0.0, abs=1e-12) + assert out_Y[0, 1] == pytest.approx(1.0) + + # node 1 at theta = pi: (0*cos, 0*sin, 5) = (0, 0, 5) regardless + assert out_X[1, 2] == pytest.approx(0.0) + assert out_Y[1, 2] == pytest.approx(0.0, abs=1e-12) + assert out_Z[1, 2] == pytest.approx(5.0) + + +def test_tensor_pure_radial_stress(): + """Pure sigma_rr at theta=0 -> sigma_xx in Cartesian; at theta=pi/2 -> sigma_yy.""" + n_nodes = 1 + # (s_rr, s_zz, s_tt, s_rz, 0, 0) = (1, 0, 0, 0, 0, 0) + field = np.array([[1.0], [0.0], [0.0], [0.0], [0.0], [0.0]]) + theta = np.array([0.0, np.pi / 2]) + out = _revolve_axi_field(field, theta, n_nodes) + assert out.shape == (6, 2) + + # at theta = 0: sigma_xx = 1, sigma_yy = 0, sigma_zz = 0, off-diag = 0 + assert out[0, 0] == pytest.approx(1.0) + assert out[1, 0] == pytest.approx(0.0, abs=1e-12) + assert out[3, 0] == pytest.approx(0.0, abs=1e-12) + + # at theta = pi/2: sigma_xx = 0, sigma_yy = 1 + assert out[0, 1] == pytest.approx(0.0, abs=1e-12) + assert out[1, 1] == pytest.approx(1.0) + assert out[3, 1] == pytest.approx(0.0, abs=1e-12) + + +def test_tensor_pure_hoop_stress(): + """Pure sigma_tt -> rotates to sigma_yy at theta=0 and sigma_xx at theta=pi/2.""" + n_nodes = 1 + field = np.array([[0.0], [0.0], [1.0], [0.0], [0.0], [0.0]]) + theta = np.array([0.0, np.pi / 2]) + out = _revolve_axi_field(field, theta, n_nodes) + + # at theta=0: e_theta = (0, 1, 0), so sigma_tt -> sigma_yy + assert out[0, 0] == pytest.approx(0.0, abs=1e-12) + assert out[1, 0] == pytest.approx(1.0) + # at theta=pi/2: e_theta = (-1, 0, 0), so sigma_tt -> sigma_xx + assert out[0, 1] == pytest.approx(1.0) + assert out[1, 1] == pytest.approx(0.0, abs=1e-12) + + +def test_tensor_pure_axial_stress_invariant(): + """sigma_zz_cyl maps directly to sigma_zz_cartesian at every theta.""" + n_nodes = 1 + field = np.array([[0.0], [3.0], [0.0], [0.0], [0.0], [0.0]]) + theta = _theta(8) + out = _revolve_axi_field(field, theta, n_nodes) + # slot 2 of the Cartesian 6-vector should be 3.0 at every theta + assert np.allclose(out[2], 3.0) + # all other slots should be zero + assert np.allclose(out[0], 0.0) + assert np.allclose(out[1], 0.0) + assert np.allclose(out[3], 0.0, atol=1e-12) + assert np.allclose(out[4], 0.0, atol=1e-12) + assert np.allclose(out[5], 0.0, atol=1e-12) + + +def test_tensor_rz_shear_rotates_to_xz_yz(): + """sigma_rz at theta=0 -> sigma_xz, at theta=pi/2 -> sigma_yz.""" + n_nodes = 1 + field = np.array([[0.0], [0.0], [0.0], [2.0], [0.0], [0.0]]) + theta = np.array([0.0, np.pi / 2, np.pi]) + out = _revolve_axi_field(field, theta, n_nodes) + + # theta = 0: e_r = (1, 0, 0), so sigma_rz -> sigma_xz (slot 4) + assert out[4, 0] == pytest.approx(2.0) + assert out[5, 0] == pytest.approx(0.0, abs=1e-12) + # theta = pi/2: e_r = (0, 1, 0), so sigma_rz -> sigma_yz (slot 5) + assert out[4, 1] == pytest.approx(0.0, abs=1e-12) + assert out[5, 1] == pytest.approx(2.0) + # theta = pi: e_r = (-1, 0, 0), so sigma_rz -> -sigma_xz + assert out[4, 2] == pytest.approx(-2.0) + assert out[5, 2] == pytest.approx(0.0, abs=1e-12) + + +def test_tensor_isotropic_pressure_invariant(): + """An isotropic cylindrical state -p*I rotates to -p*I in Cartesian.""" + n_nodes = 1 + p = 5.0 + field = np.array([[-p], [-p], [-p], [0.0], [0.0], [0.0]]) + theta = _theta(8) + out = _revolve_axi_field(field, theta, n_nodes) + # diagonals all -p, off-diagonals all 0 + assert np.allclose(out[0], -p) + assert np.allclose(out[1], -p) + assert np.allclose(out[2], -p) + assert np.allclose(out[3], 0.0, atol=1e-12) + assert np.allclose(out[4], 0.0, atol=1e-12) + assert np.allclose(out[5], 0.0, atol=1e-12) + + +def test_strain_field_uses_engineering_shear_for_xy(): + """For a Strain field, slot 3 of the output is gamma_xy = 2 * eps_xy. + + A pure radial strain eps_rr at theta=pi/4 in Cartesian gives + eps_xy_tensor = sc * (eps_rr - 0) = (1/2) * eps_rr + gamma_xy_engineering = 2 * eps_xy_tensor = eps_rr. + """ + n_nodes = 1 + eps_rr = 0.04 + field = np.array([[eps_rr], [0.0], [0.0], [0.0], [0.0], [0.0]]) + theta = np.array([np.pi / 4]) + out = _revolve_axi_field(field, theta, n_nodes, field_name="Strain") + assert out[3, 0] == pytest.approx(eps_rr, rel=1e-12) + + +def test_stress_field_uses_tensor_xy_no_factor_of_2(): + """For a stress (default), slot 3 = sigma_xy with no engineering factor.""" + n_nodes = 1 + sigma_rr = 100.0 + field = np.array([[sigma_rr], [0.0], [0.0], [0.0], [0.0], [0.0]]) + theta = np.array([np.pi / 4]) + out = _revolve_axi_field(field, theta, n_nodes, field_name="Stress") + # sigma_xy = sin*cos*(sigma_rr - 0) = (1/2)*sigma_rr at theta = pi/4 + assert out[3, 0] == pytest.approx(0.5 * sigma_rr, rel=1e-12) + + +def test_explicit_kind_overrides_field_name_heuristic(): + """When kind is explicit, it wins over any name-based detection.""" + n_nodes = 1 + val = 0.04 + field = np.array([[val], [0.0], [0.0], [0.0], [0.0], [0.0]]) + theta = np.array([np.pi / 4]) + + # name says "Strain" but kind="stress" -> no factor of 2 + out_stress = _revolve_axi_field( + field, theta, n_nodes, kind="stress", field_name="MyStrainLikeField" + ) + assert out_stress[3, 0] == pytest.approx(0.5 * val, rel=1e-12) + + # name says "Stress" but kind="strain" -> factor of 2 + out_strain = _revolve_axi_field( + field, theta, n_nodes, kind="strain", field_name="MyStressLikeField" + ) + assert out_strain[3, 0] == pytest.approx(val, rel=1e-12) + + +def test_invalid_kind_raises(): + """Unknown kind must raise.""" + n_nodes = 1 + field = np.zeros((6, 1)) + theta = np.array([0.0]) + with pytest.raises(ValueError, match="kind"): + _revolve_axi_field(field, theta, n_nodes, kind="bogus")