AI @ FElupe

[adtzlr]

This Thread is a collection of ideas on how to use AI with FElupe.

a) Segment-to-Segment contact formulation (claude 3.7)

import numpy as np
from scipy.sparse import csr_matrix, lil_matrix
import felupe as fem


class SegmentToSegmentContact:
    """A class to handle segment-to-segment contact in FEM simulations.

    This class implements a penalty-based contact formulation between two
    segment surfaces. It detects contact pairs, computes contact forces,
    and assembles the corresponding residual vector and tangent matrix.

    Parameters
    ----------
    slave_field : fem.Field
        The field containing the slave surface segments
    master_field : fem.Field
        The field containing the master surface segments
    penalty : float
        Penalty parameter for contact enforcement
    friction_coef : float, optional
        Coefficient of friction, default is 0.0 (frictionless)
    """

    def __init__(self, slave_field, master_field, penalty=1e5, friction_coef=0.0):
        self.field = slave_field
        self.slave_field = slave_field
        self.master_field = master_field
        self.penalty = penalty
        self.friction_coef = friction_coef

        # Create assemble object to handle vector and matrix assembly
        self.assemble = fem.mechanics.Assemble(vector=self._vector, matrix=self._matrix)

        # Store results
        self.results = fem.mechanics.Results()

        # Contact pairs and status
        self.contact_pairs = []
        self.active_contacts = []
        self.gap_values = []
        self.contact_forces = []

        # Initialize contact search data
        self._initialize_contact_search()

    def _initialize_contact_search(self):
        """Initialize data structures for contact search."""
        # Get node coordinates for both surfaces
        self.slave_nodes = self.slave_field.region.mesh.points
        self.master_nodes = self.master_field.region.mesh.points

        # Get segment connectivity
        self.slave_segments = self.slave_field.region.mesh.cells
        self.master_segments = self.master_field.region.mesh.cells

    def update_contact_pairs(self):
        """Detect potential contact pairs between slave and master segments."""
        self.contact_pairs = []
        self.active_contacts = []
        self.gap_values = []

        # Get current nodal positions (including displacements)
        slave_positions = self.slave_nodes + self.slave_field[0].values
        master_positions = self.master_nodes + self.master_field[0].values

        # Simple all-to-all contact search (can be optimized with spatial partitioning)
        for i, slave_segment in enumerate(self.slave_segments):
            # Get slave segment nodes
            slave_points = slave_positions[slave_segment]

            for j, master_segment in enumerate(self.master_segments):
                # Get master segment nodes
                master_points = master_positions[master_segment]

                # Check if these segments are in contact
                is_contact, gap, contact_point, local_coords = (
                    self._check_segment_contact(slave_points, master_points)
                )

                if is_contact:
                    self.contact_pairs.append((i, j))
                    self.active_contacts.append(True)
                    self.gap_values.append(gap)

    def _check_segment_contact(self, slave_points, master_points):
        """Check if a slave segment is in contact with a master segment.

        Returns
        -------
        is_contact : bool
            True if contact detected
        gap : float
            Normal gap between segments (negative indicates penetration)
        contact_point : ndarray
            Coordinates of contact point
        local_coords : ndarray
            Local coordinates on segments
        """
        # For line segments in 2D or 3D
        # This is a simplified projection-based approach

        # Project slave segment midpoint to master segment
        slave_midpoint = np.mean(slave_points, axis=0)

        # Master segment as a line
        master_vec = master_points[1] - master_points[0]
        master_length = np.linalg.norm(master_vec)
        master_dir = master_vec / master_length

        # Vector from master segment start to slave midpoint
        to_slave = slave_midpoint - master_points[0]

        # Projection onto master segment direction (normalized length along segment)
        xi = np.dot(to_slave, master_dir) / master_length

        # Clamp to segment
        xi_clamped = max(0.0, min(1.0, xi))

        # Find closest point on master segment
        closest_point = master_points[0] + xi_clamped * master_vec

        # Compute gap vector and distance
        gap_vector = slave_midpoint - closest_point
        gap_distance = np.linalg.norm(gap_vector)

        # Normal direction to master segment
        if len(master_points[0]) == 2:  # 2D
            normal = np.array([-master_dir[1], master_dir[0]])
        else:  # 3D (requires additional definition, simplified here)
            # In 3D, we'd need a better definition of normal
            # This is highly simplified
            ref_vec = np.array([0, 0, 1])
            if np.allclose(np.cross(master_dir, ref_vec), 0):
                ref_vec = np.array([0, 1, 0])
            normal_dir = np.cross(master_dir, ref_vec)
            normal = normal_dir / np.linalg.norm(normal_dir)

        # Compute signed gap (negative means penetration)
        signed_gap = np.dot(gap_vector, normal)

        # Check if in contact (penetration)
        is_contact = signed_gap < 0

        return is_contact, signed_gap, closest_point, xi_clamped

    def _vector(self, field=None, **kwargs):
        """Compute the contact residual vector.

        Returns
        -------
        residual : csr_matrix
            Sparse contact residual vector
        """
        self.slave_field.link(field)
        self.master_field.link(field)

        # Update contact pairs
        self.update_contact_pairs()

        if not self.active_contacts:
            # No active contacts, return zero vector
            total_dofs = (
                self.slave_field.region.mesh.npoints * self.slave_field.region.mesh.dim
            )
            return csr_matrix(([0.0], ([0], [0])), shape=(total_dofs, 1))

        # Initialize sparse residual vector in LIL format for efficient assembly
        total_dofs = (
            self.slave_field.region.mesh.npoints * self.slave_field.region.mesh.dim
        )
        residual = lil_matrix((total_dofs, 1))

        # Get current nodal positions
        slave_positions = self.slave_nodes + self.slave_field[0].values
        master_positions = self.master_nodes + self.master_field[0].values

        # Process each active contact
        self.contact_forces = []
        for idx, (slave_idx, master_idx) in enumerate(self.contact_pairs):
            if not self.active_contacts[idx]:
                continue

            # Get segment nodes
            slave_segment = self.slave_segments[slave_idx]
            master_segment = self.master_segments[master_idx]

            slave_points = slave_positions[slave_segment]
            master_points = master_positions[master_segment]

            # Recompute contact data
            _, gap, contact_point, _ = self._check_segment_contact(
                slave_points, master_points
            )

            # Compute contact force (penalty method)
            if gap < 0:  # Negative gap = penetration
                # Master segment direction
                master_vec = master_points[1] - master_points[0]
                master_length = np.linalg.norm(master_vec)
                master_dir = master_vec / master_length

                # Normal direction
                if len(master_points[0]) == 2:  # 2D
                    normal = np.array([-master_dir[1], master_dir[0]])
                else:  # 3D
                    ref_vec = np.array([0, 0, 1])
                    if np.allclose(np.cross(master_dir, ref_vec), 0):
                        ref_vec = np.array([0, 1, 0])
                    normal_dir = np.cross(master_dir, ref_vec)
                    normal = normal_dir / np.linalg.norm(normal_dir)

                # Contact force magnitude (penalty * penetration depth)
                force_magnitude = -self.penalty * gap

                # Contact force vector
                contact_force = force_magnitude * normal
                self.contact_forces.append(contact_force)

                # Apply force to slave nodes (distribute equally for simplicity)
                for i, node_idx in enumerate(slave_segment):
                    for d in range(len(contact_force)):
                        dof = node_idx * self.slave_field.region.mesh.dim + d
                        residual[dof, 0] += contact_force[d] / len(slave_segment)

                # Apply opposite force to master nodes
                for i, node_idx in enumerate(master_segment):
                    for d in range(len(contact_force)):
                        dof = node_idx * self.master_field.region.mesh.dim + d
                        residual[dof, 0] -= contact_force[d] / len(master_segment)

        # Return as CSR matrix for efficient operations
        return residual.tocsr()

    def _matrix(self, field=None, **kwargs):
        """Compute the contact tangent matrix.

        Returns
        -------
        tangent : csr_matrix
            Sparse contact tangent matrix
        """
        # For simplicity, use finite difference approximation
        # In a production code, you would derive and implement the analytical tangent

        self.slave_field.link(field)
        self.master_field.link(field)

        total_dofs = (
            self.slave_field.region.mesh.npoints * self.slave_field.region.mesh.dim
        )
        tangent = lil_matrix((total_dofs, total_dofs))

        if not self.active_contacts:
            return tangent.tocsr()

        # Very simple tangent approximation for active contacts
        # This is a highly simplified version
        for idx, (slave_idx, master_idx) in enumerate(self.contact_pairs):
            if not self.active_contacts[idx]:
                continue

            slave_segment = self.slave_segments[slave_idx]
            master_segment = self.master_segments[master_idx]

            # Add penalty contribution to tangent matrix
            for i, slave_node in enumerate(slave_segment):
                for j, master_node in enumerate(master_segment):
                    for d in range(self.slave_field.region.mesh.dim):
                        slave_dof = slave_node * self.slave_field.region.mesh.dim + d
                        master_dof = master_node * self.master_field.region.mesh.dim + d

                        # Simplified tangent stiffness
                        k_contact = self.penalty / (
                            len(slave_segment) * len(master_segment)
                        )

                        # Add to tangent matrix
                        tangent[slave_dof, slave_dof] += k_contact
                        tangent[master_dof, master_dof] += k_contact
                        tangent[slave_dof, master_dof] -= k_contact
                        tangent[master_dof, slave_dof] -= k_contact

        return tangent.tocsr()

    def compute_contact_status(self):
        """Compute and return contact status information.

        Returns
        -------
        dict
            Dictionary containing contact information
        """
        status = {
            "num_contacts": len(self.active_contacts),
            "active_contacts": sum(self.active_contacts),
            "contact_pairs": self.contact_pairs,
            "gap_values": self.gap_values,
            "contact_forces": self.contact_forces,
        }

        return status


mesh = fem.Cube(n=4)
region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])

boundaries, loadcase = fem.dof.uniaxial(field, clamped=True)
umat = fem.NeoHooke(mu=1, bulk=2)
solid = fem.SolidBody(umat=umat, field=field)
contact = SegmentToSegmentContact(field, field)

move = fem.math.linsteps([0, 1], num=5)
ramp = {boundaries["move"]: move}
step = fem.Step(items=[solid, contact], ramp=ramp, boundaries=boundaries)

job = fem.Job(steps=[step])
job.evaluate()