{ "cells": [ { "cell_type": "markdown", "id": "c5aff1a82352ab46", "metadata": { "collapsed": false }, "source": [ "# Traction force microscopy - an ngsolve example with Landweber iteration\n", "\n", "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](https://ngsolve.org/downloads).\n", "\n", "## Mathematical Model\n", "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.\n", "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\n", "\n", "$$ div (\\sigma(u)) = 0 \\quad \\text{in } \\Omega $$\n", "$$ \\sigma(u) n = t \\quad \\text{on } \\Gamma_{Top} $$\n", "$$ \\sigma(u) n = 0 \\quad \\text{on } \\Gamma_{Sides} $$\n", "$$ u = 0 \\quad \\text{on } \\Gamma_{Bottom} $$\n", "\n", "\n", "where $\\sigma$ is given by a material law, in this case the linear Hooke's law\n", "\n", "$$ \\sigma(u) = 2 \\mu \\varepsilon + \\lambda tr( \\varepsilon) I $$\n", "\n", "for the Lamé parameters $\\mu$, $\\lambda$ and $\\varepsilon$ is the linearized strain tensor\n", "\n", "$$ \\varepsilon = \\frac{1}{2} (\\nabla u + (\\nabla u)^T ) $$\n", "\n", "\n", "In the weak formulation we get \n", "$$ a(u,v) = l_t(v) , \\quad u,v \\in H^1_{0,\\Gamma_D}(\\Omega,\\mathbb{R}^3)$$\n", "\n", "for the bilinear form\n", "\n", "$$ a(u,v) = \\int_\\Omega 2 \\mu \\varepsilon(u) : \\varepsilon(v) + \\lambda div(u) div(v) dx $$\n", "\n", "and the linear form\n", "\n", "$$ l_t(v) = \\int_{\\Gamma_T} t v ds $$\n", "\n", "## Implementation\n", "\n", "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.\n", "\n", "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).\n", "\n", "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. )\n", "\n", "---\n", "\n", "### Parameters\n", "`domain` : `regpy.vecsps.ngsolve.NgsSpace`\n", " > 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.\n", " \n", "`codomain` : `regpy.vecsps.ngsolve.NgsSpace`\n", " > The codomain on which the operator is defined. Should be as domain, but without `bdr`defined.\n", "\n", "`mu`,`lambda`: `float`\n", "> Lamé parameters used in Hooke's law\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "e25c5301-cad2-4eab-a867-0608b2a40bfd", "metadata": { "ExecuteTime": { "end_time": "2024-04-18T06:29:57.557681978Z", "start_time": "2024-04-18T06:29:56.739288182Z" } }, "outputs": [], "source": [ "import ngsolve as ngs\n", "from regpy.operators.ngsolve import NGSolveOperator\n", "\n", "def strain(u):\n", " return 0.5 * (ngs.Grad(u) + ngs.Grad(u).trans)\n", "\n", "class TFM(NGSolveOperator):\n", "\n", " def __init__(self, domain, codomain, mu, lam):\n", " codomain = codomain\n", " # Need to know the boundary to calculate Neumann bdr condition\n", " assert domain.bdr is not None\n", " super().__init__(domain, codomain, linear=True)\n", " # From NgSolveOperator\n", " #self.gfu_read_in = ngs.GridFunction(self.domain.fes)\n", "\n", " # Lamé Parameters for substrate\n", " self.mu = mu\n", " self.lam = lam\n", "\n", " self.fes_domain = domain.fes\n", " self.fes_codomain = codomain.fes\n", "\n", "\n", " # grid functions for later use\n", " self.gfu_eval = ngs.GridFunction(self.fes_codomain) # solution, return value of _eval\n", " self.gfu_adjoint = ngs.GridFunction(self.fes_domain) # grid function return value of adjoint (trace of gfu_adjoint_sol)\n", "\n", " self.gfu_bf = ngs.GridFunction(self.fes_codomain) # grid function for defining integrator (bilinearform)\n", " self.gfu_lf = ngs.GridFunction(self.fes_codomain) # grid function for defining right hand side (linearform), f\n", "\n", "\n", "\n", " u, v = self.fes_codomain.TnT()\n", "\n", " # Define Bilinearform, will be assembled later\n", " self.a = ngs.BilinearForm(self.fes_codomain, symmetric=True)\n", " self.a += (2 * mu * ngs.InnerProduct(strain(u), strain(v)) + lam * ngs.div(u) * ngs.div(v)) * ngs.dx\n", "\n", "\n", " # Define Linearform for evaluation, will be assembled later\n", " self.b = ngs.LinearForm(self.fes_domain)\n", " self.b += self.gfu_lf * v * ngs.ds(domain.bdr)\n", "\n", "\n", " # Define linearform to trick ngsolve for computation of discrete adjoint\n", " self.b_help = ngs.LinearForm(self.fes_domain)\n", " self.b_help.Assemble()\n", "\n", "\n", " # Initialize preconditioner for solving the Dirichlet problems by ngs.BVP\n", " self.prec = ngs.Preconditioner(self.a, 'direct')\n", " self.a.Assemble()\n", "\n", "\n", " # Left term: Bilinearform self.a\n", " # Right term: Linearform self.b\n", " def _eval(self, traction, differentiate=False):\n", "\n", " # Assemble Linearform, boundary term\n", " self._read_in(traction, self.gfu_lf)\n", " self.b.Assemble()\n", "\n", " self._solve_dirichlet_problem(bf=self.a, lf=self.b, gf=self.gfu_eval, prec=self.prec)\n", "\n", " return self.gfu_eval.vec.FV().NumPy()[:].copy()\n", "\n", "\n", " def _adjoint(self, displacement):\n", " # Bilinearform already assembled in init -> initialization with 0, s.t. object exists\n", " # Diskrete Adjoint w.r.t. standard inner product\n", "\n", " self.b_help.vec.FV().NumPy()[:] = displacement.copy()\n", " self._solve_dirichlet_problem(bf=self.a, lf=self.b_help, gf=self.gfu_adjoint, prec=self.prec)\n", " self._read_in(self.gfu_adjoint.vec.FV().NumPy()[:].copy(), self.gfu_lf)\n", " self.b.Assemble()\n", "\n", " return self.b.vec.FV().NumPy()[:].copy()" ] }, { "cell_type": "markdown", "id": "fe66e281-383e-41a7-a6f2-c9ba38eff4d3", "metadata": {}, "source": [ "## Create Operators\n", "\n", "First we import logging to later get printed messages during our solution algorithm.\n", "\n", "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.\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "b44ff62f-6ef5-483a-b99e-04d574f397f4", "metadata": {}, "outputs": [], "source": [ "import logging\n", "from ngsolve.webgui import Draw \n", "from netgen.occ import *\n", "from regpy.vecsps.ngsolve import NgsSpace\n", "\n", "logging.basicConfig(\n", " level=logging.INFO,\n", " format='%(asctime)s %(levelname)s %(name)-20s :: %(message)s'\n", ")\n", "\n", "'''############################### material parameters ###############################'''\n", "\n", "# E Young's modulus, nu Poisson ratio\n", "E, nu = 10000, 0.45\n", "# Lamé parameters\n", "mu = E / 2 / (1 + nu)\n", "lam = E * nu / ((1 + nu) * (1 - 2 * nu))\n", "\n", "'''############################### mesh of substrate ################################'''\n", "\n", "# generate two triangular meshes of mesh-size h, order k, for generating and reconstructing traction forces\n", "h_gen = 0.3\n", "h_rec = 0.4\n", "k = 3\n", "\n", "\n", "box = Box(Pnt(-2, -0.3,-2), Pnt(2, 0.3,2))\n", "box.faces.Max(Y).name = 'top'\n", "box.faces.Min(Y).name = 'bottom'\n", "geo = OCCGeometry(box)\n", "\n", "mesh_gen = ngs.Mesh(geo.GenerateMesh(maxh=h_gen))\n", "mesh_gen.Curve(3)\n", "\n", "Draw(mesh_gen)\n", "\n", "mesh_rec = ngs.Mesh(geo.GenerateMesh(maxh=h_rec))\n", "mesh_rec.Curve(3)\n", "\n", "'''############################### define operator ################################'''\n", "# operator for generating data on a different discretization than reconstruction operator\n", "\n", "fes_domain_gen = ngs.VectorH1(mesh_gen, order=k, dirichlet=\"bottom\")\n", "# fes_domain = ngs.H1(mesh, order=k, dim=3)\n", "domain_gen = NgsSpace(fes_domain_gen, bdr='top')\n", "\n", "fes_codomain_gen = ngs.VectorH1(mesh_gen, order=k, dirichlet=\"bottom\")\n", "# fes_codomain = ngs.H1(mesh, order=k, dirichlet=\"bottom|side\", dim=3)\n", "codomain_gen = NgsSpace(fes_codomain_gen)\n", "\n", "op_gen = TFM(domain_gen, codomain=codomain_gen, mu=mu, lam=lam)\n", "\n", "\n", "fes_domain_rec = ngs.VectorH1(mesh_rec, order=k, dirichlet=\"bottom\")\n", "# fes_domain = ngs.H1(mesh, order=k, dim=3)\n", "domain_rec = NgsSpace(fes_domain_rec, bdr='top')\n", "\n", "fes_codomain_rec = ngs.VectorH1(mesh_rec, order=k, dirichlet=\"bottom\")\n", "# fes_codomain = ngs.H1(mesh, order=k, dirichlet=\"bottom|side\", dim=3)\n", "codomain_rec = NgsSpace(fes_codomain_rec)\n", "\n", "op_rec = TFM(domain_rec, codomain=codomain_rec, mu=mu, lam=lam)" ] }, { "cell_type": "markdown", "id": "58f65bb2-92ec-4df1-97ad-bc3338d1a4b7", "metadata": {}, "source": [ "## Generate Groundtruth\n", "\n", "In this step we generate a toy surface traction force first as a `CoefficientFunction`from 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.\n", "\n", "To change between ngsolve and regpy, the functions `to_ngs` and `from_ngs` can be used." ] }, { "cell_type": "code", "execution_count": null, "id": "fa8378b2-2d0e-484d-ae9c-c428fb2e4b9c", "metadata": { "ExecuteTime": { "end_time": "2024-04-17T07:43:33.024011907Z", "start_time": "2024-04-17T07:43:32.786030558Z" } }, "outputs": [], "source": [ "a = 4\n", "c = 1\n", "\n", "force_intensity = (-a * (ngs.sqrt(ngs.x**2 + ngs.z**2) - 0.5)**2 + c )\n", "\n", "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))))\n", "\n", "cell_xz = -(ngs.sqrt(ngs.x ** 2 + ngs.z ** 2) - 1)\n", "\n", "\n", "traction_true_cf = ngs.IfPos(c1=cell_xz, then_obj=cell_force, else_obj=ngs.CF((0, 0, 0)))\n", "\n", "traction_true_gf = ngs.GridFunction(fes_domain_gen)\n", "traction_true_gf.Set(traction_true_cf, definedon = mesh_gen.Boundaries( \"top\" ))\n", "\n", "traction_true = domain_gen.from_ngs(traction_true_gf)\n", "\n", "Draw(traction_true_gf, mesh_gen)" ] }, { "cell_type": "markdown", "id": "7656a146-8fb8-4cd4-ac2d-860f85c8941a", "metadata": {}, "source": [ "## Generate Data\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "d4c65800-9c35-4277-8341-6765367e558b", "metadata": { "ExecuteTime": { "end_time": "2024-04-17T07:43:34.055490909Z", "start_time": "2024-04-17T07:43:33.215227413Z" } }, "outputs": [], "source": [ "import numpy as np\n", "\n", "displacement_true = op_gen(traction_true)\n", "\n", "\n", "Draw(codomain_gen.to_ngs(displacement_true))\n", "\n", "\n", "\n", "np.random.seed(42)\n", "\n", "\n", "noise = 1e-6 * codomain_rec.randn()\n", "data_gf = ngs.GridFunction(fes_codomain_rec)\n", "data_gf.Set(codomain_gen.to_ngs(displacement_true))\n", "data = codomain_rec.from_ngs(data_gf) + noise\n", "\n", "Draw(codomain_rec.to_ngs(data), mesh_rec)\n" ] }, { "cell_type": "markdown", "id": "4f852088-0c80-4430-a2e9-d390b3119280", "metadata": {}, "source": [ "## Solve the Inverse Problem with regularization\n", "\n", "Now, with available data we can use, e.g., landweber iteration to reconstruct the groundtruth traction force from the noisy displacement.\n", "\n", "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 `AbstractSpace`s and the according implementation can be found in `regpy.hilbert.ngsolve`. \n", "\n", "As a next step we define the we choose an initial guess, in this case we choose zero. \n", "\n", "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`.\n", "\n", "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`. \n", "\n", "After everything is defined run the solver with the specified stopping rule using the method `run()` of the solver. " ] }, { "cell_type": "code", "execution_count": null, "id": "347b3904-6ccb-4612-bfb0-624493e8ea1f", "metadata": { "ExecuteTime": { "start_time": "2024-04-17T07:43:34.054942815Z" }, "is_executing": true }, "outputs": [], "source": [ "from regpy.solvers import RegularizationSetting\n", "from regpy.solvers.nonlinear.landweber import Landweber\n", "import regpy.stoprules as rules\n", "from regpy.hilbert import L2Boundary, Hm0\n", "\n", "setting = RegularizationSetting(op=op_rec, penalty=L2Boundary, data_fid=Hm0)\n", "init = domain_rec.from_ngs((0, 0, 0))\n", "\n", "landweber = Landweber(setting, data, init)\n", "\n", "stoprule = (\n", " rules.CountIterations(100) +\n", " rules.Discrepancy(setting.h_codomain.norm, data, noiselevel=setting.h_codomain.norm(noise), tau=2.2)\n", ")\n", "\n", "reco, reco_data = landweber.run(stoprule)\n", "\n", "reco_gf = domain_rec.to_ngs(reco)\n" ] }, { "cell_type": "markdown", "id": "1ad4e521-3639-432a-b184-5607186a4c65", "metadata": {}, "source": [ "## Compute Error\n", "Now we can calculate the relative error on $L^2(\\Gamma_{Top})$." ] }, { "cell_type": "code", "execution_count": null, "id": "cf7fb9a1-829e-4868-bb6d-6f0d6472ab45", "metadata": {}, "outputs": [], "source": [ "# calculate relative L2 error on boundary\n", "err_abs = ngs.sqrt(ngs.Integrate((reco_gf - traction_true_cf)**2, mesh_rec,definedon=mesh_rec.Boundaries('top')))\n", "norm_true = ngs.sqrt(ngs.Integrate((traction_true_cf)**2, mesh_rec,definedon=mesh_rec.Boundaries('top'))) \n", "norm_rec = ngs.sqrt(ngs.Integrate((reco_gf)**2, mesh_rec,definedon=mesh_rec.Boundaries('top')))\n", "err_rel = err_abs/norm_true * 100\n", "\n", "print('relative error:' , err_rel, '%')" ] }, { "cell_type": "markdown", "id": "533cbf0f-0a59-441c-b67e-a1157722641c", "metadata": {}, "source": [ "## Plot Reconstruction\n", "\n", "We plot the true traction forces on the reconstruction mesh, the reconstruction and the error." ] }, { "cell_type": "code", "execution_count": null, "id": "8b0157b5-8118-4bc1-b766-fa359014e0a2", "metadata": { "is_executing": true }, "outputs": [], "source": [ "Draw(traction_true_gf, domain_rec.fes.mesh)\n", "\n", "Draw(reco_gf)\n", "\n", "error = traction_true_gf - domain_rec.to_ngs(reco)\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "8c4d6043-a3e3-457f-a565-19c6ed7ab3f6", "metadata": {}, "outputs": [], "source": [ "Draw(error, domain_rec.fes.mesh)" ] }, { "cell_type": "code", "execution_count": null, "id": "caee4912-8889-474e-91a3-d4190bfc84d9", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.12" } }, "nbformat": 4, "nbformat_minor": 5 }