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) )
2026-01-29 15:10:33,231 WARNING  Discrepancy          :: Initialization of Discrepancy with two positional arguments (norm, data) deprecated. Use one positional argument (noiselevel) and specifiy norm and data via a keyword argument setting!
2026-01-29 15:10:33,231 WARNING Discrepancy          :: Initialization of Discrepancy with two positional arguments (norm, data) deprecated. Use one positional argument (noiselevel) and specifiy norm and data via a keyword argument setting!
2026-01-29 15:10:33,236 INFO     CombineRules         :: it. 0>=50 | discr./noiselevel = 6.98< 1.30
2026-01-29 15:10:33,236 INFO CombineRules         :: it. 0>=50 | discr./noiselevel = 6.98< 1.30
2026-01-29 15:10:33,328 INFO     CombineRules         :: it. 1>=50 | discr./noiselevel = 6.13< 1.30
2026-01-29 15:10:33,328 INFO CombineRules         :: it. 1>=50 | discr./noiselevel = 6.13< 1.30
2026-01-29 15:10:33,380 INFO     CombineRules         :: it. 2>=50 | discr./noiselevel = 5.50< 1.30
2026-01-29 15:10:33,380 INFO CombineRules         :: it. 2>=50 | discr./noiselevel = 5.50< 1.30
2026-01-29 15:10:33,431 INFO     CombineRules         :: it. 3>=50 | discr./noiselevel = 4.96< 1.30
2026-01-29 15:10:33,431 INFO CombineRules         :: it. 3>=50 | discr./noiselevel = 4.96< 1.30
2026-01-29 15:10:33,482 INFO     CombineRules         :: it. 4>=50 | discr./noiselevel = 4.50< 1.30
2026-01-29 15:10:33,482 INFO CombineRules         :: it. 4>=50 | discr./noiselevel = 4.50< 1.30
2026-01-29 15:10:33,534 INFO     CombineRules         :: it. 5>=50 | discr./noiselevel = 4.09< 1.30
2026-01-29 15:10:33,534 INFO CombineRules         :: it. 5>=50 | discr./noiselevel = 4.09< 1.30
2026-01-29 15:10:33,584 INFO     CombineRules         :: it. 6>=50 | discr./noiselevel = 3.74< 1.30
2026-01-29 15:10:33,584 INFO CombineRules         :: it. 6>=50 | discr./noiselevel = 3.74< 1.30
2026-01-29 15:10:33,637 INFO     CombineRules         :: it. 7>=50 | discr./noiselevel = 3.43< 1.30
2026-01-29 15:10:33,637 INFO CombineRules         :: it. 7>=50 | discr./noiselevel = 3.43< 1.30
2026-01-29 15:10:33,688 INFO     CombineRules         :: it. 8>=50 | discr./noiselevel = 3.16< 1.30
2026-01-29 15:10:33,688 INFO CombineRules         :: it. 8>=50 | discr./noiselevel = 3.16< 1.30
2026-01-29 15:10:33,738 INFO     CombineRules         :: it. 9>=50 | discr./noiselevel = 2.92< 1.30
2026-01-29 15:10:33,738 INFO CombineRules         :: it. 9>=50 | discr./noiselevel = 2.92< 1.30
2026-01-29 15:10:33,790 INFO     CombineRules         :: it. 10>=50 | discr./noiselevel = 2.71< 1.30
2026-01-29 15:10:33,790 INFO CombineRules         :: it. 10>=50 | discr./noiselevel = 2.71< 1.30
2026-01-29 15:10:33,842 INFO     CombineRules         :: it. 11>=50 | discr./noiselevel = 2.53< 1.30
2026-01-29 15:10:33,842 INFO CombineRules         :: it. 11>=50 | discr./noiselevel = 2.53< 1.30
2026-01-29 15:10:33,894 INFO     CombineRules         :: it. 12>=50 | discr./noiselevel = 2.37< 1.30
2026-01-29 15:10:33,894 INFO CombineRules         :: it. 12>=50 | discr./noiselevel = 2.37< 1.30
2026-01-29 15:10:33,946 INFO     CombineRules         :: it. 13>=50 | discr./noiselevel = 2.22< 1.30
2026-01-29 15:10:33,946 INFO CombineRules         :: it. 13>=50 | discr./noiselevel = 2.22< 1.30
2026-01-29 15:10:33,999 INFO     CombineRules         :: it. 14>=50 | discr./noiselevel = 2.10< 1.30
2026-01-29 15:10:33,999 INFO CombineRules         :: it. 14>=50 | discr./noiselevel = 2.10< 1.30
2026-01-29 15:10:34,051 INFO     CombineRules         :: it. 15>=50 | discr./noiselevel = 1.99< 1.30
2026-01-29 15:10:34,051 INFO CombineRules         :: it. 15>=50 | discr./noiselevel = 1.99< 1.30
2026-01-29 15:10:34,101 INFO     CombineRules         :: it. 16>=50 | discr./noiselevel = 1.90< 1.30
2026-01-29 15:10:34,101 INFO CombineRules         :: it. 16>=50 | discr./noiselevel = 1.90< 1.30
2026-01-29 15:10:34,151 INFO     CombineRules         :: it. 17>=50 | discr./noiselevel = 1.82< 1.30
2026-01-29 15:10:34,151 INFO CombineRules         :: it. 17>=50 | discr./noiselevel = 1.82< 1.30
2026-01-29 15:10:34,202 INFO     CombineRules         :: it. 18>=50 | discr./noiselevel = 1.74< 1.30
2026-01-29 15:10:34,202 INFO CombineRules         :: it. 18>=50 | discr./noiselevel = 1.74< 1.30
2026-01-29 15:10:34,254 INFO     CombineRules         :: it. 19>=50 | discr./noiselevel = 1.68< 1.30
2026-01-29 15:10:34,254 INFO CombineRules         :: it. 19>=50 | discr./noiselevel = 1.68< 1.30
2026-01-29 15:10:34,305 INFO     CombineRules         :: it. 20>=50 | discr./noiselevel = 1.62< 1.30
2026-01-29 15:10:34,305 INFO CombineRules         :: it. 20>=50 | discr./noiselevel = 1.62< 1.30
2026-01-29 15:10:34,354 INFO     CombineRules         :: it. 21>=50 | discr./noiselevel = 1.58< 1.30
2026-01-29 15:10:34,354 INFO CombineRules         :: it. 21>=50 | discr./noiselevel = 1.58< 1.30
2026-01-29 15:10:34,404 INFO     CombineRules         :: it. 22>=50 | discr./noiselevel = 1.53< 1.30
2026-01-29 15:10:34,404 INFO CombineRules         :: it. 22>=50 | discr./noiselevel = 1.53< 1.30
2026-01-29 15:10:34,455 INFO     CombineRules         :: it. 23>=50 | discr./noiselevel = 1.50< 1.30
2026-01-29 15:10:34,455 INFO CombineRules         :: it. 23>=50 | discr./noiselevel = 1.50< 1.30
2026-01-29 15:10:34,505 INFO     CombineRules         :: it. 24>=50 | discr./noiselevel = 1.46< 1.30
2026-01-29 15:10:34,505 INFO CombineRules         :: it. 24>=50 | discr./noiselevel = 1.46< 1.30
2026-01-29 15:10:34,555 INFO     CombineRules         :: it. 25>=50 | discr./noiselevel = 1.44< 1.30
2026-01-29 15:10:34,555 INFO CombineRules         :: it. 25>=50 | discr./noiselevel = 1.44< 1.30
2026-01-29 15:10:34,606 INFO     CombineRules         :: it. 26>=50 | discr./noiselevel = 1.41< 1.30
2026-01-29 15:10:34,606 INFO CombineRules         :: it. 26>=50 | discr./noiselevel = 1.41< 1.30
2026-01-29 15:10:34,655 INFO     CombineRules         :: it. 27>=50 | discr./noiselevel = 1.39< 1.30
2026-01-29 15:10:34,655 INFO CombineRules         :: it. 27>=50 | discr./noiselevel = 1.39< 1.30
2026-01-29 15:10:34,706 INFO     CombineRules         :: it. 28>=50 | discr./noiselevel = 1.37< 1.30
2026-01-29 15:10:34,706 INFO CombineRules         :: it. 28>=50 | discr./noiselevel = 1.37< 1.30
2026-01-29 15:10:34,756 INFO     CombineRules         :: it. 29>=50 | discr./noiselevel = 1.35< 1.30
2026-01-29 15:10:34,756 INFO CombineRules         :: it. 29>=50 | discr./noiselevel = 1.35< 1.30
2026-01-29 15:10:34,807 INFO     CombineRules         :: it. 30>=50 | discr./noiselevel = 1.33< 1.30
2026-01-29 15:10:34,807 INFO CombineRules         :: it. 30>=50 | discr./noiselevel = 1.33< 1.30
2026-01-29 15:10:34,857 INFO     CombineRules         :: it. 31>=50 | discr./noiselevel = 1.32< 1.30
2026-01-29 15:10:34,857 INFO CombineRules         :: it. 31>=50 | discr./noiselevel = 1.32< 1.30
2026-01-29 15:10:34,908 INFO     CombineRules         :: it. 32>=50 | discr./noiselevel = 1.31< 1.30
2026-01-29 15:10:34,908 INFO CombineRules         :: it. 32>=50 | discr./noiselevel = 1.31< 1.30
2026-01-29 15:10:34,961 INFO     CombineRules         :: it. 33>=50 | discr./noiselevel = 1.30< 1.30
Rule Discrepancy(noiselevel=1.8788546038475205e-05, tau=1.3) triggered.
2026-01-29 15:10:34,961 INFO CombineRules         :: it. 33>=50 | discr./noiselevel = 1.30< 1.30
Rule Discrepancy(noiselevel=1.8788546038475205e-05, tau=1.3) triggered.

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.7729625344089 %

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