# Solving linear elasticity with a near nullspace

In this tutorial, we explore a linear elasticity discretisation. As we did in Solving the Poisson problem with different preconditioning strategies, we will solve the linear system arising from the discretisation using a PETSc KSP. We begin by creating a discretisation for the weak formulation of linear elasticity with Lamé coefficients \(\mu\) and \(\lambda\), i.e. find \(\vec{u}\in [H^1_0(\Omega)]^d\) such that

\[a(u,v) := 2\mu \int_{\Omega} \epsilon(\vec{u}) : \epsilon(\vec{v}) \; d\vec{x} + \lambda \int_\Omega (\nabla \cdot \vec{u})\; d\vec{x} = L(v) := \int_{\Omega} fv\; d\vec{x}\qquad \vec{v}\in [H^1_0(\Omega)]^d.\]

We can easily discretise this problem using NGSolve:

    from ngsolve import *
    import netgen.gui
    import netgen.meshing as ngm
    from mpi4py.MPI import COMM_WORLD

    mesh = Mesh(unit_square.GenerateMesh(maxh=0.1,comm=COMM_WORLD))

    E, nu = 210, 0.2
    mu  = E / 2 / (1+nu)
    lam = E * nu / ((1+nu)*(1-2*nu))

    def Stress(strain):
       return 2*mu*strain + lam*Trace(strain)*Id(2)

    fes = VectorH1(mesh, order=1, dirichlet="left")
    u,v = fes.TnT()

    a = BilinearForm(InnerProduct(Stress(Sym(Grad(u))), Sym(Grad(v)))*dx)
    a.Assemble()

    force = CF( (0,1) )
    f = LinearForm(force*v*ds("right")).Assemble()

We begin solving the linear system using PETSc’s own implementation of an algebraic multigrid preconditioner.

    from ngsPETSc import KrylovSolver
    opts = {'ksp_type': 'cg',
            'ksp_monitor': None,
            'pc_type': 'gamg'}
    solver = KrylovSolver(a,fes, solverParameters=opts)
    gfu = GridFunction(fes, name="gfu")
    solver.solve(f.vec, gfu.vec)
    Draw (gfu)

Preconditioner performance for PETSc GAMG.
Preconditioner | Iterations
--- | ---
PETSc GAMG | 19 (6.27e-08)

To improve the performance of PETSc GAMG we pass it the near-null space composed of the rigid body motions. Using the `near` flag we tell PETSc KSP to pass the nullspace as a near nullspace to PETSc GAMG.

    from ngsPETSc import NullSpace
    rbms = []
    for val in [(1,0), (0,1), (-y,x)]:
       rbm = GridFunction(fes)
       rbm.Set(CF(val))
       rbms.append(rbm.vec)
    nullspace = NullSpace(fes, rbms, near=True)
    opts = {'ksp_type': 'cg',
            'pc_monitor': None,
            'pc_type': 'gamg'}
    solver = KrylovSolver(a,fes, solverParameters=opts, nullspace=nullspace)
    gfu = GridFunction(fes, name="gfu_near")
    solver.solve(f.vec, gfu.vec)
    Draw (gfu)

Preconditioner performance for PETSc GAMG.
Preconditioner | Iterations
--- | ---
PETSc GAMG | 19 (6.27e-08)
PETSc GAMG (with near nullspace) | 6 (4.55e-07)