diff --git a/.gitignore b/.gitignore new file mode 100644 index 00000000..e43b0f98 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +.DS_Store diff --git a/docs/Install.rst b/docs/Install.rst index c1e2bdf9..84930d29 100644 --- a/docs/Install.rst +++ b/docs/Install.rst @@ -13,38 +13,80 @@ Minimal installation with pip: $ pip install fedoo -Full installation with pip with all dependencies: - -.. code-block:: none - - $ pip install fedoo[all] - - The required dependencies that are automatically installed with fedoo are: * `Numpy `_ - - * `Scipy `_ mainly for sparse matrices. -In addition, the conda package also includes some recommanded dependencies: + * `Scipy `_ mainly for sparse matrices. - * `Simcoon `_ +In addition, the conda package also includes some recommended dependencies: + + * `Simcoon `_ brings a lot of features (finite strain, non-linear constitutive laws, ...). Simcoon can be installed using conda or pip. - * `PyVista `_ - for results visulisation and mesh utils. - + * `PyVista `_ + for results visualization and mesh utils. + * An efficient sparse matrix solver (pypardiso or petsc4py) depending on the processor as described below. -It is highly recommanded to also install a fast direct sparse matrix solver +Full pip install +---------------- + +It is also possible to install fedoo with all recommended dependencies +(sparse solver, plotting, IPC contact) in one line: + +.. code-block:: none + + $ pip install fedoo[all] + +This installs the following optional groups: ``solver``, ``plot``, +``simcoon``, ``test`` and ``ipc``. + +``pyvistaqt``, which is required for the viewer, is not included in the all +group. This allows you to choose your preferred Qt binding (``pyqt5``, +``pyqt6`` or ``pyside6``). We recommend installing only one of these to avoid +potential library conflicts. + +To enable the viewer, you can install the dependencies explicitly: + +.. code-block:: none + + $ pip install fedoo[all] pyvistaqt pyqt5 + +Alternatively, use the ``gui`` install group that includes ``pyvistaqt`` +and ``pyside6``: + +.. code-block:: none + + $ pip install fedoo[all, gui] + + +Individual optional groups +-------------------------- + +You can also install optional groups individually: + +.. code-block:: none + + $ pip install fedoo[solver] # fast sparse solver (pypardiso or umfpack) + $ pip install fedoo[plot] # matplotlib + pyvista + $ pip install fedoo[simcoon] # simcoon + $ pip install fedoo[ipc] # IPC contact (ipctk) + $ pip install fedoo[gui] # pyvistaqt + pyside6 + + +Sparse solvers +-------------- + +It is highly recommended to install a fast direct sparse matrix solver to improve performances: - * `Pypardiso `_ + * `Pypardiso `_ for intel processors (binding to the pardiso solver) - + * `Petsc4Py `_ mainly compatible with linux or macos including the MUMPS solver. @@ -53,11 +95,20 @@ to improve performances: To be able to launch the fedoo viewer, the module `pyvistaqt `_ is also required. -As mentioned earlier, a lot of features (finite strain, non-linear -constitutive laws, ...) requires the installation of simcoon. Simcoon can -be installed from pip or conda. To install simcoon using conda: + +Simcoon +------- + +Many features (such as finite strain and non-linear constitutive laws) require +Simcoon to be installed. Simcoon is available via both pip and conda. +To install Simcoon individually, use either: .. code-block:: none $ conda install -c conda-forge -c set3mah simcoon - \ No newline at end of file + +Or: + +.. code-block:: none + + $ pip install simcoon \ No newline at end of file diff --git a/docs/boundary_conditions.rst b/docs/boundary_conditions.rst index 1769848c..8d58aceb 100644 --- a/docs/boundary_conditions.rst +++ b/docs/boundary_conditions.rst @@ -135,10 +135,16 @@ fedoo. They can be created and add to the problem with the pb.bc.add method. Contact ============================== -A first implementation of the contact is proposed for 2D problems. Contact -is still in development and subject to slight change in future versions. +Fedoo provides two contact approaches: a **penalty-based** method and an +**IPC (Incremental Potential Contact)** method. Both are implemented as +assembly objects and can be combined with other assemblies using +:py:meth:`fedoo.Assembly.sum`. -For now, only frictionless contact is implemented. +Penalty-based contact +______________________ + +The penalty method uses a node-to-surface formulation. It is available for +2D problems and supports frictionless contact. .. autosummary:: :toctree: generated/ @@ -153,10 +159,70 @@ to a problem, we need first to create the contact assembly (using the class assembly with :py:meth:`fedoo.Assembly.sum`. -Example --------------- +IPC contact +______________________ + +The IPC (Incremental Potential Contact) method uses barrier potentials from +the `ipctk `_ library to guarantee intersection-free +configurations. It supports both 2D and 3D problems, friction, and optional +CCD (Continuous Collision Detection) line search. + +Unlike the penalty method, IPC does **not** require tuning a penalty +parameter. The barrier stiffness :math:`\kappa` is automatically computed +and adaptively updated to balance the elastic and contact forces. The only +physical parameter is ``dhat`` — the barrier activation distance that +controls the minimum gap between surfaces (default: 0.1% of the bounding +box diagonal). + +**Choosing dhat** --- The default ``dhat=1e-3`` (relative) means the +barrier activates when surfaces are within 0.1 % of the bounding-box +diagonal. For problems with a very small initial gap, increase ``dhat`` +(e.g. ``1e-2``) so the barrier catches contact early. For tight-fitting +assemblies where a visible gap is unacceptable, decrease it (e.g. +``1e-4``), but expect more Newton–Raphson iterations. + +**CCD line search** --- Enabling ``use_ccd=True`` is recommended for +problems where first contact occurs suddenly (e.g. a punch hitting a +plate) or where self-contact can cause rapid topology changes. + +**Energy-based backtracking** --- When ``use_ccd=True``, an +energy-based backtracking phase is automatically enabled after CCD +filtering: the step is halved until total energy (exact barrier + +quadratic elastic approximation) decreases. This matches the +reference IPC algorithm and improves convergence robustness. Set +``line_search_energy=False`` to disable (faster but may degrade +convergence for difficult contact scenarios). + +**Convergence criterion** --- The ``'Force'`` convergence criterion +is recommended for IPC contact problems. It measures the relative +decrease of the force residual, matching the gradient-norm convergence +used by reference IPC implementations:: + + pb.set_nr_criterion('Force', tol=5e-3, max_subiter=15) + +The ``'Displacement'`` criterion may become unreliable as contact +stiffness grows. + +The ``ipctk`` package is required and can be installed with: + +.. code-block:: bash + + pip install ipctk + # or + pip install fedoo[ipc] + +.. autosummary:: + :toctree: generated/ + :template: custom-class-template.rst + + fedoo.constraint.IPCContact + fedoo.constraint.IPCSelfContact -Here an example of a contact between a square and a disk. + +Penalty contact example +__________________________ + +Here an example of a contact between a square and a disk using the penalty method. .. code-block:: python @@ -195,13 +261,13 @@ Here an example of a contact between a square and a disk. contact.contact_search_once = True contact.eps_n = 5e5 - #---- Material properties -------------- + #---- Material properties -------------- props = np.array([200e3, 0.3, 1e-5, 300, 1000, 0.3]) - # E, nu, alpha (non used), Re, k, m - material_rect = fd.constitutivelaw.Simcoon("EPICP", props) + # E, nu, alpha (non used), Re, k, m + material_rect = fd.constitutivelaw.Simcoon("EPICP", props) material_disk = fd.constitutivelaw.ElasticIsotrop(50e3, 0.3) #E, nu material = fd.constitutivelaw.Heterogeneous( - (material_rect, material_disk), + (material_rect, material_disk), ('rect', 'disk') ) @@ -245,13 +311,92 @@ Here an example of a contact between a square and a disk. # Write movie with default options # ------------------------------------ results.write_movie(filename, - 'Stress', - component = 'XX', - data_type = 'Node', - framerate = 24, - quality = 5, + 'Stress', + component = 'XX', + data_type = 'Node', + framerate = 24, + quality = 5, clim = [-3e3, 3e3] ) -Video of results: +Video of results: :download:`contact video <./_static/examples/disk_rectangle_contact.mp4>` + + +IPC contact example +__________________________ + +The same disk-rectangle contact problem can be solved with the IPC method. +The IPC method does not require choosing slave/master nodes or tuning a +penalty parameter. The barrier stiffness is automatically computed and +adapted. + +.. code-block:: python + + import fedoo as fd + import numpy as np + + fd.ModelingSpace("2D") + + #---- Create geometries (same as penalty example) -------------- + mesh_rect = fd.mesh.rectangle_mesh(nx=11, ny=21, + x_min=0, x_max=1, y_min=0, y_max=1, + elm_type='quad4', name='Domain') + mesh_disk = fd.mesh.disk_mesh(radius=0.5, nr=6, nt=6, elm_type='quad4') + mesh_disk.nodes += np.array([1.5, 0.48]) + mesh = fd.Mesh.stack(mesh_rect, mesh_disk) + + #---- Define IPC contact -------------- + surf = fd.mesh.extract_surface(mesh) + ipc_contact = fd.constraint.IPCContact( + mesh, surface_mesh=surf, + friction_coefficient=0.3, # Coulomb friction coefficient + use_ccd=True, # enable CCD line search for robustness + ) + # barrier_stiffness is auto-computed; dhat defaults to 1e-3 * bbox_diag + + #---- Material and problem setup -------------- + material = fd.constitutivelaw.ElasticIsotrop(200e3, 0.3) + wf = fd.weakform.StressEquilibrium(material, nlgeom=True) + solid_assembly = fd.Assembly.create(wf, mesh) + assembly = fd.Assembly.sum(solid_assembly, ipc_contact) + + pb = fd.problem.NonLinear(assembly) + res = pb.add_output('results', solid_assembly, ['Disp', 'Stress']) + # ... add BCs ... + pb.nlsolve(dt=0.005, tmax=1) + +.. note:: + + When using ``add_output``, pass the **solid assembly** (not the sum). + ``AssemblySum`` objects cannot be used directly for output extraction. + + +IPC self-contact example +__________________________ + +For self-contact problems, use :py:class:`~fedoo.constraint.IPCSelfContact` +which automatically extracts the surface from the volumetric mesh. + +.. code-block:: python + + import fedoo as fd + import numpy as np + + fd.ModelingSpace("3D") + + mesh = fd.Mesh.read("gyroid.vtk") + material = fd.constitutivelaw.ElasticIsotrop(1e5, 0.3) + + # Self-contact: auto surface extraction, auto barrier stiffness + # Add line_search_energy=True for extra robustness (slower) + contact = fd.constraint.IPCSelfContact(mesh, use_ccd=True) + + wf = fd.weakform.StressEquilibrium(material, nlgeom="UL") + solid = fd.Assembly.create(wf, mesh) + assembly = fd.Assembly.sum(solid, contact) + + pb = fd.problem.NonLinear(assembly) + res = pb.add_output("results", solid, ["Disp", "Stress"]) + # ... add BCs ... + pb.nlsolve(dt=0.05, tmax=1, update_dt=True) diff --git a/examples/03-advanced/tube_compression.py b/examples/03-advanced/tube_compression.py index 434b8651..f472333f 100644 --- a/examples/03-advanced/tube_compression.py +++ b/examples/03-advanced/tube_compression.py @@ -55,7 +55,7 @@ # contact parameters contact.contact_search_once = True -contact.eps_n = 1e6 # contact penalty +contact.eps_n = 1e6 / (2 * np.pi * 24) # contact penalty (scaled for 2πr weighting) contact.max_dist = 1.5 # max distance for the contact search ############################################################################### diff --git a/examples/contact/benchmark/hertz_axisymmetric.py b/examples/contact/benchmark/hertz_axisymmetric.py new file mode 100644 index 00000000..b5272959 --- /dev/null +++ b/examples/contact/benchmark/hertz_axisymmetric.py @@ -0,0 +1,383 @@ +""" +Hertz axisymmetric contact benchmark: rigid sphere on elastic half-space +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Validates penalty contact in fedoo's ``2Daxi`` modeling space against the +classical Hertz (1882) analytical solution for a rigid sphere pressed +onto an elastic half-space. + +The 2D axisymmetric formulation makes the 3D problem computationally +cheap (2D mesh on the r-z plane) while producing the full 3D result +including the ``2*pi*r`` circumferential integration. + +**Two-body setup (r-z plane):** + - Elastic substrate (half-space approximation): large rectangle + - Quasi-rigid parabolic indenter: thin body with bottom surface + following z(r) = H_sub + gap + r^2/(2R) + +Analytical solution (rigid indenter on elastic body): + - Effective modulus: ``E* = E / (1 - nu^2)`` + - Force-displacement: ``F = 4/3 * E* * sqrt(R) * delta^(3/2)`` + - Contact radius: ``a = sqrt(R * delta)`` + - Max pressure: ``p0 = 2*E*/(pi) * sqrt(delta/R)`` + +.. note:: + This benchmark is penalty-only. IPCContact does not support ``2Daxi`` + (no radial ``2*pi*r`` weighting in the IPC barrier formulation). +""" + +import fedoo as fd +import numpy as np +import os +from time import time + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + +# ========================================================================= +# Parameters +# ========================================================================= +R = 10.0 # sphere radius +E = 1000.0 # elastic modulus (substrate) +nu = 0.3 # Poisson's ratio +E_star = E / (1.0 - nu**2) # effective modulus (rigid indenter) +E_rigid = 1e6 # quasi-rigid indenter +nu_rigid = 0.3 + +gap = 0.05 # initial gap at r=0 between substrate top and indenter +delta_max = 0.5 # max indentation at tmax=1 + +a_max = np.sqrt(R * delta_max) # max contact radius ~ 2.236 + +R_domain = 25.0 # substrate radial extent (~11 * a_max) +H_sub = 25.0 # substrate height +R_indenter = 5.0 # indenter radial extent (~2.2 * a_max) +H_indenter = 2.0 # indenter thickness + +# Mesh resolution (number of nodes along each axis) +nx_sub = 101 +ny_sub = 101 +nx_ind = 51 +ny_ind = 5 + + +# ========================================================================= +# Analytical solution +# ========================================================================= +def hertz_analytical(R, E_star, delta): + """Hertz force-displacement for rigid sphere on elastic half-space.""" + delta = np.asarray(delta, dtype=float) + d_pos = np.maximum(delta, 0.0) + F = 4.0 / 3.0 * E_star * np.sqrt(R) * d_pos**1.5 + a = np.sqrt(R * d_pos) + p0 = np.where(delta > 0, 2.0 * E_star / np.pi * np.sqrt(d_pos / R), 0.0) + return F, a, p0 + + +# ========================================================================= +# Mesh +# ========================================================================= +fd.ModelingSpace("2Daxi") + +# --- Substrate (elastic half-space approximation) --- +mesh_sub = fd.mesh.rectangle_mesh( + nx=nx_sub, + ny=ny_sub, + x_min=0, + x_max=1, + y_min=0, + y_max=1, + elm_type="quad4", +) +# Biased node placement: dense near r=0 and z=H_sub (contact zone) +nodes = mesh_sub.nodes +nodes[:, 0] = R_domain * nodes[:, 0] ** 1.5 # radial grading +nodes[:, 1] = H_sub * (1.0 - (1.0 - nodes[:, 1]) ** 1.5) # axial: dense at top +mesh_sub.element_sets["substrate"] = np.arange(mesh_sub.n_elements) + +# --- Indenter (quasi-rigid parabolic body) --- +mesh_ind = fd.mesh.rectangle_mesh( + nx=nx_ind, + ny=ny_ind, + x_min=0, + x_max=1, + y_min=0, + y_max=1, + elm_type="quad4", +) +ind_nodes = mesh_ind.nodes +r_frac = ind_nodes[:, 0].copy() +z_frac = ind_nodes[:, 1].copy() + +# Radial grading (dense near r=0) +r_ind = R_indenter * r_frac**1.5 + +# Bottom surface follows parabola: z(r) = H_sub + gap + r^2/(2R) +# Top surface: z = z_bottom + H_indenter +z_bottom = H_sub + gap + r_ind**2 / (2.0 * R) +ind_nodes[:, 0] = r_ind +ind_nodes[:, 1] = z_bottom + z_frac * H_indenter + +mesh_ind.element_sets["indenter"] = np.arange(mesh_ind.n_elements) + +# Store indenter bottom node indices before stacking +indenter_bottom_local = list(mesh_ind.node_sets["bottom"]) + +# --- Stack meshes --- +n_sub = mesh_sub.n_nodes +mesh = fd.Mesh.stack(mesh_sub, mesh_ind) + +# Indenter bottom/top nodes in global numbering +indenter_bottom_global = set(n + n_sub for n in indenter_bottom_local) +indenter_top_local = list(mesh_ind.node_sets["top"]) +nodes_ind_top = np.array([n + n_sub for n in indenter_top_local]) + + +# ========================================================================= +# Contact setup +# ========================================================================= +# Slave nodes: top surface of substrate +nodes_sub_top = mesh.find_nodes("Y", H_sub) +nodes_sub_top = nodes_sub_top[nodes_sub_top < n_sub] # only substrate nodes + +# Master surface: bottom of indenter only (filter out top/side edges) +surf_ind_full = fd.mesh.extract_surface(mesh.extract_elements("indenter")) +keep = [ + i + for i, elem in enumerate(surf_ind_full.elements) + if all(int(n) in indenter_bottom_global for n in elem) +] +from fedoo.core.mesh import Mesh as _Mesh + +surf_indenter = _Mesh( + mesh.nodes, + surf_ind_full.elements[keep], + surf_ind_full.elm_type, + name="indenter_bottom", +) + +print(f"Substrate top nodes: {len(nodes_sub_top)}") +print(f"Indenter bottom edges: {len(keep)}") +print(f"a_max = {a_max:.3f}, E* = {E_star:.1f}, gap = {gap}") + +contact = fd.constraint.Contact( + nodes_sub_top, + surf_indenter, + search_algorithm="bucket", +) +contact.contact_search_once = False +contact.eps_n = 1e5 # ~100*E, penalty stiffness +contact.max_dist = gap + delta_max + 0.5 + + +# ========================================================================= +# Material +# ========================================================================= +mat_sub = fd.constitutivelaw.ElasticIsotrop(E, nu) +mat_ind = fd.constitutivelaw.ElasticIsotrop(E_rigid, nu_rigid) +material = fd.constitutivelaw.Heterogeneous( + (mat_sub, mat_ind), + ("substrate", "indenter"), +) + + +# ========================================================================= +# Assembly and problem +# ========================================================================= +wf = fd.weakform.StressEquilibrium(material, nlgeom=False) +solid = fd.Assembly.create(wf, mesh) +assembly = fd.Assembly.sum(solid, contact) + +pb = fd.problem.NonLinear(assembly) + +if not os.path.isdir("results"): + os.mkdir("results") +pb.add_output("results/hertz_axi", solid, ["Disp", "Stress"]) + + +# ========================================================================= +# Boundary conditions +# ========================================================================= +# Bottom of substrate: fix axial (Y), free radial +nodes_bottom = mesh.find_nodes("Y", 0) +pb.bc.add("Dirichlet", nodes_bottom, "DispY", 0) + +# Symmetry axis (r=0): fix radial displacement +nodes_axis = mesh.find_nodes("X", 0) +pb.bc.add("Dirichlet", nodes_axis, "DispX", 0) + +# Top of indenter: prescribe downward displacement +pb.bc.add("Dirichlet", nodes_ind_top, "DispY", -delta_max) + + +# ========================================================================= +# Solver +# ========================================================================= +pb.set_nr_criterion("Displacement", tol=5e-3, max_subiter=15) + +# --- Tracking --- +r_top = mesh.nodes[nodes_sub_top, 0] +history = {"time": [], "reaction_y": [], "a_num": []} + + +def track(pb): + history["time"].append(pb.time) + F = pb.get_ext_forces("Disp") + history["reaction_y"].append(np.sum(F[1, nodes_ind_top])) + + # Estimate contact radius from displacement profile. + # In Hertz theory u_z(a) = delta/2, so the contact edge is where + # the displacement drops below half the center value. + disp = pb.get_disp() + uy = disp[1, nodes_sub_top] + uy_center = np.min(uy) # most negative = max compression + if uy_center < -1e-10: + in_contact = uy < 0.5 * uy_center + a = np.max(r_top[in_contact]) if np.any(in_contact) else 0.0 + else: + a = 0.0 + history["a_num"].append(a) + + +print("\n" + "=" * 62) +print(" HERTZ AXISYMMETRIC CONTACT BENCHMARK (Penalty)") +print("=" * 62) + +t0 = time() +pb.nlsolve( + dt=0.05, + tmax=1, + update_dt=True, + print_info=1, + interval_output=0.2, + callback=track, + exec_callback_at_each_iter=True, +) +solve_time = time() - t0 +print(f"\nSolve time: {solve_time:.2f} s") + + +# ========================================================================= +# Comparison with analytical solution +# ========================================================================= +t_num = np.array(history["time"]) +F_num = -np.array(history["reaction_y"]) # reaction is negative (compression) +delta_num = delta_max * t_num - gap # effective indentation + +# Analytical curve +delta_ana = np.linspace(0, delta_max - gap, 200) +F_ana, a_ana, p0_ana = hertz_analytical(R, E_star, delta_ana) + +# Numerical at final step +delta_final = delta_num[-1] +F_final = F_num[-1] +F_ana_final = hertz_analytical(R, E_star, delta_final)[0] + +print("\n" + "-" * 62) +print(" RESULTS SUMMARY") +print("-" * 62) +print(f"{'Quantity':<30} {'Numerical':>15} {'Analytical':>15}") +print("-" * 62) +print(f"{'Final indentation':<30} {delta_final:15.4f} {delta_max - gap:15.4f}") +print(f"{'Final force F':<30} {F_final:15.2f} {float(F_ana_final):15.2f}") +err_pct = 100 * (F_final - float(F_ana_final)) / float(F_ana_final) +print(f"{'F error':<30} {err_pct:14.2f}%") + +# Contact radius at final step +a_num_arr = np.array(history["a_num"]) +a_num_final = a_num_arr[-1] +a_ana_val = np.sqrt(R * max(delta_final, 0)) +print(f"{'Contact radius a':<30} {a_num_final:15.4f} {a_ana_val:15.4f}") +if a_ana_val > 0: + a_err = 100 * (a_num_final - a_ana_val) / a_ana_val + print(f"{'a error':<30} {a_err:14.2f}%") +print() + + +# ========================================================================= +# Per-increment detail +# ========================================================================= +print("-" * 62) +print(" Per-increment detail") +print("-" * 62) +print( + f"{'Inc':>4} {'Time':>8} {'Delta':>10} {'F_num':>12} {'F_ana':>12} {'Err%':>8} {'a_num':>8} {'a_ana':>8}" +) +for i in range(len(history["time"])): + d = delta_num[i] + f_n = F_num[i] + a_n = a_num_arr[i] + if d > 0: + f_a = 4.0 / 3.0 * E_star * np.sqrt(R) * d**1.5 + err = 100 * (f_n - f_a) / f_a if f_a > 1e-10 else 0 + a_a = np.sqrt(R * d) + else: + f_a = 0.0 + err = 0.0 + a_a = 0.0 + print( + f"{i+1:4d} {t_num[i]:8.4f} {d:10.4f} {f_n:12.2f} {f_a:12.2f} {err:8.2f} {a_n:8.4f} {a_a:8.4f}" + ) +print() + + +# ========================================================================= +# Plots +# ========================================================================= +try: + import matplotlib.pyplot as plt + + fig, axes = plt.subplots(1, 2, figsize=(14, 6)) + fig.suptitle("Hertz Axisymmetric Contact Benchmark (Penalty)", fontsize=14) + + # --- Force vs indentation --- + ax = axes[0] + ax.plot(delta_ana, F_ana, "k-", linewidth=2, label="Hertz analytical") + mask = delta_num > 0 + if np.any(mask): + ax.plot( + delta_num[mask], + F_num[mask], + "o-", + color="tab:blue", + markersize=4, + label="FEM (penalty)", + ) + ax.set_xlabel(r"Indentation $\delta$") + ax.set_ylabel(r"Contact force $F$") + ax.set_title(r"Force vs Indentation: $F = \frac{4}{3} E^* \sqrt{R}\, \delta^{3/2}$") + ax.legend() + ax.grid(True, alpha=0.3) + ax.set_xlim(left=0) + ax.set_ylim(bottom=0) + + # --- Contact radius vs indentation --- + ax = axes[1] + ax.plot( + delta_ana, a_ana, "k-", linewidth=2, label=r"Analytical $a = \sqrt{R\,\delta}$" + ) + # Plot FEM contact radius at each increment (only where delta > 0) + mask_a = delta_num > 0 + if np.any(mask_a): + ax.plot( + delta_num[mask_a], + a_num_arr[mask_a], + "o-", + color="tab:blue", + markersize=4, + label="FEM (penalty)", + ) + ax.set_xlabel(r"Indentation $\delta$") + ax.set_ylabel(r"Contact radius $a$") + ax.set_title(r"Contact Radius: $a = \sqrt{R\,\delta}$") + ax.legend() + ax.grid(True, alpha=0.3) + ax.set_xlim(left=0) + ax.set_ylim(bottom=0) + + plt.tight_layout() + plt.savefig("results/hertz_axisymmetric_benchmark.png", dpi=150) + print("Benchmark plot saved to results/hertz_axisymmetric_benchmark.png") + plt.show() + +except ImportError: + print("matplotlib not available -- skipping plots.") diff --git a/examples/contact/benchmark/westergaard_sinusoidal.py b/examples/contact/benchmark/westergaard_sinusoidal.py new file mode 100644 index 00000000..4046e13a --- /dev/null +++ b/examples/contact/benchmark/westergaard_sinusoidal.py @@ -0,0 +1,516 @@ +""" +Westergaard sinusoidal contact benchmark: IPC vs penalty +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Validates both IPC and penalty contact against the **exact** Westergaard +analytical solution for periodic sinusoidal contact (plane strain). + +A rigid flat plate is pressed onto an elastic body whose top surface has +a sinusoidal profile z(x) = A*cos(2*pi*x/lam). Periodic BCs in x make the +domain naturally finite — no infinite half-space approximation error. + +Analytical formulas (Westergaard 1939, Johnson 1985): + - Complete contact pressure: p* = pi * E* * A / lam + - Mean pressure vs contact half-width: p_bar = p* sin^2(pi*a/lam) + - Mean displacement (from ContactMechanics): + d(a) = 2A * [-0.5 + 0.5*cos^2(pi*a/lam) + + sin^2(pi*a/lam) * ln(sin(pi*a/lam))] + +.. note:: + The IPC method requires the ``ipctk`` package: + ``pip install ipctk`` or ``pip install fedoo[ipc]``. +""" + +import fedoo as fd +import numpy as np +import os +from time import time + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + +# ========================================================================= +# Parameters +# ========================================================================= +lam = 1.0 # wavelength +A = 0.02 # amplitude (2% of wavelength) +H = 2.0 # elastic block height (2*lam) +H_rigid = 0.1 # rigid plate thickness +gap = 0.005 # initial gap between peak and plate + +E = 1000.0 +nu = 0.3 +E_star = E / (1.0 - nu**2) # plane strain modulus +E_rigid = 1e6 # quasi-rigid +nu_rigid = 0.3 + +p_star = np.pi * E_star * A / lam # complete contact pressure + +# Mesh resolution +nx = 80 +ny_block = 80 +ny_plate = 2 + +# Prescribed displacement at tmax=1 +delta_max = 0.06 # total descent of rigid plate top + + +# ========================================================================= +# Analytical solution (parametric in contact half-width a) +# ========================================================================= +def westergaard_analytical(lam, A, gap, E_star, H_block=None): + """Return (delta_ana, p_bar_ana) arrays parametric in a. + + If H_block is given, adds the finite-block bulk compression correction: + delta_total = gap - d_mean + p_bar * H / E_star + (accounts for the laterally-constrained block compliance that the + semi-infinite half-space solution does not include). + """ + p_star_val = np.pi * E_star * A / lam + a = np.linspace(0.001, 0.499, 500) * lam + s = np.sin(np.pi * a / lam) + c = np.cos(np.pi * a / lam) + p_bar = p_star_val * s**2 + # Mean surface displacement (compression into the block) + d_mean = 2.0 * A * (-0.5 + 0.5 * c**2 + s**2 * np.log(s)) + # Rigid plate descent = gap closure + elastic compression + bulk + delta = gap - d_mean + if H_block is not None: + delta += p_bar * H_block / E_star + return delta, p_bar, p_star_val + + +# ========================================================================= +# Helpers +# ========================================================================= +def add_symmetry_bc(pb, mesh): + """Apply DispX=0 on left/right faces — exact for the cosine symmetry. + + The cosine surface z = A*cos(2*pi*x/lam) has mirror symmetry about + x=0 and x=lam/2, so u_x = 0 at x=0 and x=lam. This is equivalent + to periodic BC with E_xx=0 for this specific geometry. + """ + nodes_left = mesh.find_nodes("X", 0) + nodes_right = mesh.find_nodes("X", lam) + pb.bc.add("Dirichlet", nodes_left, "DispX", 0) + pb.bc.add("Dirichlet", nodes_right, "DispX", 0) + + +def filter_surface_no_lateral(surf_mesh, vol_mesh, tol=1e-8): + """Remove surface edges on left/right boundaries (x=0 or x=lam). + + This prevents IPC from detecting spurious contacts between vertical + edges of the stacked bodies. + """ + nodes = vol_mesh.nodes + keep = [] + for i, elem in enumerate(surf_mesh.elements): + xs = nodes[elem, 0] + # Drop edge if all nodes are on left OR all on right boundary + if np.all(np.abs(xs) < tol) or np.all(np.abs(xs - lam) < tol): + continue + keep.append(i) + if len(keep) == len(surf_mesh.elements): + return surf_mesh # nothing to filter + from fedoo.core.mesh import Mesh + + new_elements = surf_mesh.elements[keep] + return Mesh(nodes, new_elements, surf_mesh.elm_type, name="filtered_surface") + + +# ========================================================================= +# Mesh and material builder +# ========================================================================= +def build_mesh_and_material(): + """Build two-body mesh (elastic block + rigid plate) with sinusoidal top. + + Returns (mesh, material, nodes_bottom, nodes_top_plate, n_block). + """ + fd.ModelingSpace("2D") + + # --- Elastic block --- + mesh_block = fd.mesh.rectangle_mesh( + nx=nx + 1, + ny=ny_block + 1, + x_min=0, + x_max=lam, + y_min=0, + y_max=H, + elm_type="tri3", + ) + # Sinusoidal deformation of top surface (blended through height) + nodes = mesh_block.nodes + nodes[:, 1] += A * np.cos(2 * np.pi * nodes[:, 0] / lam) * (nodes[:, 1] / H) + mesh_block.element_sets["block"] = np.arange(mesh_block.n_elements) + + # --- Rigid flat plate --- + y_plate_bot = H + A + gap + y_plate_top = y_plate_bot + H_rigid + mesh_plate = fd.mesh.rectangle_mesh( + nx=nx + 1, + ny=ny_plate + 1, + x_min=0, + x_max=lam, + y_min=y_plate_bot, + y_max=y_plate_top, + elm_type="tri3", + ) + mesh_plate.element_sets["plate"] = np.arange(mesh_plate.n_elements) + + # --- Stack --- + n_block = mesh_block.n_nodes + mesh = fd.Mesh.stack(mesh_block, mesh_plate) + + # Node sets + nodes_bottom = mesh.find_nodes("Y", 0) + nodes_top_plate = mesh.find_nodes("Y", y_plate_top) + + # --- Material (heterogeneous) --- + mat_block = fd.constitutivelaw.ElasticIsotrop(E, nu) + mat_plate = fd.constitutivelaw.ElasticIsotrop(E_rigid, nu_rigid) + material = fd.constitutivelaw.Heterogeneous( + (mat_block, mat_plate), + ("block", "plate"), + ) + + return mesh, material, nodes_bottom, nodes_top_plate, n_block + + +# ========================================================================= +# Approach 1: IPC contact +# ========================================================================= +print("=" * 62) +print(" IPC CONTACT (Westergaard sinusoidal benchmark)") +print("=" * 62) + +mesh, material, nodes_bottom, nodes_top_plate, n_block = build_mesh_and_material() + +# Extract surface and remove left/right edges to avoid spurious contacts +surf_raw = fd.mesh.extract_surface(mesh) +surf_ipc = filter_surface_no_lateral(surf_raw, mesh) + +ipc_contact = fd.constraint.IPCContact( + mesh, + surface_mesh=surf_ipc, + dhat=0.003, + dhat_is_relative=False, + friction_coefficient=0.0, + use_ccd=True, +) + +wf = fd.weakform.StressEquilibrium(material, nlgeom=False) +solid = fd.Assembly.create(wf, mesh) +assembly = fd.Assembly.sum(solid, ipc_contact) + +pb_ipc = fd.problem.NonLinear(assembly) + +if not os.path.isdir("results"): + os.mkdir("results") +pb_ipc.add_output("results/westergaard_ipc", solid, ["Disp", "Stress"]) + +# Fix bottom of elastic block +pb_ipc.bc.add("Dirichlet", nodes_bottom, "Disp", 0) + +# Prescribe vertical displacement on top of rigid plate +pb_ipc.bc.add("Dirichlet", nodes_top_plate, "DispY", -delta_max) + +# Symmetry: DispX=0 on left/right faces (exact for cosine geometry) +add_symmetry_bc(pb_ipc, mesh) + +pb_ipc.set_nr_criterion("Displacement", tol=5e-3, max_subiter=15) + +# --- Tracking --- +history_ipc = { + "time": [], + "reaction_y": [], + "n_collisions": [], + "kappa": [], + "min_distance": [], +} + + +def track_ipc(pb): + history_ipc["time"].append(pb.time) + F = pb.get_ext_forces("Disp") + history_ipc["reaction_y"].append(np.sum(F[1, nodes_top_plate])) + history_ipc["n_collisions"].append(len(ipc_contact._collisions)) + history_ipc["kappa"].append(ipc_contact._kappa) + verts = ipc_contact._get_current_vertices(pb) + if len(ipc_contact._collisions) > 0: + min_d = ipc_contact._collisions.compute_minimum_distance( + ipc_contact._collision_mesh, + verts, + ) + else: + min_d = float("inf") + history_ipc["min_distance"].append(min_d) + + +t0 = time() +pb_ipc.nlsolve( + dt=0.05, + tmax=1, + update_dt=True, + print_info=1, + interval_output=0.2, + callback=track_ipc, + exec_callback_at_each_iter=True, +) +ipc_time = time() - t0 +print(f"\nIPC solve time: {ipc_time:.2f} s") + + +# ========================================================================= +# Approach 2: Penalty contact +# ========================================================================= +print("\n" + "=" * 62) +print(" PENALTY CONTACT (Westergaard sinusoidal benchmark)") +print("=" * 62) + +mesh2, material2, nodes_bottom2, nodes_top_plate2, n_block2 = build_mesh_and_material() + +# Slave nodes: top of elastic block (sinusoidal surface) +surf_block = fd.mesh.extract_surface(mesh2.extract_elements("block")) +block_boundary = set(np.unique(surf_block.elements).tolist()) +nodes_block_top = np.array( + [ + n + for n in block_boundary + if mesh2.nodes[n, 1] > H * 0.9 # top ~10% of block + ] +) + +# Master surface: bottom of rigid plate +surf_plate = fd.mesh.extract_surface(mesh2.extract_elements("plate")) + +penalty_contact = fd.constraint.Contact( + nodes_block_top, + surf_plate, + search_algorithm="bucket", +) +penalty_contact.contact_search_once = False +penalty_contact.eps_n = 1e5 +penalty_contact.max_dist = gap + A + 0.01 + +wf2 = fd.weakform.StressEquilibrium(material2, nlgeom=False) +solid2 = fd.Assembly.create(wf2, mesh2) +assembly2 = fd.Assembly.sum(solid2, penalty_contact) + +pb_penalty = fd.problem.NonLinear(assembly2) +pb_penalty.add_output("results/westergaard_penalty", solid2, ["Disp", "Stress"]) + +# Fix bottom of elastic block +pb_penalty.bc.add("Dirichlet", nodes_bottom2, "Disp", 0) + +# Prescribe vertical displacement on top of rigid plate +pb_penalty.bc.add("Dirichlet", nodes_top_plate2, "DispY", -delta_max) + +# Symmetry: DispX=0 on left/right faces +add_symmetry_bc(pb_penalty, mesh2) + +pb_penalty.set_nr_criterion("Displacement", tol=5e-3, max_subiter=15) + +# --- Tracking --- +history_penalty = {"time": [], "reaction_y": []} + + +def track_penalty(pb): + history_penalty["time"].append(pb.time) + F = pb.get_ext_forces("Disp") + history_penalty["reaction_y"].append(np.sum(F[1, nodes_top_plate2])) + + +t0 = time() +pb_penalty.nlsolve( + dt=0.05, + tmax=1, + update_dt=True, + print_info=1, + interval_output=0.2, + callback=track_penalty, + exec_callback_at_each_iter=True, +) +penalty_time = time() - t0 +print(f"\nPenalty solve time: {penalty_time:.2f} s") + + +# ========================================================================= +# Comparison summary +# ========================================================================= +print("\n") +print("=" * 62) +print(" PERFORMANCE COMPARISON: IPC vs Penalty") +print("=" * 62) + +rows = [ + ("Total solve time", f"{ipc_time:.2f} s", f"{penalty_time:.2f} s"), + ( + "Total increments", + str(len(history_ipc["time"])), + str(len(history_penalty["time"])), + ), + ( + "Final reaction Fy", + f"{history_ipc['reaction_y'][-1]:.4f}" if history_ipc["reaction_y"] else "N/A", + f"{history_penalty['reaction_y'][-1]:.4f}" + if history_penalty["reaction_y"] + else "N/A", + ), + ("p* (analytical)", f"{p_star:.4f}", f"{p_star:.4f}"), + ( + "IPC min gap distance", + f"{min(d for d in history_ipc['min_distance'] if np.isfinite(d)):.2e}" + if any(np.isfinite(d) for d in history_ipc["min_distance"]) + else "inf", + "N/A", + ), + ( + "IPC final kappa", + f"{history_ipc['kappa'][-1]:.2e}" if history_ipc["kappa"] else "N/A", + "N/A", + ), +] + +w_label, w_val = 28, 18 +print(f"{'Metric':<{w_label}} {'IPC':>{w_val}} {'Penalty':>{w_val}}") +print("-" * (w_label + 2 * w_val + 2)) +for label, v_ipc, v_pen in rows: + print(f"{label:<{w_label}} {v_ipc:>{w_val}} {v_pen:>{w_val}}") +print() + + +# ========================================================================= +# Per-increment detail: IPC +# ========================================================================= +print("-" * 62) +print(" Per-increment detail: IPC") +print("-" * 62) +print( + f"{'Inc':>4} {'Time':>8} {'Reaction Fy':>14} " + f"{'Collisions':>11} {'Min gap':>12} {'Kappa':>12}" +) +for i in range(len(history_ipc["time"])): + min_d = history_ipc["min_distance"][i] + min_d_str = f"{min_d:.2e}" if np.isfinite(min_d) else "inf" + print( + f"{i+1:4d} {history_ipc['time'][i]:8.4f} " + f"{history_ipc['reaction_y'][i]:14.4f} " + f"{history_ipc['n_collisions'][i]:11d} " + f"{min_d_str:>12} " + f"{history_ipc['kappa'][i]:12.2e}" + ) +print() + + +# ========================================================================= +# Per-increment detail: Penalty +# ========================================================================= +print("-" * 62) +print(" Per-increment detail: PENALTY") +print("-" * 62) +print(f"{'Inc':>4} {'Time':>8} {'Reaction Fy':>14}") +for i in range(len(history_penalty["time"])): + print( + f"{i+1:4d} {history_penalty['time'][i]:8.4f} " + f"{history_penalty['reaction_y'][i]:14.4f}" + ) +print() + + +# ========================================================================= +# Plots: numerical vs analytical +# ========================================================================= +try: + import matplotlib.pyplot as plt + + delta_ana, p_bar_ana, _ = westergaard_analytical(lam, A, gap, E_star, H_block=H) + + fig, axes = plt.subplots(1, 2, figsize=(14, 6)) + fig.suptitle("Westergaard Sinusoidal Contact Benchmark", fontsize=14) + + # --- Normalized pressure vs displacement --- + ax = axes[0] + ax.plot( + delta_ana / A, + p_bar_ana / p_star, + "k-", + linewidth=2, + label="Analytical (corrected)", + ) + + # IPC results + if history_ipc["reaction_y"]: + delta_ipc = np.array(history_ipc["time"]) * delta_max + # Mean pressure = -reaction_y / lam (reaction is negative for compression) + p_bar_ipc = -np.array(history_ipc["reaction_y"]) / lam + ax.plot( + delta_ipc / A, + p_bar_ipc / p_star, + "s-", + color="tab:blue", + markersize=4, + label="IPC", + ) + + # Penalty results + if history_penalty["reaction_y"]: + delta_pen = np.array(history_penalty["time"]) * delta_max + p_bar_pen = -np.array(history_penalty["reaction_y"]) / lam + ax.plot( + delta_pen / A, + p_bar_pen / p_star, + "o-", + color="tab:orange", + markersize=4, + label="Penalty", + ) + + ax.set_xlabel(r"$\delta / A$") + ax.set_ylabel(r"$\bar{p} / p^*$") + ax.set_title("Mean Pressure vs Displacement (normalized)") + ax.legend() + ax.grid(True, alpha=0.3) + # Clip to numerical data range + max_delta_A = delta_max / A * 1.1 + ax.set_xlim(0, max_delta_A) + ax.set_ylim(bottom=0) + + # --- IPC-specific: min gap and kappa --- + ax = axes[1] + finite_gaps = [ + (t, d) + for t, d in zip(history_ipc["time"], history_ipc["min_distance"]) + if np.isfinite(d) + ] + if finite_gaps: + t_gap, d_gap = zip(*finite_gaps) + ax.plot(t_gap, d_gap, "s-", color="tab:green", markersize=3, label="Min gap") + if hasattr(ipc_contact, "_actual_dhat") and ipc_contact._actual_dhat is not None: + ax.axhline( + y=ipc_contact._actual_dhat, + color="blue", + linestyle="--", + alpha=0.5, + label=f"dhat = {ipc_contact._actual_dhat:.4f}", + ) + ax.axhline(y=0, color="red", linestyle="--", alpha=0.5, label="Zero (penetration)") + ax.set_xlabel("Time") + ax.set_ylabel("Min gap distance") + ax.set_title("IPC Minimum Gap Distance") + ax.legend() + ax.grid(True, alpha=0.3) + + plt.tight_layout() + plt.savefig("results/westergaard_benchmark.png", dpi=150) + print("Benchmark plot saved to results/westergaard_benchmark.png") + plt.show() + +except ImportError: + print("matplotlib not available -- skipping plots.") + + +# ========================================================================= +# Post-processing (requires pyvista) +# ========================================================================= +# Uncomment the lines below to visualise the final deformed configuration. +# res_ipc.plot("Stress", "vm", "Node", show=True, scale=1, show_nodes=True) +# res_penalty.plot("Stress", "vm", "Node", show=True, scale=1, show_nodes=True) diff --git a/examples/contact/comparison/ccd_vs_ogc.py b/examples/contact/comparison/ccd_vs_ogc.py new file mode 100644 index 00000000..289e33fa --- /dev/null +++ b/examples/contact/comparison/ccd_vs_ogc.py @@ -0,0 +1,310 @@ +""" +CCD vs OGC line-search comparison +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Runs the 2D disk indentation benchmark twice — once with CCD +(Continuous Collision Detection) scalar line search and once with +OGC (Offset Geometric Contact) per-vertex trust-region filtering — +then compares convergence behaviour and accuracy. + +Metrics compared: + + - **NR iterations per converged increment** (fewer is better) + - **Total wall-clock solve time** + - **Force–indentation curve** (both should agree with each other + and with the Hertz analytical prediction) + +.. note:: + Requires ``ipctk``, ``gmsh``, and ``matplotlib``. +""" + +import fedoo as fd +import numpy as np +import os +import tempfile +from time import time + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + +# ========================================================================= +# Parameters (shared by both runs) +# ========================================================================= + +E_plate = 1e3 +E_disk = 1e5 +nu = 0.3 +R = 5.0 +plate_half = 30.0 +plate_h = 40.0 +gap = 0.1 +imposed_disp = -2.0 +dhat_abs = 0.05 + +DT = 0.05 +TMAX = 1.0 +TOL = 5e-3 +MAX_SUBITER = 8 + + +# ========================================================================= +# Helper: build the full problem (mesh + material + BCs) +# ========================================================================= + + +def build_problem(method): + """Build and return (pb, solid, nodes_disk_top, history). + + Parameters + ---------- + method : str + ``"ccd"`` or ``"ogc"``. + """ + fd.ModelingSpace("2D") + + # --- Plate mesh (gmsh) --- + import gmsh + + gmsh.initialize() + gmsh.option.setNumber("General.Verbosity", 1) + + plate_tag = gmsh.model.occ.addRectangle( + -plate_half, + 0, + 0, + 2 * plate_half, + plate_h, + ) + gmsh.model.occ.synchronize() + gmsh.model.addPhysicalGroup(2, [plate_tag], tag=1, name="plate") + + gmsh.model.mesh.field.add("Box", 1) + gmsh.model.mesh.field.setNumber(1, "VIn", 0.3) + gmsh.model.mesh.field.setNumber(1, "VOut", 3.0) + gmsh.model.mesh.field.setNumber(1, "XMin", -6) + gmsh.model.mesh.field.setNumber(1, "XMax", 6) + gmsh.model.mesh.field.setNumber(1, "YMin", plate_h - 4) + gmsh.model.mesh.field.setNumber(1, "YMax", plate_h) + gmsh.model.mesh.field.setNumber(1, "Thickness", 4) + gmsh.model.mesh.field.setAsBackgroundMesh(1) + gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) + gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) + gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) + + gmsh.model.mesh.generate(2) + plate_msh = os.path.join(tempfile.gettempdir(), "plate_cmp.msh") + gmsh.write(plate_msh) + gmsh.finalize() + + mesh_plate = fd.mesh.import_msh(plate_msh, mesh_type="surface") + if mesh_plate.nodes.shape[1] == 3: + mesh_plate.nodes = mesh_plate.nodes[:, :2] + mesh_plate.element_sets["plate"] = np.arange(mesh_plate.n_elements) + + # --- Disk mesh --- + mesh_disk = fd.mesh.disk_mesh(radius=R, nr=12, nt=24, elm_type="tri3") + mesh_disk.nodes += np.array([0, plate_h + R + gap]) + mesh_disk.element_sets["disk"] = np.arange(mesh_disk.n_elements) + + mesh = fd.Mesh.stack(mesh_plate, mesh_disk) + + # --- IPC contact --- + surf = fd.mesh.extract_surface(mesh) + ipc_contact = fd.constraint.IPCContact( + mesh, + surface_mesh=surf, + dhat=dhat_abs, + dhat_is_relative=False, + use_ccd=(method == "ccd"), + use_ogc=(method == "ogc"), + ) + + # --- Material & assembly --- + mat = fd.constitutivelaw.Heterogeneous( + ( + fd.constitutivelaw.ElasticIsotrop(E_plate, nu), + fd.constitutivelaw.ElasticIsotrop(E_disk, nu), + ), + ("plate", "disk"), + ) + wf = fd.weakform.StressEquilibrium(mat, nlgeom=False) + solid = fd.Assembly.create(wf, mesh) + assembly = fd.Assembly.sum(solid, ipc_contact) + + pb = fd.problem.NonLinear(assembly) + + nodes_bottom = mesh.find_nodes("Y", 0) + nodes_disk_top = mesh.find_nodes("Y", mesh.bounding_box.ymax) + + pb.bc.add("Dirichlet", nodes_bottom, "Disp", 0) + pb.bc.add("Dirichlet", nodes_disk_top, "Disp", [0, imposed_disp]) + pb.set_nr_criterion("Displacement", tol=TOL, max_subiter=MAX_SUBITER) + + history = {"time": [], "reaction_y": [], "nr_iters": []} + + return pb, solid, nodes_disk_top, history + + +# ========================================================================= +# Instrumented solve: capture NR iteration count per increment +# ========================================================================= + + +def solve_instrumented(pb, solid, nodes_disk_top, history, label): + """Run nlsolve while tracking per-increment NR iterations.""" + + if not os.path.isdir("results"): + os.mkdir("results") + + # Monkey-patch solve_time_increment to capture NR iteration counts + _original = pb.solve_time_increment + + def _patched(max_subiter=None, tol_nr=None): + convergence, nr_iter, err = _original(max_subiter, tol_nr) + if convergence: + history["nr_iters"].append(nr_iter) + return convergence, nr_iter, err + + pb.solve_time_increment = _patched + + def track(pb): + history["time"].append(pb.time) + F = pb.get_ext_forces("Disp") + history["reaction_y"].append(np.sum(F[1, nodes_disk_top])) + + print("=" * 60) + print(f"2D DISK INDENTATION -- {label}") + print("=" * 60) + t0 = time() + pb.nlsolve( + dt=DT, + tmax=TMAX, + update_dt=True, + print_info=1, + callback=track, + ) + wall_time = time() - t0 + print(f"{label} solve time: {wall_time:.2f} s") + + # Restore original + pb.solve_time_increment = _original + + return wall_time + + +# ========================================================================= +# Run both methods +# ========================================================================= + +pb_ccd, solid_ccd, ndt_ccd, hist_ccd = build_problem("ccd") +wt_ccd = solve_instrumented(pb_ccd, solid_ccd, ndt_ccd, hist_ccd, "CCD") + +pb_ogc, solid_ogc, ndt_ogc, hist_ogc = build_problem("ogc") +wt_ogc = solve_instrumented(pb_ogc, solid_ogc, ndt_ogc, hist_ogc, "OGC") + + +# ========================================================================= +# Summary table +# ========================================================================= + + +def summary(label, hist, wt): + iters = np.array(hist["nr_iters"]) + return { + "label": label, + "increments": len(iters), + "total_nr": int(iters.sum()), + "mean_nr": iters.mean(), + "max_nr": int(iters.max()), + "wall_time": wt, + } + + +s_ccd = summary("CCD", hist_ccd, wt_ccd) +s_ogc = summary("OGC", hist_ogc, wt_ogc) + +print("\n" + "=" * 60) +print("COMPARISON SUMMARY") +print("=" * 60) +header = f"{'':12s} {'Increments':>10s} {'Total NR':>10s} {'Mean NR':>10s} {'Max NR':>10s} {'Time (s)':>10s}" +print(header) +print("-" * len(header)) +for s in (s_ccd, s_ogc): + print( + f"{s['label']:12s} {s['increments']:10d} {s['total_nr']:10d} " + f"{s['mean_nr']:10.2f} {s['max_nr']:10d} {s['wall_time']:10.2f}" + ) + + +# ========================================================================= +# Plots +# ========================================================================= + +try: + import matplotlib.pyplot as plt + + # --- Force–indentation comparison --- + t_ccd = np.array(hist_ccd["time"]) + Fy_ccd = -np.array(hist_ccd["reaction_y"]) + delta_ccd = np.maximum(t_ccd * abs(imposed_disp) - gap, 0) + + t_ogc = np.array(hist_ogc["time"]) + Fy_ogc = -np.array(hist_ogc["reaction_y"]) + delta_ogc = np.maximum(t_ogc * abs(imposed_disp) - gap, 0) + + # Hertz + E_star = 1.0 / ((1 - nu**2) / E_plate + (1 - nu**2) / E_disk) + + def hertz_2d(a): + d = a**2 / (4 * R) * (2 * np.log(4 * R / a) - 1) + P = np.pi * E_star / (4 * R) * a**2 + return d, P + + a_vals = np.linspace(0.01, R * 0.95, 500) + delta_h, F_h = np.array([hertz_2d(a) for a in a_vals]).T + dmax = max(delta_ccd.max(), delta_ogc.max()) + mask = delta_h <= dmax * 1.05 + delta_h, F_h = delta_h[mask], F_h[mask] + + fig, axes = plt.subplots(1, 2, figsize=(13, 5)) + + # Left: Force vs indentation + ax = axes[0] + ax.plot(delta_ccd, Fy_ccd, "o-", ms=4, label="CCD") + ax.plot(delta_ogc, Fy_ogc, "s-", ms=4, label="OGC") + ax.plot(delta_h, F_h, "--", lw=2, color="gray", label="Hertz (2D)") + ax.set_xlabel("Indentation depth (mm)") + ax.set_ylabel("Force per unit thickness (N/mm)") + ax.set_title("Force--Indentation") + ax.legend() + ax.grid(True, alpha=0.3) + + # Right: NR iterations per increment + ax = axes[1] + ax.bar( + np.arange(len(hist_ccd["nr_iters"])) - 0.2, + hist_ccd["nr_iters"], + width=0.4, + label="CCD", + alpha=0.8, + ) + ax.bar( + np.arange(len(hist_ogc["nr_iters"])) + 0.2, + hist_ogc["nr_iters"], + width=0.4, + label="OGC", + alpha=0.8, + ) + ax.set_xlabel("Increment") + ax.set_ylabel("NR iterations") + ax.set_title("Newton--Raphson iterations per increment") + ax.legend() + ax.grid(True, alpha=0.3, axis="y") + + fig.suptitle("CCD vs OGC -- 2D Disk Indentation", fontsize=14, y=1.02) + fig.tight_layout() + fig.savefig("results/ccd_vs_ogc.png", dpi=150, bbox_inches="tight") + print("\nComparison plot saved to results/ccd_vs_ogc.png") + plt.show() + +except ImportError: + print("matplotlib not available -- skipping comparison plots") diff --git a/examples/contact/comparison/ccd_vs_ogc_selfcontact.py b/examples/contact/comparison/ccd_vs_ogc_selfcontact.py new file mode 100644 index 00000000..3af1bfa1 --- /dev/null +++ b/examples/contact/comparison/ccd_vs_ogc_selfcontact.py @@ -0,0 +1,210 @@ +""" +CCD vs OGC self-contact comparison (hole plate compression) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Runs the 2D hole-plate self-contact benchmark twice — once with CCD +and once with OGC — then compares convergence behaviour, timing, and +reaction-force curves. + +.. note:: + Requires ``ipctk`` and ``simcoon``. +""" + +import fedoo as fd +import numpy as np +import os +from time import time + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + +# ========================================================================= +# Parameters +# ========================================================================= + +DHAT_REL = 3e-3 +IMPOSED_DISP = -70.0 +DT = 0.01 +TMAX = 1.0 +TOL = 5e-3 +MAX_SUBITER = 10 + + +# ========================================================================= +# Helper: build problem +# ========================================================================= + + +def build_problem(method): + """Build hole-plate self-contact problem with CCD or OGC.""" + fd.ModelingSpace("2D") + + mesh = fd.mesh.hole_plate_mesh(nr=15, nt=15, length=100, height=100, radius=45) + + contact = fd.constraint.IPCSelfContact( + mesh, + dhat=DHAT_REL, + dhat_is_relative=True, + use_ccd=(method == "ccd"), + use_ogc=(method == "ogc"), + ) + + nodes_top = mesh.find_nodes("Y", mesh.bounding_box.ymax) + nodes_bottom = mesh.find_nodes("Y", mesh.bounding_box.ymin) + + E, nu = 200e3, 0.3 + props = np.array([E, nu, 1e-5, 300, 1000, 0.3]) + material = fd.constitutivelaw.Simcoon("EPICP", props) + + wf = fd.weakform.StressEquilibrium(material, nlgeom="UL") + solid = fd.Assembly.create(wf, mesh) + assembly = fd.Assembly.sum(solid, contact) + + pb = fd.problem.NonLinear(assembly) + + if not os.path.isdir("results"): + os.mkdir("results") + + pb.bc.add("Dirichlet", nodes_bottom, "Disp", 0) + pb.bc.add("Dirichlet", nodes_top, "Disp", [0, IMPOSED_DISP]) + pb.set_nr_criterion("Displacement", tol=TOL, max_subiter=MAX_SUBITER) + + history = {"time": [], "reaction_y": [], "nr_iters": []} + return pb, solid, nodes_top, history + + +# ========================================================================= +# Instrumented solve +# ========================================================================= + + +def solve_instrumented(pb, solid, nodes_top, history, label): + _original = pb.solve_time_increment + + def _patched(max_subiter=None, tol_nr=None): + convergence, nr_iter, err = _original(max_subiter, tol_nr) + if convergence: + history["nr_iters"].append(nr_iter) + return convergence, nr_iter, err + + pb.solve_time_increment = _patched + + def track(pb): + history["time"].append(pb.time) + F = pb.get_ext_forces("Disp") + history["reaction_y"].append(np.sum(F[1, nodes_top])) + + print("=" * 60) + print(f"HOLE PLATE SELF-CONTACT -- {label}") + print("=" * 60) + t0 = time() + pb.nlsolve( + dt=DT, + tmax=TMAX, + update_dt=True, + print_info=1, + callback=track, + ) + wall_time = time() - t0 + print(f"{label} solve time: {wall_time:.2f} s") + + pb.solve_time_increment = _original + return wall_time + + +# ========================================================================= +# Run both +# ========================================================================= + +pb_ccd, solid_ccd, nt_ccd, hist_ccd = build_problem("ccd") +wt_ccd = solve_instrumented(pb_ccd, solid_ccd, nt_ccd, hist_ccd, "CCD") + +pb_ogc, solid_ogc, nt_ogc, hist_ogc = build_problem("ogc") +wt_ogc = solve_instrumented(pb_ogc, solid_ogc, nt_ogc, hist_ogc, "OGC") + + +# ========================================================================= +# Summary +# ========================================================================= + + +def summary(label, hist, wt): + iters = np.array(hist["nr_iters"]) + return { + "label": label, + "increments": len(iters), + "total_nr": int(iters.sum()), + "mean_nr": iters.mean(), + "max_nr": int(iters.max()), + "wall_time": wt, + } + + +s_ccd = summary("CCD", hist_ccd, wt_ccd) +s_ogc = summary("OGC", hist_ogc, wt_ogc) + +print("\n" + "=" * 60) +print("COMPARISON SUMMARY") +print("=" * 60) +header = f"{'':12s} {'Increments':>10s} {'Total NR':>10s} {'Mean NR':>10s} {'Max NR':>10s} {'Time (s)':>10s}" +print(header) +print("-" * len(header)) +for s in (s_ccd, s_ogc): + print( + f"{s['label']:12s} {s['increments']:10d} {s['total_nr']:10d} " + f"{s['mean_nr']:10.2f} {s['max_nr']:10d} {s['wall_time']:10.2f}" + ) + + +# ========================================================================= +# Plots +# ========================================================================= + +try: + import matplotlib.pyplot as plt + + t_ccd = np.array(hist_ccd["time"]) + Fy_ccd = -np.array(hist_ccd["reaction_y"]) + + t_ogc = np.array(hist_ogc["time"]) + Fy_ogc = -np.array(hist_ogc["reaction_y"]) + + fig, axes = plt.subplots(1, 2, figsize=(13, 5)) + + ax = axes[0] + ax.plot(t_ccd, Fy_ccd, "o-", ms=3, label="CCD") + ax.plot(t_ogc, Fy_ogc, "s-", ms=3, label="OGC") + ax.set_xlabel("Normalized time") + ax.set_ylabel("Reaction force (Y)") + ax.set_title("Reaction Force vs Time") + ax.legend() + ax.grid(True, alpha=0.3) + + ax = axes[1] + ax.bar( + np.arange(len(hist_ccd["nr_iters"])) - 0.2, + hist_ccd["nr_iters"], + width=0.4, + label="CCD", + alpha=0.8, + ) + ax.bar( + np.arange(len(hist_ogc["nr_iters"])) + 0.2, + hist_ogc["nr_iters"], + width=0.4, + label="OGC", + alpha=0.8, + ) + ax.set_xlabel("Increment") + ax.set_ylabel("NR iterations") + ax.set_title("Newton-Raphson iterations per increment") + ax.legend() + ax.grid(True, alpha=0.3, axis="y") + + fig.suptitle("CCD vs OGC -- Self-Contact (Hole Plate)", fontsize=14, y=1.02) + fig.tight_layout() + fig.savefig("results/ccd_vs_ogc_selfcontact.png", dpi=150, bbox_inches="tight") + print("\nComparison plot saved to results/ccd_vs_ogc_selfcontact.png") + plt.show() + +except ImportError: + print("matplotlib not available -- skipping plots") diff --git a/examples/contact/comparison/disk_rectangle.py b/examples/contact/comparison/disk_rectangle.py new file mode 100644 index 00000000..70da2592 --- /dev/null +++ b/examples/contact/comparison/disk_rectangle.py @@ -0,0 +1,475 @@ +""" +Disk-rectangle contact: penalty vs IPC comparison +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This example compares the two contact approaches available in fedoo +on a 2D problem: a stiff disk pushed into an elastic rectangle. + + - **Penalty method** (``fd.constraint.Contact``): node-to-surface + formulation with a user-tuned penalty parameter ``eps_n``. + - **IPC method** (``fd.constraint.IPCContact``): barrier-potential + formulation from the ipctk library guaranteeing intersection-free + configurations, with optional friction and CCD line search. + +Both simulations use the same geometry and boundary conditions. +Per-increment metrics are tracked and a comparison summary is printed. + +.. note:: + The IPC method requires the ``ipctk`` package: + ``pip install ipctk`` or ``pip install fedoo[ipc]``. +""" + +import fedoo as fd +import numpy as np +import os +from time import time + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + + +def build_mesh_and_material(): + """Build the shared geometry and material (fresh ModelingSpace each call).""" + fd.ModelingSpace("2D") + + # --- Rectangle mesh --- + mesh_rect = fd.mesh.rectangle_mesh( + nx=11, + ny=21, + x_min=0, + x_max=1, + y_min=0, + y_max=1, + elm_type="quad4", + ) + mesh_rect.element_sets["rect"] = np.arange(mesh_rect.n_elements) + + # --- Disk mesh --- + mesh_disk = fd.mesh.disk_mesh(radius=0.5, nr=6, nt=6, elm_type="quad4") + mesh_disk.nodes += np.array([1.5, 0.48]) + mesh_disk.element_sets["disk"] = np.arange(mesh_disk.n_elements) + + # --- Stack into one mesh --- + mesh = fd.Mesh.stack(mesh_rect, mesh_disk) + + # --- Material --- + E, nu = 200e3, 0.3 + material_rect = fd.constitutivelaw.ElasticIsotrop(E, nu) + material_disk = fd.constitutivelaw.ElasticIsotrop(50e3, nu) + material = fd.constitutivelaw.Heterogeneous( + (material_rect, material_disk), ("rect", "disk") + ) + + return mesh, material + + +# ========================================================================= +# Approach 1 : Penalty contact +# ========================================================================= + +print("=" * 60) +print("PENALTY CONTACT") +print("=" * 60) + +mesh, material = build_mesh_and_material() + +nodes_left = mesh.find_nodes("X", 0) +nodes_bc = mesh.find_nodes("X>1.5") +nodes_bc = list(set(nodes_bc).intersection(mesh.node_sets["boundary"])) + +# Slave nodes = right face of rectangle (exclude disk nodes at X=1) +rect_nodes = set(np.unique(mesh.extract_elements("rect").elements).tolist()) +nodes_contact = np.array([n for n in mesh.find_nodes("X", 1) if n in rect_nodes]) +surf = fd.mesh.extract_surface(mesh.extract_elements("disk")) + +penalty_contact = fd.constraint.Contact(nodes_contact, surf) +penalty_contact.contact_search_once = True +penalty_contact.eps_n = 5e5 +penalty_contact.max_dist = 1.0 + +wf = fd.weakform.StressEquilibrium(material, nlgeom=True) +solid_assembly = fd.Assembly.create(wf, mesh) +assembly = fd.Assembly.sum(solid_assembly, penalty_contact) + +pb_penalty = fd.problem.NonLinear(assembly) + +if not os.path.isdir("results"): + os.mkdir("results") +res_penalty = pb_penalty.add_output( + "results/penalty_contact", solid_assembly, ["Disp", "Stress"] +) + +pb_penalty.bc.add("Dirichlet", nodes_left, "Disp", 0) +pb_penalty.bc.add("Dirichlet", nodes_bc, "Disp", [-0.05, 0.025]) +pb_penalty.set_nr_criterion("Displacement", tol=5e-3, max_subiter=5) + +# --- Tracking callback --- +history_penalty = { + "time": [], + "nr_iters": [], + "reaction_x": [], + "contact_force_norm": [], + "max_disp_x": [], +} +_penalty_prev_niter = [0] # mutable container for closure + + +def track_penalty(pb): + history_penalty["time"].append(pb.time) + # NR iterations for this increment + nr_this = pb.n_iter - _penalty_prev_niter[0] + _penalty_prev_niter[0] = pb.n_iter + history_penalty["nr_iters"].append(nr_this) + # Reaction force at left boundary (X component) + F = pb.get_ext_forces("Disp") + history_penalty["reaction_x"].append(np.sum(F[0, nodes_left])) + # Contact force norm (use .current which holds the updated state) + gv = penalty_contact.current.global_vector + history_penalty["contact_force_norm"].append( + np.linalg.norm(gv) if not np.isscalar(gv) else 0.0 + ) + # Max X-displacement at contact interface + disp = pb.get_disp() + history_penalty["max_disp_x"].append(np.max(np.abs(disp[0, nodes_contact]))) + + +t0 = time() +pb_penalty.nlsolve( + dt=0.05, + tmax=1, + update_dt=True, + print_info=1, + interval_output=0.1, + callback=track_penalty, + exec_callback_at_each_iter=True, +) +penalty_solve_time = time() - t0 +print(f"Penalty solve time: {penalty_solve_time:.2f} s") + + +# ========================================================================= +# Approach 2 : IPC contact +# ========================================================================= + +print("\n" + "=" * 60) +print("IPC CONTACT") +print("=" * 60) + +mesh2, material2 = build_mesh_and_material() + +nodes_left2 = mesh2.find_nodes("X", 0) +nodes_bc2 = mesh2.find_nodes("X>1.5") +nodes_bc2 = list(set(nodes_bc2).intersection(mesh2.node_sets["boundary"])) + +# Nodes at rectangle right face only (exclude disk nodes at X=1) +rect_nodes2 = set(np.unique(mesh2.extract_elements("rect").elements).tolist()) +nodes_contact2 = np.array([n for n in mesh2.find_nodes("X", 1) if n in rect_nodes2]) + +# IPC contact: extract the whole surface -- no slave/master distinction needed +surf_ipc = fd.mesh.extract_surface(mesh2) + +ipc_contact = fd.constraint.IPCContact( + mesh2, + surface_mesh=surf_ipc, + dhat=1e-3, # relative to bbox diagonal + dhat_is_relative=True, + friction_coefficient=0.0, # frictionless (same as penalty) + use_ccd=True, # CCD line search for robustness +) + +wf2 = fd.weakform.StressEquilibrium(material2, nlgeom=True) +solid_assembly2 = fd.Assembly.create(wf2, mesh2) +assembly2 = fd.Assembly.sum(solid_assembly2, ipc_contact) + +pb_ipc = fd.problem.NonLinear(assembly2) + +res_ipc = pb_ipc.add_output("results/ipc_contact", solid_assembly2, ["Disp", "Stress"]) + +pb_ipc.bc.add("Dirichlet", nodes_left2, "Disp", 0) +pb_ipc.bc.add("Dirichlet", nodes_bc2, "Disp", [-0.05, 0.025]) +pb_ipc.set_nr_criterion("Displacement", tol=5e-3, max_subiter=5) + +# --- Tracking callback --- +history_ipc = { + "time": [], + "nr_iters": [], + "reaction_x": [], + "contact_force_norm": [], + "max_disp_x": [], + "n_collisions": [], + "kappa": [], + "min_distance": [], +} +_ipc_prev_niter = [0] + + +def track_ipc(pb): + history_ipc["time"].append(pb.time) + # NR iterations for this increment + nr_this = pb.n_iter - _ipc_prev_niter[0] + _ipc_prev_niter[0] = pb.n_iter + history_ipc["nr_iters"].append(nr_this) + # Reaction force at left boundary (X component) + F = pb.get_ext_forces("Disp") + history_ipc["reaction_x"].append(np.sum(F[0, nodes_left2])) + # Contact force norm + history_ipc["contact_force_norm"].append(np.linalg.norm(ipc_contact.global_vector)) + # Max X-displacement at contact interface (rectangle right face) + disp = pb.get_disp() + history_ipc["max_disp_x"].append(np.max(np.abs(disp[0, nodes_contact2]))) + # IPC-specific metrics + history_ipc["n_collisions"].append(len(ipc_contact._collisions)) + history_ipc["kappa"].append(ipc_contact._kappa) + # Minimum distance (gap) + verts = ipc_contact._get_current_vertices(pb) + if len(ipc_contact._collisions) > 0: + min_d = ipc_contact._collisions.compute_minimum_distance( + ipc_contact._collision_mesh, verts + ) + else: + min_d = float("inf") + history_ipc["min_distance"].append(min_d) + + +t0 = time() +pb_ipc.nlsolve( + dt=0.05, + tmax=1, + update_dt=True, + print_info=1, + interval_output=0.1, + callback=track_ipc, + exec_callback_at_each_iter=True, +) +ipc_solve_time = time() - t0 +print(f"IPC solve time: {ipc_solve_time:.2f} s") + + +# ========================================================================= +# Comparison summary +# ========================================================================= + +print("\n") +print("=" * 62) +print(" PERFORMANCE COMPARISON: Penalty vs IPC Contact") +print("=" * 62) + +total_nr_penalty = sum(history_penalty["nr_iters"]) +total_nr_ipc = sum(history_ipc["nr_iters"]) + +# Build rows: (label, penalty_value, ipc_value) +rows = [ + ("Total solve time", f"{penalty_solve_time:.2f} s", f"{ipc_solve_time:.2f} s"), + ( + "Total increments", + str(len(history_penalty["time"])), + str(len(history_ipc["time"])), + ), + ("Total NR iterations", str(total_nr_penalty), str(total_nr_ipc)), + ( + "Avg NR iter / increment", + f"{total_nr_penalty / max(len(history_penalty['time']), 1):.1f}", + f"{total_nr_ipc / max(len(history_ipc['time']), 1):.1f}", + ), + ( + "Final reaction Fx", + f"{history_penalty['reaction_x'][-1]:.1f}" + if history_penalty["reaction_x"] + else "N/A", + f"{history_ipc['reaction_x'][-1]:.1f}" if history_ipc["reaction_x"] else "N/A", + ), + ( + "Max contact force norm", + f"{max(history_penalty['contact_force_norm']):.1f}" + if history_penalty["contact_force_norm"] + else "N/A", + f"{max(history_ipc['contact_force_norm']):.1f}" + if history_ipc["contact_force_norm"] + else "N/A", + ), + ( + "Final contact force norm", + f"{history_penalty['contact_force_norm'][-1]:.1f}" + if history_penalty["contact_force_norm"] + else "N/A", + f"{history_ipc['contact_force_norm'][-1]:.1f}" + if history_ipc["contact_force_norm"] + else "N/A", + ), + ( + "Max |disp_x|", + f"{max(history_penalty['max_disp_x']):.6f}" + if history_penalty["max_disp_x"] + else "N/A", + f"{max(history_ipc['max_disp_x']):.6f}" if history_ipc["max_disp_x"] else "N/A", + ), + ( + "IPC min gap distance", + "N/A", + f"{min(d for d in history_ipc['min_distance'] if np.isfinite(d)):.2e}" + if any(np.isfinite(d) for d in history_ipc["min_distance"]) + else "inf", + ), + ( + "IPC final barrier kappa", + "N/A", + f"{history_ipc['kappa'][-1]:.2e}" if history_ipc["kappa"] else "N/A", + ), + ( + "IPC max active collisions", + "N/A", + str(max(history_ipc["n_collisions"])) if history_ipc["n_collisions"] else "N/A", + ), +] + +# Print formatted table +w_label = 28 +w_val = 16 +print(f"{'Metric':<{w_label}} {'Penalty':>{w_val}} {'IPC':>{w_val}}") +print("-" * (w_label + 2 * w_val + 2)) +for label, v_pen, v_ipc in rows: + print(f"{label:<{w_label}} {v_pen:>{w_val}} {v_ipc:>{w_val}}") +print() + + +# ========================================================================= +# Per-increment detail tables +# ========================================================================= + +print("-" * 62) +print(" Per-increment detail: PENALTY") +print("-" * 62) +print(f"{'Inc':>4} {'Time':>8} {'NR iter':>8} {'Reaction Fx':>14} {'|F_contact|':>14}") +for i in range(len(history_penalty["time"])): + print( + f"{i+1:4d} {history_penalty['time'][i]:8.4f} " + f"{history_penalty['nr_iters'][i]:8d} " + f"{history_penalty['reaction_x'][i]:14.2f} " + f"{history_penalty['contact_force_norm'][i]:14.2f}" + ) + +print() +print("-" * 62) +print(" Per-increment detail: IPC") +print("-" * 62) +print( + f"{'Inc':>4} {'Time':>8} {'NR iter':>8} {'Reaction Fx':>14} " + f"{'|F_contact|':>14} {'Collisions':>11} {'Min gap':>12} {'Kappa':>12}" +) +for i in range(len(history_ipc["time"])): + min_d = history_ipc["min_distance"][i] + min_d_str = f"{min_d:.2e}" if np.isfinite(min_d) else "inf" + print( + f"{i+1:4d} {history_ipc['time'][i]:8.4f} " + f"{history_ipc['nr_iters'][i]:8d} " + f"{history_ipc['reaction_x'][i]:14.2f} " + f"{history_ipc['contact_force_norm'][i]:14.2f} " + f"{history_ipc['n_collisions'][i]:11d} " + f"{min_d_str:>12} " + f"{history_ipc['kappa'][i]:12.2e}" + ) + +print() + +# ========================================================================= +# Optional plots (requires matplotlib) +# ========================================================================= + +try: + import matplotlib.pyplot as plt + + fig, axes = plt.subplots(2, 2, figsize=(12, 8)) + fig.suptitle("Penalty vs IPC Contact Comparison", fontsize=14) + + # Reaction force vs time + ax = axes[0, 0] + ax.plot( + history_penalty["time"], + history_penalty["reaction_x"], + "o-", + label="Penalty", + markersize=3, + ) + ax.plot( + history_ipc["time"], history_ipc["reaction_x"], "s-", label="IPC", markersize=3 + ) + ax.set_xlabel("Time") + ax.set_ylabel("Reaction Fx") + ax.set_title("Reaction Force (X) at Left Boundary") + ax.legend() + ax.grid(True, alpha=0.3) + + # Contact force norm vs time + ax = axes[0, 1] + ax.plot( + history_penalty["time"], + history_penalty["contact_force_norm"], + "o-", + label="Penalty", + markersize=3, + ) + ax.plot( + history_ipc["time"], + history_ipc["contact_force_norm"], + "s-", + label="IPC", + markersize=3, + ) + ax.set_xlabel("Time") + ax.set_ylabel("||F_contact||") + ax.set_title("Contact Force Norm") + ax.legend() + ax.grid(True, alpha=0.3) + + # NR iterations per increment + ax = axes[1, 0] + ax.bar( + np.arange(len(history_penalty["nr_iters"])) - 0.2, + history_penalty["nr_iters"], + width=0.4, + label="Penalty", + ) + ax.bar( + np.arange(len(history_ipc["nr_iters"])) + 0.2, + history_ipc["nr_iters"], + width=0.4, + label="IPC", + ) + ax.set_xlabel("Increment") + ax.set_ylabel("NR Iterations") + ax.set_title("Newton-Raphson Iterations per Increment") + ax.legend() + ax.grid(True, alpha=0.3) + + # IPC-specific: min gap and kappa + ax = axes[1, 1] + finite_gaps = [ + (t, d) + for t, d in zip(history_ipc["time"], history_ipc["min_distance"]) + if np.isfinite(d) + ] + if finite_gaps: + t_gap, d_gap = zip(*finite_gaps) + ax.plot(t_gap, d_gap, "s-", color="tab:green", markersize=3, label="Min gap") + ax.set_xlabel("Time") + ax.set_ylabel("Min Gap Distance") + ax.set_title("IPC Minimum Gap Distance") + ax.axhline(y=0, color="red", linestyle="--", alpha=0.5, label="Zero (penetration)") + ax.legend() + ax.grid(True, alpha=0.3) + + plt.tight_layout() + plt.savefig("results/penalty_vs_ipc_comparison.png", dpi=150) + print("Comparison plot saved to results/penalty_vs_ipc_comparison.png") + plt.show() + +except ImportError: + print("matplotlib not available -- skipping plots.") + + +# ========================================================================= +# Post-processing (requires pyvista) +# ========================================================================= +# Uncomment the lines below to visualise and compare the results. + +# res_penalty.plot("Stress", "vm", "Node", show=True, scale=1, show_nodes=True) +# res_ipc.plot("Stress", "vm", "Node", show=True, scale=1, show_nodes=True) diff --git a/examples/contact/comparison/gyroid_selfcontact.py b/examples/contact/comparison/gyroid_selfcontact.py new file mode 100644 index 00000000..b6cac509 --- /dev/null +++ b/examples/contact/comparison/gyroid_selfcontact.py @@ -0,0 +1,477 @@ +""" +Gyroid self-contact under compression: penalty vs IPC comparison +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This example compares the two self-contact approaches available in fedoo +on a 3D gyroid unit cell under uniaxial compression. + + - **Penalty method** (``fd.constraint.SelfContact``): node-to-surface + formulation with a user-tuned penalty parameter. + - **IPC method** (``fd.constraint.IPCSelfContact``): barrier-potential + formulation from the ipctk library guaranteeing intersection-free + configurations. + +The gyroid mesh (30k nodes, 114k tet4 elements) is loaded from a periodic +VTK file generated by microgen. Simple compression BCs are applied +(no periodic BC) so that the gyroid walls buckle and self-contact occurs. + +**Key finding:** The penalty self-contact method struggles on this mesh +(30k surface nodes, all-to-all contact search in pure Python) while IPC +(using ipctk's C++ hash-grid broad phase) scales well. + +.. note:: + The IPC method requires the ``ipctk`` package: + ``pip install ipctk`` or ``pip install fedoo[ipc]``. +""" + +import fedoo as fd +import numpy as np +import os +import sys +from time import time + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + +MESH_PATH = os.path.join( + os.path.dirname(os.path.abspath(__file__)), "../../../util/meshes/gyroid_per.vtk" +) +COMPRESSION = 0.15 # 15% compressive strain (height = 1.0) +RUN_PENALTY = "--penalty" in sys.argv # opt-in: very slow on this mesh + + +def build_gyroid(): + """Load gyroid mesh and define elastic material.""" + fd.ModelingSpace("3D") + mesh = fd.Mesh.read(MESH_PATH) + + # Linear elastic material + material = fd.constitutivelaw.ElasticIsotrop(1e5, 0.3) + + return mesh, material + + +# ========================================================================= +# Approach 1 : Penalty self-contact (opt-in, slow on 30k surface nodes) +# ========================================================================= + +penalty_solve_time = None +history_penalty = { + "time": [], + "nr_iters": [], + "reaction_z": [], + "contact_force_norm": [], +} + +if RUN_PENALTY: + print("=" * 60) + print("PENALTY SELF-CONTACT (Gyroid compression)") + print("=" * 60) + + mesh, material = build_gyroid() + + nodes_top = mesh.find_nodes("Z", mesh.bounding_box.zmax) + nodes_bottom = mesh.find_nodes("Z", mesh.bounding_box.zmin) + + surf = fd.mesh.extract_surface(mesh) + penalty_contact = fd.constraint.SelfContact( + surf, "linear", search_algorithm="bucket" + ) + penalty_contact.contact_search_once = True + penalty_contact.eps_n = 1e5 + # max_dist must be very small for gyroid (walls ~0.05-0.1 apart) + penalty_contact.max_dist = 0.03 + + wf = fd.weakform.StressEquilibrium(material, nlgeom="UL") + solid_assembly = fd.Assembly.create(wf, mesh) + assembly = fd.Assembly.sum(solid_assembly, penalty_contact) + + pb_penalty = fd.problem.NonLinear(assembly) + + if not os.path.isdir("results"): + os.mkdir("results") + res_penalty = pb_penalty.add_output( + "results/gyroid_penalty", solid_assembly, ["Disp", "Stress"] + ) + + pb_penalty.bc.add("Dirichlet", nodes_bottom, "Disp", 0) + pb_penalty.bc.add("Dirichlet", nodes_top, "Disp", [0, 0, -COMPRESSION]) + pb_penalty.set_nr_criterion("Displacement", tol=5e-3, max_subiter=10) + + _penalty_prev_niter = [0] + + def track_penalty(pb): + history_penalty["time"].append(pb.time) + nr_this = pb.n_iter - _penalty_prev_niter[0] + _penalty_prev_niter[0] = pb.n_iter + history_penalty["nr_iters"].append(nr_this) + F = pb.get_ext_forces("Disp") + history_penalty["reaction_z"].append(np.sum(F[2, nodes_top])) + gv = penalty_contact.current.global_vector + history_penalty["contact_force_norm"].append( + np.linalg.norm(gv) if not np.isscalar(gv) else 0.0 + ) + + t0 = time() + try: + pb_penalty.nlsolve( + dt=0.01, + tmax=1, + update_dt=True, + print_info=2, + interval_output=0.1, + callback=track_penalty, + exec_callback_at_each_iter=True, + ) + penalty_solve_time = time() - t0 + print(f"Penalty solve time: {penalty_solve_time:.2f} s") + except Exception as e: + penalty_solve_time = time() - t0 + print(f"Penalty FAILED after {penalty_solve_time:.2f} s: {e}") +else: + print("=" * 60) + print("PENALTY SELF-CONTACT: SKIPPED") + print(" (30k surface nodes -- pure-Python contact search is very slow)") + print(" Run with --penalty flag to include penalty method.") + print("=" * 60) + + +# ========================================================================= +# Approach 2 : IPC self-contact +# ========================================================================= + +print("\n" + "=" * 60) +print("IPC SELF-CONTACT (Gyroid compression)") +print("=" * 60) + +mesh2, material2 = build_gyroid() + +nodes_top2 = mesh2.find_nodes("Z", mesh2.bounding_box.zmax) +nodes_bottom2 = mesh2.find_nodes("Z", mesh2.bounding_box.zmin) + +ipc_contact = fd.constraint.IPCSelfContact( + mesh2, + dhat=5e-3, + dhat_is_relative=True, + friction_coefficient=0.0, + use_ccd=True, +) + +wf2 = fd.weakform.StressEquilibrium(material2, nlgeom="UL") +solid_assembly2 = fd.Assembly.create(wf2, mesh2) +assembly2 = fd.Assembly.sum(solid_assembly2, ipc_contact) + +pb_ipc = fd.problem.NonLinear(assembly2) + +if not os.path.isdir("results"): + os.mkdir("results") +res_ipc = pb_ipc.add_output("results/gyroid_ipc", solid_assembly2, ["Disp", "Stress"]) + +pb_ipc.bc.add("Dirichlet", nodes_bottom2, "Disp", 0) +pb_ipc.bc.add("Dirichlet", nodes_top2, "Disp", [0, 0, -COMPRESSION]) +pb_ipc.set_nr_criterion("Displacement", tol=5e-3, max_subiter=10) + +# --- Tracking callback --- +history_ipc = { + "time": [], + "nr_iters": [], + "reaction_z": [], + "contact_force_norm": [], + "n_collisions": [], + "kappa": [], + "min_distance": [], +} +_ipc_prev_niter = [0] + + +def track_ipc(pb): + history_ipc["time"].append(pb.time) + nr_this = pb.n_iter - _ipc_prev_niter[0] + _ipc_prev_niter[0] = pb.n_iter + history_ipc["nr_iters"].append(nr_this) + # Reaction force at top (Z component) + F = pb.get_ext_forces("Disp") + history_ipc["reaction_z"].append(np.sum(F[2, nodes_top2])) + # Contact force norm + history_ipc["contact_force_norm"].append(np.linalg.norm(ipc_contact.global_vector)) + # IPC-specific metrics + history_ipc["n_collisions"].append(len(ipc_contact._collisions)) + history_ipc["kappa"].append(ipc_contact._kappa) + # Minimum distance (gap) + verts = ipc_contact._get_current_vertices(pb) + if len(ipc_contact._collisions) > 0: + min_d = ipc_contact._collisions.compute_minimum_distance( + ipc_contact._collision_mesh, verts + ) + else: + min_d = float("inf") + history_ipc["min_distance"].append(min_d) + + +t0 = time() +pb_ipc.nlsolve( + dt=0.05, + tmax=1, + update_dt=True, + print_info=1, + interval_output=0.1, + callback=track_ipc, + exec_callback_at_each_iter=True, +) +ipc_solve_time = time() - t0 +print(f"IPC solve time: {ipc_solve_time:.2f} s") + + +# ========================================================================= +# Comparison summary +# ========================================================================= + +print("\n") +print("=" * 62) +print(" PERFORMANCE COMPARISON: Penalty vs IPC Self-Contact") +print(" (Gyroid unit cell, {:.0f}% compression)".format(COMPRESSION * 100)) +print("=" * 62) + +total_nr_ipc = sum(history_ipc["nr_iters"]) + +if history_penalty["time"]: + total_nr_penalty = sum(history_penalty["nr_iters"]) + pen_time_str = f"{penalty_solve_time:.2f} s" if penalty_solve_time else "FAILED" + pen_inc_str = str(len(history_penalty["time"])) + pen_nr_str = str(total_nr_penalty) + pen_avg_str = f"{total_nr_penalty / max(len(history_penalty['time']), 1):.1f}" + pen_fz_str = ( + f"{history_penalty['reaction_z'][-1]:.2f}" + if history_penalty["reaction_z"] + else "N/A" + ) + pen_fc_str = ( + f"{max(history_penalty['contact_force_norm']):.2f}" + if history_penalty["contact_force_norm"] + else "N/A" + ) + pen_fc_final = ( + f"{history_penalty['contact_force_norm'][-1]:.2f}" + if history_penalty["contact_force_norm"] + else "N/A" + ) +else: + pen_time_str = "SKIPPED" + pen_inc_str = pen_nr_str = pen_avg_str = pen_fz_str = pen_fc_str = pen_fc_final = ( + "---" + ) + +rows = [ + ("Total solve time", pen_time_str, f"{ipc_solve_time:.2f} s"), + ("Total increments", pen_inc_str, str(len(history_ipc["time"]))), + ("Total NR iterations", pen_nr_str, str(total_nr_ipc)), + ( + "Avg NR iter / increment", + pen_avg_str, + f"{total_nr_ipc / max(len(history_ipc['time']), 1):.1f}", + ), + ( + "Final reaction Fz", + pen_fz_str, + f"{history_ipc['reaction_z'][-1]:.2f}" if history_ipc["reaction_z"] else "N/A", + ), + ( + "Max contact force norm", + pen_fc_str, + f"{max(history_ipc['contact_force_norm']):.2f}" + if history_ipc["contact_force_norm"] + else "N/A", + ), + ( + "Final contact force norm", + pen_fc_final, + f"{history_ipc['contact_force_norm'][-1]:.2f}" + if history_ipc["contact_force_norm"] + else "N/A", + ), + ( + "IPC min gap distance", + "N/A", + f"{min(d for d in history_ipc['min_distance'] if np.isfinite(d)):.2e}" + if any(np.isfinite(d) for d in history_ipc["min_distance"]) + else "inf", + ), + ( + "IPC final barrier kappa", + "N/A", + f"{history_ipc['kappa'][-1]:.2e}" if history_ipc["kappa"] else "N/A", + ), + ( + "IPC max active collisions", + "N/A", + str(max(history_ipc["n_collisions"])) if history_ipc["n_collisions"] else "N/A", + ), +] + +w_label = 28 +w_val = 16 +print(f"{'Metric':<{w_label}} {'Penalty':>{w_val}} {'IPC':>{w_val}}") +print("-" * (w_label + 2 * w_val + 2)) +for label, v_pen, v_ipc in rows: + print(f"{label:<{w_label}} {v_pen:>{w_val}} {v_ipc:>{w_val}}") +print() + + +# ========================================================================= +# Per-increment detail tables +# ========================================================================= + +if history_penalty["time"]: + print("-" * 62) + print(" Per-increment detail: PENALTY") + print("-" * 62) + print( + f"{'Inc':>4} {'Time':>8} {'NR iter':>8} {'Reaction Fz':>14} {'|F_contact|':>14}" + ) + for i in range(len(history_penalty["time"])): + print( + f"{i+1:4d} {history_penalty['time'][i]:8.4f} " + f"{history_penalty['nr_iters'][i]:8d} " + f"{history_penalty['reaction_z'][i]:14.2f} " + f"{history_penalty['contact_force_norm'][i]:14.2f}" + ) + print() + +print("-" * 62) +print(" Per-increment detail: IPC") +print("-" * 62) +print( + f"{'Inc':>4} {'Time':>8} {'NR iter':>8} {'Reaction Fz':>14} " + f"{'|F_contact|':>14} {'Collisions':>11} {'Min gap':>12} {'Kappa':>12}" +) +for i in range(len(history_ipc["time"])): + min_d = history_ipc["min_distance"][i] + min_d_str = f"{min_d:.2e}" if np.isfinite(min_d) else "inf" + print( + f"{i+1:4d} {history_ipc['time'][i]:8.4f} " + f"{history_ipc['nr_iters'][i]:8d} " + f"{history_ipc['reaction_z'][i]:14.2f} " + f"{history_ipc['contact_force_norm'][i]:14.2f} " + f"{history_ipc['n_collisions'][i]:11d} " + f"{min_d_str:>12} " + f"{history_ipc['kappa'][i]:12.2e}" + ) + +print() + +# ========================================================================= +# Optional plots (requires matplotlib) +# ========================================================================= + +try: + import matplotlib.pyplot as plt + + fig, axes = plt.subplots(2, 2, figsize=(12, 8)) + fig.suptitle( + f"Gyroid Self-Contact -- IPC ({COMPRESSION*100:.0f}% compression)", + fontsize=14, + ) + + # Reaction force vs time + ax = axes[0, 0] + if history_penalty["time"]: + ax.plot( + history_penalty["time"], + history_penalty["reaction_z"], + "o-", + label="Penalty", + markersize=3, + ) + ax.plot( + history_ipc["time"], history_ipc["reaction_z"], "s-", label="IPC", markersize=3 + ) + ax.set_xlabel("Time") + ax.set_ylabel("Reaction Fz") + ax.set_title("Reaction Force (Z) at Top") + ax.legend() + ax.grid(True, alpha=0.3) + + # Contact force norm vs time + ax = axes[0, 1] + if history_penalty["time"]: + ax.plot( + history_penalty["time"], + history_penalty["contact_force_norm"], + "o-", + label="Penalty", + markersize=3, + ) + ax.plot( + history_ipc["time"], + history_ipc["contact_force_norm"], + "s-", + label="IPC", + markersize=3, + ) + ax.set_xlabel("Time") + ax.set_ylabel("||F_contact||") + ax.set_title("Contact Force Norm") + ax.legend() + ax.grid(True, alpha=0.3) + + # NR iterations per increment + ax = axes[1, 0] + if history_penalty["time"]: + ax.bar( + np.arange(len(history_penalty["nr_iters"])) - 0.2, + history_penalty["nr_iters"], + width=0.4, + label="Penalty", + ) + if history_ipc["time"]: + offset = 0.2 if history_penalty["time"] else 0.0 + ax.bar( + np.arange(len(history_ipc["nr_iters"])) + offset, + history_ipc["nr_iters"], + width=0.4, + label="IPC", + ) + ax.set_xlabel("Increment") + ax.set_ylabel("NR Iterations") + ax.set_title("Newton-Raphson Iterations per Increment") + ax.legend() + ax.grid(True, alpha=0.3) + + # IPC-specific: min gap distance + ax = axes[1, 1] + finite_gaps = [ + (t, d) + for t, d in zip(history_ipc["time"], history_ipc["min_distance"]) + if np.isfinite(d) + ] + if finite_gaps: + t_gap, d_gap = zip(*finite_gaps) + ax.plot(t_gap, d_gap, "s-", color="tab:green", markersize=3, label="Min gap") + ax.set_xlabel("Time") + ax.set_ylabel("Min Gap Distance") + ax.set_title("IPC Minimum Gap Distance") + ax.axhline(y=0, color="red", linestyle="--", alpha=0.5, label="Zero (penetration)") + ax.legend() + ax.grid(True, alpha=0.3) + + plt.tight_layout() + plt.savefig("results/gyroid_penalty_vs_ipc_comparison.png", dpi=150) + print("Comparison plot saved to results/gyroid_penalty_vs_ipc_comparison.png") + plt.show() + +except ImportError: + print("matplotlib not available -- skipping plots.") + + +# ========================================================================= +# Post-processing (requires pyvista) +# ========================================================================= +# Uncomment the lines below to visualise interactively. + +# if RUN_PENALTY: +# res_penalty.plot("Stress", "vm", "Node", show=True, scale=1) +# res_ipc.plot("Stress", "vm", "Node", show=True, scale=1) + +# --- Video output (uncomment to export MP4) --- +# res_ipc.write_movie("results/gyroid_ipc", "Stress", "vm", "Node", +# framerate=10, quality=8) +# print("Movie saved to results/gyroid_ipc.mp4") diff --git a/examples/contact/comparison/ring_crush.py b/examples/contact/comparison/ring_crush.py new file mode 100644 index 00000000..25d98522 --- /dev/null +++ b/examples/contact/comparison/ring_crush.py @@ -0,0 +1,408 @@ +""" +Hole plate self-contact: penalty vs IPC comparison +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +A plate with a large circular hole (radius 45 mm in a 100 x 100 mm plate) +is compressed vertically under plane strain until the hole walls come into +self-contact during buckling. + +Two contact methods are compared: + + - **Penalty method** (``fd.constraint.SelfContact``): node-to-surface + formulation with a user-tuned penalty parameter. + - **IPC method** (``fd.constraint.IPCSelfContact``): barrier-potential + formulation from the ipctk library guaranteeing intersection-free + configurations, with CCD line search. + +Both simulations use the same geometry, material and boundary conditions. +Per-increment reaction forces and IPC-specific metrics are tracked and +a comparison summary is printed. + +Expected behaviour +------------------ +Before self-contact occurs (roughly the first half of the loading), both +methods produce identical results. Once the hole walls come into contact: + +- The **penalty method** allows some interpenetration and the contact forces + may grow rapidly if the penalty parameter is too large relative to the + stiffness. +- The **IPC method** prevents all interpenetration (minimum gap > 0) but the + adaptive barrier stiffness (kappa) may grow aggressively, especially + when sustained contacts with decreasing gap are present. + +.. note:: + This example requires the ``simcoon`` package for the elasto-plastic + material and the ``ipctk`` package for the IPC method. +""" + +import fedoo as fd +import numpy as np +import os +from time import time + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + + +# ========================================================================= +# Shared geometry and material builder +# ========================================================================= + + +def build_mesh_and_material(): + """Build the hole plate geometry and elasto-plastic material. + + Returns fresh objects each call (new ModelingSpace). + """ + fd.ModelingSpace("2D") + + mesh = fd.mesh.hole_plate_mesh( + nr=15, + nt=15, + length=100, + height=100, + radius=45, + ) + + # Material: elasto-plastic (EPICP) — same as tube_compression example + E, nu = 200e3, 0.3 + sigma_y = 300 + k, m = 1000, 0.3 + props = np.array([E, nu, 1e-5, sigma_y, k, m]) + material = fd.constitutivelaw.Simcoon("EPICP", props) + + return mesh, material + + +def get_bc_nodes(mesh): + """Return (bottom, top) node arrays.""" + nodes_bot = mesh.find_nodes("Y", mesh.bounding_box.ymin) + nodes_top = mesh.find_nodes("Y", mesh.bounding_box.ymax) + return nodes_bot, nodes_top + + +# ========================================================================= +# Approach 1 : Penalty self-contact +# ========================================================================= + +print("=" * 62) +print(" PENALTY SELF-CONTACT (hole plate crushing)") +print("=" * 62) + +mesh, material = build_mesh_and_material() +nodes_bot, nodes_top = get_bc_nodes(mesh) + +surf = fd.mesh.extract_surface(mesh) +penalty_contact = fd.constraint.SelfContact( + surf, + "linear", + search_algorithm="bucket", +) +penalty_contact.contact_search_once = True +penalty_contact.eps_n = 1e6 +penalty_contact.max_dist = 1.5 + +wf = fd.weakform.StressEquilibrium(material, nlgeom="UL") +solid_assembly = fd.Assembly.create(wf, mesh) +assembly = fd.Assembly.sum(solid_assembly, penalty_contact) + +pb_penalty = fd.problem.NonLinear(assembly) + +if not os.path.isdir("results"): + os.mkdir("results") +res_penalty = pb_penalty.add_output( + "results/ring_crush_penalty", + solid_assembly, + ["Disp", "Stress"], +) + +pb_penalty.bc.add("Dirichlet", nodes_bot, "Disp", 0) +pb_penalty.bc.add("Dirichlet", nodes_top, "Disp", [0, -40]) +pb_penalty.set_nr_criterion("Displacement", tol=5e-3, max_subiter=5) + +# --- Tracking callback --- +history_penalty = {"time": [], "reaction_y": []} + + +def track_penalty(pb): + history_penalty["time"].append(pb.time) + F = pb.get_ext_forces("Disp") + history_penalty["reaction_y"].append(np.sum(F[1, nodes_top])) + + +t0_penalty = time() +pb_penalty.nlsolve( + dt=0.01, + tmax=1, + update_dt=True, + print_info=1, + callback=track_penalty, + exec_callback_at_each_iter=True, +) +penalty_time = time() - t0_penalty +print(f"\nPenalty solve time: {penalty_time:.2f} s") + + +# ========================================================================= +# Approach 2 : IPC contact +# ========================================================================= + +print("\n" + "=" * 62) +print(" IPC CONTACT (hole plate crushing)") +print("=" * 62) + +mesh2, material2 = build_mesh_and_material() +nodes_bot2, nodes_top2 = get_bc_nodes(mesh2) + +ipc_contact = fd.constraint.IPCSelfContact( + mesh2, + dhat=5e-3, + dhat_is_relative=True, + friction_coefficient=0.0, + use_ccd=True, +) + +wf2 = fd.weakform.StressEquilibrium(material2, nlgeom="UL") +solid_assembly2 = fd.Assembly.create(wf2, mesh2) +assembly2 = fd.Assembly.sum(solid_assembly2, ipc_contact) + +pb_ipc = fd.problem.NonLinear(assembly2) + +res_ipc = pb_ipc.add_output( + "results/ring_crush_ipc", + solid_assembly2, + ["Disp", "Stress"], +) + +pb_ipc.bc.add("Dirichlet", nodes_bot2, "Disp", 0) +pb_ipc.bc.add("Dirichlet", nodes_top2, "Disp", [0, -40]) +pb_ipc.set_nr_criterion("Displacement", tol=5e-3, max_subiter=5) + +# --- Tracking callback --- +history_ipc = { + "time": [], + "reaction_y": [], + "n_collisions": [], + "kappa": [], + "min_distance": [], +} + + +def track_ipc(pb): + history_ipc["time"].append(pb.time) + F = pb.get_ext_forces("Disp") + history_ipc["reaction_y"].append(np.sum(F[1, nodes_top2])) + history_ipc["n_collisions"].append(len(ipc_contact._collisions)) + history_ipc["kappa"].append(ipc_contact._kappa) + verts = ipc_contact._get_current_vertices(pb) + if len(ipc_contact._collisions) > 0: + min_d = ipc_contact._collisions.compute_minimum_distance( + ipc_contact._collision_mesh, + verts, + ) + else: + min_d = float("inf") + history_ipc["min_distance"].append(min_d) + + +t0_ipc = time() +pb_ipc.nlsolve( + dt=0.01, + tmax=1, + update_dt=True, + print_info=1, + callback=track_ipc, + exec_callback_at_each_iter=True, +) +ipc_time = time() - t0_ipc +print(f"\nIPC solve time: {ipc_time:.2f} s") + + +# ========================================================================= +# Comparison summary +# ========================================================================= + +print("\n") +print("=" * 62) +print(" PERFORMANCE COMPARISON: Penalty vs IPC (hole plate)") +print("=" * 62) + +n_inc_pen = len(history_penalty["time"]) +n_inc_ipc = len(history_ipc["time"]) + +rows = [ + ("Total solve time", f"{penalty_time:.2f} s", f"{ipc_time:.2f} s"), + ("Total increments", str(n_inc_pen), str(n_inc_ipc)), + ( + "Final reaction Fy (top)", + f"{history_penalty['reaction_y'][-1]:.1f}" + if history_penalty["reaction_y"] + else "N/A", + f"{history_ipc['reaction_y'][-1]:.1f}" if history_ipc["reaction_y"] else "N/A", + ), + ( + "IPC min gap distance", + "N/A", + f"{min(d for d in history_ipc['min_distance'] if np.isfinite(d)):.2e}" + if any(np.isfinite(d) for d in history_ipc["min_distance"]) + else "inf", + ), + ( + "IPC final barrier kappa", + "N/A", + f"{history_ipc['kappa'][-1]:.2e}" if history_ipc["kappa"] else "N/A", + ), + ( + "IPC max active collisions", + "N/A", + str(max(history_ipc["n_collisions"])) if history_ipc["n_collisions"] else "N/A", + ), +] + +w_label, w_val = 28, 16 +print(f"{'Metric':<{w_label}} {'Penalty':>{w_val}} {'IPC':>{w_val}}") +print("-" * (w_label + 2 * w_val + 2)) +for label, v_pen, v_ipc in rows: + print(f"{label:<{w_label}} {v_pen:>{w_val}} {v_ipc:>{w_val}}") +print() + + +# ========================================================================= +# Per-increment detail tables +# ========================================================================= + +print("-" * 62) +print(" Per-increment detail: PENALTY") +print("-" * 62) +print(f"{'Inc':>4} {'Time':>8} {'Reaction Fy':>14}") +for i in range(n_inc_pen): + print( + f"{i+1:4d} {history_penalty['time'][i]:8.4f} " + f"{history_penalty['reaction_y'][i]:14.2f}" + ) + +print() +print("-" * 62) +print(" Per-increment detail: IPC") +print("-" * 62) +print( + f"{'Inc':>4} {'Time':>8} {'Reaction Fy':>14} " + f"{'Collisions':>11} {'Min gap':>12} {'Kappa':>12}" +) +for i in range(n_inc_ipc): + min_d = history_ipc["min_distance"][i] + min_d_str = f"{min_d:.2e}" if np.isfinite(min_d) else "inf" + print( + f"{i+1:4d} {history_ipc['time'][i]:8.4f} " + f"{history_ipc['reaction_y'][i]:14.2f} " + f"{history_ipc['n_collisions'][i]:11d} " + f"{min_d_str:>12} " + f"{history_ipc['kappa'][i]:12.2e}" + ) +print() + + +# ========================================================================= +# Optional plots (requires matplotlib) +# ========================================================================= + +try: + import matplotlib.pyplot as plt + + fig, axes = plt.subplots(2, 2, figsize=(12, 8)) + fig.suptitle("Hole Plate Crushing: Penalty vs IPC Comparison", fontsize=14) + + # --- Force-displacement --- + disp_pen = [t * 40 for t in history_penalty["time"]] + disp_ipc = [t * 40 for t in history_ipc["time"]] + + ax = axes[0, 0] + ax.plot( + disp_pen, + [-f for f in history_penalty["reaction_y"]], + "o-", + label="Penalty", + markersize=3, + ) + ax.plot( + disp_ipc, + [-f for f in history_ipc["reaction_y"]], + "s-", + label="IPC", + markersize=3, + ) + ax.set_xlabel("Top displacement (mm)") + ax.set_ylabel("Reaction force Fy (N/mm)") + ax.set_title("Force-Displacement Curve") + ax.legend() + ax.grid(True, alpha=0.3) + + # --- IPC collisions --- + ax = axes[0, 1] + ax.plot( + history_ipc["time"], + history_ipc["n_collisions"], + "s-", + color="tab:red", + markersize=3, + ) + ax.set_xlabel("Time") + ax.set_ylabel("Active collisions") + ax.set_title("IPC Active Collision Count") + ax.grid(True, alpha=0.3) + + # --- IPC kappa evolution --- + ax = axes[1, 0] + ax.plot( + history_ipc["time"], + history_ipc["kappa"], + "s-", + color="tab:orange", + markersize=3, + ) + ax.set_xlabel("Time") + ax.set_ylabel("Barrier stiffness kappa") + ax.set_title("IPC Barrier Stiffness") + ax.set_yscale("log") + ax.grid(True, alpha=0.3) + + # --- IPC min gap --- + ax = axes[1, 1] + finite_gaps = [ + (t, d) + for t, d in zip(history_ipc["time"], history_ipc["min_distance"]) + if np.isfinite(d) + ] + if finite_gaps: + t_gap, d_gap = zip(*finite_gaps) + ax.plot(t_gap, d_gap, "s-", color="tab:green", markersize=3, label="Min gap") + if ipc_contact._actual_dhat is not None: + ax.axhline( + y=ipc_contact._actual_dhat, + color="blue", + linestyle="--", + alpha=0.5, + label="dhat", + ) + ax.axhline(y=0, color="red", linestyle="--", alpha=0.5, label="Zero (penetration)") + ax.set_xlabel("Time") + ax.set_ylabel("Min gap distance") + ax.set_title("IPC Minimum Gap Distance") + ax.legend() + ax.grid(True, alpha=0.3) + + plt.tight_layout() + plt.savefig("results/ring_crush_comparison.png", dpi=150) + print("Comparison plot saved to results/ring_crush_comparison.png") + plt.show() + +except ImportError: + print("matplotlib not available -- skipping plots.") + + +# ========================================================================= +# Final deformed shape with von Mises stress (requires pyvista) +# ========================================================================= +# Uncomment the lines below to visualise the final deformed configuration. + +# res_penalty.plot("Stress", "vm", "Node", show=True, scale=1, show_nodes=True) +# res_ipc.plot("Stress", "vm", "Node", show=True, scale=1, show_nodes=True) diff --git a/examples/contact/comparison/self_contact.py b/examples/contact/comparison/self_contact.py new file mode 100644 index 00000000..f10aa671 --- /dev/null +++ b/examples/contact/comparison/self_contact.py @@ -0,0 +1,135 @@ +""" +Self-contact: penalty vs IPC comparison +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This example compares the two self-contact approaches available in fedoo +on a 2D hole plate being crushed: + + - **Penalty method** (``fd.constraint.SelfContact``): node-to-surface + formulation with a user-tuned penalty parameter. + - **IPC method** (``fd.constraint.IPCSelfContact``): barrier-potential + formulation from the ipctk library guaranteeing intersection-free + configurations. + +Both simulations use the same geometry (a plate with a large hole) +and boundary conditions (vertical compression until self-contact occurs). + +.. note:: + This example requires the ``simcoon`` package for the elasto-plastic + material and the ``ipctk`` package for the IPC method. +""" + +import fedoo as fd +import numpy as np +import os +from time import time + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + + +def build_mesh_and_material(): + """Build the shared geometry and material (fresh ModelingSpace each call).""" + fd.ModelingSpace("2D") + + # Plate with a large circular hole -- will self-contact under compression + mesh = fd.mesh.hole_plate_mesh(nr=15, nt=15, length=100, height=100, radius=45) + + # Material: elasto-plastic (EPICP) + E, nu = 200e3, 0.3 + alpha = 1e-5 + Re = 300 # yield stress + k = 1000 # hardening modulus + m = 0.3 # hardening exponent + props = np.array([E, nu, alpha, Re, k, m]) + material = fd.constitutivelaw.Simcoon("EPICP", props) + + return mesh, material + + +# ========================================================================= +# Approach 1 : Penalty self-contact +# ========================================================================= + +print("=" * 60) +print("PENALTY SELF-CONTACT") +print("=" * 60) + +mesh, material = build_mesh_and_material() + +nodes_top = mesh.find_nodes("Y", mesh.bounding_box.ymax) +nodes_bottom = mesh.find_nodes("Y", mesh.bounding_box.ymin) + +surf = fd.mesh.extract_surface(mesh) +penalty_contact = fd.constraint.SelfContact(surf, "linear", search_algorithm="bucket") +penalty_contact.contact_search_once = True +penalty_contact.eps_n = 1e6 + +wf = fd.weakform.StressEquilibrium(material, nlgeom="UL") +solid_assembly = fd.Assembly.create(wf, mesh) +assembly = fd.Assembly.sum(solid_assembly, penalty_contact) + +pb_penalty = fd.problem.NonLinear(assembly) + +if not os.path.isdir("results"): + os.mkdir("results") +res_penalty = pb_penalty.add_output( + "results/self_contact_penalty", solid_assembly, ["Disp", "Stress"] +) + +pb_penalty.bc.add("Dirichlet", nodes_bottom, "Disp", 0) +pb_penalty.bc.add("Dirichlet", nodes_top, "Disp", [0, -70]) +pb_penalty.set_nr_criterion("Displacement", tol=5e-3, max_subiter=5) + +t0 = time() +pb_penalty.nlsolve(dt=0.01, tmax=1, update_dt=True, print_info=1, interval_output=0.1) +print(f"Penalty self-contact solve time: {time() - t0:.2f} s") + + +# ========================================================================= +# Approach 2 : IPC self-contact +# ========================================================================= + +print("\n" + "=" * 60) +print("IPC SELF-CONTACT") +print("=" * 60) + +mesh2, material2 = build_mesh_and_material() + +nodes_top2 = mesh2.find_nodes("Y", mesh2.bounding_box.ymax) +nodes_bottom2 = mesh2.find_nodes("Y", mesh2.bounding_box.ymin) + +# IPC self-contact: auto-extracts surface from the mesh +ipc_contact = fd.constraint.IPCSelfContact( + mesh2, + dhat=5e-3, + dhat_is_relative=True, + friction_coefficient=0.0, + use_ccd=True, +) + +wf2 = fd.weakform.StressEquilibrium(material2, nlgeom="UL") +solid_assembly2 = fd.Assembly.create(wf2, mesh2) +assembly2 = fd.Assembly.sum(solid_assembly2, ipc_contact) + +pb_ipc = fd.problem.NonLinear(assembly2) + +res_ipc = pb_ipc.add_output( + "results/self_contact_ipc", solid_assembly2, ["Disp", "Stress"] +) + +pb_ipc.bc.add("Dirichlet", nodes_bottom2, "Disp", 0) +pb_ipc.bc.add("Dirichlet", nodes_top2, "Disp", [0, -70]) +pb_ipc.set_nr_criterion("Displacement", tol=5e-3, max_subiter=5) + +t0 = time() +pb_ipc.nlsolve(dt=0.01, tmax=1, update_dt=True, print_info=1, interval_output=0.1) +print(f"IPC self-contact solve time: {time() - t0:.2f} s") + + +# ========================================================================= +# Post-processing (requires pyvista) +# ========================================================================= +# Uncomment the lines below to visualise and compare the results. + +# res_penalty.plot("Stress", "vm", "Node", show=True, scale=1, show_nodes=True) +# res_ipc.plot("Stress", "vm", "Node", show=True, scale=1, show_nodes=True) diff --git a/examples/contact/disk_rectangle_contact.py b/examples/contact/disk_rectangle_contact.py deleted file mode 100644 index 1672007a..00000000 --- a/examples/contact/disk_rectangle_contact.py +++ /dev/null @@ -1,169 +0,0 @@ -import fedoo as fd -import numpy as np -from time import time -import os -import pylab as plt -from numpy import linalg - -import pyvista as pv - - -start = time() -# --------------- Pre-Treatment -------------------------------------------------------- - -fd.ModelingSpace("2D") - -NLGEOM = "UL" -# Units: N, mm, MPa -h = 1 -L = 1 -E = 200e3 -nu = 0.3 -alpha = 1e-5 - -filename = "disk_rectangle_contact" -res_dir = "results/" - -mesh_rect = fd.mesh.rectangle_mesh( - nx=11, ny=21, x_min=0, x_max=L, y_min=0, y_max=h, elm_type="quad4", name="Domain" -) -mesh_rect.element_sets["rect"] = np.arange(0, mesh_rect.n_elements) -mesh_disk = fd.mesh.disk_mesh(radius=L / 2, nr=6, nt=6, elm_type="quad4") -# surf = fd.constraint.contact.extract_surface_mesh(fd.Mesh.from_pyvista(pv.Circle(radius=0.5, resolution = 30).triangulate())) -mesh_disk.nodes += np.array([1.5, 0.48]) -mesh_disk.element_sets["disk"] = np.arange(0, mesh_disk.n_elements) - -mesh = fd.Mesh.stack(mesh_rect, mesh_disk) - -# node sets for boundary conditions -nodes_left = mesh.find_nodes("X", 0) -nodes_right = mesh.find_nodes("X", L) - -# nodes_top = mesh.find_nodes('Y',1) -nodes_bc = mesh.find_nodes("X>1.5") -nodes_bc = list(set(nodes_bc).intersection(mesh.node_sets["boundary"])) - -# if slave surface = disk -# nodes_contact = mesh.node_sets['boundary'] -# surf = fd.mesh.extract_surface(mesh.extract_elements('rect')) #extract the surface of the rectangle -# surf = surf.extract_elements(surf.get_elements_from_nodes(nodes_right)) - -# if slave surface = rectangle -nodes_contact = nodes_right -surf = fd.mesh.extract_surface( - mesh.extract_elements("disk") -) # extract the surface of the disk - -contact = fd.constraint.Contact(nodes_contact, surf) -contact.contact_search_once = True -contact.eps_n = 5e5 -contact.max_dist = 1 - -mat = 1 -if mat == 0: - props = np.array([E, nu, alpha]) - material_rect = fd.constitutivelaw.Simcoon("ELISO", props, name="ConstitutiveLaw") -elif mat == 1 or mat == 2: - Re = 300 - k = 1000 # 1500 - m = 0.3 # 0.25 - if mat == 1: - props = np.array([E, nu, alpha, Re, k, m]) - material_rect = fd.constitutivelaw.Simcoon( - "EPICP", props, name="ConstitutiveLaw" - ) - elif mat == 2: - material_rect = fd.constitutivelaw.ElastoPlasticity( - E, nu, Re, name="ConstitutiveLaw" - ) - material_rect.SetHardeningFunction("power", H=k, beta=m) -else: - material_rect = fd.constitutivelaw.ElasticIsotrop(E, nu, name="ConstitutiveLaw") - -material_disk = fd.constitutivelaw.ElasticIsotrop(50e3, nu, name="ConstitutiveLaw") - -material = fd.constitutivelaw.Heterogeneous( - (material_rect, material_disk), ("rect", "disk") -) - -wf = fd.weakform.StressEquilibrium(material, nlgeom=NLGEOM) -solid_assembly = fd.Assembly.create(wf, mesh) - -# assembly = fd.Assembly.sum(solid_assembly1, solid_assembly2, contact) -assembly = fd.Assembly.sum(solid_assembly, contact) - -pb = fd.problem.NonLinear(assembly) - -# create a 'result' folder and set the desired ouputs -if not (os.path.isdir("results")): - os.mkdir("results") -# results = pb.add_output(res_dir+filename, 'Assembling', ['Disp'], output_type='Node', file_format ='npz') -# results = pb.add_output(res_dir+filename, 'Assembling', ['Cauchy', 'PKII', 'Strain', 'Cauchy_vm', 'Statev', 'Wm'], output_type='GaussPoint', file_format ='npz') - -# results = pb.add_output(res_dir+filename, solid_assembly, ['Disp', 'Cauchy', 'PKII', 'Strain', 'Cauchy_vm', 'Statev', 'Wm']) -results = pb.add_output( - res_dir + filename, solid_assembly, ["Disp", "Stress", "Strain", "Fext"] -) - -# Problem.add_output(res_dir+filename, 'Assembling', ['cauchy', 'PKII', 'strain', 'cauchy_vm', 'statev'], output_type='Element', file_format ='vtk') - -pb.bc.add("Dirichlet", nodes_left, "Disp", 0) -pb.bc.add("Dirichlet", nodes_bc, "Disp", [-0.4, 0.2]) - -# Problem.set_solver('cg', precond = True) -# pb.set_nr_criterion("Force", err0 = None, tol = 5e-3, max_subiter = 10) -pb.set_nr_criterion("Displacement", err0=None, tol=5e-3, max_subiter=5) - -# pb.nlsolve(dt = 0.001, tmax = 1, update_dt = False, print_info = 2, interval_output = 0.001) -pb.nlsolve(dt=0.005, tmax=1, update_dt=True, print_info=1, interval_output=0.01) - - -pb.bc.remove(-1) # remove last boundary contidion -pb.bc.add("Dirichlet", nodes_bc, "Disp", [0, 0]) - -pb.nlsolve(dt=0.005, tmax=1, update_dt=True, print_info=1, interval_output=0.01) - -print(time() - start) - - -# ============================================================= -# Example of plots with pyvista - uncomment the desired plot -# ============================================================= - -# ------------------------------------ -# Simple plot with default options -# ------------------------------------ -results.plot("Stress", "vm", "Node", show=True, scale=1, show_nodes=True) -results.plot("Stress", "XX", "Node", show=True, scale=1, show_nodes=True) -# results.plot('Fext', 'X', 'Node', show = True, scale = 1, show_nodes=True) - -results.plot("Disp", 0, "Node", show=True, scale=1, show_nodes=True) - -# ------------------------------------ -# Write movie with default options -# ------------------------------------ -# results.write_movie(res_dir+filename, 'Stress', 'vm', framerate = 5, quality = 5) - -results.write_movie( - res_dir + filename, - "Stress", - "XX", - "Node", - framerate=24, - quality=5, - clim=[-3e3, 3e3], -) - -# ------------------------------------ -# Save pdf plot -# ------------------------------------ -# pl = results.plot('Stress', 'vm', show = False) -# pl.save_graphic('test.pdf', title='PyVista Export', raster=True, painter=True) - -# ------------------------------------ -# Plot time history -# ------------------------------------ -# from matplotlib import pylab -# t, sigma = results.get_history(('Time','Stress'), (0,12), component = 3) -# pylab.plot(t,sigma) -# #or results.plot_history('Stress', 12) diff --git a/examples/contact/ipc/gyroid_compression.py b/examples/contact/ipc/gyroid_compression.py new file mode 100644 index 00000000..669a5725 --- /dev/null +++ b/examples/contact/ipc/gyroid_compression.py @@ -0,0 +1,61 @@ +""" +Gyroid self-contact under compression (3D, IPC method) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Compresses a 3D gyroid unit cell with IPC self-contact to prevent +wall intersections. 50% compression with Augmented Lagrangian +outer loop (``al_max_iter=5``) for robust convergence at high strain. + +.. note:: + Requires ``ipctk``. +""" + +import fedoo as fd +import numpy as np +import os + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + +fd.ModelingSpace("3D") + +# --- Geometry --- +MESH_PATH = os.path.join( + os.path.dirname(os.path.abspath(__file__)), "../../../util/meshes/gyroid_per.vtk" +) +mesh = fd.Mesh.read(MESH_PATH) +material = fd.constitutivelaw.ElasticIsotrop(1e5, 0.3) + +# --- IPC self-contact --- +contact = fd.constraint.IPCSelfContact( + mesh, + dhat=1.0e-2, + dhat_is_relative=True, + use_ccd=True, +) + +wf = fd.weakform.StressEquilibrium(material, nlgeom="UL") +solid_assembly = fd.Assembly.create(wf, mesh) +assembly = fd.Assembly.sum(solid_assembly, contact) + +pb = fd.problem.NonLinear(assembly) + +if not os.path.isdir("results"): + os.mkdir("results") +res = pb.add_output("results/gyroid_ipc", solid_assembly, ["Disp", "Stress"]) + +# --- BCs: compression 50% --- +nodes_top = mesh.find_nodes("Z", mesh.bounding_box.zmax) +nodes_bottom = mesh.find_nodes("Z", mesh.bounding_box.zmin) +pb.bc.add("Dirichlet", nodes_bottom, "Disp", 0) +pb.bc.add("Dirichlet", nodes_top, "Disp", [0, 0, -0.5]) +pb.set_nr_criterion("Displacement", tol=5.0e-3, max_subiter=25) + +pb.nlsolve(dt=0.05, tmax=1, update_dt=True, print_info=1, interval_output=0.05) + +# --- Static plot --- +res.plot("Stress", "vm", "Node", show=False, scale=1) + +# --- Video output (uncomment to export MP4) --- +# res.write_movie("results/gyroid_ipc", "Stress", "vm", "Node", +# framerate=10, quality=8) +# print("Movie saved to results/gyroid_ipc.mp4") diff --git a/examples/contact/ipc/gyroid_compression_periodic.py b/examples/contact/ipc/gyroid_compression_periodic.py new file mode 100644 index 00000000..2b8c4996 --- /dev/null +++ b/examples/contact/ipc/gyroid_compression_periodic.py @@ -0,0 +1,91 @@ +""" +Gyroid compression with periodic BCs (3D, IPC method) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Compresses a 3D gyroid unit cell (50% strain) with IPC self-contact +and periodic boundary conditions. This is the most advanced contact +example, combining large strain (UL), periodic BCs, and IPC. + +A diagnostic callback prints the number of active collisions, the +barrier stiffness, and the contact force norm at each converged step. + +.. note:: + Requires ``ipctk``. +""" + +import fedoo as fd +import numpy as np +import os + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + +fd.ModelingSpace("3D") + +# --- Geometry --- +MESH_PATH = os.path.join( + os.path.dirname(os.path.abspath(__file__)), "../../../util/meshes/gyroid_per.vtk" +) +mesh = fd.Mesh.read(MESH_PATH) +material = fd.constitutivelaw.ElasticIsotrop(1e5, 0.3) + +# --- IPC self-contact (auto-extracts surface) --- +contact = fd.constraint.IPCSelfContact( + mesh, dhat=5e-3, dhat_is_relative=True, use_ccd=True +) + +# --- Assembly --- +wf = fd.weakform.StressEquilibrium(material, nlgeom="UL") +solid_assembly = fd.Assembly.create(wf, mesh) +assembly = fd.Assembly.sum(solid_assembly, contact) + +# --- Problem --- +pb = fd.problem.NonLinear(assembly) +if not os.path.isdir("results"): + os.mkdir("results") +res = pb.add_output("results/gyroid_ipc_periodic", solid_assembly, ["Disp", "Stress"]) + +# --- BCs: periodic, compress 50% --- +bc_periodic = fd.constraint.PeriodicBC("finite_strain", tol=1e-3) +pb.bc.add(bc_periodic) + +# block a node near the center to avoid rigid body motion +pb.bc.add("Dirichlet", mesh.nearest_node(mesh.bounding_box.center), "Disp", 0) + +# Uniaxial compression: prescribe DU_xx, fix off-diagonal terms to prevent rotation +pb.bc.add("Dirichlet", "DU_xx", -0.5) +pb.bc.add("Dirichlet", "DU_xy", 0) +pb.bc.add("Dirichlet", "DU_xz", 0) +pb.bc.add("Dirichlet", "DU_yx", 0) +pb.bc.add("Dirichlet", "DU_yz", 0) +pb.bc.add("Dirichlet", "DU_zx", 0) +pb.bc.add("Dirichlet", "DU_zy", 0) + +pb.set_nr_criterion("Displacement", tol=5e-3, max_subiter=10) + + +# --- Diagnostic callback --- +def diag(pb): + n_coll = len(contact._collisions) if contact._collisions is not None else 0 + kappa = getattr(contact, "_kappa", None) + gv_norm = ( + np.linalg.norm(contact.global_vector) + if contact.global_vector is not None + else 0 + ) + print( + f" [IPC] t={pb.time:.4f} collisions={n_coll} kappa={kappa} |Fcontact|={gv_norm:.4e}" + ) + + +# --- Solve --- +pb.nlsolve( + dt=0.05, tmax=1, update_dt=True, print_info=1, interval_output=0.1, callback=diag +) + +# --- Static plot --- +# res.plot("Stress", "vm", "Node", show=True, scale=1) + +# --- Video output (uncomment to export MP4) --- +# res.write_movie("results/gyroid_ipc_periodic", "Stress", "vm", "Node", +# framerate=10, quality=8) +# print("Movie saved to results/gyroid_ipc_periodic.mp4") diff --git a/examples/contact/ipc/gyroid_compression_plastic.py b/examples/contact/ipc/gyroid_compression_plastic.py new file mode 100644 index 00000000..038cbdae --- /dev/null +++ b/examples/contact/ipc/gyroid_compression_plastic.py @@ -0,0 +1,64 @@ +""" +Gyroid self-contact under compression (3D, IPC method) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Compresses a 3D gyroid unit cell with IPC self-contact to prevent +wall intersections. 50% compression with Augmented Lagrangian +outer loop (``al_max_iter=5``) for robust convergence at high strain. + +.. note:: + Requires ``ipctk``. +""" + +import fedoo as fd +import numpy as np +import os + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + +fd.ModelingSpace("3D") + +# --- Geometry --- +MESH_PATH = os.path.join( + os.path.dirname(os.path.abspath(__file__)), "../../../util/meshes/gyroid_per.vtk" +) +mesh = fd.Mesh.read(MESH_PATH) + +E, nu = 200e3, 0.3 +props = np.array([E, nu, 1e-5, 300, 1000, 0.3]) +material = fd.constitutivelaw.Simcoon("EPICP", props) + +# --- IPC self-contact --- +contact = fd.constraint.IPCSelfContact( + mesh, + dhat=1.0e-2, + dhat_is_relative=True, + use_ccd=True, +) + +wf = fd.weakform.StressEquilibrium(material, nlgeom="UL") +solid_assembly = fd.Assembly.create(wf, mesh) +assembly = fd.Assembly.sum(solid_assembly, contact) + +pb = fd.problem.NonLinear(assembly) + +if not os.path.isdir("results"): + os.mkdir("results") +res = pb.add_output("results/gyroid_ipc_plastic", solid_assembly, ["Disp", "Stress"]) + +# --- BCs: compression 50% --- +nodes_top = mesh.find_nodes("Z", mesh.bounding_box.zmax) +nodes_bottom = mesh.find_nodes("Z", mesh.bounding_box.zmin) +pb.bc.add("Dirichlet", nodes_bottom, "Disp", 0) +pb.bc.add("Dirichlet", nodes_top, "Disp", [0, 0, -0.5]) +pb.set_nr_criterion("Displacement", tol=5.0e-3, max_subiter=25) + +pb.nlsolve(dt=0.05, tmax=1, update_dt=True, print_info=1, interval_output=0.05) + +# --- Static plot --- +res.plot("Stress", "vm", "Node", show=False, scale=1) + +# --- Video output (uncomment to export MP4) --- +# res.write_movie("results/gyroid_ipc", "Stress", "vm", "Node", +# framerate=10, quality=8) +# print("Movie saved to results/gyroid_ipc.mp4") diff --git a/examples/contact/ipc/indentation_2d.py b/examples/contact/ipc/indentation_2d.py new file mode 100644 index 00000000..7082321b --- /dev/null +++ b/examples/contact/ipc/indentation_2d.py @@ -0,0 +1,240 @@ +""" +Disk indentation (2D plane strain) — IPC contact with Hertz comparison +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +A stiff elastic disk is pressed vertically into a softer elastic rectangle +(plane strain). This is a classical Hertz-like contact benchmark. + +The plate is meshed with gmsh, using a graded size field that concentrates +elements near the contact zone. The disk uses fedoo's built-in disk_mesh. + +The simulation uses the IPC barrier method for contact and compares the +computed force-indentation curve to the 2D Hertz analytical prediction +(Johnson, *Contact Mechanics*, Ch. 4): + +* Contact half-width: ``a = sqrt(4*R*P / (pi*E*))`` +* Indentation depth: ``delta = a^2/(4*R) * (2*ln(4*R/a) - 1)`` +* Force per unit length: ``P = pi*E*/(4*R) * a^2`` + +The FEM force is expected to be lower than the Hertz prediction because +the 2D (plane-strain) Hertz solution assumes an infinite half-space, but +the finite plate is significantly more compliant. This discrepancy is +inherent to 2D contact problems, where the half-space displacement under a +line load has a logarithmic divergence. + +.. note:: + Requires ``ipctk`` and ``gmsh``. +""" + +import fedoo as fd +import numpy as np +import os +import tempfile + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + +# ========================================================================= +# Parameters +# ========================================================================= + +fd.ModelingSpace("2D") + +# Units: N, mm, MPa +E_plate = 1e3 # soft plate +E_disk = 1e5 # stiff disk (quasi-rigid, 100x stiffer) +nu = 0.3 +R = 5.0 # disk radius +plate_half = 30.0 # half-width of the plate +plate_h = 40.0 # plate height +gap = 0.1 # initial gap between disk bottom and plate top +imposed_disp = -2.0 # total vertical displacement of disk top + +# ========================================================================= +# Plate mesh (gmsh -- graded refinement near contact) +# ========================================================================= + +import gmsh + +gmsh.initialize() +gmsh.option.setNumber("General.Verbosity", 1) + +plate_tag = gmsh.model.occ.addRectangle( + -plate_half, + 0, + 0, + 2 * plate_half, + plate_h, +) +gmsh.model.occ.synchronize() +gmsh.model.addPhysicalGroup(2, [plate_tag], tag=1, name="plate") + +# Fine elements near the top-centre (contact zone), coarse far away +gmsh.model.mesh.field.add("Box", 1) +gmsh.model.mesh.field.setNumber(1, "VIn", 0.3) +gmsh.model.mesh.field.setNumber(1, "VOut", 3.0) +gmsh.model.mesh.field.setNumber(1, "XMin", -6) +gmsh.model.mesh.field.setNumber(1, "XMax", 6) +gmsh.model.mesh.field.setNumber(1, "YMin", plate_h - 4) +gmsh.model.mesh.field.setNumber(1, "YMax", plate_h) +gmsh.model.mesh.field.setNumber(1, "Thickness", 4) +gmsh.model.mesh.field.setAsBackgroundMesh(1) +gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) + +gmsh.model.mesh.generate(2) +plate_msh = os.path.join(tempfile.gettempdir(), "plate_indent2d.msh") +gmsh.write(plate_msh) +gmsh.finalize() + +mesh_plate = fd.mesh.import_msh(plate_msh, mesh_type="surface") +# gmsh always stores 3D coordinates; strip Z for 2D modelling space +if mesh_plate.nodes.shape[1] == 3: + mesh_plate.nodes = mesh_plate.nodes[:, :2] +mesh_plate.element_sets["plate"] = np.arange(mesh_plate.n_elements) +print(f"Plate mesh: {mesh_plate.n_nodes} nodes, {mesh_plate.n_elements} elems") + +# ========================================================================= +# Disk mesh (fedoo built-in -- structured quad4 -> same type as plate) +# ========================================================================= + +# Use tri3 for the disk to be compatible with gmsh plate mesh +mesh_disk = fd.mesh.disk_mesh(radius=R, nr=12, nt=24, elm_type="tri3") +# Position: disk center at (0, plate_top + R + gap) +mesh_disk.nodes += np.array([0, plate_h + R + gap]) +mesh_disk.element_sets["disk"] = np.arange(mesh_disk.n_elements) +print(f"Disk mesh: {mesh_disk.n_nodes} nodes, {mesh_disk.n_elements} elems") + +# Stack into single mesh +mesh = fd.Mesh.stack(mesh_plate, mesh_disk) +print( + f"Total mesh: {mesh.n_nodes} nodes, {mesh.n_elements} elems, type={mesh.elm_type}" +) + +# ========================================================================= +# IPC contact +# ========================================================================= + +surf = fd.mesh.extract_surface(mesh) +ipc_contact = fd.constraint.IPCContact( + mesh, + surface_mesh=surf, + dhat=0.05, # absolute dhat (< gap) + dhat_is_relative=False, + use_ccd=True, +) + +# ========================================================================= +# Material, assembly, BCs +# ========================================================================= + +mat_plate = fd.constitutivelaw.ElasticIsotrop(E_plate, nu) +mat_disk = fd.constitutivelaw.ElasticIsotrop(E_disk, nu) +material = fd.constitutivelaw.Heterogeneous( + (mat_plate, mat_disk), + ("plate", "disk"), +) + +wf = fd.weakform.StressEquilibrium(material, nlgeom=False) +solid = fd.Assembly.create(wf, mesh) +assembly = fd.Assembly.sum(solid, ipc_contact) + +pb = fd.problem.NonLinear(assembly) + +nodes_bottom = mesh.find_nodes("Y", 0) +nodes_disk_top = mesh.find_nodes("Y", mesh.bounding_box.ymax) + +pb.bc.add("Dirichlet", nodes_bottom, "Disp", 0) +pb.bc.add("Dirichlet", nodes_disk_top, "Disp", [0, imposed_disp]) +pb.set_nr_criterion("Displacement", tol=5e-3, max_subiter=8) + +# ========================================================================= +# Output and tracking callback +# ========================================================================= + +if not os.path.isdir("results"): + os.mkdir("results") +res = pb.add_output("results/indentation_2d", solid, ["Disp", "Stress"]) + +history = {"time": [], "reaction_y": []} + + +def track(pb): + history["time"].append(pb.time) + F = pb.get_ext_forces("Disp") + history["reaction_y"].append(np.sum(F[1, nodes_disk_top])) + + +# ========================================================================= +# Solve +# ========================================================================= + +print("=" * 60) +print("2D DISK INDENTATION -- IPC CONTACT") +print("=" * 60) +pb.nlsolve( + dt=0.05, + tmax=1, + update_dt=True, + print_info=1, + callback=track, +) + +# ========================================================================= +# Hertz comparison +# ========================================================================= + +try: + import matplotlib.pyplot as plt + + t = np.array(history["time"]) + Fy = -np.array(history["reaction_y"]) # positive = pushing down + + # Approximate indentation depth: prescribed displacement minus gap + # (valid because the disk is 100x stiffer, so it barely deforms) + delta = np.maximum(t * abs(imposed_disp) - gap, 0) + + # Hertz theory for 2D plane strain (cylinder on elastic half-space) + # Effective modulus: 1/E* = (1-nu1^2)/E1 + (1-nu2^2)/E2 + E_star = 1.0 / ((1 - nu**2) / E_plate + (1 - nu**2) / E_disk) + + # 2D Hertz (Johnson, Contact Mechanics, Ch. 4): + # P = pi*E*/(4*R) * a^2 (force per unit length) + # delta = a^2/(4*R) * (2*ln(4*R/a) - 1) (indentation depth) + # These two equations are solved parametrically in a. + + def hertz_2d(a): + """Return (delta, P) for given contact half-width a.""" + d = a**2 / (4 * R) * (2 * np.log(4 * R / a) - 1) + P = np.pi * E_star / (4 * R) * a**2 + return d, P + + # Build Hertz curve by sweeping contact half-width + a_vals = np.linspace(0.01, R * 0.95, 500) + delta_hertz, F_hertz = np.array([hertz_2d(a) for a in a_vals]).T + + # Keep only the range that overlaps with FEM data + mask = delta_hertz <= delta.max() * 1.05 + delta_hertz = delta_hertz[mask] + F_hertz = F_hertz[mask] + + fig, ax = plt.subplots(figsize=(7, 5)) + ax.plot(delta, Fy, "o-", ms=4, label="FEM (IPC)") + ax.plot(delta_hertz, F_hertz, "--", lw=2, label="Hertz (2D half-space)") + ax.set_xlabel("Indentation depth (mm)") + ax.set_ylabel("Force per unit thickness (N/mm)") + ax.set_title("2D Hertz Indentation -- Disk on Rectangle") + ax.legend() + ax.grid(True, alpha=0.3) + fig.tight_layout() + fig.savefig("results/indentation_2d_hertz.png", dpi=150) + print("Hertz comparison saved to results/indentation_2d_hertz.png") + plt.show() +except ImportError: + print("matplotlib not available -- skipping Hertz comparison plot") + +# ========================================================================= +# Stress plot +# ========================================================================= + +res.plot("Stress", "vm", "Node", show=False, scale=1) diff --git a/examples/contact/ipc/indentation_3d.py b/examples/contact/ipc/indentation_3d.py new file mode 100644 index 00000000..8d7fa0ce --- /dev/null +++ b/examples/contact/ipc/indentation_3d.py @@ -0,0 +1,255 @@ +""" +Sphere indentation (3D) — IPC contact with Hertz comparison +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +A stiff elastic sphere is pressed vertically into a softer elastic block. +This is the classical 3D Hertz contact benchmark. + +Both bodies are meshed with gmsh, using graded size fields that concentrate +elements near the contact zone. + +The simulation uses the IPC barrier method for contact and compares the +computed force-indentation curve to the Hertz analytical prediction: +``F = 4/3 * E* * sqrt(R) * delta^(3/2)``. + +The FEM force is expected to be ~15-20 % higher than the Hertz prediction +because the plate has finite dimensions and a clamped base (whereas Hertz +assumes an elastic half-space). + +.. note:: + Requires ``ipctk`` and ``gmsh``. +""" + +import fedoo as fd +import numpy as np +import os +import tempfile + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + +# ========================================================================= +# Parameters +# ========================================================================= + +fd.ModelingSpace("3D") + +# Units: N, mm, MPa +E_plate = 1e3 # soft plate +E_sphere = 1e5 # stiff sphere (quasi-rigid, 100x stiffer) +nu = 0.3 +R = 5.0 # sphere radius +plate_half = 25.0 # half-width of the plate +plate_h = 25.0 # plate thickness +gap = 0.1 # initial gap between sphere bottom and plate top +imposed_disp = -2.0 # total vertical displacement of sphere top + +sphere_cz = plate_h + R + gap # sphere centre z-coordinate + +# ========================================================================= +# Sphere mesh (gmsh) +# ========================================================================= + +import gmsh + +gmsh.initialize() +gmsh.option.setNumber("General.Verbosity", 1) + +sphere_tag = gmsh.model.occ.addSphere(0, 0, sphere_cz, R) +gmsh.model.occ.synchronize() +gmsh.model.addPhysicalGroup(3, [sphere_tag], tag=2, name="sphere") + +# Fine elements near the contact tip (bottom of sphere), coarser elsewhere +gmsh.model.mesh.field.add("Box", 1) +gmsh.model.mesh.field.setNumber(1, "VIn", 0.4) +gmsh.model.mesh.field.setNumber(1, "VOut", 1.2) +gmsh.model.mesh.field.setNumber(1, "XMin", -R) +gmsh.model.mesh.field.setNumber(1, "XMax", R) +gmsh.model.mesh.field.setNumber(1, "YMin", -R) +gmsh.model.mesh.field.setNumber(1, "YMax", R) +gmsh.model.mesh.field.setNumber(1, "ZMin", sphere_cz - R) +gmsh.model.mesh.field.setNumber(1, "ZMax", sphere_cz - R + 3) +gmsh.model.mesh.field.setNumber(1, "Thickness", 2) +gmsh.model.mesh.field.setAsBackgroundMesh(1) +gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) + +gmsh.model.mesh.generate(3) +sphere_msh = os.path.join(tempfile.gettempdir(), "sphere_indent3d.msh") +gmsh.write(sphere_msh) +gmsh.finalize() + +mesh_sphere = fd.mesh.import_msh(sphere_msh) +mesh_sphere.element_sets["sphere"] = mesh_sphere.element_sets.get( + "sphere", np.arange(mesh_sphere.n_elements) +) +print(f"Sphere mesh: {mesh_sphere.n_nodes} nodes, {mesh_sphere.n_elements} tet4") + +# ========================================================================= +# Plate mesh (gmsh -- graded refinement near contact) +# ========================================================================= + +gmsh.initialize() +gmsh.option.setNumber("General.Verbosity", 1) + +plate_tag = gmsh.model.occ.addBox( + -plate_half, + -plate_half, + 0, + 2 * plate_half, + 2 * plate_half, + plate_h, +) +gmsh.model.occ.synchronize() +gmsh.model.addPhysicalGroup(3, [plate_tag], tag=1, name="plate") + +# Fine elements beneath the indenter, coarse far away +gmsh.model.mesh.field.add("Box", 1) +gmsh.model.mesh.field.setNumber(1, "VIn", 0.6) +gmsh.model.mesh.field.setNumber(1, "VOut", 5.0) +gmsh.model.mesh.field.setNumber(1, "XMin", -6) +gmsh.model.mesh.field.setNumber(1, "XMax", 6) +gmsh.model.mesh.field.setNumber(1, "YMin", -6) +gmsh.model.mesh.field.setNumber(1, "YMax", 6) +gmsh.model.mesh.field.setNumber(1, "ZMin", plate_h - 5) +gmsh.model.mesh.field.setNumber(1, "ZMax", plate_h) +gmsh.model.mesh.field.setNumber(1, "Thickness", 3) +gmsh.model.mesh.field.setAsBackgroundMesh(1) +gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) + +gmsh.model.mesh.generate(3) +plate_msh = os.path.join(tempfile.gettempdir(), "plate_indent3d.msh") +gmsh.write(plate_msh) +gmsh.finalize() + +mesh_plate = fd.mesh.import_msh(plate_msh) +mesh_plate.element_sets["plate"] = mesh_plate.element_sets.get( + "plate", np.arange(mesh_plate.n_elements) +) +print(f"Plate mesh: {mesh_plate.n_nodes} nodes, {mesh_plate.n_elements} tet4") + +# ========================================================================= +# Stack into single mesh +# ========================================================================= + +mesh = fd.Mesh.stack(mesh_plate, mesh_sphere) +print(f"Total mesh: {mesh.n_nodes} nodes, {mesh.n_elements} tet4") + +# ========================================================================= +# IPC contact +# ========================================================================= + +surf = fd.mesh.extract_surface(mesh, quad2tri=True) +ipc_contact = fd.constraint.IPCContact( + mesh, + surface_mesh=surf, + dhat=0.05, # absolute dhat (< gap to avoid initial contact) + dhat_is_relative=False, + use_ccd=True, +) + +# ========================================================================= +# Material, assembly, BCs +# ========================================================================= + +mat_plate = fd.constitutivelaw.ElasticIsotrop(E_plate, nu) +mat_sphere = fd.constitutivelaw.ElasticIsotrop(E_sphere, nu) +material = fd.constitutivelaw.Heterogeneous( + (mat_plate, mat_sphere), + ("plate", "sphere"), +) + +wf = fd.weakform.StressEquilibrium(material, nlgeom=False) +solid = fd.Assembly.create(wf, mesh) +assembly = fd.Assembly.sum(solid, ipc_contact) + +pb = fd.problem.NonLinear(assembly) + +nodes_bottom = mesh.find_nodes("Z", 0) +nodes_sphere_top = mesh.find_nodes("Z", mesh.bounding_box.zmax) + +pb.bc.add("Dirichlet", nodes_bottom, "Disp", 0) +pb.bc.add("Dirichlet", nodes_sphere_top, "Disp", [0, 0, imposed_disp]) +pb.set_nr_criterion("Displacement", tol=5e-3, max_subiter=8) + +# ========================================================================= +# Output and tracking callback +# ========================================================================= + +if not os.path.isdir("results"): + os.mkdir("results") +res = pb.add_output("results/indentation_3d", solid, ["Disp", "Stress"]) + +history = {"time": [], "reaction_z": []} + + +def track(pb): + history["time"].append(pb.time) + F = pb.get_ext_forces("Disp") + history["reaction_z"].append(np.sum(F[2, nodes_sphere_top])) + + +# ========================================================================= +# Solve +# ========================================================================= + +print("=" * 60) +print("3D SPHERE INDENTATION -- IPC CONTACT") +print("=" * 60) +pb.nlsolve( + dt=0.05, + tmax=1, + update_dt=True, + print_info=1, + callback=track, +) + +# ========================================================================= +# Hertz comparison +# ========================================================================= + +try: + import matplotlib.pyplot as plt + + t = np.array(history["time"]) + Fz = -np.array(history["reaction_z"]) # positive = pushing down + + # Approximate indentation depth (prescribed - gap) + # valid because the sphere is 100x stiffer so it barely deforms + delta = np.maximum(t * abs(imposed_disp) - gap, 0) + + # Hertz theory for 3D -- two elastic bodies + # 1/E* = (1-nu1^2)/E1 + (1-nu2^2)/E2 + E_star = 1.0 / ((1 - nu**2) / E_plate + (1 - nu**2) / E_sphere) + + # F = 4/3 * E* * sqrt(R) * delta^(3/2) + delta_hertz = np.linspace(0, delta.max(), 200) + F_hertz = (4.0 / 3.0) * E_star * np.sqrt(R) * delta_hertz**1.5 + + fig, ax = plt.subplots(figsize=(7, 5)) + ax.plot(delta, Fz, "o-", ms=4, label="FEM (IPC)") + ax.plot(delta_hertz, F_hertz, "--", lw=2, label="Hertz (3D half-space)") + ax.set_xlabel("Indentation depth (mm)") + ax.set_ylabel("Force (N)") + ax.set_title("3D Hertz Indentation -- Sphere on Plate") + ax.legend() + ax.grid(True, alpha=0.3) + fig.tight_layout() + fig.savefig("results/indentation_3d_hertz.png", dpi=150) + print("Hertz comparison saved to results/indentation_3d_hertz.png") + plt.show() +except ImportError: + print("matplotlib not available -- skipping Hertz comparison plot") + +# ========================================================================= +# Stress plot +# ========================================================================= + +res.plot("Stress", "vm", "Node", show=False, scale=1, elevation=75, azimuth=20) + +# --- Video output (uncomment to export MP4) --- +# res.write_movie("results/indentation_3d", "Stress", "vm", "Node", +# framerate=10, quality=8, elevation=75, azimuth=20) +# print("Movie saved to results/indentation_3d.mp4") diff --git a/examples/contact/ipc/punch_indentation.py b/examples/contact/ipc/punch_indentation.py new file mode 100644 index 00000000..65d6c3c9 --- /dev/null +++ b/examples/contact/ipc/punch_indentation.py @@ -0,0 +1,203 @@ +""" +Hemispherical punch indentation (3D, IPC method) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +A cylindrical punch with a hemispherical tip is pressed into a thick +elastic plate. Both bodies deform visibly at the contact interface -- +a classic contact mechanics benchmark. + +Both bodies are meshed with gmsh for proper curved geometry. + +This example uses the IPC contact method which handles the curved +punch geometry robustly. + +.. note:: + Requires ``ipctk`` and ``gmsh``. +""" + +import fedoo as fd +import numpy as np +import os +import tempfile + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + +fd.ModelingSpace("3D") + +E_plate, E_punch, nu = 1e4, 3e4, 0.3 +R_punch = 3.0 # hemisphere radius +H_cyl = 5.0 # cylinder height above hemisphere +plate_half = 10.0 # plate half-width +plate_h = 10.0 # plate thickness +gap = 0.05 # initial gap + +# Punch bottom (hemisphere apex) sits at z = plate_h + gap +punch_base_z = plate_h + gap + +# ========================================================================= +# Punch mesh (gmsh: hemisphere + cylinder) +# ========================================================================= + +import gmsh + +gmsh.initialize() +gmsh.option.setNumber("General.Verbosity", 1) + +# Hemisphere: centred at (0, 0, punch_base_z + R), only bottom half +sphere_tag = gmsh.model.occ.addSphere(0, 0, punch_base_z + R_punch, R_punch) +# Cut box removes the top half of the sphere +cut_box = gmsh.model.occ.addBox( + -R_punch - 1, + -R_punch - 1, + punch_base_z + R_punch, + 2 * (R_punch + 1), + 2 * (R_punch + 1), + R_punch + 1, +) +hemi = gmsh.model.occ.cut([(3, sphere_tag)], [(3, cut_box)])[0] +hemi_tag = hemi[0][1] + +# Cylinder on top of hemisphere +cyl_tag = gmsh.model.occ.addCylinder( + 0, + 0, + punch_base_z + R_punch, + 0, + 0, + H_cyl, + R_punch, +) + +# Fuse hemisphere + cylinder into one volume +punch_parts = gmsh.model.occ.fuse([(3, hemi_tag)], [(3, cyl_tag)])[0] +punch_vol_tag = punch_parts[0][1] +gmsh.model.occ.synchronize() +gmsh.model.addPhysicalGroup(3, [punch_vol_tag], tag=2, name="punch") + +# Mesh size: fine at the tip, coarser on the cylinder +gmsh.model.mesh.field.add("Box", 1) +gmsh.model.mesh.field.setNumber(1, "VIn", 0.4) +gmsh.model.mesh.field.setNumber(1, "VOut", 1.5) +gmsh.model.mesh.field.setNumber(1, "XMin", -R_punch) +gmsh.model.mesh.field.setNumber(1, "XMax", R_punch) +gmsh.model.mesh.field.setNumber(1, "YMin", -R_punch) +gmsh.model.mesh.field.setNumber(1, "YMax", R_punch) +gmsh.model.mesh.field.setNumber(1, "ZMin", punch_base_z) +gmsh.model.mesh.field.setNumber(1, "ZMax", punch_base_z + R_punch) +gmsh.model.mesh.field.setNumber(1, "Thickness", 1) +gmsh.model.mesh.field.setAsBackgroundMesh(1) +gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) + +gmsh.model.mesh.generate(3) +punch_msh = os.path.join(tempfile.gettempdir(), "punch_hemi.msh") +gmsh.write(punch_msh) +gmsh.finalize() + +mesh_punch = fd.mesh.import_msh(punch_msh) +mesh_punch.element_sets["punch"] = mesh_punch.element_sets.get( + "punch", np.arange(mesh_punch.n_elements) +) +print(f"Punch mesh: {mesh_punch.n_nodes} nodes, {mesh_punch.n_elements} tet4") + +# ========================================================================= +# Plate mesh (gmsh) +# ========================================================================= + +gmsh.initialize() +gmsh.option.setNumber("General.Verbosity", 1) + +plate_tag = gmsh.model.occ.addBox( + -plate_half, + -plate_half, + 0, + 2 * plate_half, + 2 * plate_half, + plate_h, +) +gmsh.model.occ.synchronize() +gmsh.model.addPhysicalGroup(3, [plate_tag], tag=1, name="plate") + +gmsh.model.mesh.field.add("Box", 1) +gmsh.model.mesh.field.setNumber(1, "VIn", 0.5) +gmsh.model.mesh.field.setNumber(1, "VOut", 4.0) +gmsh.model.mesh.field.setNumber(1, "XMin", -4) +gmsh.model.mesh.field.setNumber(1, "XMax", 4) +gmsh.model.mesh.field.setNumber(1, "YMin", -4) +gmsh.model.mesh.field.setNumber(1, "YMax", 4) +gmsh.model.mesh.field.setNumber(1, "ZMin", plate_h - 4) +gmsh.model.mesh.field.setNumber(1, "ZMax", plate_h) +gmsh.model.mesh.field.setNumber(1, "Thickness", 2) +gmsh.model.mesh.field.setAsBackgroundMesh(1) +gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) + +gmsh.model.mesh.generate(3) +plate_msh = os.path.join(tempfile.gettempdir(), "plate_punch.msh") +gmsh.write(plate_msh) +gmsh.finalize() + +mesh_plate = fd.mesh.import_msh(plate_msh) +mesh_plate.element_sets["plate"] = mesh_plate.element_sets.get( + "plate", np.arange(mesh_plate.n_elements) +) +print(f"Plate mesh: {mesh_plate.n_nodes} nodes, {mesh_plate.n_elements} tet4") + +# ========================================================================= +# Stack into single mesh +# ========================================================================= + +mesh = fd.Mesh.stack(mesh_plate, mesh_punch) +print(f"Total mesh: {mesh.n_nodes} nodes, {mesh.n_elements} tet4") + +# ========================================================================= +# IPC contact, material, BCs +# ========================================================================= + +nodes_bottom = mesh.find_nodes("Z", 0) +nodes_punch_top = mesh.find_nodes("Z", mesh.bounding_box.zmax) + +surf_ipc = fd.mesh.extract_surface(mesh, quad2tri=True) +ipc_contact = fd.constraint.IPCContact( + mesh, + surface_mesh=surf_ipc, + dhat=1e-2, + dhat_is_relative=True, + use_ccd=True, +) + +mat = fd.constitutivelaw.Heterogeneous( + ( + fd.constitutivelaw.ElasticIsotrop(E_plate, nu), + fd.constitutivelaw.ElasticIsotrop(E_punch, nu), + ), + ("plate", "punch"), +) +wf = fd.weakform.StressEquilibrium(mat, nlgeom=True) +solid = fd.Assembly.create(wf, mesh) +assembly = fd.Assembly.sum(solid, ipc_contact) + +pb = fd.problem.NonLinear(assembly) + +if not os.path.isdir("results"): + os.mkdir("results") +res = pb.add_output("results/punch_ipc", solid, ["Disp", "Stress"]) + +pb.bc.add("Dirichlet", nodes_bottom, "Disp", 0) +pb.bc.add("Dirichlet", nodes_punch_top, "Disp", [0, 0, -2.0]) +pb.set_nr_criterion("Displacement", tol=5e-3, max_subiter=8) + +print("=" * 60) +print("HEMISPHERICAL PUNCH -- IPC CONTACT") +print("=" * 60) +pb.nlsolve(dt=0.05, tmax=1, update_dt=True, print_info=1, interval_output=0.05) + +# --- Static plot --- +res.plot("Stress", "vm", "Node", show=False, scale=1, elevation=75, azimuth=20) + +# --- Video output (uncomment to export MP4) --- +# res.write_movie("results/punch_ipc", "Stress", "vm", "Node", +# framerate=10, quality=8, elevation=75, azimuth=20) +# print("Movie saved to results/punch_ipc.mp4") diff --git a/examples/contact/ipc/self_contact.py b/examples/contact/ipc/self_contact.py new file mode 100644 index 00000000..6099b349 --- /dev/null +++ b/examples/contact/ipc/self_contact.py @@ -0,0 +1,80 @@ +""" +Self-contact — hole plate compression (2D, IPC method) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +A plate with a large circular hole is compressed until self-contact +occurs on the hole's inner surface. Uses the IPC self-contact method +(``fd.constraint.IPCSelfContact``) with CCD line search. + +Large strain (Updated Lagrangian with logarithmic corotational) +and EPICP elasto-plastic material. + +.. note:: + Requires ``ipctk`` and ``simcoon``. +""" + +import fedoo as fd +import numpy as np +import os + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + +fd.ModelingSpace("2D") + +# --- Geometry --- +mesh = fd.mesh.hole_plate_mesh(nr=15, nt=15, length=100, height=100, radius=45) + +# --- IPC self-contact --- +# dhat=3e-3 relative (~0.3% of bbox diagonal) -- slightly larger than +# the default 1e-3 to give the barrier room when the hole collapses. +# use_ccd=True prevents surfaces from crossing between NR iterations. +contact = fd.constraint.IPCSelfContact( + mesh, + dhat=3e-3, + dhat_is_relative=True, + use_ccd=True, +) + +nodes_top = mesh.find_nodes("Y", mesh.bounding_box.ymax) +nodes_bottom = mesh.find_nodes("Y", mesh.bounding_box.ymin) + +# --- Material: EPICP (E, nu, alpha, sigmaY, H, beta) --- +E, nu = 200e3, 0.3 +props = np.array([E, nu, 1e-5, 300, 1000, 0.3]) +material = fd.constitutivelaw.Simcoon("EPICP", props) + +# --- Large strain: UL + logarithmic corotational (default corate="log") --- +wf = fd.weakform.StressEquilibrium(material, nlgeom="UL") +solid_assembly = fd.Assembly.create(wf, mesh) +assembly = fd.Assembly.sum(solid_assembly, contact) + +pb = fd.problem.NonLinear(assembly) + +if not os.path.isdir("results"): + os.mkdir("results") +res = pb.add_output( + "results/self_contact_ipc", solid_assembly, ["Disp", "Stress", "Strain"] +) + +pb.bc.add("Dirichlet", nodes_bottom, "Disp", 0) +pb.bc.add("Dirichlet", nodes_top, "Disp", [0, -70]) +pb.set_nr_criterion("Displacement", tol=5e-3, max_subiter=15) +# pb.set_nr_criterion("Force", tol=1e-2, max_subiter=15) + + +pb.nlsolve( + dt=0.05, + tmax=1, + update_dt=True, + dt_increase_niter=8, + print_info=1, + interval_output=0.01, +) + +# --- Static plot --- +res.plot("Stress", "vm", "Node", show=False, scale=1) + +# --- Video output (uncomment to export MP4) --- +# res.write_movie("results/self_contact_ipc", "Stress", "vm", "Node", +# framerate=24, quality=8, clim=[0, 1.5e3]) +# print("Movie saved to results/self_contact_ipc.mp4") diff --git a/examples/contact/ogc/indentation_2d.py b/examples/contact/ogc/indentation_2d.py new file mode 100644 index 00000000..c93a8753 --- /dev/null +++ b/examples/contact/ogc/indentation_2d.py @@ -0,0 +1,207 @@ +""" +Disk indentation (2D plane strain) — OGC contact with Hertz comparison +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Same benchmark as ``ipc/indentation_2d.py`` but using the OGC +(Offset Geometric Contact) trust-region method instead of CCD. +OGC filters the displacement per vertex rather than uniformly +scaling the step, which can improve convergence. + +.. note:: + Requires ``ipctk`` and ``gmsh``. +""" + +import fedoo as fd +import numpy as np +import os +import tempfile + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + +# ========================================================================= +# Parameters +# ========================================================================= + +fd.ModelingSpace("2D") + +E_plate = 1e3 +E_disk = 1e5 +nu = 0.3 +R = 5.0 +plate_half = 30.0 +plate_h = 40.0 +gap = 0.1 +imposed_disp = -2.0 + +# ========================================================================= +# Plate mesh (gmsh -- graded refinement near contact) +# ========================================================================= + +import gmsh + +gmsh.initialize() +gmsh.option.setNumber("General.Verbosity", 1) + +plate_tag = gmsh.model.occ.addRectangle( + -plate_half, + 0, + 0, + 2 * plate_half, + plate_h, +) +gmsh.model.occ.synchronize() +gmsh.model.addPhysicalGroup(2, [plate_tag], tag=1, name="plate") + +gmsh.model.mesh.field.add("Box", 1) +gmsh.model.mesh.field.setNumber(1, "VIn", 0.3) +gmsh.model.mesh.field.setNumber(1, "VOut", 3.0) +gmsh.model.mesh.field.setNumber(1, "XMin", -6) +gmsh.model.mesh.field.setNumber(1, "XMax", 6) +gmsh.model.mesh.field.setNumber(1, "YMin", plate_h - 4) +gmsh.model.mesh.field.setNumber(1, "YMax", plate_h) +gmsh.model.mesh.field.setNumber(1, "Thickness", 4) +gmsh.model.mesh.field.setAsBackgroundMesh(1) +gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) + +gmsh.model.mesh.generate(2) +plate_msh = os.path.join(tempfile.gettempdir(), "plate_indent2d.msh") +gmsh.write(plate_msh) +gmsh.finalize() + +mesh_plate = fd.mesh.import_msh(plate_msh, mesh_type="surface") +if mesh_plate.nodes.shape[1] == 3: + mesh_plate.nodes = mesh_plate.nodes[:, :2] +mesh_plate.element_sets["plate"] = np.arange(mesh_plate.n_elements) +print(f"Plate mesh: {mesh_plate.n_nodes} nodes, {mesh_plate.n_elements} elems") + +# ========================================================================= +# Disk mesh +# ========================================================================= + +mesh_disk = fd.mesh.disk_mesh(radius=R, nr=12, nt=24, elm_type="tri3") +mesh_disk.nodes += np.array([0, plate_h + R + gap]) +mesh_disk.element_sets["disk"] = np.arange(mesh_disk.n_elements) +print(f"Disk mesh: {mesh_disk.n_nodes} nodes, {mesh_disk.n_elements} elems") + +mesh = fd.Mesh.stack(mesh_plate, mesh_disk) +print( + f"Total mesh: {mesh.n_nodes} nodes, {mesh.n_elements} elems, type={mesh.elm_type}" +) + +# ========================================================================= +# IPC contact with OGC trust-region +# ========================================================================= + +surf = fd.mesh.extract_surface(mesh) +ipc_contact = fd.constraint.IPCContact( + mesh, + surface_mesh=surf, + dhat=0.05, + dhat_is_relative=False, + use_ogc=True, +) + +# ========================================================================= +# Material, assembly, BCs +# ========================================================================= + +mat_plate = fd.constitutivelaw.ElasticIsotrop(E_plate, nu) +mat_disk = fd.constitutivelaw.ElasticIsotrop(E_disk, nu) +material = fd.constitutivelaw.Heterogeneous( + (mat_plate, mat_disk), + ("plate", "disk"), +) + +wf = fd.weakform.StressEquilibrium(material, nlgeom=False) +solid = fd.Assembly.create(wf, mesh) +assembly = fd.Assembly.sum(solid, ipc_contact) + +pb = fd.problem.NonLinear(assembly) + +nodes_bottom = mesh.find_nodes("Y", 0) +nodes_disk_top = mesh.find_nodes("Y", mesh.bounding_box.ymax) + +pb.bc.add("Dirichlet", nodes_bottom, "Disp", 0) +pb.bc.add("Dirichlet", nodes_disk_top, "Disp", [0, imposed_disp]) +pb.set_nr_criterion("Displacement", tol=5e-3, max_subiter=8) + +# ========================================================================= +# Output and tracking callback +# ========================================================================= + +if not os.path.isdir("results"): + os.mkdir("results") +res = pb.add_output("results/indentation_2d", solid, ["Disp", "Stress"]) + +history = {"time": [], "reaction_y": []} + + +def track(pb): + history["time"].append(pb.time) + F = pb.get_ext_forces("Disp") + history["reaction_y"].append(np.sum(F[1, nodes_disk_top])) + + +# ========================================================================= +# Solve +# ========================================================================= + +print("=" * 60) +print("2D DISK INDENTATION -- OGC TRUST-REGION") +print("=" * 60) +pb.nlsolve( + dt=0.05, + tmax=1, + update_dt=True, + print_info=1, + callback=track, +) + +# ========================================================================= +# Hertz comparison +# ========================================================================= + +try: + import matplotlib.pyplot as plt + + t = np.array(history["time"]) + Fy = -np.array(history["reaction_y"]) + + delta = np.maximum(t * abs(imposed_disp) - gap, 0) + + E_star = 1.0 / ((1 - nu**2) / E_plate + (1 - nu**2) / E_disk) + + def hertz_2d(a): + d = a**2 / (4 * R) * (2 * np.log(4 * R / a) - 1) + P = np.pi * E_star / (4 * R) * a**2 + return d, P + + a_vals = np.linspace(0.01, R * 0.95, 500) + delta_hertz, F_hertz = np.array([hertz_2d(a) for a in a_vals]).T + + mask = delta_hertz <= delta.max() * 1.05 + delta_hertz = delta_hertz[mask] + F_hertz = F_hertz[mask] + + fig, ax = plt.subplots(figsize=(7, 5)) + ax.plot(delta, Fy, "o-", ms=4, label="FEM (OGC)") + ax.plot(delta_hertz, F_hertz, "--", lw=2, label="Hertz (2D half-space)") + ax.set_xlabel("Indentation depth (mm)") + ax.set_ylabel("Force per unit thickness (N/mm)") + ax.set_title("2D Hertz Indentation -- OGC Trust-Region") + ax.legend() + ax.grid(True, alpha=0.3) + fig.tight_layout() + fig.savefig("results/indentation_2d_hertz.png", dpi=150) + print("Hertz comparison saved to results/indentation_2d_hertz.png") + plt.show() +except ImportError: + print("matplotlib not available -- skipping Hertz comparison plot") + +# ========================================================================= +# Stress plot +# ========================================================================= + +res.plot("Stress", "vm", "Node", show=False, scale=1) diff --git a/examples/contact/ogc/indentation_3d.py b/examples/contact/ogc/indentation_3d.py new file mode 100644 index 00000000..3b4e616d --- /dev/null +++ b/examples/contact/ogc/indentation_3d.py @@ -0,0 +1,236 @@ +""" +Sphere indentation (3D) — OGC contact with Hertz comparison +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Same benchmark as ``ipc/indentation_3d.py`` but using the OGC +(Offset Geometric Contact) trust-region method instead of CCD. + +.. note:: + Requires ``ipctk`` and ``gmsh``. +""" + +import fedoo as fd +import numpy as np +import os +import tempfile + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + +# ========================================================================= +# Parameters +# ========================================================================= + +fd.ModelingSpace("3D") + +E_plate = 1e3 +E_sphere = 1e5 +nu = 0.3 +R = 5.0 +plate_half = 25.0 +plate_h = 25.0 +gap = 0.1 +imposed_disp = -2.0 + +sphere_cz = plate_h + R + gap + +# ========================================================================= +# Sphere mesh (gmsh) +# ========================================================================= + +import gmsh + +gmsh.initialize() +gmsh.option.setNumber("General.Verbosity", 1) + +sphere_tag = gmsh.model.occ.addSphere(0, 0, sphere_cz, R) +gmsh.model.occ.synchronize() +gmsh.model.addPhysicalGroup(3, [sphere_tag], tag=2, name="sphere") + +gmsh.model.mesh.field.add("Box", 1) +gmsh.model.mesh.field.setNumber(1, "VIn", 0.4) +gmsh.model.mesh.field.setNumber(1, "VOut", 1.2) +gmsh.model.mesh.field.setNumber(1, "XMin", -R) +gmsh.model.mesh.field.setNumber(1, "XMax", R) +gmsh.model.mesh.field.setNumber(1, "YMin", -R) +gmsh.model.mesh.field.setNumber(1, "YMax", R) +gmsh.model.mesh.field.setNumber(1, "ZMin", sphere_cz - R) +gmsh.model.mesh.field.setNumber(1, "ZMax", sphere_cz - R + 3) +gmsh.model.mesh.field.setNumber(1, "Thickness", 2) +gmsh.model.mesh.field.setAsBackgroundMesh(1) +gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) + +gmsh.model.mesh.generate(3) +sphere_msh = os.path.join(tempfile.gettempdir(), "sphere_indent3d.msh") +gmsh.write(sphere_msh) +gmsh.finalize() + +mesh_sphere = fd.mesh.import_msh(sphere_msh) +mesh_sphere.element_sets["sphere"] = mesh_sphere.element_sets.get( + "sphere", np.arange(mesh_sphere.n_elements) +) +print(f"Sphere mesh: {mesh_sphere.n_nodes} nodes, {mesh_sphere.n_elements} tet4") + +# ========================================================================= +# Plate mesh (gmsh -- graded refinement near contact) +# ========================================================================= + +gmsh.initialize() +gmsh.option.setNumber("General.Verbosity", 1) + +plate_tag = gmsh.model.occ.addBox( + -plate_half, + -plate_half, + 0, + 2 * plate_half, + 2 * plate_half, + plate_h, +) +gmsh.model.occ.synchronize() +gmsh.model.addPhysicalGroup(3, [plate_tag], tag=1, name="plate") + +gmsh.model.mesh.field.add("Box", 1) +gmsh.model.mesh.field.setNumber(1, "VIn", 0.6) +gmsh.model.mesh.field.setNumber(1, "VOut", 5.0) +gmsh.model.mesh.field.setNumber(1, "XMin", -6) +gmsh.model.mesh.field.setNumber(1, "XMax", 6) +gmsh.model.mesh.field.setNumber(1, "YMin", -6) +gmsh.model.mesh.field.setNumber(1, "YMax", 6) +gmsh.model.mesh.field.setNumber(1, "ZMin", plate_h - 5) +gmsh.model.mesh.field.setNumber(1, "ZMax", plate_h) +gmsh.model.mesh.field.setNumber(1, "Thickness", 3) +gmsh.model.mesh.field.setAsBackgroundMesh(1) +gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) + +gmsh.model.mesh.generate(3) +plate_msh = os.path.join(tempfile.gettempdir(), "plate_indent3d.msh") +gmsh.write(plate_msh) +gmsh.finalize() + +mesh_plate = fd.mesh.import_msh(plate_msh) +mesh_plate.element_sets["plate"] = mesh_plate.element_sets.get( + "plate", np.arange(mesh_plate.n_elements) +) +print(f"Plate mesh: {mesh_plate.n_nodes} nodes, {mesh_plate.n_elements} tet4") + +# ========================================================================= +# Stack into single mesh +# ========================================================================= + +mesh = fd.Mesh.stack(mesh_plate, mesh_sphere) +print(f"Total mesh: {mesh.n_nodes} nodes, {mesh.n_elements} tet4") + +# ========================================================================= +# IPC contact with OGC trust-region +# ========================================================================= + +surf = fd.mesh.extract_surface(mesh, quad2tri=True) +ipc_contact = fd.constraint.IPCContact( + mesh, + surface_mesh=surf, + dhat=0.05, + dhat_is_relative=False, + use_ogc=True, +) + +# ========================================================================= +# Material, assembly, BCs +# ========================================================================= + +mat_plate = fd.constitutivelaw.ElasticIsotrop(E_plate, nu) +mat_sphere = fd.constitutivelaw.ElasticIsotrop(E_sphere, nu) +material = fd.constitutivelaw.Heterogeneous( + (mat_plate, mat_sphere), + ("plate", "sphere"), +) + +wf = fd.weakform.StressEquilibrium(material, nlgeom=False) +solid = fd.Assembly.create(wf, mesh) +assembly = fd.Assembly.sum(solid, ipc_contact) + +pb = fd.problem.NonLinear(assembly) + +nodes_bottom = mesh.find_nodes("Z", 0) +nodes_sphere_top = mesh.find_nodes("Z", mesh.bounding_box.zmax) + +pb.bc.add("Dirichlet", nodes_bottom, "Disp", 0) +pb.bc.add("Dirichlet", nodes_sphere_top, "Disp", [0, 0, imposed_disp]) +pb.set_nr_criterion("Displacement", tol=5e-3, max_subiter=8) + +# ========================================================================= +# Output and tracking callback +# ========================================================================= + +if not os.path.isdir("results"): + os.mkdir("results") +res = pb.add_output("results/indentation_3d", solid, ["Disp", "Stress"]) + +history = {"time": [], "reaction_z": []} + + +def track(pb): + history["time"].append(pb.time) + F = pb.get_ext_forces("Disp") + history["reaction_z"].append(np.sum(F[2, nodes_sphere_top])) + + +# ========================================================================= +# Solve +# ========================================================================= + +print("=" * 60) +print("3D SPHERE INDENTATION -- OGC TRUST-REGION") +print("=" * 60) +pb.nlsolve( + dt=0.05, + tmax=1, + update_dt=True, + print_info=1, + callback=track, +) + +# ========================================================================= +# Hertz comparison +# ========================================================================= + +try: + import matplotlib.pyplot as plt + + t = np.array(history["time"]) + Fz = -np.array(history["reaction_z"]) + + delta = np.maximum(t * abs(imposed_disp) - gap, 0) + + E_star = 1.0 / ((1 - nu**2) / E_plate + (1 - nu**2) / E_sphere) + + delta_hertz = np.linspace(0, delta.max(), 200) + F_hertz = (4.0 / 3.0) * E_star * np.sqrt(R) * delta_hertz**1.5 + + fig, ax = plt.subplots(figsize=(7, 5)) + ax.plot(delta, Fz, "o-", ms=4, label="FEM (OGC)") + ax.plot(delta_hertz, F_hertz, "--", lw=2, label="Hertz (3D half-space)") + ax.set_xlabel("Indentation depth (mm)") + ax.set_ylabel("Force (N)") + ax.set_title("3D Hertz Indentation -- OGC Trust-Region") + ax.legend() + ax.grid(True, alpha=0.3) + fig.tight_layout() + fig.savefig("results/indentation_3d_hertz.png", dpi=150) + print("Hertz comparison saved to results/indentation_3d_hertz.png") + plt.show() +except ImportError: + print("matplotlib not available -- skipping Hertz comparison plot") + +# ========================================================================= +# Stress plot +# ========================================================================= + +res.plot("Stress", "vm", "Node", show=False, scale=1, elevation=75, azimuth=20) + +# --- Video output (uncomment to export MP4) --- +# res.write_movie("results/indentation_3d", "Stress", "vm", "Node", +# framerate=10, quality=8, elevation=75, azimuth=20) +# print("Movie saved to results/indentation_3d.mp4") diff --git a/examples/contact/ogc/punch_indentation.py b/examples/contact/ogc/punch_indentation.py new file mode 100644 index 00000000..45d5fc52 --- /dev/null +++ b/examples/contact/ogc/punch_indentation.py @@ -0,0 +1,193 @@ +""" +Hemispherical punch indentation (3D, OGC method) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Same benchmark as ``ipc/punch_indentation.py`` but using the OGC +(Offset Geometric Contact) trust-region method instead of CCD. + +Both bodies are meshed with gmsh for proper curved geometry. + +.. note:: + Requires ``ipctk`` and ``gmsh``. +""" + +import fedoo as fd +import numpy as np +import os +import tempfile + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + +fd.ModelingSpace("3D") + +E_plate, E_punch, nu = 1e4, 3e4, 0.3 +R_punch = 3.0 +H_cyl = 5.0 +plate_half = 10.0 +plate_h = 10.0 +gap = 0.05 + +punch_base_z = plate_h + gap + +# ========================================================================= +# Punch mesh (gmsh: hemisphere + cylinder) +# ========================================================================= + +import gmsh + +gmsh.initialize() +gmsh.option.setNumber("General.Verbosity", 1) + +sphere_tag = gmsh.model.occ.addSphere(0, 0, punch_base_z + R_punch, R_punch) +cut_box = gmsh.model.occ.addBox( + -R_punch - 1, + -R_punch - 1, + punch_base_z + R_punch, + 2 * (R_punch + 1), + 2 * (R_punch + 1), + R_punch + 1, +) +hemi = gmsh.model.occ.cut([(3, sphere_tag)], [(3, cut_box)])[0] +hemi_tag = hemi[0][1] + +cyl_tag = gmsh.model.occ.addCylinder( + 0, + 0, + punch_base_z + R_punch, + 0, + 0, + H_cyl, + R_punch, +) + +punch_parts = gmsh.model.occ.fuse([(3, hemi_tag)], [(3, cyl_tag)])[0] +punch_vol_tag = punch_parts[0][1] +gmsh.model.occ.synchronize() +gmsh.model.addPhysicalGroup(3, [punch_vol_tag], tag=2, name="punch") + +gmsh.model.mesh.field.add("Box", 1) +gmsh.model.mesh.field.setNumber(1, "VIn", 0.4) +gmsh.model.mesh.field.setNumber(1, "VOut", 1.5) +gmsh.model.mesh.field.setNumber(1, "XMin", -R_punch) +gmsh.model.mesh.field.setNumber(1, "XMax", R_punch) +gmsh.model.mesh.field.setNumber(1, "YMin", -R_punch) +gmsh.model.mesh.field.setNumber(1, "YMax", R_punch) +gmsh.model.mesh.field.setNumber(1, "ZMin", punch_base_z) +gmsh.model.mesh.field.setNumber(1, "ZMax", punch_base_z + R_punch) +gmsh.model.mesh.field.setNumber(1, "Thickness", 1) +gmsh.model.mesh.field.setAsBackgroundMesh(1) +gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) + +gmsh.model.mesh.generate(3) +punch_msh = os.path.join(tempfile.gettempdir(), "punch_hemi.msh") +gmsh.write(punch_msh) +gmsh.finalize() + +mesh_punch = fd.mesh.import_msh(punch_msh) +mesh_punch.element_sets["punch"] = mesh_punch.element_sets.get( + "punch", np.arange(mesh_punch.n_elements) +) +print(f"Punch mesh: {mesh_punch.n_nodes} nodes, {mesh_punch.n_elements} tet4") + +# ========================================================================= +# Plate mesh (gmsh) +# ========================================================================= + +gmsh.initialize() +gmsh.option.setNumber("General.Verbosity", 1) + +plate_tag = gmsh.model.occ.addBox( + -plate_half, + -plate_half, + 0, + 2 * plate_half, + 2 * plate_half, + plate_h, +) +gmsh.model.occ.synchronize() +gmsh.model.addPhysicalGroup(3, [plate_tag], tag=1, name="plate") + +gmsh.model.mesh.field.add("Box", 1) +gmsh.model.mesh.field.setNumber(1, "VIn", 0.5) +gmsh.model.mesh.field.setNumber(1, "VOut", 4.0) +gmsh.model.mesh.field.setNumber(1, "XMin", -4) +gmsh.model.mesh.field.setNumber(1, "XMax", 4) +gmsh.model.mesh.field.setNumber(1, "YMin", -4) +gmsh.model.mesh.field.setNumber(1, "YMax", 4) +gmsh.model.mesh.field.setNumber(1, "ZMin", plate_h - 4) +gmsh.model.mesh.field.setNumber(1, "ZMax", plate_h) +gmsh.model.mesh.field.setNumber(1, "Thickness", 2) +gmsh.model.mesh.field.setAsBackgroundMesh(1) +gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) +gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) + +gmsh.model.mesh.generate(3) +plate_msh = os.path.join(tempfile.gettempdir(), "plate_punch.msh") +gmsh.write(plate_msh) +gmsh.finalize() + +mesh_plate = fd.mesh.import_msh(plate_msh) +mesh_plate.element_sets["plate"] = mesh_plate.element_sets.get( + "plate", np.arange(mesh_plate.n_elements) +) +print(f"Plate mesh: {mesh_plate.n_nodes} nodes, {mesh_plate.n_elements} tet4") + +# ========================================================================= +# Stack into single mesh +# ========================================================================= + +mesh = fd.Mesh.stack(mesh_plate, mesh_punch) +print(f"Total mesh: {mesh.n_nodes} nodes, {mesh.n_elements} tet4") + +# ========================================================================= +# IPC contact with OGC trust-region +# ========================================================================= + +nodes_bottom = mesh.find_nodes("Z", 0) +nodes_punch_top = mesh.find_nodes("Z", mesh.bounding_box.zmax) + +surf_ipc = fd.mesh.extract_surface(mesh, quad2tri=True) +ipc_contact = fd.constraint.IPCContact( + mesh, + surface_mesh=surf_ipc, + dhat=1e-2, + dhat_is_relative=True, + use_ogc=True, +) + +mat = fd.constitutivelaw.Heterogeneous( + ( + fd.constitutivelaw.ElasticIsotrop(E_plate, nu), + fd.constitutivelaw.ElasticIsotrop(E_punch, nu), + ), + ("plate", "punch"), +) +wf = fd.weakform.StressEquilibrium(mat, nlgeom=True) +solid = fd.Assembly.create(wf, mesh) +assembly = fd.Assembly.sum(solid, ipc_contact) + +pb = fd.problem.NonLinear(assembly) + +if not os.path.isdir("results"): + os.mkdir("results") +res = pb.add_output("results/punch_ogc", solid, ["Disp", "Stress"]) + +pb.bc.add("Dirichlet", nodes_bottom, "Disp", 0) +pb.bc.add("Dirichlet", nodes_punch_top, "Disp", [0, 0, -2.0]) +pb.set_nr_criterion("Displacement", tol=5e-3, max_subiter=8) + +print("=" * 60) +print("HEMISPHERICAL PUNCH -- OGC TRUST-REGION") +print("=" * 60) +pb.nlsolve(dt=0.05, tmax=1, update_dt=True, print_info=1, interval_output=0.05) + +# --- Static plot --- +res.plot("Stress", "vm", "Node", show=False, scale=1, elevation=75, azimuth=20) + +# --- Video output (uncomment to export MP4) --- +# res.write_movie("results/punch_ogc", "Stress", "vm", "Node", +# framerate=10, quality=8, elevation=75, azimuth=20) +# print("Movie saved to results/punch_ogc.mp4") diff --git a/examples/contact/ogc/self_contact.py b/examples/contact/ogc/self_contact.py new file mode 100644 index 00000000..1f73534d --- /dev/null +++ b/examples/contact/ogc/self_contact.py @@ -0,0 +1,67 @@ +""" +Self-contact — hole plate compression (2D, OGC method) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Same benchmark as ``ipc/self_contact.py`` but using the OGC +(Offset Geometric Contact) trust-region method instead of CCD. + +Large strain (Updated Lagrangian with logarithmic corotational) +and EPICP elasto-plastic material. + +.. note:: + Requires ``ipctk`` and ``simcoon``. +""" + +import fedoo as fd +import numpy as np +import os + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + +fd.ModelingSpace("2D") + +# --- Geometry --- +mesh = fd.mesh.hole_plate_mesh(nr=15, nt=15, length=100, height=100, radius=45) + +# --- IPC self-contact with OGC trust-region --- +contact = fd.constraint.IPCSelfContact( + mesh, + dhat=3e-3, + dhat_is_relative=True, + use_ogc=True, +) + +nodes_top = mesh.find_nodes("Y", mesh.bounding_box.ymax) +nodes_bottom = mesh.find_nodes("Y", mesh.bounding_box.ymin) + +# --- Material: EPICP (E, nu, alpha, sigmaY, H, beta) --- +E, nu = 200e3, 0.3 +props = np.array([E, nu, 1e-5, 300, 1000, 0.3]) +material = fd.constitutivelaw.Simcoon("EPICP", props) + +# --- Large strain: UL + logarithmic corotational (default corate="log") --- +wf = fd.weakform.StressEquilibrium(material, nlgeom="UL") +solid_assembly = fd.Assembly.create(wf, mesh) +assembly = fd.Assembly.sum(solid_assembly, contact) + +pb = fd.problem.NonLinear(assembly) + +if not os.path.isdir("results"): + os.mkdir("results") +res = pb.add_output( + "results/self_contact_ogc", solid_assembly, ["Disp", "Stress", "Strain"] +) + +pb.bc.add("Dirichlet", nodes_bottom, "Disp", 0) +pb.bc.add("Dirichlet", nodes_top, "Disp", [0, -70]) +pb.set_nr_criterion("Displacement", tol=5e-3, max_subiter=10) + +pb.nlsolve(dt=0.01, tmax=1, update_dt=True, print_info=1, interval_output=0.01) + +# --- Static plot --- +res.plot("Stress", "vm", "Node", show=False, scale=1) + +# --- Video output (uncomment to export MP4) --- +# res.write_movie("results/self_contact_ogc", "Stress", "vm", "Node", +# framerate=24, quality=8, clim=[0, 1.5e3]) +# print("Movie saved to results/self_contact_ogc.mp4") diff --git a/examples/contact/penalty/disk_rectangle.py b/examples/contact/penalty/disk_rectangle.py new file mode 100644 index 00000000..1181306f --- /dev/null +++ b/examples/contact/penalty/disk_rectangle.py @@ -0,0 +1,91 @@ +""" +Disk-rectangle contact (2D, penalty method) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +A stiff elastic disk is pushed into an elasto-plastic rectangle (quad4, +plane strain) using the penalty contact method, then released. + +Demonstrates the legacy penalty-based contact approach +(``fd.constraint.Contact``) with node-to-surface formulation. + +.. note:: + Requires ``simcoon`` for the EPICP elasto-plastic material. +""" + +import fedoo as fd +import numpy as np +import os + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + +fd.ModelingSpace("2D") + +# --- Geometry --- +mesh_rect = fd.mesh.rectangle_mesh( + nx=11, + ny=21, + x_min=0, + x_max=1, + y_min=0, + y_max=1, + elm_type="quad4", +) +mesh_rect.element_sets["rect"] = np.arange(mesh_rect.n_elements) + +mesh_disk = fd.mesh.disk_mesh(radius=0.5, nr=6, nt=6, elm_type="quad4") +mesh_disk.nodes += np.array([1.5, 0.48]) +mesh_disk.element_sets["disk"] = np.arange(mesh_disk.n_elements) + +mesh = fd.Mesh.stack(mesh_rect, mesh_disk) + +nodes_left = mesh.find_nodes("X", 0) +nodes_right = mesh.find_nodes("X", 1) +nodes_bc = mesh.find_nodes("X>1.5") +nodes_bc = list(set(nodes_bc).intersection(mesh.node_sets["boundary"])) + +# --- Contact (penalty) --- +surf = fd.mesh.extract_surface(mesh.extract_elements("disk")) +contact = fd.constraint.Contact(nodes_right, surf) +contact.contact_search_once = True +contact.eps_n = 5e5 +contact.max_dist = 1 + +# --- Material --- +E, nu = 200e3, 0.3 +props = np.array([E, nu, 1e-5, 300, 1000, 0.3]) +material_rect = fd.constitutivelaw.Simcoon("EPICP", props) +material_disk = fd.constitutivelaw.ElasticIsotrop(50e3, nu) +material = fd.constitutivelaw.Heterogeneous( + (material_rect, material_disk), ("rect", "disk") +) + +wf = fd.weakform.StressEquilibrium(material, nlgeom="UL") +solid_assembly = fd.Assembly.create(wf, mesh) +assembly = fd.Assembly.sum(solid_assembly, contact) + +pb = fd.problem.NonLinear(assembly) + +if not os.path.isdir("results"): + os.mkdir("results") +res = pb.add_output( + "results/disk_rectangle_contact", solid_assembly, ["Disp", "Stress", "Strain"] +) + +# --- Step 1: push --- +pb.bc.add("Dirichlet", nodes_left, "Disp", 0) +pb.bc.add("Dirichlet", nodes_bc, "Disp", [-0.4, 0.2]) +pb.set_nr_criterion("Displacement", tol=5e-3, max_subiter=5) +pb.nlsolve(dt=0.005, tmax=1, update_dt=True, print_info=1, interval_output=0.01) + +# --- Step 2: release --- +pb.bc.remove(-1) +pb.bc.add("Dirichlet", nodes_bc, "Disp", [0, 0]) +pb.nlsolve(dt=0.005, tmax=1, update_dt=True, print_info=1, interval_output=0.01) + +# --- Static plot --- +res.plot("Stress", "XX", "Node", show=False, scale=1) + +# --- Video output (uncomment to export MP4) --- +# res.write_movie("results/disk_rectangle_contact", "Stress", "XX", "Node", +# framerate=24, quality=5, clim=[-3e3, 3e3]) +# print("Movie saved to results/disk_rectangle_contact.mp4") diff --git a/examples/contact/penalty/self_contact.py b/examples/contact/penalty/self_contact.py new file mode 100644 index 00000000..d2529996 --- /dev/null +++ b/examples/contact/penalty/self_contact.py @@ -0,0 +1,65 @@ +""" +Self-contact — hole plate compression (2D, penalty method) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +A plate with a large circular hole is compressed until self-contact +occurs on the hole's inner surface. Uses the penalty self-contact +method (``fd.constraint.SelfContact``) with an elasto-plastic material. + +Large strain (Updated Lagrangian with logarithmic corotational). + +.. note:: + Requires ``simcoon`` for the EPICP elasto-plastic material. +""" + +import fedoo as fd +import numpy as np +import os + +os.chdir(os.path.dirname(os.path.abspath(__file__)) or ".") + +fd.ModelingSpace("2D") + +# --- Geometry --- +mesh = fd.mesh.hole_plate_mesh(nr=15, nt=15, length=100, height=100, radius=45) +surf = fd.mesh.extract_surface(mesh) + +# --- Contact (penalty self-contact) --- +contact = fd.constraint.SelfContact(surf, "linear", search_algorithm="bucket") +contact.contact_search_once = True +contact.eps_n = 1e6 + +nodes_top = mesh.find_nodes("Y", mesh.bounding_box.ymax) +nodes_bottom = mesh.find_nodes("Y", mesh.bounding_box.ymin) + +# --- Material: EPICP (E, nu, alpha, sigmaY, H, beta) --- +E, nu = 200e3, 0.3 +props = np.array([E, nu, 1e-5, 300, 1000, 0.3]) +material = fd.constitutivelaw.Simcoon("EPICP", props) + +# --- Large strain --- +wf = fd.weakform.StressEquilibrium(material, nlgeom="UL") +solid_assembly = fd.Assembly.create(wf, mesh) +assembly = fd.Assembly.sum(solid_assembly, contact) + +pb = fd.problem.NonLinear(assembly) + +if not os.path.isdir("results"): + os.mkdir("results") +res = pb.add_output( + "results/self_contact_penalty", solid_assembly, ["Disp", "Stress", "Strain"] +) + +pb.bc.add("Dirichlet", nodes_bottom, "Disp", 0) +pb.bc.add("Dirichlet", nodes_top, "Disp", [0, -70]) +pb.set_nr_criterion("Displacement", tol=5e-3, max_subiter=5) + +pb.nlsolve(dt=0.01, tmax=1, update_dt=True, print_info=1, interval_output=0.01) + +# --- Static plot --- +res.plot("Stress", "vm", "Node", show=False, scale=1) + +# --- Video output (uncomment to export MP4) --- +# res.write_movie("results/self_contact_penalty", "Stress", "vm", "Node", +# framerate=24, quality=8, clim=[0, 1.5e3]) +# print("Movie saved to results/self_contact_penalty.mp4") diff --git a/examples/contact/self_contact.py b/examples/contact/self_contact.py deleted file mode 100644 index a6c3290c..00000000 --- a/examples/contact/self_contact.py +++ /dev/null @@ -1,170 +0,0 @@ -import fedoo as fd -import numpy as np -from time import time -import os -import pylab as plt -from numpy import linalg - -import pyvista as pv - - -start = time() -# --------------- Pre-Treatment -------------------------------------------------------- - -fd.ModelingSpace("2D") - -NLGEOM = "UL" -# Units: N, mm, MPa -h = 100 -w = 1 -L = 100 -E = 200e3 -nu = 0.3 -alpha = 1e-5 # ??? -meshname = "Domain" -uimp = 1 - - -filename = "self_contact" -res_dir = "results/" - - -# mesh = fd.mesh.box_mesh(nx=11, ny=11, nz=11, x_min=0, x_max=L, y_min=0, y_max=h, z_min = 0, z_max = w, elm_type = 'hex8', name = meshname) - -mesh = fd.mesh.hole_plate_mesh(nr=15, nt=15, length=100, height=100, radius=45) -surf = fd.mesh.extract_surface(mesh) - -# surf.plot() #to plot the mesh -# surf.plot_normals(5) #to check normals orientation - -# contact = fd.constraint.contact.SelfContact(surf, 'biliear') -contact = fd.constraint.contact.SelfContact(surf, "linear", search_algorithm="bucket") -contact.contact_search_once = True -contact.eps_n = 1e6 -# contact.eps_a = 1e4 -# contact.limit_soft_contact = 0.01 - -# node sets for boundary conditions -nodes_top = mesh.find_nodes("Y", mesh.bounding_box.ymax) -nodes_bottom = mesh.find_nodes("Y", mesh.bounding_box.ymin) - -mat = 1 -if mat == 0: - props = np.array([E, nu, alpha]) - material = fd.constitutivelaw.Simcoon("ELISO", props, name="ConstitutiveLaw") -elif mat == 1 or mat == 2: - Re = 300 - k = 1000 # 1500 - m = 0.3 # 0.25 - if mat == 1: - props = np.array([E, nu, alpha, Re, k, m]) - material = fd.constitutivelaw.Simcoon("EPICP", props, name="ConstitutiveLaw") - elif mat == 2: - material = fd.constitutivelaw.ElastoPlasticity( - E, nu, Re, name="ConstitutiveLaw" - ) - material.SetHardeningFunction("power", H=k, beta=m) -else: - material = fd.constitutivelaw.ElasticIsotrop(E, nu, name="ConstitutiveLaw") - - -#### trouver pourquoi les deux fonctions suivantes ne donnent pas la même chose !!!! -wf = fd.weakform.StressEquilibrium("ConstitutiveLaw", nlgeom=NLGEOM) - -solid_assembly = fd.Assembly.create( - wf, mesh, name="Assembling" -) # uses MeshChange=True when the mesh change during the time - -assembly = fd.Assembly.sum(solid_assembly, contact) -# assembly = solid_assembly - -pb = fd.problem.NonLinear(assembly) - -# create a 'result' folder and set the desired ouputs -if not (os.path.isdir("results")): - os.mkdir("results") -# results = pb.add_output(res_dir+filename, 'Assembling', ['Disp'], output_type='Node', file_format ='npz') -# results = pb.add_output(res_dir+filename, 'Assembling', ['Cauchy', 'PKII', 'Strain', 'Cauchy_vm', 'Statev', 'Wm'], output_type='GaussPoint', file_format ='npz') - -# results = pb.add_output(res_dir+filename, solid_assembly, ['Disp', 'Cauchy', 'PKII', 'Strain', 'Cauchy_vm', 'Statev', 'Wm']) -results = pb.add_output( - res_dir + filename, solid_assembly, ["Disp", "Stress", "Strain", "Fext"] -) - -# Problem.add_output(res_dir+filename, 'Assembling', ['cauchy', 'PKII', 'strain', 'cauchy_vm', 'statev'], output_type='Element', file_format ='vtk') - - -################### step 1 ################################ -pb.bc.add("Dirichlet", nodes_bottom, "Disp", 0) -pb.bc.add("Dirichlet", nodes_top, "Disp", [0, -70]) - -# Problem.set_solver('cg', precond = True) -# pb.set_nr_criterion("Force", err0 = None, tol = 5e-3, max_subiter = 10) -pb.set_nr_criterion("Displacement", err0=None, tol=5e-3, max_subiter=5) - -# pb.nlsolve(dt = 0.001, tmax = 1, update_dt = False, print_info = 2, interval_output = 0.001) -pb.nlsolve(dt=0.01, tmax=1, update_dt=True, print_info=1, interval_output=0.01) - -E = np.array( - fd.Assembly.get_all()["Assembling"].get_strain( - pb.get_dof_solution(), "GaussPoint", False - ) -).T - -# ################### step 2 ################################ -# bc.Remove() -# #We set initial condition to the applied force to relax the load -# F_app = Problem.get_ext_forces('DispY')[nodes_topCenter] -# bc = Problem.bc.add('Neumann','DispY', 0, nodes_topCenter, initialValue=F_app)#face_center) - -# Problem.nlsolve(dt = 1., update_dt = True, ToleranceNR = 0.01) - -print(time() - start) - - -# ============================================================= -# Example of plots with pyvista - uncomment the desired plot -# ============================================================= - -# ------------------------------------ -# Simple plot with default options -# ------------------------------------ -results.plot("Stress", "vm", "Node", show=True, scale=1) -# results.plot('Stress', 'XX', 'Node', show = True, scale = 1, show_nodes=True) -# results.plot('Stress', 'XX', 'Node', show = True, scale = 1, show_nodes=True, node_labels =True) -# results.plot('Fext', 'X', 'Node', show = True, scale = 1, show_nodes=True) - -# ------------------------------------ -# Simple plot with default options and save to png -# ------------------------------------ -# pl = results.plot('Stress', 0, show = False) -# pl.show(screenshot = "test.png") - -# ------------------------------------ -# Write movie with default options -# ------------------------------------ -# results.write_movie(res_dir+filename, 'Stress', 'vm', framerate = 5, quality = 5) - -results.write_movie( - res_dir + filename, - "Stress", - "vm", - "Node", - framerate=24, - quality=10, - clim=[0, 1.5e3], -) - -# ------------------------------------ -# Save pdf plot -# ------------------------------------ -# pl = results.plot('Stress', 'vm', show = False) -# pl.save_graphic('test.pdf', title='PyVista Export', raster=True, painter=True) - -# ------------------------------------ -# Plot time history -# ------------------------------------ -# from matplotlib import pylab -# t, sigma = results.get_history(('Time','Cauchy'), (0,12), component = 3) -# pylab.plot(t,sigma) -# #or results.plot_history('Cauchy', 12) diff --git a/fedoo/constraint/__init__.py b/fedoo/constraint/__init__.py index fe0fb808..f198ad81 100644 --- a/fedoo/constraint/__init__.py +++ b/fedoo/constraint/__init__.py @@ -2,3 +2,4 @@ from .periodic_bc import PeriodicBC from .contact import Contact, SelfContact from .distributed_load import Pressure, DistributedForce, SurfaceForce +from .ipc_contact import IPCContact, IPCSelfContact diff --git a/fedoo/constraint/contact.py b/fedoo/constraint/contact.py index 7235c479..d17c9f1f 100644 --- a/fedoo/constraint/contact.py +++ b/fedoo/constraint/contact.py @@ -1,20 +1,15 @@ from __future__ import annotations -# from fedoo.core.base import ProblemBase import numpy as np from fedoo.core.base import AssemblyBase -# from fedoo.core.boundary_conditions import BCBase, MPC, ListBC -from fedoo.core.base import ProblemBase -from fedoo.core.mesh import MeshBase, Mesh +from fedoo.core.mesh import Mesh from fedoo.lib_elements.element_list import get_element from scipy import sparse, spatial from fedoo.core.modelingspace import ModelingSpace from fedoo.mesh import extract_surface as extract_surface_mesh from copy import copy -import time - class Contact(AssemblyBase): """Contact Assembly based on a node 2 surface formulation""" @@ -24,43 +19,46 @@ def __init__( slave_nodes: list[int] | Mesh, surface_mesh: Mesh, normal_law: str = "linear", - search_algorithm: str = "bucket", #'bucket', #'search_nearest' + search_algorithm: str = "bucket", # 'bucket', 'search_nearest', 'ipctk' space: ModelingSpace | None = None, name: str = "Contact nodes 2 surface", ): - """ - Assembly related to surface 2 surface contact, using a node 2 surface formulation - with a penality method. + """Contact Assembly based on a node 2 surface penality formulation. Parameters ---------- - slave_nodes: numpy array of int or Mesh List of nodes indices: nodes that belong to the slave surface - or Mesh of the slave surface (in this case, the slave node indices are - extracted). + or Mesh of the slave surface (in this case, the slave node indices + are extracted). surface_mesh: fedoo.Mesh Mesh of the master surface normal_law: str in {'linear', 'bilinear'}, default = 'linear' Type of contact law for the normal contact. space: ModelingSpace - Modeling space associated to the weakform. If None is specified, the active ModelingSpace is considered. + Modeling space associated to the weakform. If None is specified, + the active ModelingSpace is considered. name: str The name of contact assembly Notes - -------------------------- + ----- Several attributes can be modified to change the contact behavior: - * eps_n: contact penalty parameter. Need to be adjusted depending of the rigidity of the material in contact. + * eps_n: contact penalty parameter. Need to be adjusted depending of + the rigidity of the material in contact. * clearance: Contact clearance to adjust the two surfaces. - * contact_search_once: If True (default) the elmement in contact are only searched at the begining of iteration. - It avoid oscilations between contact and non contact state during the NR iterations. + * contact_search_once: If True (default), the elmement in contact + are only searched at the begining of iteration. + It avoid oscilations between contact and non contact state during + the NR iterations. * tol: Tolerance for possible slide of a node outside an element. * max_dist: Max distance from nodes at which contact is considered. - * eps_a: Penalty parameter for soft contact in bilinear contact law (only used if bilinear law is choosen). - * limit_soft_contact: For bilinear contact law, define the penetration limit between soft contact - (stiffness eps_a), and hard contact (stiffness eps_n) - only used if bilinear law is choosen. - + * eps_a: Penalty parameter for soft contact in bilinear contact law + (only used if bilinear law is choosen). + * limit_soft_contact: For bilinear contact law, define the + penetration limit between soft contact (stiffness eps_a), + and hard contact (stiffness eps_n) - only used if bilinear law is + choosen. """ if space is None: space = ModelingSpace.get_active() @@ -76,33 +74,10 @@ def __init__( self.search_algorithm = search_algorithm.lower() if self.search_algorithm == "bucket": - # by default bucket_size should be set to max sqrt(2)*max(edge_size) - if self.mesh.elements.shape[1] == 2: - max_edge_size = np.linalg.norm( - self.mesh.nodes[self.mesh.elements[:, 1]] - - self.mesh.nodes[self.mesh.elements[:, 0]], - axis=1, - ).max() - elif self.mesh.elements.shape[1] == 3: - max_edge_size = np.max( - [ - np.linalg.norm( - self.mesh.nodes[self.mesh.elements[:, ind[1]]] - - self.mesh.nodes[self.mesh.elements[:, ind[0]]], - axis=1, - ).max() - for ind in [[0, 1], [1, 2], [2, 0]] - ] - ) - else: - raise NotImplementedError( - f"'{self.mesh.elm_type}' elements are not " - "supported in contact. Use either 'tri3' " - "or 'lin2' elements." - ) - - self.bucket_size = np.sqrt(2) * max_edge_size + self.bucket_size = np.sqrt(2) * self._compute_max_edge_size() # self.bucket_size = self.mesh.bounding_box.size.max()/10 #bucket size #bounding box includes all nodes + elif self.search_algorithm == "ipctk": + self._setup_ipctk_broad_phase(space.ndim) elif search_algorithm.lower() != "search_nearest": raise (NameError, f"Search alogithm {search_algorithm} not unknown.") @@ -135,9 +110,6 @@ def __init__( self.sv = {} """ Dictionary of state variables associated to the associated for the current problem.""" self.sv_start = {} - # self.bc_type = 'Contact' - # BCBase.__init__(self, name) - # self._update_during_inc = 1 self.normal_law = normal_law.lower() """Name of the contact law. The normal law can be run using the "run_normal_law" method.""" @@ -151,9 +123,134 @@ def __init__( def global_search(self): if self.search_algorithm == "bucket": return self._nearest_node_bucket_sort() + elif self.search_algorithm == "ipctk": + return self._nearest_element_ipctk() else: #'search_nearest' return self._nearest_node() + def _compute_max_edge_size(self): + """Maximum edge length across all master surface elements.""" + elm_crd = self.mesh.nodes[self.mesh.elements] + n_per_elm = elm_crd.shape[1] + if n_per_elm == 2: # lin2 + return np.linalg.norm(elm_crd[:, 1] - elm_crd[:, 0], axis=1).max() + elif n_per_elm == 3: # tri3 + return max( + np.linalg.norm(elm_crd[:, j] - elm_crd[:, i], axis=1).max() + for i, j in [(0, 1), (1, 2), (2, 0)] + ) + else: + raise NotImplementedError( + f"'{self.mesh.elm_type}' elements are not supported in " + "contact. Use either 'tri3' or 'lin2' elements." + ) + + def _setup_ipctk_broad_phase(self, ndim): + """Build ipctk CollisionMesh for broad-phase contact search. + + Slave nodes already in master elements get codimensional + duplicates so ipctk still generates EV/FV candidates for them. + """ + try: + import ipctk + + self._ipctk = ipctk + except ImportError: + raise ImportError( + "The 'ipctk' package is required for search_algorithm='ipctk'. " + "Install it with: pip install ipctk" + ) + + # Compact vertex set: master ∪ slave + all_relevant = np.union1d(self.master_nodes, self.slave_nodes) + n_base = len(all_relevant) + + # Shared slaves need codimensional duplicates for EV/FV candidates + shared_mask = np.isin(self.slave_nodes, self.master_nodes) + shared_slaves = self.slave_nodes[shared_mask] + + # Global → compact index mapping + g2c = np.full(self.mesh.n_nodes, -1, dtype=int) + g2c[all_relevant] = np.arange(n_base) + + # Compact vertex ID for each slave (shared ones get duplicate IDs) + slave_cvids = g2c[self.slave_nodes].copy() + slave_cvids[shared_mask] = n_base + np.arange(shared_mask.sum()) + + # Reverse lookup: compact vertex ID → slave index (-1 = not slave) + n_verts = n_base + len(shared_slaves) + cvid_to_slave = np.full(n_verts, -1, dtype=int) + cvid_to_slave[slave_cvids] = np.arange(len(self.slave_nodes)) + + # For self-contact: map duplicate IDs back to element-space IDs + vid_to_orig = np.arange(n_verts, dtype=int) + if len(shared_slaves): + vid_to_orig[n_base:] = g2c[shared_slaves] + + self._ipctk_cvid_to_slave = cvid_to_slave + self._ipctk_vid_to_orig = vid_to_orig + self._ipctk_compact_elements = g2c[self.mesh.elements] + self._ipctk_all_relevant = all_relevant + self._ipctk_shared_slaves = shared_slaves + self._ipctk_ndim = ndim + + # Build CollisionMesh + pos = self._ipctk_build_vertices() + elems = self._ipctk_compact_elements + if ndim == 2: + edges, faces = elems, np.empty((0, 3), dtype=int) + else: + edges, faces = ipctk.edges(elems), elems + self._ipctk_cmesh = ipctk.CollisionMesh(pos, edges, faces) + self._ipctk_inflation = np.sqrt(2) * self._compute_max_edge_size() + + def _ipctk_build_vertices(self): + """Current vertex positions for ipctk (compact + shared dups, 3D, F-contiguous).""" + pos = self.mesh.nodes[self._ipctk_all_relevant] + if len(self._ipctk_shared_slaves): + pos = np.vstack([pos, self.mesh.nodes[self._ipctk_shared_slaves]]) + if pos.shape[1] == 2: + pos = np.hstack([pos, np.zeros((len(pos), 1))]) + return np.asfortranarray(pos, dtype=float) + + def _nearest_element_ipctk(self): + """Broad-phase search using ipctk HashGrid. + + Returns list of candidate master element indices per slave node. + """ + cands = self._ipctk.Candidates() + cands.build( + self._ipctk_cmesh, + self._ipctk_build_vertices(), + min(self.max_dist, self._ipctk_inflation), + broad_phase=self._ipctk.HashGrid(), + ) + + elems = self._ipctk_compact_elements + cvid_to_slave = self._ipctk_cvid_to_slave + vid_to_orig = self._ipctk_vid_to_orig + is_self = isinstance(self, SelfContact) + result = [[] for _ in range(len(self.slave_nodes))] + + if self._ipctk_ndim == 2: + for c in cands.ev_candidates: + si = cvid_to_slave[c.vertex_id] + if si < 0: + continue + if is_self and vid_to_orig[c.vertex_id] in elems[c.edge_id]: + continue + result[si].append(c.edge_id) + else: + for c in cands.fv_candidates: + si = cvid_to_slave[c.vertex_id] + if si < 0: + continue + if is_self and vid_to_orig[c.vertex_id] in elems[c.face_id]: + continue + result[si].append(c.face_id) + + return result + def assemble_global_mat(self, compute="all"): pass @@ -200,6 +297,10 @@ def contact_search(self, contact_list={}, update_contact=True): contact_elements = [] contact_g = [] + is_axi = self.space._dimension == "2Daxi" + if is_axi: + contact_r = [] # radial coordinate for 2πr weighting + data_Ns = [] if dim == 1: data_N0s = [] @@ -218,12 +319,8 @@ def contact_search(self, contact_list={}, update_contact=True): first_indices_disp ) # number of dof involved in a contact point - list_nodes = np.unique(surf.elements) - if update_contact: - t0 = time.time() list_possible_elements = self.global_search() - # print('temps: ', time.time()-t0) for i_nd, slave_node in enumerate(nodes): # if self.no_slip and slave_node in contact_list: @@ -405,6 +502,8 @@ def contact_search(self, contact_list={}, update_contact=True): contact_g.append(g) contact_elements.append(el) + if is_axi: + contact_r.append(surf.nodes[slave_node, 0]) # if vec_xi > 1: # shape_func_val = np.array([0.,1.]) @@ -525,7 +624,15 @@ def contact_search(self, contact_list={}, update_contact=True): # print('Warning, contact have been loosing') Fcontact, eps = self.run_normal_law(contact_g - self.clearance) - # print(Fcontact) + + # Apply 2πr circumferential weighting for axisymmetric problems + if is_axi: + axi_weight = 2 * np.pi * np.array(contact_r) + Fcontact = Fcontact * axi_weight + if np.isscalar(eps): + eps = eps * axi_weight + else: + eps = eps * axi_weight if dim == 1: F_div_l = sparse.diags( @@ -729,13 +836,9 @@ def update(self, pb, compute="all"): self.sv["contact_list"] = self.current.contact_search( self.sv["contact_list"], self.update_contact ) - # if self.update_contact: - # print(self.sv['contact_list']) if self.contact_search_once: self.update_contact = False - # self.current.assemble_global_mat(compute) - def run_normal_law(self, g): if self.normal_law == "linear": return self._linear_law(g) @@ -988,48 +1091,6 @@ def __init__( super().__init__(nodes, mesh, normal_law, search_algorithm, space, name) - # def global_search(self): - # """Find a list of elements that may be in contact for all slave nodes. - - # This method find the nearest neighbor of all slave nodes and return the - # associated elements if the nearest_neighbor is near enough - # (criterion givien by the max_dist attribute). - - # Returns - # ------- - # possible_elements: list[list[int]] - # possible_elements[i] is the list of indices of the master surface elements - # that may be in contact with the ith slave nodes - # are given in possible_elements[i]. If no element can be - # in contact, possible_elements[i] is an empty list. - # """ - # return self.nearest_node_self() - - # possible_elements = [] - # for id_nd,slave_node in enumerate(self.slave_nodes): - # dist_slave_nodes = np.linalg.norm(self.mesh.nodes[slave_node]-self.mesh.nodes[self.master_nodes], axis=1) - - # #get the indice with the minimal distance without considering the distance of slave_node between itself - # trial_nodes_indices = [] - # if id_nd != 0: - # trial_nodes_indices.append(dist_slave_nodes[:id_nd].argmin()) - - # if id_nd != len(self.master_nodes)-1: - # trial_nodes_indices.append(dist_slave_nodes[id_nd+1:].argmin() + id_nd+1) - - # trial_node_indice = trial_nodes_indices[np.argmin(dist_slave_nodes[trial_nodes_indices])] - - # if dist_slave_nodes[trial_node_indice] > self.max_dist: - # #to improve performance, ignore contact if distance to the closest node is to high - # possible_elements.append([]) - # else: - - # nearest_neighbors = self.master_nodes[trial_node_indice] - # possible_elements.append([el for el in range(self.mesh.n_elements) if nearest_neighbors in self.mesh.elements[el] and slave_node not in self.mesh.elements[el]]) - # # possible_elements.append([el for el in range(self.mesh.n_elements) if nearest_neighbors in self.mesh.elements[el]]) - - # return possible_elements - def _nearest_node(self): """Find a list of elements that may be in contact for all slave nodes. @@ -1235,177 +1296,3 @@ def _nearest_node_bucket_sort(self): ).nonzero()[0] for slave_nd, nn in zip(slave_nodes, nearest_neighbors) ] - - -# Have not been tested for now -class NodeContact(AssemblyBase): - """Class that define a node 2 node contact""" - - def __init__( - self, mesh, nodes1, nodes2, space=None, name="Contact nodes 2 surface" - ): - """ - In development. - """ - if space is None: - space = ModelingSpace.get_active() - AssemblyBase.__init__(self, name, space) - - self.nodes1 = nodes1 - self.nodes2 = nodes2 - self.mesh = mesh - - self.current = self - - self.eps_n = 1e7 - """ Contact penalty parameter. """ - - self.max_dist = 0.5 - """ Max distance from nodes at which contact is considered""" - - self.sv = {} - """ Dictionary of state variables associated to the associated for the current problem.""" - self.sv_start = {} - - def assemble_global_mat(self, compute="all"): - pass - - def contact_search(self): - nodes1 = self.nodes1 - nodes2 = self.nodes2 - if update_contact: - # update contact connection - new_contact_list = {} - - # look for contact - contact_nodes = [] - contact_g = [] - - # matrix for special treatment -> node 2 node if two nodes are close - # should be move in a node_2_node contact assembly - Xs = [] - indices_n2n = [] - data_n2n = [] - # end - - list_nodes = np.unique(surf.elements) - - for nd in nodes1: - dist_nodes = np.linalg.norm(mesh.nodes[nd] - mesh.nodes[nodes2], axis=1) - trial_node_indice = dist_nodes.argmin() - - if dist_nodes[trial_node_indice] >= self.max_dist: - # to improve performance, ignore contact if distance to the closest node is to high - continue - - # asses which elements are in contact - list_contact_nodes = np.where(dist_nodes < self.max_dist)[0] - - for nd2 in list_contact_nodes: - t = mesh.nodes[nd] - mesh.nodes[nd2] # x1-x2 - Xs.extend(t) - - # n2n_global_vector[[slave_node,slave_node+surf.n_nodes]] = -self.eps_nn*t - # n2n_global_vector[[master_node,master_node+surf.n_nodes]] = self.eps_nn*t - normal = t / np.linalg.norm(t) - - contact_nodes = np.array([nd, nd2]) - - sorted_indices = contact_nodes.argsort() - indices_n2n.extend( - list( - ( - contact_nodes[sorted_indices] - + np.array([[0], [mesh.n_nodes]]) - ).ravel() - ) - ) - data_n2n.extend(list(np.tile(np.array([-1, 1])[sorted_indices], 2))) - - M_n2n = sparse.csr_array( - (data_n2n, indices_n2n, np.arange(0, len(indices_n2n) + 1, 2)), - shape=(len(indices_n2n) // 2, self.space.nvar * surf.n_nodes), - ) - - # contact law -> put it in function - # Fc0 = self.eps_a*self.clearance - # eps = (contact_g <= 0) * self.eps_n + (contact_g > 0) * self.eps_a - # Fcontact = (-self.eps_n * contact_g + Fc0) * (contact_g <= 0) \ - # +(self.eps_a * (self.clearance - contact_g)) * (contact_g > 0) - - self.global_matrix = (-self.eps_n) * M_n2n.T @ M_n2n - - self.global_vector = self.eps_n * M_n2n.T @ np.array(Xs) - - # voir eq 9.35 et 9.36 (page 241 du pdf) avec def 9.18 et 9.19 page 239 - - def set_disp(self, disp): - if np.isscalar(disp) and disp == 0: - self.current = self - else: - new_crd = self.mesh.nodes + disp.T - if self.current == self: - # initialize a new - new_mesh = copy(self.mesh) - new_mesh.nodes = new_crd - new_assembly = copy(self) - new_assembly.mesh = new_mesh - self.current = new_assembly - else: - self.current.mesh.nodes = new_crd - - def initialize(self, pb): - # self.update(problem) - # initialize the contact list - self.current.contact_search() # initialize contact state - self.sv_start = dict(self.sv) - - def set_start(self, problem): - self.sv_start = dict( - self.sv - ) # create a new dict with alias inside (not deep copy) - - def to_start(self, pb): - self.sv = dict(self.sv_start) - self.set_disp(pb.get_disp()) - self.current.contact_search() # initialize global_matrix and global_vector - - def update(self, pb, compute="all"): - self.set_disp(pb.get_disp()) - self.current.contact_search() - - -# def global_search_bucket_sort(self): -# possible_elements = [] - -# #bucket_sort -# w = self.mesh.bounding_box.size.max()/10 #bucket size -# n_buckets = (self.mesh.bounding_box.size // w).astype(int) - -# temp = ((self.mesh.nodes[self.slave_nodes] - self.mesh.bounding_box[0])/w).astype(int) -# bucket_id = temp[:,0] + temp[:,1] * n_buckets[0] + temp[:,2] * n_buckets[0]*n_buckets[1] #only work in 3D -> improve - -# sorted_indices = bucket_id.argsort() - -# temp = [0] + [i for i in range(1,len(bucket_id)) if bucket_id[sorted_indices[i]] != bucket_id[sorted_indices[i-1]]] - -# bucket = {bucket_id[sorted_indices[temp[i]]]: sorted_indices[temp[i]:temp[i+1]] for i in range(0, len(temp)-1)} #bucket[2] contain the ith bucket -> perhaps best to include all bucket coordinates (list of list for instance) - -# neighbors_buckets = {i: (i+np.array([-1,0,1])+n_buckets[0]*np.array([-1,0,1]).reshape(-1,1) + n_buckets[0]*n_buckets[1]*np.array([-1,0,1]).reshape(-1,1,1)).flatten() for i in bucket} - -# for id_b in bucket: -# #closest_node technique: -# slave_node = bucket[id_b] -# master_nodes = sum([list(bucket.get(i,[])) for i in neighbors_buckets[id_b]],[]) - -# # dist_slave_nodes = np.linalg.norm(self.mesh.nodes[slave_node]-self.mesh.nodes[master_nodes], axis=1) -# dist_slave_nodes = self.mesh.nodes[slave_node] @ self.mesh.nodes[master_nodes].T -# trial_node_indice = dist_slave_nodes.argmin() -# if dist_slave_nodes[trial_node_indice] > self.max_dist: -# #to improve performance, ignore contact if distance to the closest node is to high -# possible_elements.append([]) -# else: - -# nearest_neighbors = self.master_nodes[trial_node_indice] -# # possible_elements.append([el for el in range(self.mesh.n_elements) if nearest_neighbors in self.mesh.elements[el] and slave_node not in self.mesh.elements[el]]) -# possible_elements.append([el for el in range(self.mesh.n_elements) if nearest_neighbors in self.mesh.elements[el]]) diff --git a/fedoo/constraint/ipc_contact.py b/fedoo/constraint/ipc_contact.py new file mode 100644 index 00000000..5c01c79f --- /dev/null +++ b/fedoo/constraint/ipc_contact.py @@ -0,0 +1,1201 @@ +from __future__ import annotations + +import numpy as np +from scipy import sparse +from fedoo.core.base import AssemblyBase +from fedoo.core.modelingspace import ModelingSpace + + +def _import_ipctk(): + """Lazy import of ipctk with a clear error message.""" + try: + import ipctk + + return ipctk + except ImportError: + raise ImportError( + "The 'ipctk' package is required for IPC contact. " + "Install it with: pip install ipctk\n" + "Or install fedoo with IPC support: pip install fedoo[ipc]" + ) + + +class IPCContact(AssemblyBase): + r"""Contact Assembly based on the IPC (Incremental Potential Contact) method. + + Uses the `ipctk `_ library to guarantee + intersection-free configurations through barrier potentials and optional + CCD (Continuous Collision Detection) line search. + + The IPC method adds a barrier energy term :math:`\kappa\,B(\mathbf{x})` + to the total potential energy. The barrier :math:`B` grows to infinity + as the distance between surfaces approaches zero and vanishes when the + distance exceeds ``dhat``. This guarantees that the deformed + configuration remains intersection-free at every Newton–Raphson + iteration. + + **Barrier stiffness auto-tuning** — When ``barrier_stiffness`` is left + to ``None`` (recommended), the stiffness :math:`\kappa` is automatically + initialised by ``ipctk.initial_barrier_stiffness`` so that it balances + the elastic and barrier energy gradients. If no contact exists at the + initial configuration (common for self-contact problems), :math:`\kappa` + is re-initialised the first time collisions appear, using the actual + elastic forces at that stage. Afterwards, :math:`\kappa` is updated + adaptively at each converged time increment via + ``ipctk.update_barrier_stiffness``. + + **Convergence safeguards (SDI)** — When the set of active collisions + changes during a Newton–Raphson iteration (new contacts appear or + existing ones separate), or when the minimum gap drops below + 0.1 × d_hat, the solver is forced to perform at least one additional + NR iteration even if the displacement criterion is already satisfied. + When surfaces are dangerously close (gap < 0.01 × d_hat), at least + three extra iterations are forced. This prevents premature + convergence when the contact state is still evolving. + + **CCD line search** — When ``use_ccd=True``, the Newton–Raphson step + size is clamped to the largest value that keeps all surfaces + intersection-free. For the elastic prediction of an increment where + no contacts existed at the start, a conservative minimum distance of + 0.1 × d_hat is enforced to prevent surfaces from jumping deep into + the barrier zone. If this conservative bound yields a zero step + (e.g. due to shared mesh edges at zero distance in self-contact), + the CCD falls back to the standard ``min_distance=0`` computation. + + **OGC trust-region** — When ``use_ogc=True``, the Offset Geometric + Contact (OGC) trust-region method from ``ipctk.ogc`` replaces CCD. + Instead of uniformly scaling the displacement step by a single + scalar, OGC filters the displacement per vertex, clamping each + vertex's move to lie within a trust region that guarantees no + intersection. This can provide better convergence by not + penalising distant vertices for a tight contact elsewhere. + ``use_ogc`` and ``use_ccd`` are mutually exclusive. + + .. warning:: + OGC is only supported for **shell / surface meshes** where every + node lies on the contact surface. Solid meshes with interior + nodes must use ``use_ccd=True`` instead. An error is raised at + initialization if ``use_ogc=True`` is used with a solid mesh. + + Parameters + ---------- + mesh : fedoo.Mesh + Full volumetric mesh (needed for DOF compatibility). + surface_mesh : fedoo.Mesh, optional + Pre-extracted surface mesh (``lin2`` in 2D, ``tri3`` in 3D). + extract_surface : bool, default=False + If ``True``, automatically extract the surface from the volumetric + mesh. + dhat : float, default=1e-3 + Barrier activation distance. Pairs of primitives closer than + ``dhat`` receive a repulsive barrier force. A smaller value gives + a tighter contact (less visible gap) but requires finer time steps. + dhat_is_relative : bool, default=True + If ``True``, ``dhat`` is interpreted as a fraction of the bounding + box diagonal of the surface mesh. + barrier_stiffness : float, optional + Barrier stiffness :math:`\kappa`. If ``None`` (default), it is + automatically computed and adaptively updated — **this is the + recommended setting**. When specified explicitly, the value is + kept fixed unless ``adaptive_barrier_stiffness`` is ``True``, in + which case it is used as the initial value and may be updated. + friction_coefficient : float, default=0.0 + Coulomb friction coefficient :math:`\mu`. Set to ``0`` for + frictionless contact. + eps_v : float, default=1e-3 + Friction smoothing velocity (only relevant when + ``friction_coefficient > 0``). + broad_phase : str, default="hash_grid" + Broad-phase collision detection method. One of ``"hash_grid"``, + ``"brute_force"``, ``"spatial_hash"`` or ``"bvh"``. + adaptive_barrier_stiffness : bool, default=True + If ``True``, :math:`\kappa` is updated adaptively at each converged + time increment using ``ipctk.update_barrier_stiffness``. + use_ccd : bool, default=False + Enable CCD (Continuous Collision Detection) line search. When + enabled, the Newton–Raphson step size is limited so that no + intersection can occur between iterations. Mutually exclusive + with ``use_ogc``. + line_search_energy : bool or None, default=None + Controls energy-based backtracking after CCD line search. The + step is halved until total energy (exact barrier + quadratic + elastic approximation) decreases. When ``None`` (default), + automatically enabled when ``use_ccd=True`` — this matches the + reference IPC algorithm. Set to ``False`` to disable (faster + but may degrade convergence). Has no effect when ``use_ogc`` + is enabled. + use_ogc : bool, default=False + Enable OGC (Offset Geometric Contact) trust-region step + filtering. Per-vertex displacement clamping replaces the + uniform scalar step of CCD. Mutually exclusive with + ``use_ccd``. **Only supported for shell / surface meshes** + where all nodes are on the surface; raises ``ValueError`` + for solid meshes with interior nodes. + space : ModelingSpace, optional + Modeling space. If ``None``, the active ``ModelingSpace`` is used. + name : str, default="IPC Contact" + Name of the contact assembly. + + Notes + ----- + The IPC assembly is used by adding it to the structural assembly with + :py:meth:`fedoo.Assembly.sum`: + + .. code-block:: python + + solid = fd.Assembly.create(wf, mesh) + ipc = fd.constraint.IPCContact(mesh, surface_mesh=surf) + assembly = fd.Assembly.sum(solid, ipc) + + When using ``add_output``, pass the *solid* assembly (not the sum) + to avoid errors. + + See Also + -------- + IPCSelfContact : Convenience wrapper for self-contact problems. + fedoo.constraint.Contact : Penalty-based contact method. + """ + + def __init__( + self, + mesh, + surface_mesh=None, + extract_surface=False, + dhat=1e-3, + dhat_is_relative=True, + barrier_stiffness=None, + friction_coefficient=0.0, + eps_v=1e-3, + broad_phase="hash_grid", + adaptive_barrier_stiffness=True, + use_ccd=False, + line_search_energy=None, + use_ogc=False, + space=None, + name="IPC Contact", + ): + if use_ccd and use_ogc: + raise ValueError( + "use_ccd and use_ogc are mutually exclusive. " + "Choose one line-search / step-filtering strategy." + ) + + if space is None: + space = ModelingSpace.get_active() + AssemblyBase.__init__(self, name, space) + + self.mesh = mesh + self._surface_mesh = surface_mesh + self._extract_surface = extract_surface + self._dhat = dhat + self._dhat_is_relative = dhat_is_relative + self._barrier_stiffness_init = barrier_stiffness + self.friction_coefficient = friction_coefficient + self._eps_v = eps_v + self._broad_phase_str = broad_phase + self._adaptive_barrier_stiffness = adaptive_barrier_stiffness + self._use_ccd = use_ccd + # Auto-enable energy backtracking when CCD is active (reference IPC + # algorithm always uses energy-based line search after CCD filtering) + if line_search_energy is None: + self._line_search_energy = use_ccd + else: + self._line_search_energy = line_search_energy + self._use_ogc = use_ogc + + self.current = self + + # Will be initialized in initialize() + self._collision_mesh = None + self._barrier_potential = None + self._friction_potential = None + self._collisions = None + self._friction_collisions = None + self._kappa = barrier_stiffness + self._max_kappa = None + self._max_kappa_set = False + self._actual_dhat = None + self._rest_positions = None + self._prev_min_distance = None + self._surface_node_indices = None + self._scatter_matrix = None # maps ipctk surface DOFs to fedoo full DOFs + self._broad_phase = None + + self._n_global_dof = 0 # extra DOFs from e.g. PeriodicBC + self._last_vertices = None # cached for assemble_global_mat + self._pb = None # reference to problem (for accessing elastic gradient) + self._n_collisions_at_start = 0 # collision count at start of increment + self._kappa_boosted_this_increment = False # prevent repeated within-NR boosts + self._bbox_diag = None # bounding box diagonal (cached) + + # OGC trust-region state + self._ogc_trust_region = None + self._ogc_vertices_at_start = None + + self.sv = {} + self.sv_start = {} + + def _create_broad_phase(self): + """Create an ipctk BroadPhase instance from the string name.""" + ipctk = _import_ipctk() + mapping = { + "hash_grid": ipctk.HashGrid, + "brute_force": ipctk.BruteForce, + "spatial_hash": ipctk.SpatialHash, + "bvh": ipctk.BVH, + } + cls = mapping.get(self._broad_phase_str.lower()) + if cls is None: + raise ValueError( + f"Unknown broad phase method '{self._broad_phase_str}'. " + f"Choose from: {list(mapping.keys())}" + ) + return cls() + + def _build_scatter_matrix(self, surface_node_indices, n_nodes, ndim): + """Build sparse matrix mapping ipctk surface DOFs to fedoo full DOFs. + + ipctk uses interleaved layout for surface nodes: + [x0, y0, z0, x1, y1, z1, ...] + fedoo uses blocked layout for all nodes: + [ux_0..ux_N, uy_0..uy_N, uz_0..uz_N] + + Parameters + ---------- + surface_node_indices : array + Global indices of surface nodes in the full mesh. + n_nodes : int + Total number of nodes in the full mesh. + ndim : int + Number of spatial dimensions (2 or 3). + + Returns + ------- + P : sparse csc_matrix + Shape (nvar * n_nodes, n_surf * ndim). + """ + disp_ranks = np.array(self.space.get_rank_vector("Disp")) + n_surf = len(surface_node_indices) + nvar = self.space.nvar + + rows = np.empty(n_surf * ndim, dtype=int) + cols = np.empty(n_surf * ndim, dtype=int) + + for d in range(ndim): + start = d * n_surf + end = (d + 1) * n_surf + # fedoo blocked row: disp_ranks[d] * n_nodes + global_node + rows[start:end] = disp_ranks[d] * n_nodes + surface_node_indices + # ipctk interleaved col: local_i * ndim + d + cols[start:end] = np.arange(n_surf) * ndim + d + + data = np.ones(n_surf * ndim) + P = sparse.csc_matrix( + (data, (rows, cols)), + shape=(nvar * n_nodes, n_surf * ndim), + ) + return P + + def _get_current_vertices(self, pb): + """Extract current surface vertex positions from problem displacement. + + Returns + ------- + vertices : ndarray, shape (n_surf_verts, ndim) + Current positions of surface vertices. + """ + disp = pb.get_disp() # shape (ndim, n_nodes) + + if np.isscalar(disp) and disp == 0: + return self._rest_positions.copy() + + # disp has shape (ndim, n_nodes) + surf_disp = disp[:, self._surface_node_indices].T # (n_surf, ndim) + return self._rest_positions + surf_disp + + def _get_elastic_gradient_on_surface(self): + """Extract elastic energy gradient projected onto surface DOFs. + + Uses the global vector from the parent assembly (excluding IPC) + to get the elastic force contribution, then projects it to + ipctk's interleaved surface DOF layout via P^T. + """ + if self._pb is None: + return None + + # Get the global vector from the parent assembly + # The parent assembly is the AssemblySum containing this IPC assembly + # We need the elastic part only (without IPC contribution) + assembly = self._pb.assembly + if not hasattr(assembly, "_list_assembly"): + return None + + # Sum global vectors from all assemblies except this one + grad_energy_full = None + for a in assembly._list_assembly: + if a is self: + continue + gv = a.current.get_global_vector() + if gv is not None and not (np.isscalar(gv) and gv == 0): + if grad_energy_full is None: + grad_energy_full = np.array(gv, dtype=float) + else: + grad_energy_full += gv + + if grad_energy_full is None: + return None + + # Truncate to mesh DOFs (remove global DOFs like PeriodicBC) + n_mesh_dof = self.space.nvar * self.mesh.n_nodes + grad_energy_full = grad_energy_full[:n_mesh_dof] + + # Project to surface DOFs: P^T @ full_grad -> surface_grad + return self._scatter_matrix.T @ grad_energy_full + + def _initialize_kappa(self, vertices): + """Compute barrier stiffness from gradient balance. + + Called each time step from ``set_start`` (when auto-kappa is + enabled and contacts exist) and from ``update`` (when contacts + first appear during an increment that started contactless). + + The MAX guard ensures kappa never decreases: only a larger + value from gradient balance is accepted. + """ + ipctk = _import_ipctk() + + bbox_diag = self._bbox_diag + + # Barrier gradient on surface DOFs + n_surf_dof = len(self._surface_node_indices) * self.space.ndim + if len(self._collisions) > 0: + grad_barrier = self._barrier_potential.gradient( + self._collisions, self._collision_mesh, vertices + ) + else: + grad_barrier = np.zeros(n_surf_dof) + + # Elastic gradient projected to surface DOFs + grad_energy = self._get_elastic_gradient_on_surface() + if grad_energy is None: + grad_energy = np.zeros(n_surf_dof) + + new_kappa, new_max_kappa = ipctk.initial_barrier_stiffness( + bbox_diag, + self._barrier_potential.barrier, + self._actual_dhat, + 1.0, # average_mass (1.0 for quasi-static) + grad_energy, + grad_barrier, + ) + # MAX guard: kappa can only increase to prevent oscillations + if self._kappa is None or new_kappa > self._kappa: + self._kappa = new_kappa + # max_kappa is only updated until _max_kappa_set is True + # (standard IPC: set once, then frozen as the adaptive cap) + if not self._max_kappa_set: + if self._max_kappa is None or new_max_kappa > self._max_kappa: + self._max_kappa = new_max_kappa + self.sv["kappa"] = self._kappa + self.sv["max_kappa"] = self._max_kappa + self.sv["prev_min_distance"] = self._prev_min_distance + + def _compute_ipc_contributions(self, vertices, compute="all"): + """Compute barrier (and friction) gradient and hessian. + + ipctk returns gradient/hessian in surface interleaved DOF space. + The scatter matrix P maps them to fedoo's full blocked DOF space: + global_vector = P @ ipctk_gradient + global_matrix = P @ ipctk_hessian @ P.T + + Sets self.global_matrix and self.global_vector. + """ + n_mesh_dof = self.space.nvar * self.mesh.n_nodes + n_dof = n_mesh_dof + self._n_global_dof + if len(self._collisions) == 0: + self.global_matrix = sparse.csr_array((n_dof, n_dof)) + self.global_vector = np.zeros(n_dof) + return + + ipctk = _import_ipctk() + P = self._scatter_matrix + + # Barrier contributions + if compute != "matrix": + grad_surf = self._barrier_potential.gradient( + self._collisions, self._collision_mesh, vertices + ) + # 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. + self.global_vector = -self._kappa * (P @ grad_surf) + + if compute != "vector": + try: + hess_surf = self._barrier_potential.hessian( + self._collisions, + self._collision_mesh, + vertices, + project_hessian_to_psd=ipctk.PSDProjectionMethod.CLAMP, + ) + except RuntimeError: + # CLAMP/ABS projection can fail on degenerate contact + # pairs (zero-distance); skip projection entirely. + hess_surf = self._barrier_potential.hessian( + self._collisions, + self._collision_mesh, + vertices, + project_hessian_to_psd=ipctk.PSDProjectionMethod.NONE, + ) + self.global_matrix = self._kappa * (P @ hess_surf @ P.T) + + # Friction contributions + if self.friction_coefficient > 0 and self._friction_collisions is not None: + if len(self._friction_collisions) > 0: + # Friction gradient/hessian from ipctk already include + # the effect of kappa and mu through the TangentialCollisions + # that were built with normal_stiffness and mu. + if compute != "matrix": + fric_grad_surf = self._friction_potential.gradient( + self._friction_collisions, + self._collision_mesh, + vertices, + ) + self.global_vector += -(P @ fric_grad_surf) + + if compute != "vector": + fric_hess_surf = self._friction_potential.hessian( + self._friction_collisions, + self._collision_mesh, + vertices, + project_hessian_to_psd=ipctk.PSDProjectionMethod.CLAMP, + ) + self.global_matrix += P @ fric_hess_surf @ P.T + + # Pad with zeros to account for extra global DOFs (e.g. PeriodicBC) + if self._n_global_dof > 0: + if compute != "matrix": + self.global_vector = np.pad(self.global_vector, (0, self._n_global_dof)) + if compute != "vector": + self.global_matrix = sparse.block_diag( + [ + self.global_matrix, + sparse.csr_array((self._n_global_dof, self._n_global_dof)), + ], + ) + + def _ccd_line_search(self, pb, dX): + """Compute step size using CCD + energy-based backtracking. + + **Phase 1 — CCD**: limits alpha to the largest collision-free + step. + + **Phase 2 — Energy backtracking**: starting from the CCD alpha, + halves the step until the total energy decreases. Barrier + energy is evaluated exactly from ipctk; elastic energy change + is approximated via a quadratic model (Lt * deps at the global + level) using the already-computed elastic residual and tangent + stiffness — no ``umat()`` call needed. + + Parameters + ---------- + pb : Problem + The current problem. + dX : ndarray + The displacement increment (full DOF vector). + + Returns + ------- + alpha : float + Step size in (0, 1] that is collision-free and decreases + total energy. + """ + ipctk = _import_ipctk() + ndim = self.space.ndim + + # Current vertices + vertices_current = self._get_current_vertices(pb) + + # Compute displacement of surface nodes from dX (fedoo blocked layout) + disp_ranks = np.array(self.space.get_rank_vector("Disp")) + n_nodes = self.mesh.n_nodes + surf_disp = np.zeros((len(self._surface_node_indices), ndim)) + for d in range(ndim): + surf_disp[:, d] = dX[disp_ranks[d] * n_nodes + self._surface_node_indices] + + vertices_next = vertices_current + surf_disp + + # --- Phase 1: CCD (collision-free step size) --- + # When no contacts at step start, use conservative min_distance + # to prevent jumping deep into the barrier zone. + # When contacts already exist, use min_distance=0 to avoid + # ipctk "initial distance <= d_min" warnings and double CCD. + if self._n_collisions_at_start == 0: + min_distance = 0.1 * self._actual_dhat + else: + min_distance = 0.0 + + alpha = ipctk.compute_collision_free_stepsize( + self._collision_mesh, + vertices_current, + vertices_next, + min_distance=min_distance, + broad_phase=self._broad_phase, + ) + + # Fallback: if min_distance caused alpha=0 (due to pre-existing + # zero-distance pairs from shared mesh edges), recompute with + # standard CCD (min_distance=0). + if alpha <= 0 and min_distance > 0: + alpha = ipctk.compute_collision_free_stepsize( + self._collision_mesh, + vertices_current, + vertices_next, + min_distance=0.0, + broad_phase=self._broad_phase, + ) + + if alpha < 1.0: + alpha *= 0.9 + + if alpha <= 0: + return alpha + + # --- Phase 2: Energy-based backtracking (opt-in) --- + if not self._line_search_energy: + return alpha + + # Skip during elastic prediction: the elastic residual D_e is + # near zero (previous step converged) and doesn't account for + # external work from the BC increment, so the quadratic model + # would incorrectly predict energy increase. CCD alone + # suffices for elastic prediction safety. + is_elastic_prediction = np.isscalar(pb._dU) and pb._dU == 0 + if is_elastic_prediction: + return alpha + + # Skip if no contacts exist and none expected at trial position + if len(self._collisions) == 0 and self._n_collisions_at_start == 0: + return alpha + + # Barrier energy at current position + E_barrier_current = self._kappa * self._barrier_potential( + self._collisions, self._collision_mesh, vertices_current + ) + + # Quadratic elastic energy approximation: + # ΔΠ_elastic ≈ α*c1 + 0.5*α²*c2 + # where c1 = -D_elastic·dX (gradient·step), + # c2 = dX^T · K_elastic · dX (curvature) + # Uses already-computed elastic residual and tangent from the + # current NR iteration — no umat() needed at trial points. + c1, c2 = 0.0, 0.0 + parent = self.associated_assembly_sum + if parent is not None: + for asm in parent._list_assembly: + if asm is not self: + D_e = asm.current.get_global_vector() + K_e = asm.current.get_global_matrix() + if D_e is not None and K_e is not None: + c1 = -np.dot(D_e, dX) + c2 = np.dot(dX, K_e @ dX) + break + + for _ls_iter in range(12): + vertices_trial = vertices_current + alpha * surf_disp + + # Rebuild collisions at trial position + trial_collisions = ipctk.NormalCollisions() + trial_collisions.build( + self._collision_mesh, + vertices_trial, + self._actual_dhat, + broad_phase=self._broad_phase, + ) + E_barrier_trial = self._kappa * self._barrier_potential( + trial_collisions, self._collision_mesh, vertices_trial + ) + + dE_elastic = alpha * c1 + 0.5 * alpha**2 * c2 + dE_total = dE_elastic + (E_barrier_trial - E_barrier_current) + + if dE_total <= 0: + break + alpha *= 0.5 + + return alpha + + # ------------------------------------------------------------------ + # OGC trust-region methods + # ------------------------------------------------------------------ + + def _ogc_warm_start(self, pb, dX): + """Filter the elastic prediction step using OGC warm-start. + + Converts the solver's displacement increment *dX* to surface + vertex positions, calls ``warm_start_time_step`` which moves + vertices toward the predicted position while respecting trust + regions, then writes the filtered displacement back into *dX*. + """ + ipctk = _import_ipctk() + ndim = self.space.ndim + disp_ranks = np.array(self.space.get_rank_vector("Disp")) + n_nodes = self.mesh.n_nodes + n_surf = len(self._surface_node_indices) + + # Surface displacement extracted from dX (fedoo blocked layout) + surf_disp = np.zeros((n_surf, ndim)) + for d in range(ndim): + surf_disp[:, d] = dX[disp_ranks[d] * n_nodes + self._surface_node_indices] + + # Predicted position = start-of-increment position + displacement + # ipctk requires F-contiguous arrays; x_start must also be writeable + x_start = np.asfortranarray(self._ogc_vertices_at_start.copy()) + pred_x = np.asfortranarray(self._ogc_vertices_at_start + surf_disp) + + # warm_start_time_step modifies x_start in-place toward pred_x + self._ogc_trust_region.warm_start_time_step( + self._collision_mesh, + x_start, + pred_x, + self._collisions, + self._actual_dhat, + broad_phase=self._broad_phase, + ) + + # Per-vertex OGC filtering (all nodes are on the surface) + filtered_disp = x_start - self._ogc_vertices_at_start + for d in range(ndim): + dX[disp_ranks[d] * n_nodes + self._surface_node_indices] = filtered_disp[ + :, d + ] + + def _ogc_filter_step(self, pb, dX): + """Filter a Newton–Raphson displacement step using OGC. + + Extracts the surface displacement from *dX*, calls + ``filter_step`` which clamps per-vertex moves in-place, + then writes the filtered values back into *dX*. + """ + ipctk = _import_ipctk() + ndim = self.space.ndim + disp_ranks = np.array(self.space.get_rank_vector("Disp")) + n_nodes = self.mesh.n_nodes + n_surf = len(self._surface_node_indices) + + # Current surface vertex positions (F-contiguous for ipctk) + vertices_current = np.asfortranarray(self._get_current_vertices(pb)) + + # Surface displacement from dX (F-contiguous + writeable for ipctk) + surf_dx = np.zeros((n_surf, ndim), order="F") + for d in range(ndim): + surf_dx[:, d] = dX[disp_ranks[d] * n_nodes + self._surface_node_indices] + + # Per-vertex OGC filtering (all nodes are on the surface) + self._ogc_trust_region.filter_step( + self._collision_mesh, vertices_current, surf_dx + ) + for d in range(ndim): + dX[disp_ranks[d] * n_nodes + self._surface_node_indices] = surf_dx[:, d] + + def _ogc_step_filter_callback(self, pb, dX, is_elastic_prediction=False): + """Dispatch to warm-start or NR filter depending on context.""" + if is_elastic_prediction: + self._ogc_warm_start(pb, dX) + else: + self._ogc_filter_step(pb, dX) + + def assemble_global_mat(self, compute="all"): + """Recompute barrier contributions from cached vertex positions. + + Called by the solver when the tangent matrix needs reassembly + without a full ``update()`` cycle. Uses vertex positions + cached by the most recent ``update()`` or ``set_start()`` call. + """ + if self._last_vertices is not None: + self._compute_ipc_contributions(self._last_vertices, compute) + + def initialize(self, pb): + """Initialize the IPC contact assembly. + + Called once by the solver at the start of the problem. Builds + the collision mesh, scatter matrix, and broad-phase instance. + Computes the initial barrier stiffness (or uses the user-provided + value). When ``use_ccd`` is ``True``, registers a CCD + line-search callback on the solver. + """ + ipctk = _import_ipctk() + ndim = self.space.ndim + + self._pb = pb + + # Track extra global DOFs (e.g. from PeriodicBC) + self._n_global_dof = pb.n_global_dof + + # Extract surface mesh if needed + if self._extract_surface and self._surface_mesh is None: + from fedoo.mesh import extract_surface as extract_surface_mesh + + self._surface_mesh = extract_surface_mesh(self.mesh) + + if self._surface_mesh is None: + raise ValueError( + "A surface mesh must be provided or extract_surface must be True." + ) + + surf_mesh = self._surface_mesh + + # Get surface node indices (global indices in the full mesh) + self._surface_node_indices = np.unique(surf_mesh.elements) + + # OGC per-vertex filtering is only valid for shell/surface meshes + # where every node is on the surface. Solid meshes with interior + # DOFs cannot be handled correctly by per-vertex clamping. + if self._use_ogc and len(self._surface_node_indices) != self.mesh.n_nodes: + raise ValueError( + "use_ogc=True requires a shell/surface mesh where all nodes " + "are on the surface. This mesh has interior nodes " + f"({self.mesh.n_nodes - len(self._surface_node_indices)} " + f"interior, {len(self._surface_node_indices)} surface). " + "Use use_ccd=True instead for solid meshes." + ) + + # Build mapping from global node indices to local surface ordering + global_to_local = np.full(self.mesh.n_nodes, -1, dtype=int) + global_to_local[self._surface_node_indices] = np.arange( + len(self._surface_node_indices) + ) + + # Remap surface elements to local ordering + local_elements = global_to_local[surf_mesh.elements] + + # Rest positions of surface vertices + self._rest_positions = self.mesh.nodes[self._surface_node_indices] + + # Build scatter matrix for DOF mapping (ipctk surface -> fedoo full) + self._scatter_matrix = self._build_scatter_matrix( + self._surface_node_indices, self.mesh.n_nodes, ndim + ) + + # Create broad phase instance + self._broad_phase = self._create_broad_phase() + + # Build edges and faces for the collision mesh + if ndim == 2: + edges = local_elements # lin2 elements are edges + faces = np.empty((0, 3), dtype=int) + else: # 3D + edges = ipctk.edges(local_elements) + faces = local_elements # tri3 elements are faces + + # Create CollisionMesh without displacement_map + # (displacement_map causes segfault in ipctk 1.5.0; + # DOF reordering is handled manually via scatter matrix) + self._collision_mesh = ipctk.CollisionMesh( + self._rest_positions, + edges, + faces, + ) + + # Cache bounding-box diagonal (used for dhat, kappa, energy line search) + self._bbox_diag = np.linalg.norm( + self._rest_positions.max(axis=0) - self._rest_positions.min(axis=0) + ) + + # Compute actual dhat + if self._dhat_is_relative: + self._actual_dhat = self._dhat * self._bbox_diag + else: + self._actual_dhat = self._dhat + + # Create barrier potential + self._barrier_potential = ipctk.BarrierPotential(self._actual_dhat) + + # Create friction potential if needed + if self.friction_coefficient > 0: + self._friction_potential = ipctk.FrictionPotential(self._eps_v) + + # Build initial collisions + self._collisions = ipctk.NormalCollisions() + vertices = self._get_current_vertices(pb) + self._collisions.build( + self._collision_mesh, + vertices, + self._actual_dhat, + broad_phase=self._broad_phase, + ) + + # Auto-compute or set barrier stiffness + if self._kappa is None: + self._initialize_kappa(vertices) + else: + self._max_kappa = 100.0 * self._kappa + + # Build friction collisions if needed + if self.friction_coefficient > 0: + self._friction_collisions = ipctk.TangentialCollisions() + self._friction_collisions.build( + self._collision_mesh, + vertices, + self._collisions, + self._barrier_potential, + self._kappa, + self.friction_coefficient, + ) + + # Store minimum distance + if len(self._collisions) > 0: + self._prev_min_distance = self._collisions.compute_minimum_distance( + self._collision_mesh, vertices + ) + else: + self._prev_min_distance = np.inf + + # Compute initial contributions + self._compute_ipc_contributions(vertices) + + # Store initial state + self.sv["kappa"] = self._kappa + self.sv["max_kappa"] = self._max_kappa + self.sv["prev_min_distance"] = self._prev_min_distance + self.sv["max_kappa_set"] = self._max_kappa_set + self.sv_start = dict(self.sv) + + # Register line-search / step-filter callback + if self._use_ccd: + pb._step_size_callback = self._ccd_line_search + elif self._use_ogc: + self._ogc_trust_region = ipctk.ogc.TrustRegion(self._actual_dhat) + pb._step_filter_callback = self._ogc_step_filter_callback + + def _update_kappa_adaptive(self, vertices): + """Double kappa when the gap is small and decreasing. + + Uses ``ipctk.update_barrier_stiffness`` with + ``dhat_epsilon_scale = dhat / bbox_diag`` so the doubling + triggers within the barrier zone. Only called between time + steps (in ``set_start``), never inside the NR loop. + """ + ipctk = _import_ipctk() + bbox_diag = self._bbox_diag + + if len(self._collisions) > 0: + min_dist = self._collisions.compute_minimum_distance( + self._collision_mesh, vertices + ) + else: + min_dist = np.inf + + eps_scale = self._actual_dhat / bbox_diag + new_kappa = ipctk.update_barrier_stiffness( + self._prev_min_distance, + min_dist, + self._max_kappa, + self._kappa, + bbox_diag, + dhat_epsilon_scale=eps_scale, + ) + if new_kappa > self._kappa: + self._kappa = new_kappa + self.sv["kappa"] = self._kappa + + self._prev_min_distance = min_dist + self.sv["max_kappa"] = self._max_kappa + self.sv["prev_min_distance"] = self._prev_min_distance + + def set_start(self, pb): + """Begin a new time increment. + + Kappa management follows the reference IPC algorithm: + + 1. **Re-initialization** from gradient balance at each time + step when contacts exist. The MAX guard in + ``_initialize_kappa`` prevents kappa from decreasing. + + 2. **Adaptive doubling** — ``ipctk.update_barrier_stiffness`` + doubles kappa when the gap is small and decreasing. This + runs between time steps only (never within the NR loop, + as changing kappa mid-solve prevents NR convergence). + """ + self.sv_start = dict(self.sv) + + vertices = self._get_current_vertices(pb) + + # Rebuild collisions + self._collisions.build( + self._collision_mesh, + vertices, + self._actual_dhat, + broad_phase=self._broad_phase, + ) + + # Re-initialize kappa from gradient balance each time step + # when contacts exist. The MAX guard in _initialize_kappa + # prevents kappa from decreasing. max_kappa is set once + # (the first time contacts appear) and kept fixed so that it + # actually caps the adaptive doubling. + if self._barrier_stiffness_init is None and len(self._collisions) > 0: + self._initialize_kappa(vertices) + if not self._max_kappa_set: + self._max_kappa_set = True + + # Adaptive doubling when gap is small and decreasing + if self._adaptive_barrier_stiffness and self._max_kappa is not None: + self._update_kappa_adaptive(vertices) + + # Track collision count for detecting new contacts during NR + self._n_collisions_at_start = len(self._collisions) + self._kappa_boosted_this_increment = False + + # Store last vertices for next increment + self._last_vertices = vertices + + # Store start-of-increment vertices for OGC warm-start + if self._use_ogc: + self._ogc_vertices_at_start = vertices.copy() + + # Compute IPC contributions for the elastic prediction + self._compute_ipc_contributions(vertices) + + def update(self, pb, compute="all"): + """Update the IPC assembly for the current Newton–Raphson iteration. + + Rebuilds the collision set from the current vertex positions, + re-initialises kappa if contacts appeared for the first time in + an increment that started contact-free, enforces SDI (minimum + NR sub-iterations) when the contact state changed or surfaces + are dangerously close, and recomputes the barrier contributions. + + Parameters + ---------- + pb : NonLinear + The current problem instance. + compute : str, default "all" + What to compute: ``"all"``, ``"matrix"``, or ``"vector"``. + """ + vertices = self._get_current_vertices(pb) + self._last_vertices = vertices + + # Rebuild collision set + self._collisions.build( + self._collision_mesh, + vertices, + self._actual_dhat, + broad_phase=self._broad_phase, + ) + + # Initialize kappa when contacts first appear during an + # increment that started with zero contacts. MAX guard in + # _initialize_kappa ensures kappa never decreases. + if ( + not self._kappa_boosted_this_increment + and self._n_collisions_at_start == 0 + and self._barrier_stiffness_init is None + and len(self._collisions) > 0 + ): + self._initialize_kappa(vertices) + self._kappa_boosted_this_increment = True + # Set generous max_kappa ceiling + self._max_kappa = max(self._max_kappa, 1e4 * max(self._kappa, 1.0)) + + # SDI: force at least one NR correction when contact state + # changed (collision count differs from start) or when surfaces + # are dangerously close (min_d/dhat < 0.1). + n_collisions_now = len(self._collisions) + need_sdi = n_collisions_now != self._n_collisions_at_start + min_d_val = None + if not need_sdi and n_collisions_now > 0: + min_d_val = self._collisions.compute_minimum_distance( + self._collision_mesh, vertices + ) + if min_d_val < 0.1 * self._actual_dhat: + need_sdi = True + elif n_collisions_now > 0: + min_d_val = self._collisions.compute_minimum_distance( + self._collision_mesh, vertices + ) + if need_sdi: + # Scale min_subiter by proximity: more iterations when closer + if min_d_val is not None and min_d_val < 0.01 * self._actual_dhat: + self._pb._nr_min_subiter = max(self._pb._nr_min_subiter, 3) + + # Always force at least 1 NR iteration when IPC is active. + # The Displacement criterion's reference error grows with + # accumulated displacement, making the tolerance progressively + # looser. Without this floor, many steps are accepted at + # iter 0 (elastic prediction only), accumulating equilibrium + # errors that cause stress oscillations. + self._pb._nr_min_subiter = max(self._pb._nr_min_subiter, 1) + + # Build friction collisions if enabled + if self.friction_coefficient > 0: + self._friction_collisions.build( + self._collision_mesh, + vertices, + self._collisions, + self._barrier_potential, + self._kappa, + self.friction_coefficient, + ) + + # Update OGC trust regions after rebuilding collisions + if self._use_ogc and self._ogc_trust_region is not None: + self._ogc_trust_region.update_if_needed( + self._collision_mesh, + np.asfortranarray(vertices), + self._collisions, + broad_phase=self._broad_phase, + ) + + # Compute contributions + self._compute_ipc_contributions(vertices, compute) + + def to_start(self, pb): + """Restart the current time increment after NR failure. + + Restores kappa, max_kappa and prev_min_distance from saved + state, rebuilds collisions from the restored displacement, + and recomputes barrier contributions. + """ + self.sv = dict(self.sv_start) + self._kappa = self.sv_start.get("kappa", self._kappa) + self._max_kappa = self.sv_start.get("max_kappa", self._max_kappa) + self._prev_min_distance = self.sv_start.get( + "prev_min_distance", self._prev_min_distance + ) + self._max_kappa_set = self.sv_start.get("max_kappa_set", self._max_kappa_set) + self._kappa_boosted_this_increment = False + + # Recompute from current displacement state + vertices = self._get_current_vertices(pb) + self._last_vertices = vertices + + self._collisions.build( + self._collision_mesh, + vertices, + self._actual_dhat, + broad_phase=self._broad_phase, + ) + + if self.friction_coefficient > 0: + self._friction_collisions.build( + self._collision_mesh, + vertices, + self._collisions, + self._barrier_potential, + self._kappa, + self.friction_coefficient, + ) + + self._compute_ipc_contributions(vertices) + + # Reset OGC trust region on NR failure + if self._use_ogc: + self._ogc_vertices_at_start = vertices.copy() + self._ogc_trust_region = _import_ipctk().ogc.TrustRegion(self._actual_dhat) + + +class IPCSelfContact(IPCContact): + r"""Self-contact Assembly using the IPC method. + + Convenience wrapper around :py:class:`IPCContact` that automatically + extracts the surface mesh from the volumetric mesh. Useful for + problems where a single body may come into contact with itself (e.g. + buckling, folding, compression of porous structures). + + All barrier stiffness tuning is automatic (see :py:class:`IPCContact` + for details). + + Parameters + ---------- + mesh : fedoo.Mesh + Full volumetric mesh. + dhat : float, default=1e-3 + Barrier activation distance (relative to bounding box diagonal by + default, see ``dhat_is_relative``). + dhat_is_relative : bool, default=True + If ``True``, ``dhat`` is a fraction of the bounding box diagonal. + barrier_stiffness : float, optional + Barrier stiffness :math:`\kappa`. ``None`` (default) for automatic + computation and adaptive update — **recommended**. + friction_coefficient : float, default=0.0 + Coulomb friction coefficient :math:`\mu`. + eps_v : float, default=1e-3 + Friction smoothing velocity. + broad_phase : str, default="hash_grid" + Broad-phase collision detection method. + adaptive_barrier_stiffness : bool, default=True + Adaptively update :math:`\kappa` at each converged time increment. + use_ccd : bool, default=False + Enable CCD line search for robustness. Mutually exclusive + with ``use_ogc``. + line_search_energy : bool or None, default=None + Energy-based backtracking after CCD. ``None`` auto-enables + when ``use_ccd=True`` (see :py:class:`IPCContact`). + use_ogc : bool, default=False + Enable OGC trust-region step filtering. Mutually exclusive + with ``use_ccd``. **Only supported for shell / surface + meshes**; raises ``ValueError`` for solid meshes. + space : ModelingSpace, optional + Modeling space. + name : str, default="IPC Self Contact" + Name of the contact assembly. + + Examples + -------- + .. code-block:: python + + import fedoo as fd + + fd.ModelingSpace("3D") + mesh = fd.Mesh.read("gyroid.vtk") + material = fd.constitutivelaw.ElasticIsotrop(1e5, 0.3) + + contact = fd.constraint.IPCSelfContact(mesh, use_ccd=True) + + wf = fd.weakform.StressEquilibrium(material, nlgeom="UL") + solid = fd.Assembly.create(wf, mesh) + assembly = fd.Assembly.sum(solid, contact) + + pb = fd.problem.NonLinear(assembly) + res = pb.add_output("results", solid, ["Disp", "Stress"]) + # ... add BCs and solve ... + + See Also + -------- + IPCContact : Base class with full parameter documentation. + """ + + def __init__( + self, + mesh, + dhat=1e-3, + dhat_is_relative=True, + barrier_stiffness=None, + friction_coefficient=0.0, + eps_v=1e-3, + broad_phase="hash_grid", + adaptive_barrier_stiffness=True, + use_ccd=False, + line_search_energy=None, + use_ogc=False, + space=None, + name="IPC Self Contact", + ): + super().__init__( + mesh=mesh, + surface_mesh=None, + extract_surface=True, + dhat=dhat, + dhat_is_relative=dhat_is_relative, + barrier_stiffness=barrier_stiffness, + friction_coefficient=friction_coefficient, + eps_v=eps_v, + broad_phase=broad_phase, + adaptive_barrier_stiffness=adaptive_barrier_stiffness, + use_ccd=use_ccd, + line_search_energy=line_search_energy, + use_ogc=use_ogc, + space=space, + name=name, + ) diff --git a/fedoo/core/problem.py b/fedoo/core/problem.py index af91111f..56da0d6e 100644 --- a/fedoo/core/problem.py +++ b/fedoo/core/problem.py @@ -62,7 +62,7 @@ def __init__(self, A=None, B=0, D=0, mesh=None, name="MainProblem", space=None): self.__X = 0 # np.ndarray( self.n_dof ) #empty array self._Xbc = 0 - self._dof_slave = np.array([]) + self._dof_slave = np.array([]) # dof from dirichlet bc or mpc self._dof_free = np.array([]) # prepering output demand to export results @@ -340,8 +340,8 @@ def apply_boundary_conditions(self, t_fact=1, t_fact_old=None): F = np.zeros(n_dof) build_mpc = False - dof_blocked = set() # only dirichlet bc - dof_slave = set() + dof_blocked = set() # only dirichlet bc dof + dof_slave = set() # dirichlet + mpc elminated dof data = [] row = [] col = [] @@ -399,7 +399,7 @@ def apply_boundary_conditions(self, t_fact=1, t_fact_old=None): # modification col numbering from dof_free to np.arange(len(dof_free)) changeInd = np.full( n_dof, np.nan - ) # mettre des nan plutôt que des zeros pour générer une erreur si pb + ) # nan used instead of 0 to generate an error changeInd[dof_free] = np.arange(len(dof_free)) col = changeInd[np.hstack(col)] @@ -409,13 +409,12 @@ def apply_boundary_conditions(self, t_fact=1, t_fact_old=None): row = row[mask] data = data[mask] - # self._MFext = M.tocsr().T self._MFext = M self._dof_blocked = dof_blocked else: self._MFext = None - # #adding identity for free nodes + # adding identity for free nodes col = np.hstack( (col, np.arange(len(dof_free))) ) # np.hstack((col,dof_free)) #col.append(dof_free) @@ -502,9 +501,8 @@ def get_ext_forces(self, name="all", include_mpc=True): col = np.hstack((M.col, dof_idt)) row = np.hstack((M.row, dof_idt)) data = np.hstack((M.data, np.ones(len(dof_idt)))) - M = M.tocsr().T - # need to be checked - # M = self._MFext + scipy.sparse.identity(self.n_dof, dtype='d') + M = sparse.coo_matrix((data, (row, col)), shape=M.shape).tocsr().T + # M = M.tocsr().T if np.isscalar(self.get_D()) and self.get_D() == 0: return self._get_vect_component(M @ self.get_A() @ self.get_X(), name) else: diff --git a/fedoo/mesh/importmesh.py b/fedoo/mesh/importmesh.py index 322690d5..4cd517ea 100644 --- a/fedoo/mesh/importmesh.py +++ b/fedoo/mesh/importmesh.py @@ -263,9 +263,9 @@ def import_msh( for i in range(NbEntityBlocks): l = msh.pop(0).split() - entityDim = l[ - 0 - ] # 0 for point, 1 for curve, 2 for surface, 3 for volume + entityDim = int( + l[0] + ) # 0 for point, 1 for curve, 2 for surface, 3 for volume entityTag = l[1] # id of the entity elementType = l[2] numElementsInBlock = int(l[3]) @@ -306,7 +306,7 @@ def import_msh( if entityTag in Entities: for physicalTag in Entities[entityTag]: if physicalTag in PhysicalNames: - el_name = PhysicalNames["physicalTag"] + el_name = PhysicalNames[physicalTag].strip('"') else: el_name = "PhysicalPart" + physicalTag diff --git a/fedoo/problem/non_linear.py b/fedoo/problem/non_linear.py index 9b1ad3b9..23f26b9b 100644 --- a/fedoo/problem/non_linear.py +++ b/fedoo/problem/non_linear.py @@ -25,7 +25,8 @@ def __init__(self, assembly, nlgeom=False, name="MainProblem"): "err0": None, # default error for NR error estimation "criterion": "Displacement", "tol": 1e-3, - "max_subiter": 5, + "max_subiter": 10, + "dt_increase_niter": None, "norm_type": 2, } """ @@ -36,12 +37,24 @@ def __init__(self, assembly, nlgeom=False, name="MainProblem"): ['Displacement', 'Force', 'Work']. Default is 'Displacement'. * 'tol': Error tolerance for convergence. Default is 1e-3. - * 'max_subiter': Number of nr iteration before returning a - convergence error. Default is 5. + * 'max_subiter': Number of nr iterations before returning a + convergence error. Default is 6. + * 'dt_increase_niter': number of nr iterations threshold that + define an easy convergence. In problem allowing automatic + convergence, if the Newton–Raphson loop converges in strictly + fewer iterations, the time step is increased. + If None, the value is initialized to ``max_subiter // 3``. * 'norm_type': define the norm used to test the criterion Use numpy.inf for the max value. Default is 2. """ + self._step_size_callback = None + self._step_filter_callback = None # OGC per-vertex filter + self._nr_min_subiter = ( + 0 # SDI: minimum NR sub-iterations before accepting convergence + ) + self._t_fact_inc = None # frozen t_fact for the current increment + self.__assembly = assembly super().__init__(A, B, D, assembly.mesh, name, assembly.space) self.nlgeom = nlgeom @@ -131,6 +144,17 @@ def elastic_prediction(self): ) # not modified in principle if dt is not modified, except the very first iteration. May be optimized by testing the change of dt self.solve() + # OGC per-vertex filter or CCD scalar line search + if self._step_filter_callback is not None: + dX = self.get_X() + self._step_filter_callback(self, dX, is_elastic_prediction=True) + elif self._step_size_callback is not None: + dX = self.get_X() + alpha = self._step_size_callback(self, dX) + if alpha < 1.0: + # Scale only free DOFs; preserve prescribed Dirichlet values + self.set_X(dX * alpha + self._Xbc * (1 - alpha)) + # set the increment Dirichlet boundray conditions to 0 (i.e. will not change during the NR interations) try: self._Xbc *= 0 @@ -142,6 +166,7 @@ def elastic_prediction(self): def set_start(self, save_results=False, callback=None): # dt not used for static problem + self._nr_min_subiter = 0 # reset SDI for new increment if not (np.isscalar(self._dU) and self._dU == 0): self._U += self._dU self._dU = 0 @@ -166,12 +191,22 @@ def set_start(self, save_results=False, callback=None): def to_start(self): self._dU = 0 + self._nr_min_subiter = 0 + self._t_fact_inc = None self._err0 = self.nr_parameters["err0"] # initial error for NR error estimation self.__assembly.to_start(self) def NewtonRaphsonIncrement(self): # solve and update total displacement. A and D should up to date self.solve() + if self._step_filter_callback is not None: + dX = self.get_X() + self._step_filter_callback(self, dX, is_elastic_prediction=False) + elif self._step_size_callback is not None: + dX = self.get_X() + alpha = self._step_size_callback(self, dX) + if alpha < 1.0: + self.set_X(dX * alpha) self._dU += self.get_X() def update(self, compute="all", updateWeakForm=True): @@ -199,11 +234,11 @@ def reset(self): # self.set_A(self.__assembly.current.get_global_matrix()) #tangent stiffness # self.set_D(self.__assembly.current.get_global_vector()) - B = 0 self._U = 0 self._dU = 0 self._err0 = self.nr_parameters["err0"] # initial error for NR error estimation + self._t_fact_inc = None self.t0 = 0 self.tmax = 1 self.__iter = 0 @@ -229,34 +264,46 @@ def NewtonRaphsonError(self): dof_free = self._dof_free if len(dof_free) == 0: return 0 - if self._err0 is None: # if self._err0 is None -> initialize the value of err0 - # if self.nr_parameters["criterion"] == "Displacement": - # self._err0 = np.linalg.norm( - # (self._U + self._dU)[dof_free], norm_type - # ) # Displacement criterion - # if self._err0 == 0: - # self._err0 = 1 - # return 1 - # return np.max(np.abs(self.get_X()[dof_free])) / self._err0 - # else: - # self._err0 = 1 - # self._err0 = self.NewtonRaphsonError() - # return 1 + if self._err0 is None: # assess err0 from current state if self.nr_parameters["criterion"] == "Displacement": - err0 = np.linalg.norm((self._U + self._dU), norm_type) + # Normalize by the current increment (not the total + # accumulated displacement) so the criterion does not + # become progressively looser as loading proceeds. + err0 = np.linalg.norm(self._dU, norm_type) + # if not np.array_equal(self._U, 0): + # err0 += 0.01 * np.linalg.norm(self._U, norm_type) + # err0 += 1e-8*np.max(self.mesh.bounding_box.size) if err0 == 0: err0 = 1 - return np.max(np.abs(self.get_X()[dof_free])) / err0 - else: + return 1 + return np.linalg.norm(self.get_X()[dof_free], norm_type) / err0 + elif self.nr_parameters["criterion"] == "Force": + # Normalize by external force + err0 = np.linalg.norm( + self.get_ext_forces(include_mpc=False), + norm_type, + ) + # err0 += 1e-8 # to avoid divizion by 0 + if err0 == 0: + err0 = 1 + return 1 + if np.isscalar(self.get_D()) and self.get_D() == 0: + return np.linalg.norm(self.get_B()[dof_free], norm_type) / err0 + else: + return ( + np.linalg.norm( + self.get_B()[dof_free] + self.get_D()[dof_free], + norm_type, + ) + / err0 + ) + else: # initialize the value of self._err0 self._err0 = 1 self._err0 = self.NewtonRaphsonError() return 1 else: if self.nr_parameters["criterion"] == "Displacement": - # return np.max(np.abs(self.get_X()[dof_free]))/self._err0 #Displacement criterion - return ( - np.linalg.norm(self.get_X()[dof_free], norm_type) / self._err0 - ) # Displacement criterion + return np.linalg.norm(self.get_X()[dof_free], norm_type) / self._err0 elif self.nr_parameters["criterion"] == "Force": # Force criterion if np.isscalar(self.get_D()) and self.get_D() == 0: return ( @@ -306,8 +353,13 @@ def set_nr_criterion(self, criterion="Displacement", **kargs): If None (default), err0 is automatically computed. * 'tol': float, default is 1e-3. Error tolerance for convergence. - * 'max_subiter': int, default = 5. + * 'max_subiter': int, default = 10. Number of nr iteration before returning a convergence error. + * 'dt_increase_niter': number of nr iterations threshold that + define an easy convergence. In problem allowing automatic + convergence, if the Newton–Raphson loop converges in strictly + fewer iterations, the time step is increased. + Default is set to max_subiter//3. * 'norm_type': int or numpy.inf, default = 2. Define the norm used to test the criterion """ @@ -318,10 +370,17 @@ def set_nr_criterion(self, criterion="Displacement", **kargs): self.nr_parameters["criterion"] = criterion for key in kargs: - if key not in ["err0", "tol", "max_subiter", "norm_type"]: + if key not in [ + "err0", + "tol", + "max_subiter", + "dt_increase_niter", + "norm_type", + ]: raise NameError( "Newton Raphson parameters should be in " - "['err0', 'tol', 'max_subiter', 'norm_type']" + "['err0', 'tol', 'max_subiter', 'dt_increase_niter', " + "'norm_type']" ) self.nr_parameters[key] = kargs[key] @@ -332,6 +391,12 @@ def solve_time_increment(self, max_subiter=None, tol_nr=None): if tol_nr is None: tol_nr = self.nr_parameters["tol"] + # Freeze t_fact for the duration of this increment so that + # any modification of self.dtime (e.g. by a CCD callback) + # does not change the target time seen by boundary conditions + # or subiter print output. + self._t_fact_inc = self.t_fact + self.elastic_prediction() for subiter in range(max_subiter): # newton-raphson iterations # update Stress and initial displacement and Update stiffness matrix @@ -348,13 +413,14 @@ def solve_time_increment(self, max_subiter=None, tol_nr=None): ) ) - if normRes < tol_nr: # convergence of the NR algorithm + if ( + normRes < tol_nr and subiter >= self._nr_min_subiter + ): # convergence of the NR algorithm # Initialize the next increment + self._t_fact_inc = None return 1, subiter, normRes # --------------- Solve -------------------------------------------------------- - # self.__Assembly.current.assemble_global_mat(compute = 'matrix') - # self.set_A(self.__Assembly.current.get_global_matrix()) self.update( compute="matrix", updateWeakForm=False ) # assemble the tangeant matrix @@ -362,6 +428,7 @@ def solve_time_increment(self, max_subiter=None, tol_nr=None): self.NewtonRaphsonIncrement() + self._t_fact_inc = None return 0, subiter, normRes def nlsolve( @@ -372,6 +439,7 @@ def nlsolve( t0: float | None = None, dt_min: float = 1e-6, max_subiter: int | None = None, + dt_increase_niter: int | None = None, tol_nr: float | None = None, print_info: int | None = None, save_at_exact_time: bool | None = None, @@ -388,7 +456,8 @@ def nlsolve( update_dt: bool, default = True If True, the time increment may be modified during resolution: * decrease if the solver has not converged - * increase if the solver has converged in one iteration. + * increase if the solver has converged quickly (see + ``dt_increase_niter``). tmax: float, optional. Time at the end of the time step. If omitted, the attribute tmax is considered (default = 1.) @@ -400,9 +469,21 @@ def nlsolve( dt_min: float, default = 1e-6 Minimal time increment max_subiter: int, optional - Maximal number of newton raphson iteration at for each time increment. - If omitted, the 'max_subiter' field in the nr_parameters attribute - (ie nr_parameters['max_subiter']) is considered (default = 5). + Maximal number of newton raphson iteration allowed for each time + increment, after the initial linear guess. + If omitted, the 'max_subiter' field in the nr_parameters + attribute (ie nr_parameters['max_subiter']) is considered + (default = 10). + dt_increase_niter: int, optional + When ``update_dt`` is ``True``, the time increment is multiplied + by 1.25 if the Newton–Raphson loop converges in strictly fewer + than ``dt_increase_niter`` iterations. If omitted, the + 'dt_increase_niter' field in the nr_parameters attribute + (ie nr_parameters['dt_increase_niter']) is considered + (default = ``max_subiter // 3``). + For contact problems where NR typically needs several iterations, + a higher value (e.g. ``max_subiter // 2``) allows dt to recover + after an earlier reduction. tol_nr: float, optional Tolerance of the newton-raphson algorithm. If omitted, the 'tol' field in the nr_parameters attribute @@ -414,23 +495,27 @@ def nlsolve( If 2, iterations and newton-raphson sub iterations info are printed. If omitted, the print_info attribute is considered (default = 1). save_at_exact_time: bool, optional - If True, the time increment is modified to stop at times defined by interval_output and allow to save results. - If omitted, the save_at_exact_time attribute is considered (default = True). + If True, the time increment is modified to stop at times defined by + interval_output and allow to save results. If omitted, the + save_at_exact_time attribute is considered (default = True). The given value is stored in the save_at_exact_time attribute. interval_output: int|float, optional - Time step for output if save_at_exact_time is True (default) else number of iter increments between 2 output - If interval_output == -1, the results is saved at each initial time_step intervals or each increment depending on the save_at_exact_time value. - If omitted, the interval_output attribute is considred (default -1) + Time step for output if save_at_exact_time is True (default) else + number of iter increments between 2 output. If + interval_output == -1, the results is saved at each initial + time_step intervals or each increment depending on the + save_at_exact_time value. If omitted, the interval_output attribute + is considred (default -1) callback: function, optional - The callback function is executed automatically during the non linear resolution. - By default, the callback function is executed when output is requested - (defined by the interval_output argument). If - exec_callback_at_each_iter is True, the callback function is excuted - at each time iteration. + The callback function is executed automatically during the non + linear resolution. By default, the callback function is executed + when output is requested (defined by the interval_output argument). + If exec_callback_at_each_iter is True, the callback function is + excuted at each time iteration. exec_callback_at_each_iter, bool, default = False - If True, the callback function is executed after each time iteration. + If True, the callback function is executed after each time + iteration. """ - # parameters if tmax is not None: self.tmax = tmax @@ -438,6 +523,10 @@ def nlsolve( self.t0 = t0 # time at the start of the time step if max_subiter is None: max_subiter = self.nr_parameters["max_subiter"] + if dt_increase_niter is None: + dt_increase_niter = self.nr_parameters["dt_increase_niter"] + if dt_increase_niter is None: + dt_increase_niter = max_subiter // 3 if tol_nr is None: tol_nr = self.nr_parameters["tol"] if print_info is not None: @@ -515,9 +604,8 @@ def nlsolve( ) # Check if dt can be increased - if update_dt and nbNRiter < 2 and dt == self.dtime: + if update_dt and nbNRiter < dt_increase_niter and dt == self.dtime: dt *= 1.25 - # print('Increase the time increment to {:.5f}'.format(dt)) else: if update_dt: @@ -569,7 +657,14 @@ def assembly(self): @property def t_fact(self): - """Adimensional time used for boundary conditions.""" + """Adimensional time used for boundary conditions. + + Frozen during ``solve_time_increment`` so that modifications to + ``self.dtime`` (e.g. by a CCD step-size callback) do not alter + the target time factor mid-increment. + """ + if self._t_fact_inc is not None: + return self._t_fact_inc return (self.time + self.dtime - self.t0) / (self.tmax - self.t0) @property diff --git a/pyproject.toml b/pyproject.toml index 5bbf10aa..3fcb1e4b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -30,7 +30,8 @@ dependencies = [ ] [project.optional-dependencies] -all = ['fedoo[solver,plot,simcoon,test]'] +all = ['fedoo[solver,plot,simcoon,test,ipc]'] +ipc = ['ipctk'] solver = [ 'pypardiso ; platform_machine!="arm64" and platform_machine!="aarch64"', 'scikit-umfpack ; platform_machine=="arm64" or platform_machine=="aarch64"' @@ -50,21 +51,16 @@ gui = [ 'pyvistaqt', 'pyside6', ] -dataio = [ - 'pandas' -] other = [ - 'microgen' + 'pandas', + 'microgen', ] [project.urls] -Documentation = 'https://fedoo.readthedocs.io/en/latest/' +Documentation = 'https://3mah.github.io/fedoo-docs/' "Bug Tracker" = 'https://github.com/3MAH/fedoo/issues' "Source Code" = 'https://github.com/3MAH/fedoo' -[tool.setuptools.dynamic] -version = {attr = 'fedoo._version.__version__'} - [tool.setuptools.packages.find] include = [ 'fedoo', diff --git a/tests/test_2DDynamicPlasticBending.py b/tests/test_2DDynamicPlasticBending.py index e6616f8e..5ed06481 100644 --- a/tests/test_2DDynamicPlasticBending.py +++ b/tests/test_2DDynamicPlasticBending.py @@ -90,7 +90,7 @@ def test_2DDynamicPlasticBending(): bc1 = pb.bc.add("Dirichlet", nodes_top1, "DispY", uimp) bc2 = pb.bc.add("Dirichlet", nodes_top2, "DispY", uimp) - pb.nlsolve(dt=0.2, tmax=1, update_dt=False, tol_nr=0.005) + pb.nlsolve(dt=0.2, tmax=1, update_dt=False, tol_nr=0.01) ################### step 2 ################################ # bc.Remove() @@ -108,3 +108,9 @@ def test_2DDynamicPlasticBending(): # assert np.abs(res.node_data["Stress"][3][234] + 70.30052080276131) < 1e-4 # REMOVE ASSERT until simcoon bug is resolved + + +if __name__ == "__main__": + import pytest + + pytest.main([__file__]) diff --git a/tests/test_2DDynamicPlasticBending_v2.py b/tests/test_2DDynamicPlasticBending_v2.py index 5fae48cd..1fcfb3a8 100644 --- a/tests/test_2DDynamicPlasticBending_v2.py +++ b/tests/test_2DDynamicPlasticBending_v2.py @@ -90,7 +90,7 @@ def test_2DDynamicPlasticBending_v2(): bc1 = pb.bc.add("Dirichlet", nodes_top1, "DispY", uimp) bc2 = pb.bc.add("Dirichlet", nodes_top2, "DispY", uimp) - pb.nlsolve(dt=0.2, tmax=1, print_info=1, update_dt=False, tol_nr=0.005) + pb.nlsolve(dt=0.2, tmax=1, print_info=1, update_dt=False, tol_nr=0.01) ################### step 2 ################################ # bc.Remove() @@ -115,3 +115,9 @@ def test_2DDynamicPlasticBending_v2(): # assert np.abs(res.node_data["Stress"][3][234] + 67.82318305757613) < 1e-4 # REMOVE ASSERT until simcoon bug is resolved + + +if __name__ == "__main__": + import pytest + + pytest.main([__file__])