conspire/domain/partition/feti/
mod.rs

1#![allow(dead_code)]
2
3#[cfg(test)]
4mod test;
5
6pub(crate) mod block;
7pub(crate) mod dual;
8pub(crate) mod dual_primal;
9pub(crate) mod interface;
10pub(crate) mod parallel;
11pub(crate) mod pcg;
12pub(crate) mod solid;
13pub(crate) mod subdomain;
14#[cfg(feature = "fem")]
15pub(crate) mod thermal;
16
17pub use block::element::{DecomposableElements, ElementSystems};
18#[cfg(feature = "fem")]
19pub use block::solve::SolveError;
20#[cfg(feature = "fem")]
21pub use dual_primal::BoundaryConditions;
22#[cfg(feature = "fem")]
23pub use pcg::Preconditioner;
24
25#[cfg(feature = "fem")]
26use crate::domain::NodalCoordinates;
27#[cfg(feature = "fem")]
28use crate::{
29    geometry::mesh::Partition,
30    math::{
31        Scalar, Vector,
32        optimize::{Krylov, KrylovMethod, LinearSolver},
33    },
34};
35#[cfg(feature = "fem")]
36use block::solve::solve_local_systems;
37
38pub(crate) const THREADS: usize = 1;
39
40/// The dual solve by GMRES, for a nonsymmetric dual operator, restarting
41/// only after a generous 100 iterations since the solve converges in tens.
42///
43/// Not the default: conjugate gradients also refuses a singular subdomain
44/// through its positive-definiteness check, which GMRES has no way to do.
45#[cfg(feature = "fem")]
46pub const GMRES: KrylovMethod = KrylovMethod::Gmres(100);
47
48/// How the subdomains are tied together.
49///
50/// Classical FETI ties every interface DOF with a Lagrange multiplier, so
51/// there are no corners and every subdomain is floating unless boundary
52/// conditions pin it. A floating subdomain's stiffness is singular along its
53/// rigid-body modes, so its local solve is a generalized inverse and the dual
54/// solve is projected against those modes.
55///
56/// A subdomain's rigid-body modes are an exact kernel of its tangent only
57/// where it carries no stress, which a subdomain cut out of a stressed body
58/// does at its interface. Classical FETI takes them as the kernel anyway, so
59/// unlike FETI-DP it is inexact for a geometrically nonlinear tangent, by an
60/// error that grows with the strain: 2e-3 of the solution at the strains of
61/// the tests, and 1e-2 of the strain at most. Inside Newton's method the
62/// solution is still the right one, but each step is approximate. The tangent
63/// must also be symmetric.
64///
65/// Its Dirichlet preconditioner is scaled by multiplicity, which it needs:
66/// unscaled it takes many times more iterations, growing with the number of
67/// subdomains. FETI-DP takes about the same either way, so it is not scaled.
68#[cfg(feature = "fem")]
69#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
70pub enum Formulation {
71    /// No corner nodes.
72    Classical,
73    /// Corner nodes are primal.
74    #[default]
75    DualPrimal,
76}
77
78/// FETI-DP solver for the linearized systems of a decomposable block.
79///
80/// The block is split by a [`Partition`], and each subdomain is solved
81/// independently, tied together through the corner DOFs and Lagrange
82/// multipliers on the interface. Any elastic block is supported. The tangent
83/// need not be symmetric, but then the dual solve must be [`GMRES`], since
84/// conjugate gradients, the default, needs a symmetric positive definite dual
85/// operator, and is what refuses a tangent that is not.
86///
87/// Only zero-displacement boundary conditions are supported. At least enough
88/// DOFs must be pinned to remove every rigid-body mode of the whole block.
89///
90/// As the linear solver of a [`NewtonRaphson`](crate::math::optimize::NewtonRaphson),
91/// it works with a fixed equality constraint only, whose fixed DOFs are the
92/// pinned ones.
93#[cfg(feature = "fem")]
94#[derive(Clone, Debug)]
95pub struct Feti {
96    pub formulation: Formulation,
97    pub method: KrylovMethod,
98    pub partition: Partition,
99    pub preconditioner: Preconditioner,
100    pub rel_tol: Scalar,
101}
102
103#[cfg(feature = "fem")]
104impl Default for Feti {
105    fn default() -> Self {
106        Self {
107            method: KrylovMethod::ConjugateGradients,
108            formulation: Formulation::DualPrimal,
109            partition: Partition::default(),
110            preconditioner: Preconditioner::Dirichlet,
111            rel_tol: Krylov::default().rel_tol,
112        }
113    }
114}
115
116#[cfg(feature = "fem")]
117impl Feti {
118    pub fn solve<B>(
119        &self,
120        block: &B,
121        nodal_coordinates: &NodalCoordinates<3>,
122        boundary_conditions: &BoundaryConditions,
123    ) -> Result<Vector, SolveError>
124    where
125        B: DecomposableElements,
126    {
127        let systems = block.element_systems(nodal_coordinates)?;
128        self.solve_systems(&systems, boundary_conditions)
129    }
130    fn solve_systems(
131        &self,
132        systems: &ElementSystems,
133        boundary_conditions: &BoundaryConditions,
134    ) -> Result<Vector, SolveError> {
135        let (stiffnesses, forces) = systems
136            .subdomains(&self.partition)
137            .map_err(SolveError::Partition)?;
138        solve_local_systems(
139            &self.partition,
140            boundary_conditions,
141            stiffnesses,
142            forces,
143            systems.positions(),
144            self.preconditioner,
145            self.rel_tol,
146            self.method,
147            self.formulation,
148        )
149    }
150}
151
152#[cfg(feature = "fem")]
153impl LinearSolver for Feti {
154    type Tangent = ElementSystems;
155    fn solve(
156        &self,
157        tangent: ElementSystems,
158        retained: &[usize],
159        _residual: &Vector,
160    ) -> Result<Vector, String> {
161        let mut fixed = vec![true; 3 * tangent.positions().len()];
162        retained.iter().for_each(|&dof| fixed[dof] = false);
163        let boundary_conditions = BoundaryConditions::new(
164            fixed
165                .iter()
166                .enumerate()
167                .filter(|&(_, &pinned)| pinned)
168                .map(|(dof, _)| (dof / 3, dof % 3))
169                .collect(),
170        );
171        let solution = self
172            .solve_systems(&tangent, &boundary_conditions)
173            .map_err(|error| error.to_string())?;
174        Ok(retained.iter().map(|&dof| solution[dof]).collect())
175    }
176}
conspire/domain/partition/feti/mod.rs

conspire/domain/partition/feti/
mod.rs
