Traction force microscopy - an ngsolve example with Landweber iteration

To run this example, regpy and ngsolve (including the extension for jupyter notebooks) have to be installed. How to install ngsolve can be found here.

Mathematical Model

Traction force microscopy (TFM) is a method from biophysics to measure cell forces. The mathematical model relies on the 3D boundary value problem from elasticity. The forward operator from TFM maps the traction forces \(t\) defined on part of the boundary to the displacement \(u\) defined on the whole domain by solving the following PDE

\[div (\sigma(u)) = 0 \quad \text{in } \Omega\]

\[\sigma(u) n = t \quad \text{on } \Gamma_{Top}\]

\[\sigma(u) n = 0 \quad \text{on } \Gamma_{Sides}\]

\[u = 0 \quad \text{on } \Gamma_{Bottom}\]

where \(\sigma\) is given by a material law, in this case the linear Hooke’s law

\[\sigma(u) = 2 \mu \varepsilon + \lambda tr( \varepsilon) I\]

for the Lamé parameters \(\mu\), \(\lambda\) and \(\varepsilon\) is the linearized strain tensor

\[\varepsilon = \frac{1}{2} (\nabla u + (\nabla u)^T )\]

In the weak formulation we get

\[a(u,v) = l_t(v) , \quad u,v \in H^1_{0,\Gamma_D}(\Omega,\mathbb{R}^3)\]

for the bilinear form

\[a(u,v) = \int_\Omega 2 \mu \varepsilon(u) : \varepsilon(v) + \lambda div(u) div(v) dx\]

and the linear form

\[l_t(v) = \int_{\Gamma_T} t v ds\]

import logging
from ngsolve.webgui import Draw
from netgen.occ import *

from regpy.util import set_rng_seed

set_rng_seed(42)

logging.basicConfig(
    level=logging.INFO,
    format='%(asctime)s %(levelname)s %(name)-20s :: %(message)s'
)

Implementation

We can use the weak formulation to define the forward operator creating a subclass `regpy.operators.ngsolve.NGSolveOperator in the next code block. In the following creating new operators of that type can simply be done by initiating with that class.

Therefore, we create a subclass of regpy.operators.ngsolve.NGSolveOperator type in the next code block. Since the operator is linear, we only need to implement the functions __init__ (for initialization when a new instance of the operator is created), _eval (for evaluating the forward operator for a given traction force), and _adjoint (for evaluating the adjoint operator for a given displacement).

In the following create new operators of that type can simply be done by initiating with that class. This class will have the following parameters. (Side remark you may also do the definition in a separate file for this example check out tfm.py in this example folder. )

---

Parameters

domain : regpy.vecsps.ngsolve.NgsSpace > The domain on which the operator is defined. Should be two- or threedimensional and boundary bdr should be set to \(\Gamma_{Top}\), the part of the boundary where the traction force applies to.

codomain : regpy.vecsps.ngsolve.NgsSpace > The codomain on which the operator is defined. Should be as domain, but without bdrdefined.

mu,lambda: float > Lamé parameters used in Hooke’s law

import ngsolve as ngs
from regpy.operators import NgsOperator

def strain(u):
    return 0.5 * (ngs.Grad(u) + ngs.Grad(u).trans)

class TFM(NgsOperator):

    def __init__(self, domain, codomain, mu, lam):
        codomain = codomain
        # Need to know the boundary to calculate Neumann bdr condition
        assert domain.bdr is not None
        super().__init__(domain, codomain, linear=True)

        # Lamé Parameters for substrate
        self.mu = mu
        self.lam = lam

        self.fes_domain = domain.fes
        self.fes_codomain = codomain.fes


        # grid functions for later use
        self.gfu_eval = ngs.GridFunction(self.fes_codomain)  # solution, return value of _eval
        self.gfu_adjoint = ngs.GridFunction(self.fes_domain) # grid function return value of adjoint (trace of gfu_adjoint_sol)

        self.gfu_bf = ngs.GridFunction(self.fes_codomain)  # grid function for defining integrator (bilinearform)
        self.gfu_lf = ngs.GridFunction(self.fes_codomain) # grid function for defining right hand side (linearform), f



        u, v = self.fes_codomain.TnT()

        # Define Bilinearform, will be assembled later
        self.a = ngs.BilinearForm(self.fes_codomain, symmetric=True)
        self.a += (2 * mu * ngs.InnerProduct(strain(u), strain(v)) + lam * ngs.div(u) * ngs.div(v)) * ngs.dx


        # Define Linearform for evaluation, will be assembled later
        self.b = ngs.LinearForm(self.fes_domain)
        self.b += self.gfu_lf * v * ngs.ds(domain.bdr)


        # Define linearform to trick ngsolve for computation of discrete adjoint
        self.b_help = ngs.LinearForm(self.fes_domain)
        self.b_help.Assemble()


        # Initialize preconditioner for solving the Dirichlet problems by ngs.BVP
        self.prec = ngs.Preconditioner(self.a, 'direct')
        self.a.Assemble()


    # Left term: Bilinearform self.a
    # Right term: Linearform self.b
    def _eval(self, traction):

        # Assemble Linearform, boundary term
        self.gfu_lf.vec.data = traction.vec
        self.b.Assemble()
        self._solve_dirichlet_problem(bf=self.a, lf=self.b, gf=self.gfu_eval,pre=self.prec, needsassembling=False, print=False)

        return self.codomain.from_ngs(self.gfu_eval, copy=True)


    def _adjoint(self, displacement):
        # Bilinearform already assembled in init -> initialization with 0, s.t. object exists
        # Diskrete Adjoint w.r.t. standard inner product

        self.b_help.vec.data = displacement.vec
        self._solve_dirichlet_problem(bf=self.a, lf=self.b_help, gf=self.gfu_adjoint, pre=self.prec, needsassembling=False, print=False)
        self.gfu_lf.vec.data = self.gfu_adjoint.vec
        self.b.Assemble()

        return self.domain.from_ngs(self.b.vec, copy=True)

Create Operators

First we import logging to later get printed messages during our solution algorithm.

To define domain and codomain of our operator, we need regpy.vecsps.ngsolve.NgsSpace. An NgsSpace in regpy consists of a finite element space from ngsolve and possibly a boundary. First we generate a cuboid as the mesh and label the top and bottom boundary. Then we can define the three dimensional H1 finite element space with dirichlet boundary conditions at the bottom. Finally we can define domain and codomain of our operator.

To prevent commiting an inverse crime, we use to different discretizations for op_gen - used to generate simulated data (the displacement) - and op_rec - used to reconstruct the traction forces.

from regpy.vecsps import NgsVectorSpace

'''############################### material parameters ###############################'''

# E Young's modulus, nu Poisson ratio
E, nu = 10000, 0.45
# Lamé parameters
mu = E / 2 / (1 + nu)
lam = E * nu / ((1 + nu) * (1 - 2 * nu))

'''############################### mesh of substrate ################################'''

# generate two triangular meshes of mesh-size h, order k, for generating and reconstructing traction forces
h_gen = 0.3
h_rec = 0.4
k = 3


box = Box(Pnt(-2, -0.3,-2), Pnt(2, 0.3,2))
box.faces.Max(Y).name = 'top'
box.faces.Min(Y).name = 'bottom'
geo = OCCGeometry(box)

mesh_gen = ngs.Mesh(geo.GenerateMesh(maxh=h_gen))
mesh_gen.Curve(3)

mesh_rec = ngs.Mesh(geo.GenerateMesh(maxh=h_rec))
mesh_rec.Curve(3)

'''############################### define operator ################################'''
# operator for generating data on a different discretization than reconstruction operator

fes_domain_gen = ngs.VectorH1(mesh_gen, order=k, dirichlet="bottom")
domain_gen = NgsVectorSpace(fes_domain_gen, bdr='top')

fes_codomain_gen = ngs.VectorH1(mesh_gen, order=k, dirichlet="bottom")
codomain_gen = NgsVectorSpace(fes_codomain_gen)

op_gen = TFM(domain_gen, codomain=codomain_gen, mu=mu, lam=lam)


fes_domain_rec = ngs.VectorH1(mesh_rec, order=k, dirichlet="bottom")
domain_rec = NgsVectorSpace(fes_domain_rec, bdr='top')

fes_codomain_rec = ngs.VectorH1(mesh_rec, order=k, dirichlet="bottom")
codomain_rec = NgsVectorSpace(fes_codomain_rec)

op_rec = TFM(domain_rec, codomain=codomain_rec, mu=mu, lam=lam)

Generate Groundtruth

In this step we generate a toy surface traction force first as a CoefficientFunctionfrom ngsolve. Because the surface traction force should just act on the boundary, we interpolate traction_true_cf to the top boundary of the generated mesh. The resulting gridfunction traction_true_gf acts then only on the top boundary.

To change between ngsolve and regpy, the functions to_ngs and from_ngs can be used.

a = 4
c = 1

force_intensity = (-a * (ngs.sqrt(ngs.x**2 + ngs.z**2) - 0.5)**2 + c )

cell_force = ngs.CF(force_intensity * (-ngs.x/(ngs.sqrt(ngs.x**2 + ngs.z**2)),0,-ngs.z/(ngs.sqrt(ngs.x**2 + ngs.z**2))))

cell_xz = -(ngs.sqrt(ngs.x ** 2 + ngs.z ** 2) - 1)


traction_true_cf = ngs.IfPos(c1=cell_xz, then_obj=cell_force, else_obj=ngs.CF((0, 0, 0)))


traction_true = domain_gen.from_ngs(traction_true_cf,definedon=mesh_gen.Boundaries( "top" ))

Draw(domain_gen.to_gf(traction_true), mesh_gen,settings={"camera": {"transformations": [{"type": "rotateY", "angle": 45},{"type": "rotateX", "angle": 45}]}})

BaseWebGuiScene

Generate Data

Next, we generate the displacement by applying the forward operator op_gen to the true traction force. Then we add random noise to the displacement and plot noiseless displacement on the generating mesh and noisy displacement on the reconstructing mesh.

displacement_true = op_gen(traction_true)

gf = codomain_gen.to_gf(displacement_true)

Draw(gf,settings={"camera": {"transformations": [{"type": "rotateY", "angle": 45},{"type": "rotateX", "angle": 45}]}})

BaseWebGuiScene

noise = 3e-6 * codomain_rec.randn()
data = codomain_rec.from_ngs(gf) + noise


Draw(codomain_rec.to_gf(data),settings={"camera": {"transformations": [{"type": "rotateY", "angle": 45},{"type": "rotateX", "angle": 45}]}})

BaseWebGuiScene

Solve the Inverse Problem with regularization

Now, with available data we can use, e.g., landweber iteration to reconstruct the groundtruth traction force from the noisy displacement.

First, we can define the regularization setting using Setting. To measure the error, we have to choose a penalty Hilbert space structure. Here we choose penalty=L2Boundary because the traction forces are elements of \(L^2(\Gamma_{Top})\). Similar to measure the data misfit, we have to choose a data fidelity Hilbert space. Here we choose our Hilbert space data_fid=Hm0 because the displacement is an element of \(H^1_{0,\Gamma_{Bottom}}(\Omega)\). L2Boundary and Hm0 are AbstractSpaces and the according implementation can be found in regpy.hilbert.ngsolve.

As a next step we define the we choose an initial guess, in this case we choose zero.

With this we can use Landweber iteration implemented in regpy.solvers.nonlinear.landweber to reconstruct the exact solution from the above constructed noisy data data.

We stop the iterative algorithm after at most \(100\) steps and have as early stopping criterion the discrepancy rule implemented. This can be easily done by summing the two instances of the regpy.stoprules.

After everything is defined run the solver with the specified stopping rule using the method run() of the solver.

from regpy.solvers import Setting
from regpy.solvers.linear import Landweber
import regpy.stoprules as rules
from regpy.hilbert import L2Boundary, Hm0

setting = Setting(op=op_rec, penalty=L2Boundary, data_fid=Hm0, data=data)
init = domain_rec.from_ngs((0, 0, 0))

landweber = Landweber(setting, init)

stoprule = (
        rules.CountIterations(50) +
        rules.Discrepancy(setting.h_codomain.norm, data, noiselevel=setting.h_codomain.norm(noise), tau=1.3)
)

reco, reco_data = landweber.run(stoprule)

reco_gf = domain_rec.to_gf(reco)

used dof inconsistency
(silence this warning by setting BilinearForm(...check_unused=False) )

Compute Error

Now we can calculate the relative error on \(L^2(\Gamma_{Top})\).

# calculate relative L2 error on boundary
err_abs = ngs.sqrt(ngs.Integrate((reco_gf - traction_true_cf)**2, mesh_rec,definedon=mesh_rec.Boundaries('top')))
norm_true = ngs.sqrt(ngs.Integrate((traction_true_cf)**2, mesh_rec,definedon=mesh_rec.Boundaries('top')))
norm_rec = ngs.sqrt(ngs.Integrate((reco_gf)**2, mesh_rec,definedon=mesh_rec.Boundaries('top')))
err_rel = err_abs/norm_true * 100

print('relative error:' , err_rel, '%')

relative error: 19.773045926237202 %

Plot Reconstruction

We plot the true traction forces on the reconstruction mesh, the reconstruction and the error.

Draw(traction_true_cf, domain_rec.fes.mesh,settings={"camera": {"transformations": [{"type": "rotateY", "angle": 45},{"type": "rotateX", "angle": 45}]}})

Draw(reco_gf,settings={"camera": {"transformations": [{"type": "rotateY", "angle": 45},{"type": "rotateX", "angle": 45}]}})

error = traction_true_cf - reco_gf

Draw(error, domain_rec.fes.mesh,settings={"camera": {"transformations": [{"type": "rotateY", "angle": 45},{"type": "rotateX", "angle": 45}]}})

BaseWebGuiScene