Add check for singular mass matrix to dynamic and eigenvalue solvers

docs/theory-manual/latex/solution.tex

@@ -344,23 +344,28 @@ \section{Dynamic Solution}
 
 Early prototypes of VirtualBow followed this approach and solved the resulting first order problem with methods provided by the C++ library Odeint \cite{bib:ahnert2011}.
 Especially the error controlled Runge-Kutta methods worked reasonably well.
-In addition to those general methods, there are also numerous methods that have been developed with finite elements/structural dynamics in mind.
-Those operate directly on the second order problem~(\ref{eq:dynamics:system_equation}) and are usually more effective for this type of problem.
+However, in addition to such general methods for first order systems, there are also numerous methods that have been developed with finite elements/structural dynamics in mind.
+They operate directly on the second order problem~(\ref{eq:dynamics:system_equation}) and are usually more effective for this type of problem.
 
-All numerical methods for solving ODEs can be grouped into either explicit and implicit ones.
+Apart from the order of the ODE, all numerical methods for solving ODEs can be grouped into either explicit and implicit ones.
 Explicit methods use the current state of the system at time~$t$ to calculate the next state at~$t + \Delta t$.
 Implicit methods instead find the next state by solving an equation that involves the current as well as the next state.
-In the context of finite element analysis, implicit methods require equilibrium iterations at every time step involving the system's tangent stiffness matrix, while explicit methods only require simple vector operations.
-Implicit methods therefore have a higher computational cost per step, but can make up for it by requiring less steps overall, as their better stability allows them to use much larger step sizes than explicit methods.
+In the context of finite element analysis, implicit methods require equilibrium iterations at every time step involving the system's tangent stiffness and damping matrices, while explicit methods only require simple vector operations to compute the next solution state.
+Implicit methods therefore have a higher computational cost per step, but can make up for it by requiring less steps overall, since their better stability properties allow them to use much larger step sizes than explicit methods.
 
 Explicit methods are easier to implement and are a good choice for small time scales, like impact and wave propagation problems.
 They are most effective when the computation of the time steps is fast, so they work well with constant, diagonal (lumped) mass matrices and simple elements with low-order shape functions.
-
 Implicit methods in contrast are well suited for larger time scales, where explicit algorithms would need an excessive number of steps.
 Examples are structural response problems or numerically stiff systems, i.e. systems whose eigenfrequencies range from very small to very large.
 
-From version 0.1 to 0.9, VirtualBow used the explicit central difference method \cite{bib:bathe2006}, which is probably one of the most popular explicit methods in finite element analysis.
-Later versions use the Generalized-$\alpha$ method instead, which is an implicit method and closely related to the also very popular Newmark-$\beta$ method.
+Based on these criteria, implicit methods should be the better fit for simulating the dynamics of a bow.
+However, from version 0.1 to 0.9, VirtualBow used the explicit central difference method \cite{bib:bathe2006}, which is one of the standard explicit methods in finite element analysis.
+The reason for this choice was the relatively simple implementation and a sufficiently good performance.
+
+Starting with version 0.10, VirtualBow uses the Newmark-$\beta$ method instead, which is a popular implicit method that is commonly used in structural dynamics.
+The initial plan was to use the Generalized-$\alpha$ method, which is closely related to the Newmark method and improves some of its properties.
+It was however decided to release VirtualBow with the simpler Newmark method first to gain some experience and feedback and only switch to the Generalized-$\alpha$ method later if a clear advantage for the bow simulation can be shown.
+That's why both of these methods are documented below even though only one is used.
 
 \subsection{Newmark-$\boldsymbol{\beta}$ Method}
 

solver/src/fem/solvers/dynamics.rs

@@ -25,6 +25,7 @@ impl Default for Settings {
 
 #[derive(PartialEq, Debug)]
 pub enum DynamicSolverError {
+    SingularMassMatrix,       // The mass matrix is singular, i.e. cannot be inverted
     LinearSolutionFailed,     // Decomposition of the tangent stiffness matrix or solution of the linear system failed
     NonFiniteStateIncrement,  // The displacement delta is not finite, i.e. contains NaN or Inf values
     MaxIterationsReached,     // Maximum number of iterations was reached without convergence
@@ -35,6 +36,7 @@ pub enum DynamicSolverError {
 impl Display for DynamicSolverError {
     fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
         match self {
+            DynamicSolverError::SingularMassMatrix      => write!(f, "The system has a singular mass matrix.")?,
             DynamicSolverError::LinearSolutionFailed    => write!(f, "Decomposition/solution of the tangent stiffness matrix failed.")?,
             DynamicSolverError::NonFiniteStateIncrement => write!(f, "Encountered a non-finite displacement increment.")?,
             DynamicSolverError::MaxIterationsReached    => write!(f, "Maximum number of iterations exceeded.")?,
@@ -76,8 +78,12 @@ impl<'a> DynamicSolver<'a> {
         // Solver state
         let mut t = self.system.get_time();
 
-        // First evaluation at start of the simulated interval to make acceleration available
-        // and provide callback information at t = t0.
+        // Check if the mass matrix is positive definite
+        if eval.get_mass_matrix().amax() < 0.0 {
+            return Err(DynamicSolverError::SingularMassMatrix);
+        }
+
+        // First evaluation at start of the simulated interval to make acceleration available and provide callback information at t = t0.
         self.system.eval_dynamics(&mut eval);
         if !callback(&self.system, &eval) {
             return Ok(());

solver/src/fem/solvers/eigen.rs

@@ -5,13 +5,15 @@ use crate::fem::system::system::System;
 
 #[derive(PartialEq, Debug)]
 pub enum EigenSolverError {
+    SingularMassMatrix,
     MatrixInversionFailed,
     NonFiniteEigenvalues,
 }
 
 impl Display for EigenSolverError {
     fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
         match self {
+            EigenSolverError::SingularMassMatrix    => write!(f, "The system has a singular mass matrix.")?,
             EigenSolverError::MatrixInversionFailed => write!(f, "Inversion of the system matrix failed.")?,
             EigenSolverError::NonFiniteEigenvalues  => write!(f, "At least one of the computed Eigenvalues is non-finite.")?,
         }
@@ -69,9 +71,13 @@ pub fn natural_frequencies_from_matrices(M: &DVector<f64>, D: &DMatrix<f64>, K:
     ];
     */
 
+    // Check if the mass matrix is positive definite
+    if M.amax() < 0.0 {
+        return Err(EigenSolverError::SingularMassMatrix);
+    }
+
     // Compute the complex eigenvalues for A*v = lambda*B*v by inverting B and solving B^(-1)*A*v = lambda*v.
     // Using a generalized eigenvalue solver would have been better, but nalgebra currently doesn't have one.
-
     let M_inv = DMatrix::from_diagonal(&M.map(|m| 1.0/m));
     let A = stack![
         0, DMatrix::identity(M.len(), M.len());