Direct image-analysis methods for surgical simulation and mixed - - PowerPoint PPT Presentation

direct image analysis methods for surgical simulation and
SMART_READER_LITE
LIVE PREVIEW

Direct image-analysis methods for surgical simulation and mixed - - PowerPoint PPT Presentation

Direct image-analysis methods for surgical simulation and mixed meshfree methods Jack S. Hale Marie Curie Postdoctoral Fellow Research Unit in Engineering Science Universit du Luxembourg Collaborators: Stephane P. A. Bordas, Universite du


slide-1
SLIDE 1

Direct image-analysis methods for surgical simulation and mixed meshfree methods

Jack S. Hale Marie Curie Postdoctoral Fellow Research Unit in Engineering Science Université du Luxembourg

Collaborators: Stephane P. A. Bordas, Universite du Luxembourg/Cardiff University Pierre Kerfriden, Cardiff University Juan J. Rodenas Garcia, Universidad Politecnica de Valencia

  • Pedro M. Baiz, Imperial College London

Alejandro Ortiz-Bernardin, Universidad de Chile Christian J. Cyron, Yale University

Université Pierre-et-Marie-Curie Institut Jean le Rond d’Alembert 26th May 1400

slide-2
SLIDE 2

An Overview

slide-3
SLIDE 3

The problems and my research

Mesh distortion in hyperelastic problems

Numerical locking in incompressible elasticity, plates and shells, Cosserat elasticity

Enriched implicit boundary methods Meshfree methods Mixed variational forms Patch projection techniques Direct image- analysis transition for surgical simulators.

slide-4
SLIDE 4

Project 1: Direct image to mesh analysis for surgical simulation

slide-5
SLIDE 5

Surgical Simulation

Model-order reduction Enriched and implicit boundary methods High Performance Computing Real-time simulation

slide-6
SLIDE 6

RealTCut and SHACRA Inria

slide-7
SLIDE 7
slide-8
SLIDE 8

COLONIX, OSIRIX

How can we move from an image…

slide-9
SLIDE 9

…or perhaps a series of images…

Source: COLONIX, OSIRIX

slide-10
SLIDE 10

Pipelines to analysis

Acquire images Segment images Mesh Surfaces Mesh Volume Perform analysis

Traditional

Acquire images Segment images Perform analysis

Implicit Boundary

NURBS

Implicit Explicit

}

}

