Applying boundary conditions to the stress stiffness matrix

When solving buckling constrained optimization, one common formulation is to enforce K + c * Ksigma to be positive definite for some load multiplier c > 1. Another formulation is to define a function that outputs the maximum eigenvalue of the system:

-Ksigma * v = lambda * K * v

and then use the constraint 1 / lambda >= c.

K is assembled from the element pseudo-densities and element stiffness matrices and Ksigma is assembled from the element stress stiffness matrices. In Ferrite and TopOpt.jl now, the boundary conditions are applied to K by making the row and column corresponding to the i^th Dirichlet bounded degree of freedom zeros and then instead of eliminating the row and column, the average diagonal value is assigned to the diagonal element of K and the RHS f is adjusted such that the i^th degree of freedom in u = K \ f has the intended Dirichlet boundary condition value. So if K_ii = 5.0 and u_i should be 2.0 then f_i = 10.0.

But what should be done with the i^th column and row of Ksigma??? Well, turns out the correct behavior in this case is to zero out the row and column WITHOUT assigning anything to the diagonal entry. Let's see why. If the i^th diagonal entry of Ksigma is 0, we can pull out an eigenvector there and then which is i^th basis vector e_i with 0 everywhere except at the i^th degree of freedom. It's obvious this is an eigenvector associated with the 0 eigenvalue because:

-Ksigma * e_i = 0 * K * e_i

In other words, the 0 eigenvalue's multiplicity was increased by 1. Since we are only interested in the maximum eigenvalue of the system:

-Ksigma * v = lambda * K * v

the 0 eigenvalue is irrelevant. Now we know that all the other eigenvectors will be K-orthonormal to e_i, i.e. v' * K * e_i = v_i * K_ii = 0. Since K_ii was assigned a positive value earlier, v_i is guaranteed to be 0 in any eigenvector associated with a non-zero eigenvalue! This is the intended behavior. Why you might ask? Great question! Keep reading.

In buckling analysis, the buckling modes represent "perturbation" directions along which there is 0 energy change in the nonlinear strain model. However by definition if a degree of freedom is fixed to a value in a Dirichlet boundary condition, it is not allowed to change and is therefore irrelevant to the buckling analysis. This means that we need to ensure our eigenvector associated with the maximum eigenvalue doesn't have a non-zero value in ith element because that would violate our boundary conditions and would therefore not be a failure mode to the constrained system.

When using the positive semidefinite constraint formulation when Ksigma_ii is 0, this degree of freedom will never cause a violation to the constraint because c > 0 and K_ii > 0. So the approach is correct in this case as well.