Solving for equilibrium in nonlinear elasticity - application to a surface matching problem


165
views
0
5 months ago by
Sorry, this is a long post, so bare with me ...
I'm trying to implement a surface matching algorithm driven by nonlinear elasticity PDEs (conservation of linear momentum). The paper is here: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3174067/
We are given two surfaces, YAS and OAS that are in vertex-wise correspondence, and have been parametrized to the sphere (OSS). The original correspondences are not optimal as they might introduce anatomically unlikely distortions. The aim of the algorithm is to find a vector v on OSS that assigns new positions to the vertices corresponding to OAS. These positions are called "Relaxed", and the resulting surface ROSS. If we interpolate the original vertices of OAS, that correspond to OSS, at locations ROSS, we obtain ROAS (Relaxed Older Anatomical Surface), and the deformation gradient between YAS and ROAS is better (less strain) than that of the original mapping.
The whole thing can be expressed as finding v in

\[
\nu \frac{\partial v}{\partial t} = \nabla \cdot \{ \mathbf{H}^{'\top} \cdot \mathbf{P} \cdot \mathbf{F}_0 \cdot \mathbf{H}^{-\top} \}  + \mathbf{H}^{-1} \cdot \mathbf{f}
\\
W = \frac{1}{2} ( \mu (I_1 I_3^{-\frac{1}{3}} + \kappa(I_3^{\frac{1}{2}} - 1)^2)
\\
\mathbf{P} = \frac{\partial W}{\partial \mathbf{F}} = \mu ( \mathbf{F} - \mu \mathbf{F}^{-\top}) + \lambda \log ( \det (\mathbf{F}))\mathbf{F}^{-\top}
\\
\mathbf{F}=\mathbf{I} + ( \nabla \mathbf{v} )^\top
\]

where \( \mathbf{F}_0 \) is the original deformation gradient between YAS and OAS, \( \mathbf{H} \) is the mapping to the sphere and \( \mathbf{H}^{'\top} \) is the mapping back at \( \mathbf{F} \). There can also be a body force \( \mathbf{f} \), which contains some information on the surface.

I implemented this in fenics, and it seems to work quite nicely. Here's a MWE of a toy example where all surfaces (YAS, OAS, OSS) are actually spheres, but OAS is distorted by a mathematical formula - the solution should reverse that distortion.

from dolfin import *
from mshr import *
import numpy as np
import sys
from ufl import grad as ufl_grad
from scipy.spatial import ConvexHull
from scipy.interpolate import LinearNDInterpolator
from numpy.core.umath_tests import inner1d
from euclid import *
from math import *


# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["quadrature_degree"] = 5
ffc_options = {"optimize": True}

# step size
dt=Constant(1)
theta=0.5 # (Crank-Nicholson)

# helper functions for deformation gradient
def Grad(v):
    return ufl_grad(v)

# Second order identity tensor
def SecondOrderIdentity(u):
    return variable(Identity(u.geometric_dimension()))

# Deformation gradient
def DeformationGradient(u):
    I = SecondOrderIdentity(u)
    return variable(I + Grad(u))

# helper functions for spherical geometry
def cart2sph(x,y,z):
    azimuth = np.arctan2(y,x)
    elevation = np.arctan2(z,np.sqrt(x**2 + y**2))
    r = np.sqrt(x**2 + y**2 + z**2)
    return azimuth, elevation, r

def sph2cart(azimuth,elevation,r):
    x = r * np.cos(elevation) * np.cos(azimuth)
    y = r * np.cos(elevation) * np.sin(azimuth)
    z = r * np.sin(elevation)
    return x, y, z

# delaunay interpolation on the sphere
def SphericalInterpolator(v,valuesOnSphere):
    f=ConvexHull(v).simplices
    
    # compute triangle areas to determine fudge factor
    def triAreas(v,f):
        a=np.sqrt(np.sum( (v[f[:,0],:] - v[f[:,1],:])**2,1));
        b=np.sqrt(np.sum( (v[f[:,1],:] - v[f[:,2],:])**2,1));
        c=np.sqrt(np.sum( (v[f[:,0],:] - v[f[:,2],:])**2,1));
        s = (a+b+c)/2;
        A = np.real(np.sqrt(s*(s-a)*(s-b)*(s-c)));
        return A

    As=triAreas(v,f)
    bigInd=np.argmax(As)

    biggest=v[f[bigInd,:],:];

    # compute max. distance between this triangle and the sphere
    ds=1-np.sin(np.arccos(inner1d((biggest[(1,2,1),:]-biggest[(0,0,2),:]),biggest[(1,2,1),:])/np.sqrt(np.sum((biggest[(0,0,2),:]-biggest[(1,2,1),:])**2,1))));
    fact=2*max(ds);

    F=LinearNDInterpolator(np.vstack((v*(1+fact),v*(1-fact))), 
                           np.vstack((valuesOnSphere,valuesOnSphere)), 
                           fill_value=np.nan, rescale=False)

    return F

# create a subdivision sphere
def SubdivisionSphereMesh(num_subdivisions):

    def icosahedron():
        """Construct a 20-sided polyhedron"""
        faces = [ \
            (0,1,2),
            (0,2,3),
            (0,3,4),
            (0,4,5),
            (0,5,1),
            (11,6,7),
            (11,7,8),
            (11,8,9),
            (11,9,10),
            (11,10,6),
            (1,2,6),
            (2,3,7),
            (3,4,8),
            (4,5,9),
            (5,1,10),
            (6,7,2),
            (7,8,3),
            (8,9,4),
            (9,10,5),
            (10,6,1) ]
        verts = [ \
            ( 0.000,  0.000,  1.000 ),
            ( 0.894,  0.000,  0.447 ),
            ( 0.276,  0.851,  0.447 ),
            (-0.724,  0.526,  0.447 ),
            (-0.724, -0.526,  0.447 ),
            ( 0.276, -0.851,  0.447 ),
            ( 0.724,  0.526, -0.447 ),
            (-0.276,  0.851, -0.447 ),
            (-0.894,  0.000, -0.447 ),
            (-0.276, -0.851, -0.447 ),
            ( 0.724, -0.526, -0.447 ),
            ( 0.000,  0.000, -1.000 ) ]
        return verts, faces

    def subdivide(verts, faces):
        """Subdivide each triangle into four triangles, pushing verts to the unit sphere"""
        triangles = len(faces)
        for faceIndex in xrange(triangles):

            # Create three new verts at the midpoints of each edge:
            face = faces[faceIndex]
            a,b,c = (Vector3(*verts[vertIndex]) for vertIndex in face)
            verts.append((a + b).normalized()[:])
            verts.append((b + c).normalized()[:])
            verts.append((a + c).normalized()[:])

            # Split the current triangle into four smaller triangles:
            i = len(verts) - 3
            j, k = i+1, i+2
            faces.append((i, j, k))
            faces.append((face[0], i, k))
            faces.append((i, face[1], j))
            faces[faceIndex] = (k, j, face[2])

        return verts, faces

    verts, faces = icosahedron()
    for x in xrange(num_subdivisions):
        verts, faces = subdivide(verts, faces)

    mesh = Mesh()
    editor = MeshEditor()
    editor.open(mesh,'triangle',2, 3)
    editor.init_vertices(np.shape(verts)[0])
    editor.init_cells(np.shape(faces)[0])
    [editor.add_vertex(i,n) for i,n in enumerate(np.array(verts))]
    [editor.add_cell(i,n) for i,n in enumerate(np.array(faces,dtype=np.uintp))]
    editor.close()

    return mesh




# ***** DATA. A 2D SPHERICAL SHELL EMBEDDING IN 3D *****
sphere_subdivisions=3
YAS = SubdivisionSphereMesh(sphere_subdivisions)
OAS = SubdivisionSphereMesh(sphere_subdivisions)
OSS = SubdivisionSphereMesh(sphere_subdivisions)
ROSS = SubdivisionSphereMesh(sphere_subdivisions)

# apply a deformation
(TH,PHI,R)=cart2sph(YAS.coordinates()[:,0],YAS.coordinates()[:,1],YAS.coordinates()[:,2]);
th = TH + (.1/np.pi) * (np.sin(4*TH) * np.cos(PHI)**2);
(OAS.coordinates()[:,0],OAS.coordinates()[:,1],OAS.coordinates()[:,2])=sph2cart(th,PHI,R);



# ******** DOMAINS ********
T1 = TensorFunctionSpace(OSS, "Lagrange", 1);  ttmp = Function(T1)
V1 = VectorFunctionSpace(OSS, "Lagrange", 1);  vtmp = Function(V1)
S1 = FunctionSpace(OSS,"Lagrange", 1);         stmp = Function(S1);  

T = TensorFunctionSpace(OSS, "Lagrange", 2); 
V = VectorFunctionSpace(OSS, "Lagrange", 2); 
S = FunctionSpace(OSS,"Lagrange", 2);

n = Expression (("x[0]", "x[1]", "x[2]"),element=V.ufl_element())
n1= Expression (("x[0]", "x[1]", "x[2]"),element=V1.ufl_element())
OSS.init_cell_orientations ( n1 )

## ******** DEFORMATION GRADIENTS BETWEEN THE TWO SURFACES YAS AND OAS, YAS/OAS AND OSS *********

# deformation gradient of transformation between YAS and OAS
vtmp.vector()[vertex_to_dof_map(V1)]=np.reshape(OAS.coordinates()-YAS.coordinates(),(-1,1))
u = interpolate(vtmp,V)
F0=DeformationGradient(u).T

# deformation gradient of transformation between YAS/OAS and OSS
vtmp.vector()[vertex_to_dof_map(V1)]=np.reshape(OAS.coordinates()-OSS.coordinates(),(-1,1))
h = interpolate(vtmp,V)
H = DeformationGradient(h)

# deformation gradient of transformation between OSS and YAS/OAS
Hp_0 = Function(T)
Hp_0.vector()[:]=project(H,T).vector().get_local()[:]

# deformation gradient of transformation between ROSS and YAS/OAS (changing during solve)
Hp = Function(T)
Hp.vector()[:]=project(H,T).vector().get_local()[:]


# ****** body force - (0, can be a function to be matched on YAS and OAS) ******
f1_=np.zeros((np.shape(OSS.coordinates())[0],1))
f2_=np.zeros((np.shape(OSS.coordinates())[0],1))

f1=Function(S1); f1.vector()[vertex_to_dof_map(S1)] = f1_
f2=Function(S1); f2.vector()[vertex_to_dof_map(S1)] = f2_

# f2 at last time-stepping iteration
f2k=Function(S1);
f2k.vector()[:] = f2.vector()[:]

# f2 on ROSS
f2_0=interpolate(f2,S1)

# current body force
B = (f1-f2)*grad(f2)

# body force at last time-stepping iteration
Bk = (f1-f2k)*grad(f2k)

## ************ CONSTITUTIVE EQUATIONS ************ 
# material properties
mu=Constant(0.01)
kappa=Constant(10)
nu=Constant(100)

def Q(v):
    # Deformation gradient
    I = Identity(v.geometric_dimension())
    F1 = I + (I-outer(n1,n1))*grad(v)

    # Deformation gradient times the spherical mappings
    iH=inv(H)
    F = Hp.T*F1*F0.T*iH.T

    # Invariants of deformation tensors
    C = F.T*F
    I1 = tr(C)
    I3 = det(C)
    J = I3**(1/2)
    I1star=I1*(I3**(-1/3))
    
    # Piola-Kirchhoff stress of neo-Hookean material
    FTi=inv(F.T)
    P = 0.5*mu*(F-FTi) + 0.25*kappa*ln(J)*FTi

    # product with 'baseline' deformation gradient
    Q = Hp.T*P*F0.T*iH.T

    return Q

g = inv(H)*B;
gk = inv(H)*Bk;


w = TestFunction(V)
vk = Function(V)
v = Function(V)

Qk = Q(v)
Q_ = Q(vk)

# initalization of displacement between OSS and ROSS
v.assign(Constant((0.0, 0.0, 0.0)))
vk.assign(Constant((0.0, 0.0, 0.0)))

## ************ VARIATIONAL PROBLEM ************ 
a = (inner((nu*v/dt),w) + inner(theta*Qk.T,grad(w)))*dx
L = (inner(theta*g + (1-theta)*gk + nu*vk/dt,w) - inner((1-theta)*Q_.T,grad(w)))*dx

FF = a - L

du = TrialFunction(V)
JJ  = derivative(FF, v, du)  # Gateaux derivative in dir. of du

problem = NonlinearVariationalProblem(FF, v, [], JJ);
solver  = NonlinearVariationalSolver(problem)

## ************ SOLVER DEFINITION ************ 
prm = solver.parameters
prm['nonlinear_solver'] = 'snes'
prm["snes_solver"]["maximum_iterations"] = 50
prm["snes_solver"]["report"] = True
prm["snes_solver"]["error_on_nonconvergence"] = False

prm['snes_solver']['linear_solver']= 'gmres'
prm['snes_solver']['preconditioner']= 'ilu'


currMeshCoords=OSS.coordinates()+np.reshape(project(v,V1).vector().get_local(),(-1,3));
currShapeCoords=OAS.coordinates()+np.reshape(project(v,V1).vector().get_local(),(-1,3));

# interpolation function for f2
f2_interp = SphericalInterpolator(currMeshCoords,np.reshape(f2_,(-1,1)))
# interpolation function for H
H_interp =  SphericalInterpolator(currMeshCoords,np.reshape(project(Hp_0,T1).vector().get_local(),(-1,9)))
# interpolation function for OAS
OAS_interp =SphericalInterpolator(currMeshCoords,currShapeCoords)


## ************ TIME STEPPING ************ 
dt=Constant(1)
maxIts=1e3
minTimeStep=1e-2
minTimeStepChange=1e-6

t = dt
it = 0
while it < maxIts:
    
        
    (no_of_iterations,converged) = solver.solve();
    
    if not(converged):
        if dt.values()[0] <= minTimeStep:
            print "COULD NOT LINEARIZE PROBLEM. SORRY."
            sys.exit(1)
        
        dt.assign(0.5*dt.values()[0]);
    
        print "time step did not converge. reducing time-step to %0.4f" % dt.values()[0]
        sys.stdout.flush()    
        
        continue
        
    it += 1
        
    changeToLast=sum((v.vector().get_local()-vk.vector().get_local())**2)

    tmp=np.reshape(v.vector().get_local()[:],(-1,3))
    print 'it %d . dt %0.3f change %0.6e (max: %0.4f)' % (it,dt,changeToLast,np.max(tmp))
    
    v_fixed=np.nan_to_num(v.vector().get_local())
    v.vector()[:]=v_fixed;
    vk.assign(v)
    f2k.assign(f2)
    
    # 'relax' spherical mapping to fit solution
    newMeshCoords=OSS.coordinates()[:]+np.reshape(project(v,V1).vector().get_local()[vertex_to_dof_map(V1)],(-1,3));
    lons,lats,rs=cart2sph(newMeshCoords[:,0],newMeshCoords[:,1],newMeshCoords[:,2])
    xx,yy,zz=sph2cart(lons,lats,np.ones_like(rs))
    newMeshCoords=np.transpose([xx,yy,zz])

    f2.vector()[vertex_to_dof_map(S1)]=f2_interp(newMeshCoords)
    ttmp.vector()[vertex_to_dof_map(T1)]=np.reshape(H_interp(newMeshCoords),(-1,1))
    Hp=interpolate(ttmp,T)
    
    if changeToLast<=minTimeStepChange:
        print 'converged (%0.6e <= %0.6e)' % (changeToLast,minTimeStepChange)
        break
    ​

now, I have a couple of questions.

Did I do anything terribly stupid? This is the first time I'm working in continuum mechanics and FEM, so this is definitely a possibility. I'm using quadratic elements for the deformation (because nonlinear elasticity) and linear elements for the body force (because it's faster??). Does that make sense?
Are there ways to increase convergence speed? I'm using Crank-Nicholson - does that make sense? Is there a better way for this type of problem? What is the optimal time-step? I did try picking a very large one and checking for convergence of the linear solver, reducing the time-step on divergence, but that does not seem optimal, and if the initial step is too large I end up with really bad solutions as well. Is there literature on time-step selection and adaptive time-stepping for this type of problem?
Even though I don't think it is the main bottleneck, is it possible to avoid reassembly in this type of problem as in the Nonlinear Poisson or Navier-Stokes tutorials?
The authors originally used COMSOL to solve the PDE and they report better optima than I was able to achieve. Also, their results seem to indicate some discrete events in the optimization (discontinuities in the functions used for convergence analysis for instance). Does anybody have an idea about what tricks comsol (or other 'big' packages) pull when solving this type of problems?

any help is greatly appreciated - it's been fun implementing all of this, but I see that there is a whole lot more to learn!
Community: FEniCS Project
Please login to add an answer/comment or follow this question.

Similar posts:
Search »