9
$\begingroup$

I want to solve $K u = b$ where $K$ is my stiffness matrix. However some constraints may be missing an therefore some rigid body motion may be still present in the system (due to eigenvalue zero). Since I'm using CG for solving the linear system this is not acceptable since sometimes CG doesn't converge on semi-positive problems (but I may sometimes converges).

Actually I'm using a penalized displacement approach in the sense that I'm adding a penalty of the form $ \alpha ||u||^2$ to the elastic energy. So the energy reads \begin{equation} \mathcal W(u) := \frac{1}{2} u^T (K + \alpha I) u - b^t u \end{equation} where $\alpha$ taken as a proportional to some diagonal entry of the stiffness matrix. But actually this has the effect to damp some deformation mode that I would sometime like to have.

Some my question is:

a) could I transform the original system so has to make it free of singularity and positive definite (such as coordinate transformation or congruence transformation or whatever) ? My idea is to use such transformation to still use CG on the transformed problem

b) Is there any standard way to deal with those singularities ?

Thank you very much !

Kind regards,

Tom

|cite|improve this question
$\endgroup$
6
$\begingroup$

The standard way is to add the constraint $u(x_0)=0$ for an arbitrarily chosen node $x_0$. This makes sure that your body can't translate or rotate and therefore takes away the zero eigenvalue. The resulting system with this constraint is positive definite even without your penalty term.

|cite|improve this answer
$\endgroup$
  • 1
    $\begingroup$ Thank you! Yep, but I'm my case I have several floating substructures and I can't tell which nodes (3 non colinear nodes in 3D) to fix. This is why I'm wondering if there is no higher-level solution since in my case the null space is well known. $\endgroup$ – Tom Oct 22 '12 at 12:15
  • $\begingroup$ If you have several structures then you need to fix one node for each structure. It doesn't matter which one, just pick one per structure. $\endgroup$ – Wolfgang Bangerth Oct 22 '12 at 13:44
  • 3
    $\begingroup$ @WolfgangBangerth This is 3D elasticity, therefore you would need to pin three non-colinear points to control the null space of dimension 6. Pinning those three displacements is a rank 9 perturbation and it's not easy to ensure that the rank-3 modification beyond the null space does not change the solution. For any choice of points and values to pin, there is a 3-dimensional family of right hand sides in which your pinned problem gives the correct answer for only one member. $\endgroup$ – Jed Brown Oct 22 '12 at 13:47
  • $\begingroup$ No, you can't pin 3 points for 9 constraints because you would then also fix their relative distances. If your boundary conditions really don't provide any other constraints (e.g. if they are no normal displacement on a circle) then you need to fix 1 point + various rotation angles on two other points for a total of 6 constraints. $\endgroup$ – Wolfgang Bangerth Oct 22 '12 at 16:31
6
$\begingroup$

If you know the null space, you can make the right hand side compatible and have the Krylov method prevent the preconditioner from causing pollution, see Why is pinning a point to remove a null space bad? for further discussion. In PETSc, this is done using the MatNullSpace object. Note that you can provide your own function to project out the null space, which would be useful to reduce projection cost when you have many floating structures.

If you do not know the null space and cannot avoid an incompatible right hand side, there are specialized Krylov methods like MINRES-QLP that can find the minimum norm solution despite. This approach may be useful if you have hinges and single-point connections that only couple some modes. Note that you must still be careful about the preconditioner causing pollution (e.g., due to LU factorization finding zero pivots, perhaps on a coarse level of multigrid).

|cite|improve this answer
$\endgroup$
  • $\begingroup$ Thank you Jed ! I thought about removing by projection the null space directly in my iterative method. But I was wondering if it was not too costly (I can create a operator that projects out the null space since it is really trivial in elasticity). Also I think that the residual should be projected as well ? $\endgroup$ – Tom Oct 22 '12 at 14:03
  • 2
    $\begingroup$ Make the right hand side compatible and project out the null space after each preconditioner application (because many preconditioners will contaminate the null space). This yields a Krylov operator $K = (I -N)P^{-1}A$ such that $\{ b, K b, K^2 b, \dotsc \}$ is orthogonal to the null space. Since your problem is symmetric, you don't need a different routine for the left and right null spaces. This is what PETSc does if you call MatSetNullSpace(). $\endgroup$ – Jed Brown Oct 22 '12 at 14:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.