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\]
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
[1]:
import ngsolve as ngs
from regpy.operators.ngsolve import NGSolveOperator
def strain(u):
return 0.5 * (ngs.Grad(u) + ngs.Grad(u).trans)
class TFM(NGSolveOperator):
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)
# From NgSolveOperator
#self.gfu_read_in = ngs.GridFunction(self.domain.fes)
# 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, differentiate=False):
# Assemble Linearform, boundary term
self._read_in(traction, self.gfu_lf)
self.b.Assemble()
self._solve_dirichlet_problem(bf=self.a, lf=self.b, gf=self.gfu_eval, prec=self.prec)
return self.gfu_eval.vec.FV().NumPy()[:].copy()
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.FV().NumPy()[:] = displacement.copy()
self._solve_dirichlet_problem(bf=self.a, lf=self.b_help, gf=self.gfu_adjoint, prec=self.prec)
self._read_in(self.gfu_adjoint.vec.FV().NumPy()[:].copy(), self.gfu_lf)
self.b.Assemble()
return self.b.vec.FV().NumPy()[:].copy()
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.
[2]:
import logging
from ngsolve.webgui import Draw
from netgen.occ import *
from regpy.vecsps.ngsolve import NgsSpace
logging.basicConfig(
level=logging.INFO,
format='%(asctime)s %(levelname)s %(name)-20s :: %(message)s'
)
'''############################### 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)
Draw(mesh_gen)
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")
# fes_domain = ngs.H1(mesh, order=k, dim=3)
domain_gen = NgsSpace(fes_domain_gen, bdr='top')
fes_codomain_gen = ngs.VectorH1(mesh_gen, order=k, dirichlet="bottom")
# fes_codomain = ngs.H1(mesh, order=k, dirichlet="bottom|side", dim=3)
codomain_gen = NgsSpace(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")
# fes_domain = ngs.H1(mesh, order=k, dim=3)
domain_rec = NgsSpace(fes_domain_rec, bdr='top')
fes_codomain_rec = ngs.VectorH1(mesh_rec, order=k, dirichlet="bottom")
# fes_codomain = ngs.H1(mesh, order=k, dirichlet="bottom|side", dim=3)
codomain_rec = NgsSpace(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.
[3]:
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_gf = ngs.GridFunction(fes_domain_gen)
traction_true_gf.Set(traction_true_cf, definedon = mesh_gen.Boundaries( "top" ))
traction_true = domain_gen.from_ngs(traction_true_gf)
Draw(traction_true_gf, mesh_gen)
[3]:
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.
[4]:
import numpy as np
displacement_true = op_gen(traction_true)
Draw(codomain_gen.to_ngs(displacement_true))
np.random.seed(42)
noise = 1e-6 * codomain_rec.randn()
data_gf = ngs.GridFunction(fes_codomain_rec)
data_gf.Set(codomain_gen.to_ngs(displacement_true))
data = codomain_rec.from_ngs(data_gf) + noise
Draw(codomain_rec.to_ngs(data), mesh_rec)
[4]:
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 RegularizationSetting. 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.
[5]:
from regpy.solvers import RegularizationSetting
from regpy.solvers.nonlinear.landweber import Landweber
import regpy.stoprules as rules
from regpy.hilbert import L2Boundary, Hm0
setting = RegularizationSetting(op=op_rec, penalty=L2Boundary, data_fid=Hm0)
init = domain_rec.from_ngs((0, 0, 0))
landweber = Landweber(setting, data, init)
stoprule = (
rules.CountIterations(100) +
rules.Discrepancy(setting.h_codomain.norm, data, noiselevel=setting.h_codomain.norm(noise), tau=2.2)
)
reco, reco_data = landweber.run(stoprule)
reco_gf = domain_rec.to_ngs(reco)
used dof inconsistency
(silence this warning by setting BilinearForm(...check_unused=False) )
2025-08-13 09:04:01,379 INFO CountIterations :: iteration = 1 / 100
2025-08-13 09:04:01,391 INFO Discrepancy :: relative discrepancy = 15.21, tolerance = 2.20
2025-08-13 09:04:02,307 INFO Landweber :: |residual| = 0.00012965419945885036
2025-08-13 09:04:02,307 INFO CountIterations :: iteration = 2 / 100
2025-08-13 09:04:02,331 INFO Discrepancy :: relative discrepancy = 13.53, tolerance = 2.20
2025-08-13 09:04:02,631 INFO Landweber :: |residual| = 0.00011537303511922831
2025-08-13 09:04:02,631 INFO CountIterations :: iteration = 3 / 100
2025-08-13 09:04:02,655 INFO Discrepancy :: relative discrepancy = 12.22, tolerance = 2.20
2025-08-13 09:04:02,967 INFO Landweber :: |residual| = 0.000104226581952904
2025-08-13 09:04:02,967 INFO CountIterations :: iteration = 4 / 100
2025-08-13 09:04:02,991 INFO Discrepancy :: relative discrepancy = 11.11, tolerance = 2.20
2025-08-13 09:04:03,303 INFO Landweber :: |residual| = 9.474028475703115e-05
2025-08-13 09:04:03,303 INFO CountIterations :: iteration = 5 / 100
2025-08-13 09:04:03,327 INFO Discrepancy :: relative discrepancy = 10.14, tolerance = 2.20
2025-08-13 09:04:03,627 INFO Landweber :: |residual| = 8.64346752992881e-05
2025-08-13 09:04:03,627 INFO CountIterations :: iteration = 6 / 100
2025-08-13 09:04:03,651 INFO Discrepancy :: relative discrepancy = 9.28, tolerance = 2.20
2025-08-13 09:04:03,955 INFO Landweber :: |residual| = 7.909324848147849e-05
2025-08-13 09:04:03,955 INFO CountIterations :: iteration = 7 / 100
2025-08-13 09:04:03,979 INFO Discrepancy :: relative discrepancy = 8.51, tolerance = 2.20
2025-08-13 09:04:04,283 INFO Landweber :: |residual| = 7.257931332128824e-05
2025-08-13 09:04:04,283 INFO CountIterations :: iteration = 8 / 100
2025-08-13 09:04:04,307 INFO Discrepancy :: relative discrepancy = 7.83, tolerance = 2.20
2025-08-13 09:04:04,607 INFO Landweber :: |residual| = 6.678819381075419e-05
2025-08-13 09:04:04,607 INFO CountIterations :: iteration = 9 / 100
2025-08-13 09:04:04,631 INFO Discrepancy :: relative discrepancy = 7.23, tolerance = 2.20
2025-08-13 09:04:04,931 INFO Landweber :: |residual| = 6.16331925001455e-05
2025-08-13 09:04:04,932 INFO CountIterations :: iteration = 10 / 100
2025-08-13 09:04:04,951 INFO Discrepancy :: relative discrepancy = 6.69, tolerance = 2.20
2025-08-13 09:04:05,259 INFO Landweber :: |residual| = 5.704036255450052e-05
2025-08-13 09:04:05,259 INFO CountIterations :: iteration = 11 / 100
2025-08-13 09:04:05,283 INFO Discrepancy :: relative discrepancy = 6.21, tolerance = 2.20
2025-08-13 09:04:05,587 INFO Landweber :: |residual| = 5.294582345509228e-05
2025-08-13 09:04:05,588 INFO CountIterations :: iteration = 12 / 100
2025-08-13 09:04:05,611 INFO Discrepancy :: relative discrepancy = 5.78, tolerance = 2.20
2025-08-13 09:04:05,915 INFO Landweber :: |residual| = 4.9293978726e-05
2025-08-13 09:04:05,915 INFO CountIterations :: iteration = 13 / 100
2025-08-13 09:04:05,939 INFO Discrepancy :: relative discrepancy = 5.40, tolerance = 2.20
2025-08-13 09:04:06,239 INFO Landweber :: |residual| = 4.6036166345875785e-05
2025-08-13 09:04:06,239 INFO CountIterations :: iteration = 14 / 100
2025-08-13 09:04:06,263 INFO Discrepancy :: relative discrepancy = 5.06, tolerance = 2.20
2025-08-13 09:04:06,571 INFO Landweber :: |residual| = 4.3129582034981416e-05
2025-08-13 09:04:06,571 INFO CountIterations :: iteration = 15 / 100
2025-08-13 09:04:06,595 INFO Discrepancy :: relative discrepancy = 4.75, tolerance = 2.20
2025-08-13 09:04:06,891 INFO Landweber :: |residual| = 4.053640381821346e-05
2025-08-13 09:04:06,892 INFO CountIterations :: iteration = 16 / 100
2025-08-13 09:04:06,915 INFO Discrepancy :: relative discrepancy = 4.48, tolerance = 2.20
2025-08-13 09:04:07,227 INFO Landweber :: |residual| = 3.8223075381269e-05
2025-08-13 09:04:07,227 INFO CountIterations :: iteration = 17 / 100
2025-08-13 09:04:07,251 INFO Discrepancy :: relative discrepancy = 4.24, tolerance = 2.20
2025-08-13 09:04:07,551 INFO Landweber :: |residual| = 3.615971788198906e-05
2025-08-13 09:04:07,551 INFO CountIterations :: iteration = 18 / 100
2025-08-13 09:04:07,575 INFO Discrepancy :: relative discrepancy = 4.03, tolerance = 2.20
2025-08-13 09:04:07,871 INFO Landweber :: |residual| = 3.4319646447579665e-05
2025-08-13 09:04:07,871 INFO CountIterations :: iteration = 19 / 100
2025-08-13 09:04:07,895 INFO Discrepancy :: relative discrepancy = 3.83, tolerance = 2.20
2025-08-13 09:04:08,187 INFO Landweber :: |residual| = 3.267897191090638e-05
2025-08-13 09:04:08,187 INFO CountIterations :: iteration = 20 / 100
2025-08-13 09:04:08,211 INFO Discrepancy :: relative discrepancy = 3.66, tolerance = 2.20
2025-08-13 09:04:08,503 INFO Landweber :: |residual| = 3.121627154517119e-05
2025-08-13 09:04:08,503 INFO CountIterations :: iteration = 21 / 100
2025-08-13 09:04:08,527 INFO Discrepancy :: relative discrepancy = 3.51, tolerance = 2.20
2025-08-13 09:04:08,847 INFO Landweber :: |residual| = 2.9912315126957238e-05
2025-08-13 09:04:08,847 INFO CountIterations :: iteration = 22 / 100
2025-08-13 09:04:08,871 INFO Discrepancy :: relative discrepancy = 3.37, tolerance = 2.20
2025-08-13 09:04:09,163 INFO Landweber :: |residual| = 2.8749834844634353e-05
2025-08-13 09:04:09,163 INFO CountIterations :: iteration = 23 / 100
2025-08-13 09:04:09,187 INFO Discrepancy :: relative discrepancy = 3.25, tolerance = 2.20
2025-08-13 09:04:09,483 INFO Landweber :: |residual| = 2.7713329513157385e-05
2025-08-13 09:04:09,483 INFO CountIterations :: iteration = 24 / 100
2025-08-13 09:04:09,507 INFO Discrepancy :: relative discrepancy = 3.14, tolerance = 2.20
2025-08-13 09:04:09,803 INFO Landweber :: |residual| = 2.6788895326422256e-05
2025-08-13 09:04:09,803 INFO CountIterations :: iteration = 25 / 100
2025-08-13 09:04:09,827 INFO Discrepancy :: relative discrepancy = 3.05, tolerance = 2.20
2025-08-13 09:04:10,147 INFO Landweber :: |residual| = 2.596407699375806e-05
2025-08-13 09:04:10,147 INFO CountIterations :: iteration = 26 / 100
2025-08-13 09:04:10,171 INFO Discrepancy :: relative discrepancy = 2.96, tolerance = 2.20
2025-08-13 09:04:10,471 INFO Landweber :: |residual| = 2.5227734556307305e-05
2025-08-13 09:04:10,471 INFO CountIterations :: iteration = 27 / 100
2025-08-13 09:04:10,495 INFO Discrepancy :: relative discrepancy = 2.88, tolerance = 2.20
2025-08-13 09:04:10,795 INFO Landweber :: |residual| = 2.456992243794405e-05
2025-08-13 09:04:10,795 INFO CountIterations :: iteration = 28 / 100
2025-08-13 09:04:10,819 INFO Discrepancy :: relative discrepancy = 2.81, tolerance = 2.20
2025-08-13 09:04:11,127 INFO Landweber :: |residual| = 2.3981778332408313e-05
2025-08-13 09:04:11,127 INFO CountIterations :: iteration = 29 / 100
2025-08-13 09:04:11,151 INFO Discrepancy :: relative discrepancy = 2.75, tolerance = 2.20
2025-08-13 09:04:11,447 INFO Landweber :: |residual| = 2.345542035405972e-05
2025-08-13 09:04:11,447 INFO CountIterations :: iteration = 30 / 100
2025-08-13 09:04:11,471 INFO Discrepancy :: relative discrepancy = 2.70, tolerance = 2.20
2025-08-13 09:04:11,771 INFO Landweber :: |residual| = 2.2983851490596843e-05
2025-08-13 09:04:11,771 INFO CountIterations :: iteration = 31 / 100
2025-08-13 09:04:11,799 INFO Discrepancy :: relative discrepancy = 2.65, tolerance = 2.20
2025-08-13 09:04:12,103 INFO Landweber :: |residual| = 2.256087081354306e-05
2025-08-13 09:04:12,103 INFO CountIterations :: iteration = 32 / 100
2025-08-13 09:04:12,119 INFO Discrepancy :: relative discrepancy = 2.60, tolerance = 2.20
2025-08-13 09:04:12,419 INFO Landweber :: |residual| = 2.2180991157961516e-05
2025-08-13 09:04:12,419 INFO CountIterations :: iteration = 33 / 100
2025-08-13 09:04:12,443 INFO Discrepancy :: relative discrepancy = 2.56, tolerance = 2.20
2025-08-13 09:04:12,747 INFO Landweber :: |residual| = 2.1839363113275534e-05
2025-08-13 09:04:12,747 INFO CountIterations :: iteration = 34 / 100
2025-08-13 09:04:12,771 INFO Discrepancy :: relative discrepancy = 2.53, tolerance = 2.20
2025-08-13 09:04:13,079 INFO Landweber :: |residual| = 2.1531705208553877e-05
2025-08-13 09:04:13,079 INFO CountIterations :: iteration = 35 / 100
2025-08-13 09:04:13,099 INFO Discrepancy :: relative discrepancy = 2.49, tolerance = 2.20
2025-08-13 09:04:13,419 INFO Landweber :: |residual| = 2.1254240160466028e-05
2025-08-13 09:04:13,419 INFO CountIterations :: iteration = 36 / 100
2025-08-13 09:04:13,443 INFO Discrepancy :: relative discrepancy = 2.46, tolerance = 2.20
2025-08-13 09:04:13,747 INFO Landweber :: |residual| = 2.100363700650539e-05
2025-08-13 09:04:13,747 INFO CountIterations :: iteration = 37 / 100
2025-08-13 09:04:13,751 INFO Discrepancy :: relative discrepancy = 2.44, tolerance = 2.20
2025-08-13 09:04:14,059 INFO Landweber :: |residual| = 2.07769588894161e-05
2025-08-13 09:04:14,059 INFO CountIterations :: iteration = 38 / 100
2025-08-13 09:04:14,083 INFO Discrepancy :: relative discrepancy = 2.41, tolerance = 2.20
2025-08-13 09:04:14,391 INFO Landweber :: |residual| = 2.057161620406855e-05
2025-08-13 09:04:14,391 INFO CountIterations :: iteration = 39 / 100
2025-08-13 09:04:14,415 INFO Discrepancy :: relative discrepancy = 2.39, tolerance = 2.20
2025-08-13 09:04:14,711 INFO Landweber :: |residual| = 2.0385324772990896e-05
2025-08-13 09:04:14,723 INFO CountIterations :: iteration = 40 / 100
2025-08-13 09:04:14,747 INFO Discrepancy :: relative discrepancy = 2.37, tolerance = 2.20
2025-08-13 09:04:15,095 INFO Landweber :: |residual| = 2.0216068684970764e-05
2025-08-13 09:04:15,095 INFO CountIterations :: iteration = 41 / 100
2025-08-13 09:04:15,115 INFO Discrepancy :: relative discrepancy = 2.35, tolerance = 2.20
2025-08-13 09:04:15,407 INFO Landweber :: |residual| = 2.0062067413330274e-05
2025-08-13 09:04:15,407 INFO CountIterations :: iteration = 42 / 100
2025-08-13 09:04:15,431 INFO Discrepancy :: relative discrepancy = 2.34, tolerance = 2.20
2025-08-13 09:04:15,723 INFO Landweber :: |residual| = 1.9921746825604232e-05
2025-08-13 09:04:15,723 INFO CountIterations :: iteration = 43 / 100
2025-08-13 09:04:15,747 INFO Discrepancy :: relative discrepancy = 2.32, tolerance = 2.20
2025-08-13 09:04:16,047 INFO Landweber :: |residual| = 1.9793713702398604e-05
2025-08-13 09:04:16,047 INFO CountIterations :: iteration = 44 / 100
2025-08-13 09:04:16,071 INFO Discrepancy :: relative discrepancy = 2.31, tolerance = 2.20
2025-08-13 09:04:16,367 INFO Landweber :: |residual| = 1.9676733397786297e-05
2025-08-13 09:04:16,367 INFO CountIterations :: iteration = 45 / 100
2025-08-13 09:04:16,391 INFO Discrepancy :: relative discrepancy = 2.30, tolerance = 2.20
2025-08-13 09:04:16,679 INFO Landweber :: |residual| = 1.956971029433068e-05
2025-08-13 09:04:16,679 INFO CountIterations :: iteration = 46 / 100
2025-08-13 09:04:16,707 INFO Discrepancy :: relative discrepancy = 2.28, tolerance = 2.20
2025-08-13 09:04:16,991 INFO Landweber :: |residual| = 1.9471670730564507e-05
2025-08-13 09:04:16,992 INFO CountIterations :: iteration = 47 / 100
2025-08-13 09:04:17,015 INFO Discrepancy :: relative discrepancy = 2.27, tolerance = 2.20
2025-08-13 09:04:17,303 INFO Landweber :: |residual| = 1.9381748105689125e-05
2025-08-13 09:04:17,303 INFO CountIterations :: iteration = 48 / 100
2025-08-13 09:04:17,327 INFO Discrepancy :: relative discrepancy = 2.26, tolerance = 2.20
2025-08-13 09:04:17,623 INFO Landweber :: |residual| = 1.9299169893954363e-05
2025-08-13 09:04:17,623 INFO CountIterations :: iteration = 49 / 100
2025-08-13 09:04:17,647 INFO Discrepancy :: relative discrepancy = 2.25, tolerance = 2.20
2025-08-13 09:04:17,947 INFO Landweber :: |residual| = 1.9223246328555225e-05
2025-08-13 09:04:17,947 INFO CountIterations :: iteration = 50 / 100
2025-08-13 09:04:17,971 INFO Discrepancy :: relative discrepancy = 2.25, tolerance = 2.20
2025-08-13 09:04:18,263 INFO Landweber :: |residual| = 1.915336054116416e-05
2025-08-13 09:04:18,263 INFO CountIterations :: iteration = 51 / 100
2025-08-13 09:04:18,287 INFO Discrepancy :: relative discrepancy = 2.24, tolerance = 2.20
2025-08-13 09:04:18,591 INFO Landweber :: |residual| = 1.908895996789645e-05
2025-08-13 09:04:18,591 INFO CountIterations :: iteration = 52 / 100
2025-08-13 09:04:18,615 INFO Discrepancy :: relative discrepancy = 2.23, tolerance = 2.20
2025-08-13 09:04:18,919 INFO Landweber :: |residual| = 1.902954885528042e-05
2025-08-13 09:04:18,919 INFO CountIterations :: iteration = 53 / 100
2025-08-13 09:04:18,943 INFO Discrepancy :: relative discrepancy = 2.23, tolerance = 2.20
2025-08-13 09:04:19,255 INFO Landweber :: |residual| = 1.8974681720525396e-05
2025-08-13 09:04:19,255 INFO CountIterations :: iteration = 54 / 100
2025-08-13 09:04:19,279 INFO Discrepancy :: relative discrepancy = 2.22, tolerance = 2.20
2025-08-13 09:04:19,571 INFO Landweber :: |residual| = 1.892395763902337e-05
2025-08-13 09:04:19,571 INFO CountIterations :: iteration = 55 / 100
2025-08-13 09:04:19,595 INFO Discrepancy :: relative discrepancy = 2.21, tolerance = 2.20
2025-08-13 09:04:19,891 INFO Landweber :: |residual| = 1.8877015248633938e-05
2025-08-13 09:04:19,891 INFO CountIterations :: iteration = 56 / 100
2025-08-13 09:04:19,915 INFO Discrepancy :: relative discrepancy = 2.21, tolerance = 2.20
2025-08-13 09:04:20,219 INFO Landweber :: |residual| = 1.883352837499616e-05
2025-08-13 09:04:20,219 INFO CountIterations :: iteration = 57 / 100
2025-08-13 09:04:20,243 INFO Discrepancy :: relative discrepancy = 2.20, tolerance = 2.20
2025-08-13 09:04:20,547 INFO Landweber :: |residual| = 1.8793202195020483e-05
2025-08-13 09:04:20,547 INFO CountIterations :: iteration = 58 / 100
2025-08-13 09:04:20,571 INFO Discrepancy :: relative discrepancy = 2.20, tolerance = 2.20
2025-08-13 09:04:20,571 INFO CombineRules :: Rule Discrepancy(noiselevel=8.525711012111121e-06, tau=2.2) triggered.
2025-08-13 09:04:20,572 INFO Landweber :: Solver converged after 57 iteration.
Compute Error¶
Now we can calculate the relative error on \(L^2(\Gamma_{Top})\).
[6]:
# 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: 13.266079090567173 %
Plot Reconstruction¶
We plot the true traction forces on the reconstruction mesh, the reconstruction and the error.
[7]:
Draw(traction_true_gf, domain_rec.fes.mesh)
Draw(reco_gf)
error = traction_true_gf - domain_rec.to_ngs(reco)
[8]:
Draw(error, domain_rec.fes.mesh)
[8]:
BaseWebGuiScene
[ ]: