Precision in hyperelastic model

[marcomusy]

Hi! I’m trying to model the tissue growth of an embryo using the hyperelastic example. I’m interested in the internal displacements.
I show here an example script which is somehow similar to the real thing: an initial circle is stretched by means of a force into its final shape. The force is variable with time (first points downthen up). This is done in 40 separate time steps:

Now imagine you want to track a “cell”. I do it by applying the computed sequence of displacements to a whatever point on the plane: these are the colored lines. The end point of these are the final position of the “cell”. Now I invert the displacement vectors, and in principle they should be walking back to where they started: these are the black lines.

While overall it seems to work well, there is a small discrepancy between the 2 lines (colored and black). I wonder if you have any idea why is that, maybe it’s something to be expected? Or can be further improved? Or am I doing something wrong?

I also tested with doing a remeshing at each step - instead of keeping the same mesh - and it doesn’t seem to change much (which is very good!):

Thanks a lot for your help!

The code:

from dolfin import *
from mshr import *
import numpy as np
import vedo

set_log_level(30)


class AllBoundaries(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0]<-10

def solve_problem(mesh, mfunc, force):
    V = VectorFunctionSpace(mesh, "CG", 1)
    u = TrialFunction(V)
    v = TestFunction(V)
    displacement = Function(V)

    bc = [DirichletBC(V, Constant((0,0)), mfunc, 1)]

    F = Constant(force)
    E = Constant(5000)
    nu = Constant(0.3)
    mu = E / (2.0*(1+nu))
    lmbda = E * nu / ((1.0+nu) * (1-2*nu))
    sigma = 2.0 * mu * sym(grad(u)) + lmbda * tr( sym(grad(u)) ) * Identity(2)
    solve(inner(sigma, grad(v)) * dx == inner(F, v) * dx, displacement, bc)
    displacement.set_allow_extrapolation(True)
    return displacement

def update(mesh, displacement):
    new_mesh = Mesh(mesh)
    ALE.move(new_mesh, displacement)
    return new_mesh

def remesh(mesh, res=50):
    if isinstance(mesh, vedo.Mesh):
        vmesh = mesh
    else:
        vmesh = vedo.Mesh([mesh.coordinates(), mesh.cells()])
    bpts = vmesh.computeNormals(cells=True).boundaries().join(reset=1) #extract boundary
    vz = vmesh.celldata["Normals"][0][2] # check z component of normals at first point
    bpts.tomesh(invert=vz<0).smooth().write('tmpmesh.xml') #vedo REMESHING + smoothing
    return Mesh("tmpmesh.xml")

#################################################################################
N = 40             # number of iterations of stretching
res = 15           # resolution of meshes
do_remesh = False  # grab the boundary and remesh the interior at each iteration
vedo.settings.useParallelProjection = True  # avoid perspective parallax

circle = Circle(Point(0, 0), 50)
mesh = generate_mesh(circle, res)

meshes = [mesh]
displacements = []
for i in range(N):
    mfunc = MeshFunction('size_t', mesh, 1, mesh.domains())
    mfunc.set_all(0)
    allb = AllBoundaries()
    allb.mark(mfunc, 1)

    F = np.array([2, (i-N/2)/N])

    displacement = solve_problem(mesh, mfunc, F)
    new_mesh = update(mesh, displacement)

    if do_remesh:
        mesh = remesh(new_mesh)
    else:
        mesh = new_mesh
    meshes.append(mesh)
    displacements.append(displacement)
    # plot things:
    txt = vedo.Text2D(f"step{i}")
    arrow = vedo.Arrow2D([0,0], F*20).z(1)
    vedo.dolfin.plot(mesh, arrow, txt, c='grey5', at=i, N=N, zoom=1.1) # PRESS q to continue

dmesh_i = meshes[0]  # initial mesh
dmesh_f = meshes[-1] # final mesh

vmesh_i = vedo.Mesh([dmesh_i.coordinates(), dmesh_i.cells()], c='grey5').z(-1)
vmesh_f = vedo.Mesh([dmesh_f.coordinates(), dmesh_f.cells()], c='grey3').wireframe()

plt = vedo.Plotter()

# move a few points along the deformation of the circle
seeds = vedo.Circle(r=50, res=res).points()
endpoints = []
for i, p in enumerate(seeds):
    line = [p]
    for disp in displacements:
        dp   = disp([p[0], p[1]])
        newp = p  + [dp[0],dp[1], 0.0]
        line.append(newp)
        p = newp
    plt += vedo.Line(line, c=i, lw=4).z(1)
    endpoints.append(newp)

# invert everything and move the end points back in place
seeds = list(reversed(endpoints))
displacements = list(reversed(displacements))
for i, p in enumerate(seeds):
    line = [p]
    for disp in displacements:
        dp   = disp([p[0], p[1]])
        newp = p  - [dp[0],dp[1], 0.0]
        line.append(newp)
        p = newp
    plt += vedo.Line(line, c='k', lw=2).z(2)

plt += [vmesh_i, vmesh_f]
plt.show(axes=1)

[kamensky]

I think the basic problem here is an off-by-one error regarding which point the n^\text{th} displacement is evaluated at. Consider a displacement function

\mathbf{u}^n : \mathbb{R}^2\to\mathbb{R}^2\text{ ,}

and, for some given point \mathbf{x}^n\in\mathbb{R}^2, define

\mathbf{x}^{n+1} := \mathbf{x}^n + \mathbf{u}^n(\mathbf{x}^n)\text{ .}

Then

\mathbf{x}^n \neq \mathbf{x}^{n+1} - \mathbf{u}^n(\mathbf{x}^{n+1})\text{ ,}

because

\mathbf{u}^n(\mathbf{x}^{n+1}) \neq \mathbf{u}^n(\mathbf{x}^n)\text{ .}

It looks like these invalid equations are being used in the attempt to trace the trajectories backward, which would explain the slight deviation.

[marcomusy]

Thanks for the answer, that indeed explains it.
So is there some meaningful way of “inverting” a set of solutions and back-trace the points or, in other words, can I solve for a in:

\mathbf{b} = \mathbf{a} + \mathbf{u}^n(\mathbf{a})\text{ .}

for a known generic b ?

[kamensky]

Assuming the deforming body never self-intersects, a unique inverse deformation exists. But, actually solving for \mathbf{a} in practice may be hard, especially since \mathbf{u}^n is not smooth. You could formally set up a Newton iteration using residual

\mathbf{R}(\mathbf{a}) = \mathbf{a} + \mathbf{u}^n(\mathbf{a}) - \mathbf{b}

and Jacobian

\mathbf{J}(\mathbf{a}) = \mathbf{I} + \nabla\mathbf{u}^n(\mathbf{a})\text{ ,}

but the fact that \nabla\mathbf{u}^n is discontinuous would make convergence unreliable. It might still work, though, if the initial guess is close to \mathbf{a} (e.g., \mathbf{b} in this case), jumps in \nabla\mathbf{u}^n are small, and some sort of line search or damping is used.