{

Implicit

slide-11
SLIDE 11

Nested Octree

Og Od M Discretisation Geometry

slide-12
SLIDE 12

How to transfer geometric information back to the discretisation?

V pd

hd (Od)

V pg

hg (Og)

pd > pg hd > hg pg = 1 Enrichment

V pd

hd (Od)

  • E[V pg

hg (Og)]

slide-13
SLIDE 13

Octree Level 5/Level 3

slide-14
SLIDE 14

Surface

6a2e86c

slide-15
SLIDE 15

Features

  • 2D and 3D problems on the same code-path.
  • parallel hybrid MPI/TBB assembly and solution.
  • fast and robust computational geometry using CGAL.
  • automatic Delaunay tessellation of integration subdomains.
  • completely separate representation of discretisation and geometry via

nested octree data structures.

  • constructs to represent soft and hard segmentations of image data.
  • implicit representation via level-sets, inside-outside functions.
  • independent hpe-type adaptivity on discretisation and geometry.
  • fast refinement and coarsening operations.
slide-16
SLIDE 16

Outlook

  • We are developing a cartesian grid implicit

boundary/enriched finite element method toolkit within deal.ii.

  • By uncoupling discretisation and geometry we

hope to produce a method particularly suited to image-based analysis.

  • Many challenges ahead, particularly with

imposition of Dirichlet boundary conditions on hyperelastic soft-tissue model.

slide-17
SLIDE 17

Project 2: Mesh-based and Meshfree Partitions

  • f Unity
slide-18
SLIDE 18

Partition of Unity

uh(x) =

N

  • i=1

φiui

N

  • i=1

φi = 1

slide-19
SLIDE 19

Finite elements

0.0 0.2 0.4 0.6 0.8 1.0 x 0.0 0.2 0.4 0.6 0.8 1.0 φ

slide-20
SLIDE 20

Mesh-based Partition of Unity

slide-21
SLIDE 21

Finite Element Method

  • Reference

Mesh

slide-22
SLIDE 22

Meshfree-based Partition of Unity

slide-23
SLIDE 23

Meshfree partition of unity is constructed in global space

Good:

  • improved approximation

properties.

  • high-order continuity.
  • less dofs for given accuracy.
  • less sensitive to poor quality

node distributions.

  • eases mesh generation.
  • Not so good:
  • expensive.
  • extra flexibility also means extra

complexity:

  • integration.
  • adaptivity/error measures.
  • mathematical proof.
  • stability.
  • boundary conditions.
slide-24
SLIDE 24

Maximum-entropy basis functions

0.0 0.2 0.4 0.6 0.8 1.0 x 0.0 0.2 0.4 0.6 0.8 1.0 φ

slide-25
SLIDE 25

0.0 0.2 0.4 0.6 0.8 1.0 x −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 0.4 u

uh u ui

0.0 0.2 0.4 0.6 0.8 1.0 x −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 0.4 u

uh u ui

uh(x) =

N

  • i=1

φiui

RBF Maximum-Entropy

slide-26
SLIDE 26

Maximum-entropy basis functions

slide-27
SLIDE 27

Maximum-entropy basis functions

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.00 0.06 0.12 0.18 0.24 0.30 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.00 0.06 0.12 0.18 0.24 0.30 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

slide-28
SLIDE 28

Maximum-entropy basis functions

uh ∈ H1

0(Ω)

All basis functions associated with interior nodes vanish on the convex hull. Vital for easy imposition of Dirichlet boundary conditions.

slide-29
SLIDE 29

The problem.

slide-30
SLIDE 30

The problem: numerical locking

Problems with small parameters crop up nearly everywhere! Hyperelasticity Incompressible fluid flow Plates and Shells Cosserat elasticity and probably many more…

slide-31
SLIDE 31

Nearly-incompressible elasticity

λ = νE (1 + ν)(1 − 2ν) Find uh ∈ Uh such that

C(uh) : (v) dΩ =

  • Γ

f · v dΩ ∀v ∈ H1

0(Ω)

slide-32
SLIDE 32

Locking

102 103 dofs 10−1 100 eL2

t = 0.1 t = 0.01 t = 0.001

slide-33
SLIDE 33

Locking

Inability of the basis functions to satisfy the constraint imposed whilst still having adequate approximation properties

||u − uh|| ≤ C(λ)hp||u||

slide-34
SLIDE 34

p-refinement

slide-35
SLIDE 35

p-refinement

102 103 104 dofs 10−6 10−5 10−4 10−3 10−2 10−1 100 eH1(θ)

p = 1 p = 2 p = 3 p = 4 p = 5

slide-36
SLIDE 36

Mixed variational formulation

Find uh ∈ Ph and ph ∈ Ph such that: Looking better

(uh) · (v) dΩ +

ph · v dΩ =

f · v dΩ v Uh

· uh q dΩ 1

phq dΩ = 0 q Ph

slide-37
SLIDE 37

Great, so the problem is solved right?

slide-38
SLIDE 38

LBB Stability

inf

q∈Ph sup v∈Uh

  • Ω q∇ · vdΩ

||v||1||q||0/R ≥ βh > 0

A question of balance

Ph Uh

slide-39
SLIDE 39

Stability

(a) Bad approximation (b) Good approximation

Image: Marie E. Rognes

slide-40
SLIDE 40

MINI element

Arnold, Brezzi, Fortin 1984 Uh Ph Linear + Bubble Linear

slide-41
SLIDE 41

The Volume Averaged Nodal Pressure Method (VANP)

slide-42
SLIDE 42

Questions I will answer…

  • 1. Can we produce a stable pair of spaces for the mixed

formulation using meshfree approximation schemes?

  • 2. Can we eliminate the pressure space to produce a generalised

displacement method?

  • 3. Can we produce a general scheme, which works for arbitrary

spaces of meshfree basis functions and even finite element basis functions?

  • 4. How does enrichment, in the manner of the MINI bubble, affect

the convergence and stability?

  • 5. Does it work in 3D?
slide-43
SLIDE 43

Questions I can answer….

  • Exactly how to construct meshfree basis functions,

including high-order reproducing radial basis functions and maximum-entropy basis functions.

  • How to accurately integrate the weak form.
  • How to implement the method in algorithmic form.
  • And of course any other questions you may have…
slide-44
SLIDE 44
  • 1. Can we produce a stable pair of spaces for the

mixed formulation using meshfree approximation schemes?

slide-45
SLIDE 45

C T (Ω) N s N b C T (Ω) N s N b∗

  • Fig. 1: Schematic representation of a two-dimensional simplicial tessellation for the VANP method.

(a) Enhanced node set N +, and (b) enhanced node set N ∗.

slide-46
SLIDE 46

Spaces

Uh := [ME(Ω; Nh, ρ)]2 Ph := CG1(Ω; Th)

For simplicity of the exposition, not required

uh(x) =

N

  • i=1

φiui ph =

M

  • i=1

Nipi

slide-47
SLIDE 47

Saddle-point problem

BT DB dΩ u +

BT mN p dΩ p =

ΦT

uf dΩ

N T

p mT B dΩ u − 1

λe

N T

p N p dΩ p = 0

B =   

∂φux ∂x1 ∂φuy ∂x2 ∂φux ∂x2 ∂φuy ∂x1

   m = 1 1 0T D =   2µ 2µ µ  

slide-48
SLIDE 48
  • 2. Can we eliminate the pressure space to produce a

generalised displacement method?

slide-49
SLIDE 49

Saddle-point problem

A B BT −C u p

  • =

f

  • How can we get rid of p?
slide-50
SLIDE 50
slide-51
SLIDE 51

Projection

slide-52
SLIDE 52

Many methods include a projection or ‘softening’

Enhanced Assumed Strains (EAS) Reduced Integration Mixed Interpolation of Tensorial Components Smoothed Finite Element Method (SFEM)

slide-53
SLIDE 53

N

  • b=1

NpamT Bb dΩ ub − 1 λ

M

  • a=1

Npa dΩ pa = 0

N

  • b=1
  • Ωa

NpamT Bb dΩ ub − 1 λ

M

  • a=1
  • Ωa

Npa dΩ pa = 0 For every pressure node: Restrict integration domain to local domain: Re-arrange to get equation for every pressure dof: pa = −λ

N

  • b=1
  • Ωa NpamT Bb dΩ
  • Ωa Npa dΩ
  • ub
slide-54
SLIDE 54

Results

slide-55
SLIDE 55
  • 1. Can we produce a stable pair of spaces for the

mixed formulation using meshfree approximation schemes?

slide-56
SLIDE 56

Timoshenko Beam

y x P

100 mm 200 mm

slide-57
SLIDE 57

Timoshenko Beam

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 −7 −6 −5 −4 −3 −2 −1

ln(1/h) ln(energy norm of the error)

MINI VANP−T1+/T1−m VANP−T1*/T1−m 1 1

slide-58
SLIDE 58

Leaky-lid cavity flow

slide-59
SLIDE 59

0.0 0.2 0.4 0.6 0.8 1.0 y −0.4 −0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 ux

LPP MaxEnt 11 × 11 MINI 10 × 10 MINI 30 × 30

Velocity across P − P

slide-60
SLIDE 60

0.0 0.2 0.4 0.6 0.8 1.0 x −0.20 −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 0.20 uy

LPP MaxEnt 11 × 11 MINI 10 × 10 MINI 30 × 30

Velocity across Q − Q

slide-61
SLIDE 61
  • 3. Can we produce a general scheme, which works

for arbitrary spaces of meshfree basis functions?

slide-62
SLIDE 62

Timoshenko Beam

0.2 0.6 1 1.4 1.8 −7 −6 −5 −4 −3 −2 −1

ln(1/h) ln(energy norm of the error)

VANP−T2+/T1−m VANP−T2*/T1−m VANP−T2/T1−m VANP−T2+/T1−r VANP−T2*/T1−r VANP−T2/T1−r 1 2

slide-63
SLIDE 63
  • 4. How does enrichment, in the manner of the MINI

bubble, affect the convergence and stability?

slide-64
SLIDE 64

MINI* element

Kim and Lee, 2000, 10.1023/A:1018973303935 Uh Ph

slide-65
SLIDE 65

C T (Ω) N s N b C T (Ω) N s N b∗

  • Fig. 1: Schematic representation of a two-dimensional simplicial tessellation for the VANP method.

(a) Enhanced node set N +, and (b) enhanced node set N ∗.

slide-66
SLIDE 66

Constrained Block

MINI First-order MaxEnt Full Bubbles

slide-67
SLIDE 67

Constrained Block

MINI First-order MaxEnt Half Bubbles

slide-68
SLIDE 68

Constrained Block

MINI No bubbles

slide-69
SLIDE 69
  • 5. Does it work in 3D?
slide-70
SLIDE 70

(a) (b) (c) (d) (e) (f) (g) (h) (i)

  • Fig. 17: Compression of a constrained block: Nodal pressure variable for (a) MINI element, (b)

VANP-T +

1 /T1-m, (c) VANP-T ∗ 1 /T1-m, (d) VANP-T + 2 /T1-m, (e) VANP-T ∗ 2 /T1-m, (f) VANP-T2/T1-

m, (g) VANP-T +

2 /T1-r, (h) VANP-T ∗ 2 /T1-r, and (i) VANP-T2/T1-r.

slide-71
SLIDE 71

Summary

  • Stable meshfree methods are developed by

mimicking existing inf-sup stable finite element methods.

  • The auxiliary pressure variable is eliminated using

a volume-averaged nodal pressure technique.

  • Currently extending to hyperelasticity, where the

robust nature of the meshfree shape functions will be a great advantage.

slide-72
SLIDE 72

Acknowledgements

  • Jack S. Hale would like to thank:
  • The support of the FNR through the AFR Marie

Curie COFUND scheme.

  • The partial support of the Université du

Luxembourg.

  • FONDECYT Chile Foreign Scholar Grant.
slide-73
SLIDE 73

Questions?