IGA_for_bsplyne.Dirichlet

  1from typing import Union
  2import numpy as np
  3from numpy.typing import NDArray
  4import scipy.sparse as sps
  5from scipy.sparse import linalg as sps_linalg
  6import numba as nb
  7from numba.typed import List as nb_List
  8import copy
  9from IGA_for_bsplyne.solvers import solve_sparse, qr_sparse
 10
 11
 12class Dirichlet:
 13    """
 14    Affine representation of linear Dirichlet constraints.
 15
 16    This class represents Dirichlet boundary conditions through an affine
 17    transformation between a reduced vector of free parameters (`dof`)
 18    and the full physical displacement vector (`u`):
 19
 20        u = C @ dof + k
 21
 22    where:
 23
 24    - `u` is the full displacement vector satisfying all imposed constraints.
 25    - `dof` is the reduced vector of unconstrained degrees of freedom.
 26    - `C` is a sparse matrix whose columns form a basis of admissible variations.
 27    - `k` is a particular displacement vector satisfying the constraints.
 28
 29    This representation allows:
 30    - elimination of constrained DOFs,
 31    - support for general linear multi-point constraints.
 32
 33    Attributes
 34    ----------
 35    C : sps.csc_matrix of shape (n_full_dofs, n_free_dofs)
 36        Sparse matrix defining the linear mapping from reduced DOFs to full DOFs.
 37    k : NDArray of float of shape (n_full_dofs,)
 38        Particular solution ensuring that `u` satisfies the Dirichlet constraints.
 39    """
 40
 41    C: sps.csc_matrix
 42    k: NDArray[np.floating]
 43
 44    def __init__(self, C: sps.spmatrix, k: NDArray[np.floating]):
 45        """
 46        Initialize an affine Dirichlet constraint representation.
 47
 48        Parameters
 49        ----------
 50        C : sps.spmatrix
 51            Sparse matrix of shape (n_full_dofs, n_free_dofs) defining
 52            the linear mapping from reduced DOFs to full DOFs.
 53        k : NDArray[np.floating]
 54            Vector of shape (n_full_dofs,) defining the affine offset.
 55
 56        Notes
 57        -----
 58        The matrix `C` is internally converted to CSC format for efficient
 59        matrix-vector products.
 60        """
 61        self.C = C.tocsc()  # type: ignore
 62        self.k = k
 63
 64    @classmethod
 65    def eye(cls, size: int):
 66        """
 67        Create an unconstrained Dirichlet mapping.
 68
 69        This method returns a `Dirichlet` object corresponding to the identity
 70        transformation:
 71
 72            u = dof
 73
 74        i.e. no Dirichlet constraints are applied.
 75
 76        Parameters
 77        ----------
 78        size : int
 79            Number of degrees of freedom.
 80
 81        Returns
 82        -------
 83        dirichlet : Dirichlet
 84            Identity `Dirichlet` object with:
 85            - C = identity(size)
 86            - k = 0
 87        """
 88        C = sps.eye(size)
 89        k = np.zeros(size, dtype="float")
 90        return cls(C, k)  # type: ignore
 91
 92    @classmethod
 93    def lock_disp_inds(cls, inds: NDArray[np.integer], k: NDArray[np.floating]):
 94        """
 95        Create Dirichlet constraints by prescribing displacement values at selected DOFs.
 96
 97        This method enforces pointwise Dirichlet conditions of the form:
 98
 99            u[i] = k[i]    for i in inds
100
101        All other DOFs remain unconstrained.
102
103        Parameters
104        ----------
105        inds : NDArray[np.integer]
106            Indices of the displacement vector `u` to be constrained.
107        k : NDArray[np.floating]
108            Target displacement vector of shape (n_full_dofs,). Only the
109            values at `inds` are enforced.
110
111        Returns
112        -------
113        dirichlet : Dirichlet
114            `Dirichlet` instance representing the prescribed displacement constraints.
115
116        Notes
117        -----
118        Constrained DOFs are removed from the reduced DOF vector. The resulting
119        number of reduced DOFs will therefore be smaller than `k.size`.
120        """
121        C = sps.eye(k.size)
122        dirichlet = cls(C, k.copy())  # type: ignore
123        dirichlet.set_u_inds_vals(inds, dirichlet.k[inds])
124        return dirichlet
125
126    def set_u_inds_vals(self, inds: NDArray[np.integer], vals: NDArray[np.floating]):
127        """
128        Impose pointwise Dirichlet conditions by prescribing displacement values
129        at selected full DOFs.
130
131        This method enforces constraints of the form:
132
133            u[i] = vals[j]     for i = inds[j]
134
135        by modifying the affine mapping:
136
137            u = C @ dof + k
138
139        The modification removes the corresponding admissible variations in `C`
140        and updates `k` so that the prescribed values are always satisfied.
141        As a result, constrained DOFs are eliminated from the reduced DOF vector.
142
143        Parameters
144        ----------
145        inds : NDArray[np.integer]
146            Indices of the full displacement vector `u` to constrain.
147        vals : NDArray[np.floating]
148            Prescribed displacement values associated with `inds`.
149            Must have the same length as `inds`.
150
151        Notes
152        -----
153        - Rows of `C` corresponding to constrained DOFs are zeroed.
154        - Columns of `C` that become unused are removed, which may reduce
155          the number of reduced DOFs.
156        - The affine offset `k` is updated to enforce the prescribed values.
157        - This operation modifies the current Dirichlet object in place.
158        """
159
160        def zero_rows(C: sps.coo_matrix, rows: NDArray[np.integer]) -> sps.coo_matrix:
161            rm = np.hstack([(C.row == row).nonzero()[0] for row in rows])
162            row = np.delete(C.row, rm)
163            col = np.delete(C.col, rm)
164            data = np.delete(C.data, rm)
165            C = sps.coo_matrix((data, (row, col)), shape=(C.shape))
166            return C
167
168        def remove_zero_cols(C: sps.coo_matrix):
169            unique, inverse = np.unique(C.col, return_inverse=True)
170            C.col = np.arange(unique.size)[inverse]
171            C._shape = (C.shape[0], unique.size)
172
173        if inds.size != 0:
174            C_coo = zero_rows(self.C.tocoo(), inds)  # type: ignore
175            remove_zero_cols(C_coo)
176            self.C = C_coo.tocsc()  # type: ignore
177            self.k[inds] = vals
178
179    def slave_reference_linear_relation(
180        self,
181        slaves: NDArray[np.integer],
182        references: NDArray[np.integer],
183        coefs: Union[NDArray[np.floating], None] = None,
184    ):
185        """
186        Impose linear multi-point constraints between full DOFs.
187
188        This method enforces relations of the form:
189
190            u[slave_i] = sum_j coefs[i, j] * u[references[i, j]]
191
192        by modifying the affine mapping:
193
194            u = C @ dof + k
195
196        Slave DOFs are expressed as linear combinations of reference DOFs.
197        The corresponding admissible variations are propagated into `C`
198        and `k`, effectively eliminating slave DOFs from the reduced parameter
199        vector while preserving the imposed constraints.
200
201        Parameters
202        ----------
203        slaves : NDArray[np.integer] of shape (n_slaves,)
204            Indices of DOFs that are constrained (slave DOFs).
205        references : NDArray[np.integer] of shape (n_slaves, n_refs)
206            Reference DOF indices controlling each slave DOF.
207            Row `i` contains the references associated with `slaves[i]`.
208        coefs : NDArray[np.floating], optional, shape (n_slaves, n_refs)
209            Linear combination coefficients linking references to slaves.
210            If None, slaves are defined as the average of their references.
211
212        Notes
213        -----
214        - Slave DOFs are removed from the reduced DOF vector.
215        - The constraint propagation accounts for hierarchical dependencies
216          between slave DOFs using a topological ordering.
217        - This operation is typically faster than using `DirichletConstraintHandler`,
218          as it directly modifies the affine mapping without constructing
219          intermediate constraint objects.
220        - Cyclic dependencies between slaves are not supported and will raise an error.
221        - This operation modifies the current Dirichlet object in place.
222
223        Examples
224        --------
225        To impose u[i] = 0.5 * (u[j] + u[k]):
226
227        >>> slaves = np.array([i])
228        >>> references = np.array([[j, k]])
229        >>> coefs = np.array([[0.5, 0.5]])
230        >>> dirichlet.slave_reference_linear_relation(slaves, references, coefs)
231        """
232
233        if coefs is None:
234            coefs = np.ones(references.shape, dtype="float") / references.shape[1]
235
236        sorted_slaves = slave_reference_linear_relation_sort(slaves, references)
237
238        C = self.C.tocsr()
239        rows, cols, data, self.k = slave_reference_linear_relation_inner(
240            C.indices,
241            C.indptr,
242            C.data,
243            self.k,
244            slaves,
245            references,
246            coefs,
247            sorted_slaves,
248        )
249        C = sps.coo_matrix((data, (rows, cols)), shape=C.shape)
250
251        unique, inverse = np.unique(C.col, return_inverse=True)
252        C.col = np.arange(unique.size)[inverse]
253        C._shape = (C.shape[0], unique.size)
254        self.C = C.tocsc()  # type: ignore
255
256    def u_du_ddof(
257        self, dof: NDArray[np.floating]
258    ) -> tuple[NDArray[np.floating], sps.csc_matrix]:
259        """
260        Evaluate the displacement field and its derivative with respect to the reduced DOFs.
261
262        The displacement field is obtained from the affine mapping:
263
264            u = C @ dof + k
265
266        Since the mapping is affine, the derivative of `u` with respect to `dof`
267        is constant and equal to `C`.
268
269        Parameters
270        ----------
271        dof : NDArray[np.floating] of shape (n_dof,)
272            Reduced degrees of freedom.
273
274        Returns
275        -------
276        u : NDArray[np.floating] of shape (n_u,)
277            Displacement field.
278        du_ddof : sps.csc_matrix of shape (n_u, n_dof)
279            Jacobian of `u` with respect to `dof` : the matrix `C`.
280
281        Notes
282        -----
283        The returned derivative is independent of `dof` because the mapping is affine.
284        """
285        u = self.C @ dof + self.k
286        du_ddof = self.C
287        return u, du_ddof
288
289    def u(self, dof: NDArray[np.floating]) -> NDArray[np.floating]:
290        """
291        Evaluate the displacement field from reduced DOFs.
292
293        The displacement field is computed using the affine mapping:
294
295            u = C @ dof + k
296
297        Parameters
298        ----------
299        dof : NDArray[np.floating]
300            Reduced degrees of freedom. Can be either:
301
302            - shape (n_dof,)
303            - shape (..., n_dof) for batched evaluation
304
305        Returns
306        -------
307        u : NDArray[np.floating]
308            Displacement field with shape:
309
310            - (n_u,) if `dof` is 1D
311            - (..., n_u) if `dof` is batched
312
313        Notes
314        -----
315        This method supports vectorized evaluation over multiple DOF vectors.
316        """
317        if dof.ndim == 1:
318            u = self.C @ dof + self.k
319        else:
320            m, n = self.C.shape
321            dof_shape = dof.shape[:-1]
322            u = (self.C @ dof.reshape((-1, n)).T + self.k.reshape((m, 1))).T.reshape(
323                (*dof_shape, m)
324            )
325        return u
326
327    def dof_lsq(self, u: NDArray[np.floating]) -> NDArray[np.floating]:
328        """
329        Compute reduced DOFs from a displacement field using a least-squares projection.
330
331        This method computes `dof` such that:
332
333            u ≈ C @ dof + k
334
335        by solving the normal equations:
336
337            (Cᵀ C) dof = Cᵀ (u - k)
338
339        Parameters
340        ----------
341        u : NDArray[np.floating]
342            Displacement field. Can be either:
343
344            - shape (n_u,)
345            - shape (..., n_u) for batched projection
346
347        Returns
348        -------
349        dof : NDArray[np.floating]
350            Reduced degrees of freedom with shape:
351
352            - (n_dof,) if `u` is 1D
353            - (..., n_dof) if `u` is batched
354
355        Notes
356        -----
357        This operation performs a least-squares inversion of the affine mapping.
358        If `C` does not have full column rank, the solution corresponds to a
359        minimum-norm least-squares solution.
360
361        The system `(Cᵀ C)` is solved using a sparse linear solver.
362        """
363        if u.ndim == 1:
364            dof = solve_sparse(self.C.T @ self.C, self.C.T @ (u - self.k))
365        else:
366            m, n = self.C.shape
367            u_shape = u.shape[:-1]
368            dof = solve_sparse(
369                self.C.T @ self.C, (u.reshape((-1, m)) - self.k[None]) @ self.C
370            ).reshape((*u_shape, n))
371        return dof
372
373
374@nb.njit(cache=True)
375def slave_reference_linear_relation_sort(
376    slaves: NDArray[np.integer], references: NDArray[np.integer]
377) -> NDArray[np.integer]:
378    """
379    Sorts the slave nodes based on reference indices to respect
380    hierarchical dependencies (each slave is processed after its references).
381
382    Parameters
383    ----------
384    slaves : NDArray[np.integer]
385        Array of slave indices.
386    references : NDArray[np.integer]
387        2D array where each row contains the reference indices controlling a slave.
388
389    Returns
390    -------
391    sorted_slaves : NDArray[np.integer]
392        Array of slave indices sorted based on dependencies.
393    """
394    slaves_set = set(slaves)
395    graph = {s: nb_List.empty_list(nb.int64) for s in slaves}
396    indegree = {s: 0 for s in slaves}
397
398    # Building the graph
399    for i in range(len(slaves)):
400        for reference in references[i]:
401            if reference in slaves_set:
402                indegree[reference] += 1
403                graph[slaves[i]].append(reference)
404
405    # Queue for slaves with no dependencies (in-degree 0)
406    queue = [s for s in slaves if indegree[s] == 0]
407
408    # Topological sorting via BFS
409    sorted_slaves = np.empty(len(slaves), dtype="int")
410    i = 0
411    while queue:
412        slave = queue.pop(0)
413        sorted_slaves[i] = slave
414        i += 1
415        for dependent_slave in graph[slave]:
416            indegree[dependent_slave] -= 1
417            if indegree[dependent_slave] == 0:
418                queue.append(dependent_slave)
419
420    if i != len(slaves):
421        raise ValueError("Cyclic dependency detected in slaves.")
422
423    sorted_slaves = sorted_slaves[::-1]
424    return sorted_slaves
425
426
427@nb.njit(cache=True)
428def slave_reference_linear_relation_inner(
429    indices: NDArray[np.integer],
430    indptr: NDArray[np.integer],
431    data: NDArray[np.floating],
432    k: NDArray[np.floating],
433    slaves: NDArray[np.integer],
434    references: NDArray[np.integer],
435    coefs: NDArray[np.floating],
436    sorted_slaves: NDArray[np.integer],
437) -> tuple[
438    NDArray[np.integer], NDArray[np.integer], NDArray[np.floating], NDArray[np.floating]
439]:
440    """
441    Applies slave-reference relations directly to CSR matrix arrays.
442
443    Parameters
444    ----------
445    indices : NDArray[np.integer]
446        Column indices of CSR matrix.
447    indptr : NDArray[np.integer]
448        Row pointers of CSR matrix.
449    data : NDArray[np.floating]
450        Nonzero values of CSR matrix.
451    k : NDArray[np.floating]
452        Vector to be updated.
453    slaves : NDArray[np.integer]
454        Array of slave indices.
455    references : NDArray[np.integer]
456        2D array where each row contains the reference indices controlling a slave.
457    coefs : NDArray[np.floating]
458        2D array of coefficients defining the linear relationship between references and
459        slaves.
460    sorted_slaves : NDArray[np.integer]
461        Array of slave indices sorted in topological order.
462
463    Returns
464    -------
465    rows : NDArray[np.integer]
466        Updated row indices of COO matrix.
467    cols : NDArray[np.integer]
468        Updated column indices of COO matrix.
469    data : NDArray[np.floating]
470        Updated nonzero values of COO matrix.
471    k : NDArray[np.floating]
472        Updated vector.
473    """
474
475    # Convert CSR to list of dicts
476    dict_list = [
477        {indices[j]: data[j] for j in range(indptr[i], indptr[i + 1])}
478        for i in range(indptr.size - 1)
479    ]
480
481    # Apply linear relation
482    slaves_to_index = {s: i for i, s in enumerate(slaves)}
483    for slave in sorted_slaves:
484        i = slaves_to_index[slave]
485        new_row = {}
486        for reference_ind, coef in zip(references[i], coefs[i]):
487            reference_row = dict_list[reference_ind]
488            for ind, val in reference_row.items():
489                if ind in new_row:
490                    new_row[ind] += coef * val
491                else:
492                    new_row[ind] = coef * val
493        new_row = {ind: val for ind, val in new_row.items() if val != 0}
494        dict_list[slave] = new_row
495        k[slave] = np.dot(k[references[i]], coefs[i])
496
497    # Convert list of dicts to COO
498    nnz = sum([len(row) for row in dict_list])
499    rows = np.empty(nnz, dtype=np.int32)
500    cols = np.empty(nnz, dtype=np.int32)
501    data = np.empty(nnz, dtype=np.float64)
502    pos = 0
503    for i, row in enumerate(dict_list):
504        for j, val in row.items():
505            rows[pos] = i
506            cols[pos] = j
507            data[pos] = val
508            pos += 1
509
510    return rows, cols, data, k
511
512
513# %%
514class DirichletConstraintHandler:
515    """
516    Accumulate affine linear constraints and construct the associated Dirichlet mapping.
517
518    This class is designed to impose linear affine constraints of the form:
519
520        D @ u = c
521
522    where `u` is the full displacement (or state) vector. Constraints are progressively
523    accumulated and, once fully specified, converted into an affine parametrization of
524    the admissible solution space:
525
526        u = C @ dof + k
527
528    where:
529
530    - `C` is a basis of the nullspace of `D` (i.e. `D @ C = 0`) obtained by
531    QR decomposition,
532    - `k` is a particular solution satisfying the constraints,
533    obtained through the QR decomposition which helps solve:
534
535        D k = c
536
537    - `dof` is a reduced vector of unconstrained degrees of freedom.
538
539    The main purpose of this class is to provide a flexible and robust interface
540    to define constraints before constructing a `Dirichlet` object representing
541    the reduced parametrization.
542
543    Typical workflow
544    ----------------
545    1. Create a handler with the number of physical DOFs (`u`).
546    2. Add constraint equations using helper methods.
547    3. **Ensure that all reference DOFs (e.g., translations or rotations introduced for
548    rigid-body relations) are fully constrained before computing `C` and `k`.**
549    4. Build the reduced representation `(C, k)` or directly create a `Dirichlet` object.
550
551    This separation allows constraints to be assembled incrementally and validated
552    before generating the final affine mapping.
553
554    Attributes
555    ----------
556    nb_dofs_init : int
557        Number of DOFs in the original unconstrained system (physical DOFs `u`).
558    lhs : sps.spmatrix
559        Accumulated constraint matrix `D`.
560    rhs : NDArray[np.floating]
561        Accumulated right-hand side vector `c`.
562
563    Notes
564    -----
565    - The constraint system may include additional reference DOFs introduced
566    to express kinematic relations or rigid-body behaviors.
567    - **All reference DOFs must be fully constrained before computing `C` and `k`;**
568    otherwise, DOFs that lie in the kernel of `D` cannot be controlled and imposed values
569    (e.g., prescribed translations) may not appear in the resulting solution `k`.
570    - The resulting affine mapping guarantees that any generated vector `u`
571    satisfies all imposed constraints.
572    """
573
574    nb_dofs_init: int
575    lhs: sps.spmatrix
576    rhs: NDArray[np.floating]
577
578    def __init__(self, nb_dofs_init: int):
579        """
580        Initialize a Dirichlet constraint handler.
581
582        Parameters
583        ----------
584        nb_dofs_init : int
585            The number of initial degrees of freedom in the unconstrained system.
586            This value is used to size the initial constraint matrix and manage
587            later extensions with reference DOFs.
588        """
589        self.nb_dofs_init = nb_dofs_init
590        self.lhs = sps.coo_matrix(np.empty((0, nb_dofs_init), dtype="float"))
591        self.rhs = np.empty(0, dtype="float")
592
593    def copy(self) -> "DirichletConstraintHandler":
594        """
595        Create a deep copy of this DirichletConstraintHandler instance.
596
597        Returns
598        -------
599        DirichletConstraintHandler
600            A new instance with the same initial number of DOFs, constraint matrix, and right-hand side vector.
601            All internal data is copied, so modifications to the returned handler do not affect the original.
602        """
603        return copy.deepcopy(self)
604
605    def add_eqs(self, lhs: sps.spmatrix, rhs: NDArray[np.floating]):
606        """
607        Append linear affine constraint equations to the global constraint system.
608
609        Adds equations of the form:
610
611            lhs @ u = rhs
612
613        to the accumulated constraint system `D @ u = c`.
614
615        If `lhs` is expressed only in terms of the initial physical DOFs
616        (`nb_dofs_init` columns), it is automatically extended with zero-padding
617        to match the current number of DOFs (e.g. after reference DOFs have been added).
618
619        Parameters
620        ----------
621        lhs : sps.spmatrix of shape (n_eqs, n_dofs)
622            Constraint matrix defining the left-hand side of the equations.
623            The number of columns must match either:
624            - the initial number of DOFs (`nb_dofs_init`), or
625            - the current number of DOFs in the handler.
626
627        rhs : NDArray[np.floating] of shape (n_eqs,)
628            Right-hand side values associated with the constraint equations.
629
630        Raises
631        ------
632        ValueError
633            If the number of columns of `lhs` is incompatible with the current
634            constraint system.
635
636        Notes
637        -----
638        Constraints are appended to the existing system and are not simplified
639        or checked for redundancy.
640        """
641        nb_eqs, nb_dofs = lhs.shape
642        assert nb_eqs == rhs.size
643        if nb_dofs == self.lhs.shape[1]:
644            self.lhs = sps.vstack((self.lhs, lhs))  # type: ignore
645        elif nb_dofs == self.nb_dofs_init:
646            zero = sps.coo_matrix((nb_eqs, self.lhs.shape[1] - self.nb_dofs_init))
647            new_lhs = sps.hstack((lhs, zero))
648            self.lhs = sps.vstack((self.lhs, new_lhs))  # type: ignore
649        else:
650            raise ValueError(
651                "lhs.shape[1] must match either the initial number of dofs or the current one."
652            )
653
654        self.rhs = np.hstack((self.rhs, rhs))
655
656    def add_ref_dofs(self, nb_dofs: int):
657        """
658        Append additional reference DOFs to the constraint system.
659
660        This method increases the size of the global DOF vector by adding
661        `nb_dofs` new reference DOFs. These DOFs are introduced without
662        imposing any constraint, i.e. they appear with zero coefficients
663        in all existing equations.
664
665        Parameters
666        ----------
667        nb_dofs : int
668            Number of reference DOFs to append to the global DOF vector.
669
670        Notes
671        -----
672        Reference DOFs are typically used to express kinematic relations,
673        rigid body motions, or other auxiliary constraint parametrizations.
674        """
675        zero = sps.coo_matrix((self.lhs.shape[0], nb_dofs))
676        self.lhs = sps.hstack((self.lhs, zero))  # type: ignore
677
678    def add_ref_dofs_with_behavior(
679        self, behavior_mat: sps.spmatrix, slave_inds: NDArray[np.integer]
680    ):
681        """
682        Introduce reference DOFs and constrain slave DOFs through a linear relation.
683
684        This method appends new reference DOFs and enforces a linear behavioral
685        relation linking these reference DOFs to existing DOFs (called *slave DOFs*).
686        The imposed constraints take the form:
687
688            behavior_mat @ ref_dofs - u[slave_inds] = 0
689
690
691        Parameters
692        ----------
693        behavior_mat : sps.spmatrix of shape (n_slaves, n_ref_dofs)
694            Linear operator defining how each reference DOF contributes to the
695            corresponding slave DOFs.
696
697        slave_inds : NDArray[np.integer] of shape (n_slaves,)
698            Indices of slave DOFs that are controlled by the reference DOFs.
699            The ordering must match the rows of `behavior_mat`.
700
701        Raises
702        ------
703        AssertionError
704            If the number of slave DOFs is inconsistent with the number of
705            rows of `behavior_mat`.
706
707        Notes
708        -----
709        This method first adds the reference DOFs to the global system and then
710        appends the corresponding constraint equations.
711        """
712        nb_slaves, nb_ref_dofs = behavior_mat.shape
713        assert nb_slaves == slave_inds.size
714        data = -np.ones(nb_slaves, dtype="float")
715        rows = np.arange(nb_slaves)
716        slaves_counterpart = sps.coo_matrix(
717            (data, (rows, slave_inds)), shape=(nb_slaves, self.lhs.shape[1])
718        )
719        lhs_to_add = sps.hstack((slaves_counterpart, behavior_mat))
720        rhs_to_add = np.zeros(nb_slaves, dtype="float")
721        self.add_ref_dofs(nb_ref_dofs)
722        self.add_eqs(lhs_to_add, rhs_to_add)  # type: ignore
723
724    def add_rigid_body_constraint(
725        self,
726        ref_point: NDArray[np.floating],
727        slaves_inds: NDArray[np.integer],
728        slaves_positions: NDArray[np.floating],
729    ):
730        """
731        Constrain slave nodes to follow a rigid body motion defined by a reference point.
732
733        This method introduces six reference DOFs representing a rigid body motion:
734        three rotations and three translations (in this order) around a reference point.
735        The displacement of each slave node is constrained to follow the rigid body motion:
736
737            u(X) = θ × (X - X_ref) + t
738
739        where `θ` is the rotation vector and `t` is the translation vector.
740
741        Parameters
742        ----------
743        ref_point : NDArray[np.floating] of shape (3,)
744            Reference point defining the center of rotation.
745
746        slaves_inds : NDArray[np.integer] of shape (3, n_nodes)
747            DOF indices of the slave nodes. Each column contains the
748            x, y, z DOF indices of a slave node.
749
750        slaves_positions : NDArray[np.floating] of shape (3, n_nodes)
751            Physical coordinates of the slave nodes.
752
753        Notes
754        -----
755        - Six reference DOFs (θx, θy, θz, tx, ty, tz) are added to represent
756          the rigid body motion.
757        - The constraint is expressed as a linear relation between the
758          reference DOFs and the slave displacements.
759        - This method is commonly used to impose rigid connections,
760          master node formulations, or reference frame constraints.
761        """
762        x, y, z = (
763            slaves_positions[:, :, None] - ref_point[:, None, None]
764        )  # x is of shape (n, 1)
765        ones = np.ones_like(x)
766        behavior = sps.bmat(
767            [
768                [None, z, -y, ones, None, None],
769                [-z, None, x, None, ones, None],
770                [y, -x, None, None, None, ones],
771            ]
772        )  # behavior matrix for U(X; theta, t) = theta \carat X + t
773        inds = slaves_inds.ravel()
774        self.add_ref_dofs_with_behavior(behavior, inds)  # type: ignore
775
776    def add_eqs_from_inds_vals(
777        self, inds: NDArray[np.integer], vals: Union[NDArray[np.floating], None] = None
778    ):
779        """
780        Impose pointwise Dirichlet conditions on selected DOFs.
781
782        This is a convenience method that adds constraint equations of the form:
783
784            u[inds] = vals
785
786
787        Parameters
788        ----------
789        inds : NDArray[np.integer] of shape (n_eqs,)
790            Indices of DOFs to constrain.
791
792        vals : NDArray[np.floating] of shape (n_eqs,), optional
793            Prescribed values associated with each constrained DOF.
794            If None, zero values are imposed.
795
796        Raises
797        ------
798        AssertionError
799            If `vals` is provided and its size differs from `inds`.
800
801        Notes
802        -----
803        This method is equivalent to adding rows of the identity matrix
804        to the constraint operator.
805        """
806        nb_eqs = inds.size
807        vals = np.zeros(nb_eqs, dtype="float") if vals is None else vals
808        assert vals.size == nb_eqs
809        data = np.ones(nb_eqs, dtype="float")
810        rows = np.arange(nb_eqs)
811        lhs = sps.coo_matrix((data, (rows, inds)), shape=(nb_eqs, self.lhs.shape[1]))
812        self.add_eqs(lhs, vals)
813
814    def make_C_k(self) -> tuple[sps.csc_matrix, NDArray[np.floating]]:
815        """
816        Construct the affine transformation (C, k) enforcing all Dirichlet constraints.
817
818        This method solves the linear constraint system:
819
820            D @ u = c
821
822        by computing:
823
824        - a matrix `C` forming a basis of the nullspace of `D`
825        (i.e., D @ C = 0),
826        - a particular solution `k` such that D @ k = c.
827
828        Any admissible vector `u` satisfying the constraints can then be written as:
829
830            u = C @ dof + k
831
832        where `dof` is the reduced vector of unconstrained degrees of freedom.
833
834        The nullspace basis and the particular solution are obtained through
835        a sparse QR factorization of Dᵀ.
836
837        Returns
838        -------
839        C : sps.csc_matrix
840            Sparse matrix of shape (n_full_dofs, n_free_dofs) whose columns
841            form a basis of the nullspace of `D`.
842
843        k : NDArray[np.floating]
844            Vector of shape (n_full_dofs,) representing a particular solution
845            satisfying the constraints.
846
847        Notes
848        -----
849        - The transformation (C, k) defines an affine parametrization of the
850        admissible displacement space.
851        - The reduced DOF vector `dof` corresponds to the coordinates of `u`
852        in the nullspace basis.
853        - The QR factorization is performed on Dᵀ to efficiently extract both
854        the nullspace basis and the particular solution.
855        """
856        # Step 1: Perform sparse QR factorization of D^T
857        # D^T = Q @ R @ P^T, where D = self.lhs and P is a column permutation
858        Q, R, perm, rank = qr_sparse(self.lhs.T)  # type: ignore
859        Q = Q.tocsc()  # type: ignore
860        R = R.tocsc()  # type: ignore
861
862        # Step 2: Extract a basis C of the nullspace of D (i.e., such that D @ C = 0)
863        # These correspond to the last n - rank columns of Q
864        C = Q[:, rank:]
865
866        # Step 3: Compute a particular solution k such that D @ k = c
867        # (c = self.rhs contains the prescribed values)
868
869        # 3.1: Apply the column permutation to the RHS vector
870        c_tilde = self.rhs[np.array(perm)]
871
872        # 3.2: Extract the leading square block R1 of R (size rank × rank)
873        R1 = R[:rank, :rank].tocsc()
874
875        # 3.3: Take the first 'rank' entries of the permuted RHS
876        c1 = c_tilde[:rank]
877
878        # 3.4: Solve the triangular system R1^T @ y1 = c1
879        y1 = sps_linalg.spsolve(R1.T, c1)
880
881        # 3.5: Complete the vector y by padding with zeros (for the nullspace part)
882        y = np.zeros(Q.shape[0])
883        y[:rank] = y1
884
885        # 3.6: Recover the particular solution k = Q @ y
886        k = Q @ y
887
888        # Step 4: Return the nullspace basis C and the particular solution k
889        return C, k  # type: ignore
890
891    def get_reduced_Ck(self) -> tuple[sps.spmatrix, NDArray[np.floating]]:
892        """
893        Extract the affine constraint transformation restricted to physical DOFs.
894
895        This method computes the full transformation (C, k) using `self.make_C_k()`
896        and returns only the rows associated with the initial (physical) degrees
897        of freedom. The resulting pair (C_u, k_u) defines:
898
899            u_phys = C_u @ dof + k_u
900
901        where `u_phys` satisfies all imposed Dirichlet constraints.
902
903        Returns
904        -------
905        C_u : sps.spmatrix
906            Reduced nullspace basis matrix of shape (nb_dofs_init, n_free_dofs).
907
908        k_u : NDArray[np.floating]
909            Reduced particular solution vector of shape (nb_dofs_init,).
910
911        Notes
912        -----
913        - The full transformation (C, k) may include auxiliary reference DOFs
914        introduced by multi-point or hierarchical constraints.
915        - Only the rows corresponding to physical DOFs are returned.
916        - A warning is emitted if reference DOFs remain unconstrained,
917        which may indicate missing boundary specifications.
918        """
919        C_ext, k_ext = self.make_C_k()
920        C_ref = C_ext[self.nb_dofs_init :, :]
921        k_ref = k_ext[self.nb_dofs_init :]
922        C_u = C_ext[: self.nb_dofs_init, :]
923        k_u = k_ext[: self.nb_dofs_init]
924        if sps_linalg.norm(C_ref) != 0:
925            print("Warning : not all reference point DoFs are specified.")
926            print(
927                f"Try specifying reference dofs number {np.abs(C_ref).sum(axis=1).A.ravel().nonzero()[0].tolist()}."  # type: ignore
928            )
929        return C_u, k_u
930
931    def create_dirichlet(self):
932        """
933        Build a `Dirichlet` object representing the reduced constraint transformation.
934
935        This is a convenience wrapper around `get_reduced_Ck()` that constructs
936        a `Dirichlet` object directly from the reduced affine mapping:
937
938            u_phys = C_u @ dof + k_u
939
940        Returns
941        -------
942        dirichlet : Dirichlet
943            A `Dirichlet` instance encapsulating the reduced transformation matrices.
944        """
945        C, k = self.get_reduced_Ck()
946        return Dirichlet(C, k)
947
948    def get_ref_multipliers_from_internal_residual(
949        self, K_u_minus_f: NDArray[np.floating]
950    ) -> NDArray[np.floating]:
951        """
952        Recover Lagrange multipliers associated with reference DOF constraints.
953
954        This method reconstructs the Lagrange multipliers λ corresponding to
955        reference DOFs using the internal residual vector of the mechanical system:
956
957            r = K @ u - f
958
959        The derivation relies on the partition of the transformation matrix:
960
961            C = [C_u;
962                C_ref]
963
964        where `C_u` maps reduced DOFs to physical DOFs and `C_ref` maps them
965        to reference DOFs.
966
967        The multipliers satisfy:
968
969            C_ref.T @ λ = - C_u.T @ r
970
971        and are obtained via the least-squares solution:
972
973            λ = - (C_ref @ C_ref.T)⁻¹ @ C_ref @ C_u.T @ r
974
975        Parameters
976        ----------
977        K_u_minus_f : NDArray[np.floating]
978            Internal residual vector of shape (nb_dofs_init,),
979            corresponding to K @ u - f restricted to physical DOFs.
980
981        Returns
982        -------
983        lamb : NDArray[np.floating]
984            Vector of Lagrange multipliers associated with reference DOF constraints.
985
986        Notes
987        -----
988        - The multipliers correspond to reaction forces or generalized constraint
989        forces transmitted through reference DOFs.
990        - The residual must be assembled consistently with the stiffness matrix
991        and load vector used to compute the displacement field.
992        - The matrix C_ref @ C_ref.T is assumed to be invertible, which holds
993        when reference constraints are linearly independent.
994        """
995        C, _ = self.make_C_k()
996        C_ref = C[self.nb_dofs_init :, :]
997        C_u = C[: self.nb_dofs_init, :]
998        lamb = -solve_sparse(C_ref @ C_ref.T, C_ref @ C_u.T @ K_u_minus_f)
999        return lamb