conspire/domain/partition/feti/block/solve/
mod.rs

1#[cfg(test)]
2mod test;
3
4use super::super::{
5    Formulation, THREADS,
6    dual::{coupling, coupling_transpose, rhs_from_forces},
7    dual_primal::{
8        BoundaryConditions, CornerSelection, build_splits,
9        coarse::{Coarse, CoarseSystem},
10        condense::Condensed,
11        rigid::{kernel, kernel_pins, removed_modes},
12        rigid_projector::{RigidProjector, add_rigid_motion, rigid_rhs},
13    },
14    interface::build_interfaces,
15    parallel::parallel_map,
16    pcg::{Preconditioner, primal_recovery, projected_pcg_with, rigid_projected_pcg},
17    subdomain::{DirichletLocal, Subdomain},
18};
19use crate::{
20    domain::ElementModelError,
21    geometry::mesh::Partition,
22    math::{
23        Scalar, SquareMatrix, Style, StyledError, Vector,
24        optimize::{KrylovError, KrylovMethod},
25        styled_error,
26    },
27};
28use std::collections::HashSet;
29
30/// Possible errors encountered when solving with FETI.
31pub enum SolveError {
32    /// Downstream error from an element.
33    Element(ElementModelError),
34    /// The partition does not describe the model it is meant to decompose.
35    Partition(String),
36    /// Downstream error from the dual PCG.
37    Krylov(KrylovError),
38    /// A subdomain left with some rigid-body modes free.
39    FloatingSubdomain { part: usize, removed: usize },
40    /// A subdomain with some part of it still free to move.
41    SingularSubdomain(usize),
42    /// The interior of a subdomain, away from the interface, is singular.
43    SingularInterior(usize),
44    /// The assembled corner problem is singular.
45    SingularCoarseProblem,
46}
47
48impl StyledError for SolveError {
49    fn message(&self, style: &Style) -> String {
50        match self {
51            Self::Element(error) => error.message(style),
52            Self::Partition(reason) => {
53                let (h, c) = (style.headline, style.frame);
54                format!("{h}The partition does not fit the model.{c}\n{reason}")
55            }
56            Self::Krylov(error) => error.message(style),
57            Self::FloatingSubdomain { part, removed } => {
58                let (h, c) = (style.headline, style.frame);
59                format!(
60                    "{h}Subdomain {part} is left floating.{c}\n\
61                    Its corner nodes and pinned degrees of freedom remove only {removed} of its \
62                    6 rigid-body modes, so its stiffness is singular. A subdomain needs three \
63                    non-collinear corner nodes (nodes shared by three or more subdomains) or \
64                    enough pinned degrees of freedom, so the partition has to bring more \
65                    subdomains together at a node."
66                )
67            }
68            Self::SingularSubdomain(part) => {
69                let (h, c) = (style.headline, style.frame);
70                format!(
71                    "{h}The stiffness of subdomain {part} is singular.{c}\n\
72                    Its corners and pinned degrees of freedom hold the subdomain as a whole, but \
73                    some part of it is still free to move: a piece attached to the rest through \
74                    only a node or an edge, or cut off from it altogether. Change the partition \
75                    so that every part of a subdomain is attached through faces. Otherwise the \
76                    model itself has a mechanism or a collapsed element."
77                )
78            }
79            Self::SingularInterior(part) => {
80                let (h, c) = (style.headline, style.frame);
81                format!(
82                    "{h}The interior of subdomain {part} is singular.{c}\n\
83                    The degrees of freedom of the subdomain away from the interface form a \
84                    singular block, though the subdomain as a whole does not, which the \
85                    Dirichlet preconditioner cannot handle. The tangent is likely not \
86                    positive definite there, as under severe compression, or the partition \
87                    leaves part of the interior loosely attached."
88                )
89            }
90            Self::SingularCoarseProblem => {
91                let (h, c) = (style.headline, style.frame);
92                format!(
93                    "{h}The problem coupling the corners is singular.{c}\n\
94                    Pin enough degrees of freedom to remove every rigid-body mode of the \
95                    whole block."
96                )
97            }
98        }
99    }
100}
101
102styled_error!(SolveError);
103
104impl From<ElementModelError> for SolveError {
105    fn from(error: ElementModelError) -> Self {
106        Self::Element(error)
107    }
108}
109
110impl From<KrylovError> for SolveError {
111    fn from(error: KrylovError) -> Self {
112        Self::Krylov(error)
113    }
114}
115
116const RELATIVE_PIVOT: Scalar = 1e-10;
117
118enum LocalError {
119    Singular,
120    Interior,
121}
122
123#[allow(clippy::too_many_arguments)]
124pub(crate) fn solve_local_systems<const D: usize>(
125    partition: &Partition,
126    boundary_conditions: &BoundaryConditions,
127    local_stiffnesses: Vec<SquareMatrix>,
128    local_forces: Vec<Vector>,
129    positions: &[[f64; D]],
130    preconditioner: Preconditioner,
131    rel_tol: Scalar,
132    method: KrylovMethod,
133    formulation: Formulation,
134) -> Result<Vector, SolveError> {
135    let corners = match formulation {
136        Formulation::Classical => CornerSelection::new(Vec::new()),
137        Formulation::DualPrimal => CornerSelection::from_partition(partition),
138    };
139    let (interfaces, num_multipliers) = build_interfaces(partition, &corners, D);
140    let (splits, corner_dofs) = build_splits(partition, &corners, boundary_conditions, D);
141    let subdomain_nodes = partition.parts_nodes();
142    let removable = D + D * (D - 1) / 2;
143    subdomain_nodes
144        .iter()
145        .zip(&splits)
146        .enumerate()
147        .try_for_each(|(part, (nodes, split))| {
148            if nodes.is_empty() || formulation == Formulation::Classical {
149                return Ok(());
150            }
151            let free: HashSet<usize> = split.dual().iter().copied().collect();
152            let constrained: Vec<usize> = (0..D * nodes.len())
153                .filter(|dof| !free.contains(dof))
154                .collect();
155            let local: Vec<[f64; D]> = nodes.iter().map(|&node| positions[node]).collect();
156            let removed = removed_modes(&local, &constrained);
157            if removed < removable {
158                Err(SolveError::FloatingSubdomain { part, removed })
159            } else {
160                Ok(())
161            }
162        })?;
163    let indices: Vec<usize> = (0..subdomain_nodes.len()).collect();
164    let condensed = parallel_map(&indices, THREADS, |&index| {
165        if formulation == Formulation::Classical {
166            return Some(Condensed::without_corners(splits[index].dual().len()));
167        }
168        Condensed::try_condense(
169            &local_stiffnesses[index],
170            &local_forces[index],
171            splits[index].primal(),
172            splits[index].dual(),
173        )
174    })
175    .into_iter()
176    .enumerate()
177    .map(|(part, condensed)| condensed.ok_or(SolveError::SingularSubdomain(part)))
178    .collect::<Result<Vec<_>, _>>()?;
179    let (schur, reduced_force) = CoarseSystem::assemble(&condensed, &splits, &corner_dofs);
180    let coarse_problem = Coarse::try_from(schur).map_err(|_| SolveError::SingularCoarseProblem)?;
181    let locals = parallel_map(&indices, THREADS, |&index| {
182        let stiffness = &local_stiffnesses[index];
183        let dual_dofs = splits[index].dual().to_vec();
184        let dual_stiffness: SquareMatrix = dual_dofs
185            .iter()
186            .map(|&row| dual_dofs.iter().map(|&col| stiffness[row][col]).collect())
187            .collect();
188        let (dual_factor, floating) = if formulation == Formulation::Classical {
189            let free: HashSet<usize> = dual_dofs.iter().copied().collect();
190            let constrained: Vec<usize> = (0..D * subdomain_nodes[index].len())
191                .filter(|dof| !free.contains(dof))
192                .collect();
193            let local: Vec<[f64; D]> = subdomain_nodes[index]
194                .iter()
195                .map(|&node| positions[node])
196                .collect();
197            let kernel = kernel(&local, &constrained);
198            let pins = kernel_pins(&kernel, &dual_dofs);
199            let keep: Vec<usize> = (0..dual_dofs.len())
200                .filter(|position| pins.binary_search(position).is_err())
201                .collect();
202            let reduced: SquareMatrix = keep
203                .iter()
204                .map(|&row| keep.iter().map(|&col| dual_stiffness[row][col]).collect())
205                .collect();
206            let factor = reduced
207                .factorize_lu()
208                .ok()
209                .filter(|factor| factor.near_zero_pivots(RELATIVE_PIVOT) == 0)
210                .ok_or(LocalError::Singular)?;
211            if kernel.is_empty() {
212                (factor, None)
213            } else {
214                (factor, Some((kernel, keep)))
215            }
216        } else {
217            let factor = dual_stiffness
218                .factorize_lu()
219                .expect("K_dd is singular, but corners should make every subdomain non-singular");
220            (factor, None)
221        };
222        DirichletLocal::try_build(stiffness, &dual_dofs, interfaces[index].dofs())
223            .map(|dirichlet| (dual_dofs, dual_stiffness, dual_factor, dirichlet, floating))
224            .ok_or(LocalError::Interior)
225    })
226    .into_iter()
227    .enumerate()
228    .map(|(part, local)| {
229        local.map_err(|error| match error {
230            LocalError::Singular => SolveError::SingularSubdomain(part),
231            LocalError::Interior => SolveError::SingularInterior(part),
232        })
233    })
234    .collect::<Result<Vec<_>, _>>()?;
235    let subdomains: Vec<Subdomain<()>> = interfaces
236        .into_iter()
237        .zip(locals)
238        .zip(splits.iter())
239        .zip(condensed.iter())
240        .zip(subdomain_nodes.iter())
241        .map(|((((interface, local), split), condensed), nodes)| {
242            let (dual_dofs, dual_stiffness, dual_factor, dirichlet, floating) = local;
243            let subdomain = Subdomain::new(
244                (),
245                interface,
246                dual_stiffness,
247                dual_factor,
248                dual_dofs,
249                nodes.len() * D,
250                condensed.dual_map.clone(),
251                condensed.primal_map.clone(),
252                split.primal().to_vec(),
253                split.primal_global().to_vec(),
254                dirichlet,
255            );
256            match floating {
257                Some((kernel, keep)) => subdomain.with_kernel(kernel, keep),
258                None => subdomain,
259            }
260        })
261        .collect();
262    let rhs = rhs_from_forces(&subdomains, &local_forces, num_multipliers)
263        - coupling(
264            &subdomains,
265            &coarse_problem.solve(&reduced_force),
266            num_multipliers,
267        );
268    let (lambda, alpha) = match formulation {
269        Formulation::DualPrimal => (
270            projected_pcg_with(
271                &subdomains,
272                &coarse_problem,
273                &rhs,
274                preconditioner,
275                rel_tol,
276                method,
277            )?,
278            None,
279        ),
280        Formulation::Classical => {
281            let projector =
282                RigidProjector::try_new(&subdomains).ok_or(SolveError::SingularCoarseProblem)?;
283            let (lambda, alpha) = rigid_projected_pcg(
284                &subdomains,
285                &projector,
286                &rhs,
287                &rigid_rhs(&subdomains, &local_forces),
288                preconditioner,
289                rel_tol,
290                method,
291            )?;
292            (lambda, Some(alpha))
293        }
294    };
295    let ct_lambda = coupling_transpose(&subdomains, &lambda, coarse_problem.len());
296    let corner_solution = coarse_problem.solve(&(reduced_force + ct_lambda));
297    let mut recovered = primal_recovery(&subdomains, &local_forces, &corner_solution, &lambda);
298    if let Some(alpha) = alpha {
299        add_rigid_motion(&subdomains, &alpha, &mut recovered);
300    }
301    let mut global = Vector::zero(positions.len() * D);
302    subdomain_nodes
303        .iter()
304        .zip(recovered.iter())
305        .for_each(|(nodes, local_solution)| {
306            nodes.iter().enumerate().for_each(|(local, &node)| {
307                (0..D).for_each(|component| {
308                    global[D * node + component] = local_solution[D * local + component]
309                })
310            })
311        });
312    Ok(global)
313}
conspire/domain/partition/feti/block/solve/mod.rs

conspire/domain/partition/feti/block/solve/
mod.rs
