Creating smooth surface models with ElasticSurfaceEmbedding.jl

Creating smooth surface models with ElasticSurfaceEmbedding.jl Yuto Horikawa (@Hyrodium) 2023-07-28

View slides online https://www.docswell.com/s/hyrodium/5JL8EQ-JuliaCon2023 1

• https://www.docswell.com/s/hyrodium/5JL8EQ-JuliaCon2023

About me - ➀ Name • Yuto Horikawa • @Hyrodium (GitHub, Twitter, bsky, etc.) I'm a maintainer of: • Rotations.jl • StaticArrays.jl • Quaternions.jl • FastGaussQuadrature.jl • IntervalSets.jl • ImageClipboard.jl • BasicBSpline.jl • ElasticSurfaceEmbedding.jl 2

• https://github.com/JuliaGeometry/Rotations.jl
• https://github.com/JuliaArrays/StaticArrays.jl
• https://github.com/JuliaGeometry/Quaternions.jl
• https://github.com/JuliaApproximation/FastGaussQuadrature.jl
• https://github.com/JuliaMath/IntervalSets.jl
• https://github.com/hyrodium/ImageClipboard.jl
• https://github.com/hyrodium/BasicBSpline.jl
• https://github.com/hyrodium/ElasticSurfaceEmbedding.jl

About me - ➁ Creating mathematical handicrafts is my lifework! 3

Contents 1. Overview 2. Our formulation framework 3. More topics on the weaving theory 4. Results 4

Contents 1. Overview 2. Our formulation framework 3. More topics on the weaving theory 4. Results 5

Overview (1) Divide the surface into strips (2) Embed the strips onto a flat plane (3) Cut out the embedded strips from a paper (4) Weave the strips together 6

Overview (1) Divide the surface into strips (2) Embed the strips onto a flat plane (3) Cut out the embedded strips from a paper (4) Weave the strips together 7

Overview (1) Divide the surface into strips (2) Embed the strips onto a flat plane (3) Cut out the embedded strips from a paper (4) Weave the strips together 8

Overview (1) Divide the surface into strips (2) Embed the strips onto a flat plane (3) Cut out the embedded strips from a paper (4) Weave the strips together 9

Overview (1) Divide the surface into strips (2) Embed the strips onto a flat plane (3) Cut out the embedded strips from a paper (4) Weave the strips together 10

Overview (1) Divide the surface into strips (2) Embed the strips onto a flat plane (3) Cut out the embedded strips from a paper (4) Weave the strips together 11

Overview (1) Divide the surface into strips (2) Embed the strips onto a flat plane (3) Cut out the embedded strips from a paper (4) Weave the strips together Questions: • How can we characterize the embedded shapes? • How can we calculate the shapes? 12

Overview (1) Divide the surface into strips (2) Embed the strips onto a flat plane (3) Cut out the embedded strips from a paper (4) Weave the strips together Questions: • How can we characterize the embedded shapes? • How can we calculate the shapes with Julia? 13

Let's start with code!

1
2
3
4
5
6
7
8
9
10
11
12

using ElasticSurfaceEmbedding
using IntervalSets
using StaticArrays
# Overload the shape definition
ElasticSurfaceEmbedding.surface(x,y) =
SVector(x, y, x^2+y^2)
# (1) split the surface into strips
dom = [(-1..1, (i-1)/10..i/10)
for i in 1:10]
# (2) Embed the strips onto a plane
res = auto_allsteps(dom)
export_pinned_steps("paraboloid", res)

(2) Embed the
strips onto a
flat plane

14

Let's start with code!

1
2
3
4
5
6
7
8
9
10
11
12

using ElasticSurfaceEmbedding
using IntervalSets
using StaticArrays
# Overload the shape definition
ElasticSurfaceEmbedding.surface(x,y) =
SVector(x, y, x^2+y^2)
# (1) split the surface into strips
dom = [(-1..1, (i-1)/10..i/10)
for i in 1:10]
# (2) Embed the strips onto a plane
res = auto_allsteps(dom)
export_pinned_steps("paraboloid", res)

(2) Embed the
strips onto a
flat plane

Define the shape of the target surface

15

Let's start with code!

1
2
3
4
5
6
7
8
9
10
11
12

using ElasticSurfaceEmbedding
using IntervalSets
using StaticArrays
# Overload the shape definition
ElasticSurfaceEmbedding.surface(x,y) =
SVector(x, y, x^2+y^2)
# (1) split the surface into strips
dom = [(-1..1, (i-1)/10..i/10)
for i in 1:10]
# (2) Embed the strips onto a plane
res = auto_allsteps(dom)
export_pinned_steps("paraboloid", res)

(2) Embed the
strips onto a
flat plane

Define the domains of the strips

16

Let's start with code!

1
2
3
4
5
6
7
8
9
10
11
12

using ElasticSurfaceEmbedding
using IntervalSets
using StaticArrays
# Overload the shape definition
ElasticSurfaceEmbedding.surface(x,y) =
SVector(x, y, x^2+y^2)
# (1) split the surface into strips
dom = [(-1..1, (i-1)/10..i/10)
for i in 1:10]
# (2) Embed the strips onto a plane
res = auto_allsteps(dom)
export_pinned_steps("paraboloid", res)

(2) Embed the
strips onto a
flat plane

Calculate the embeddings of strips

17

Let's start with code!

1
2
3
4
5
6
7
8
9
10
11
12

using ElasticSurfaceEmbedding
using IntervalSets
using StaticArrays
# Overload the shape definition
ElasticSurfaceEmbedding.surface(x,y) =
SVector(x, y, x^2+y^2)
# (1) split the surface into strips
dom = [(-1..1, (i-1)/10..i/10)
for i in 1:10]
# (2) Embed the strips onto a plane
res = auto_allsteps(dom)
export_pinned_steps("paraboloid", res)

(2) Embed the
strips onto a
flat plane

Generates SVG files in paraboloid/pinned/

18

Contents 1. Overview 2. Our formulation framework 3. More topics on the weaving theory 4. Results 19

Our formulation framework • Mathematical & physical modeling • 2D Riemannian manifold • Geometric nonlinear elasticity • Numerical methods • B-spline manifold • Galerkin method • Newton-Raphson method 20

Mathematical & physical modeling - ➀ • Geometry • Approximate 3D medium with 2D surface • Reference state M[0] and current state M[t] are diffeomorphic with Φ • M[0] = (M, g[0] ), M[t] = (M, g[t] ) are Riemannian manifolds Curved paper strip • Elasticity • Hooke's law (material linear model) • Stress is proportional to strain • S ij = C ijkl Ekl (C : stiffness tensor) • Two types of deformations • in-plane (dominant) • out-plane (ignored) 21

Mathematical & physical modeling - ➀ • Geometry • Approximate 3D medium with 2D surface • Reference state M[0] and current state M[t] are diffeomorphic with Φ • M[0] = (M, g[0] ), M[t] = (M, g[t] ) are Riemannian manifolds Curved paper strip • Elasticity • Hooke's law (material linear model) • Stress is proportional to strain • S ij = C ijkl Ekl (C : stiffness tensor) • Two types of deformations • in-plane (dominant) • out-plane (ignored) 22

Mathematical & physical modeling - ➀ • Geometry • Approximate 3D medium with 2D surface • Reference state M[0] and current state M[t] are diffeomorphic with Φ • M[0] = (M, g[0] ), M[t] = (M, g[t] ) are Riemannian manifolds • Elasticity • Hooke's law (material linear model) • Stress is proportional to strain • S ij = C ijkl Ekl (C : stiffness tensor) • Two types of deformations • in-plane (dominant) • out-plane (ignored) 23

Mathematical & physical modeling - ➀ • Geometry • Approximate 3D medium with 2D surface • Reference state M[0] and current state M[t] are diffeomorphic with Φ • M[0] = (M, g[0] ), M[t] = (M, g[t] ) are Riemannian manifolds • Elasticity • Hooke's law (material linear model) • Stress is proportional to strain • S ij = C ijkl Ekl (C : stiffness tensor) • Two types of deformations • in-plane (dominant) • out-plane (ignored) 24

Mathematical & physical modeling - ➁ • How can we characterize the embedded shape M[t] ? • Minimize the elastic energy W of the embedding Φ ∫ √ 1 W = C ijkl Eij Ekl det g[0]ij du1 ∧ du2 i,j 2 M • Energy minimization problem leads a weak from PDE: • xi[t] : unknown functions • ξ i : arbitary test funcitons ) ( ∫ s ∂xp[t] ∂xq[t] ∂ξ r ∂x[t] ijkl δrs − g[0]ij υ[0] = 0 C δpq i j k ∂u ∂u ∂u ∂ul M 25

Numerical methods • B-spline manifold • Approximate the embedded shape with piecewise polynomial xi[t] ≈ aiI N I • N I : 2D B-spline basis function • Galerkin method • Discretize PDE into nonlinear equations 0 = ajJ akK alL δij δkl AIJKL − ajJ δij B IJ =: FiI • AIJKL , B IJ : some constants • Newton-Raphson method • Solve the nonlinear equations numerically aiI = aiI − ν+1 ν ∂aiI J Fj ∂FjJ ν ν • represents ν-th iteration ν 26

Contents 1. Overview 2. Our formulation framework 3. More topics on the weaving theory 4. Results 27

Details of the shape computation The output of auto_allsteps(dom) was like this: 1 15: Initial state - domain: [-1, 1]×[0.2, 0.3] 2 └─16: Newton onestep - residual norm: 4.7580e-02, Δa norm: 2.2337e-02, computation 3 └─17: Newton onestep - residual norm: 2.7644e-03, Δa norm: 1.3057e-02, computati 4 └─18: Refinement - p₊:(0, 1), k₊:([-0.9166666666666667, -0.75, -0.583333333333 5 └─19: Newton onestep - residual norm: 1.2488e-03, Δa norm: 1.4208e-02, compu 6 └─20: Newton onestep - residual norm: 5.3362e-06, Δa norm: 3.5955e-04, com 7 └─21: 📌 Newton onestep - residual norm: 2.5164e-09, Δa norm: 1.9453e-06 Each step has its embedded shape and strain distribution: 15: Initial state - Approximate the embedded shape (somehow) 28

Details of the shape computation The output of auto_allsteps(dom) was like this: 15: Initial state - domain: [-1, 1]×[0.2, 0.3] └─16: Newton onestep - residual norm: 4.7580e-02, Δa norm: 2.2337e-02, computation └─17: Newton onestep - residual norm: 2.7644e-03, Δa norm: 1.3057e-02, computati └─18: Refinement - p₊:(0, 1), k₊:([-0.9166666666666667, -0.75, -0.583333333333 └─19: Newton onestep - residual norm: 1.2488e-03, Δa norm: 1.4208e-02, compu └─20: Newton onestep - residual norm: 5.3362e-06, Δa norm: 3.5955e-04, com └─21: 📌 Newton onestep - residual norm: 2.5164e-09, Δa norm: 1.9453e-06 Each step has its embedded shape and strain distribution: 16: Newton onestep - Relax the shape to decrease the strain energy 29

Details of the shape computation The output of auto_allsteps(dom) was like this: 15: Initial state - domain: [-1, 1]×[0.2, 0.3] └─16: Newton onestep - residual norm: 4.7580e-02, Δa norm: 2.2337e-02, computation └─17: Newton onestep - residual norm: 2.7644e-03, Δa norm: 1.3057e-02, computati └─18: Refinement - p₊:(0, 1), k₊:([-0.9166666666666667, -0.75, -0.583333333333 └─19: Newton onestep - residual norm: 1.2488e-03, Δa norm: 1.4208e-02, compu └─20: Newton onestep - residual norm: 5.3362e-06, Δa norm: 3.5955e-04, com └─21: 📌 Newton onestep - residual norm: 2.5164e-09, Δa norm: 1.9453e-06 Each step has its embedded shape and strain distribution: 17: Newton onestep - Relax the shape to decrease the strain energy 30

Details of the shape computation The output of auto_allsteps(dom) was like this: 15: Initial state - domain: [-1, 1]×[0.2, 0.3] └─16: Newton onestep - residual norm: 4.7580e-02, Δa norm: 2.2337e-02, computation └─17: Newton onestep - residual norm: 2.7644e-03, Δa norm: 1.3057e-02, computati └─18: Refinement - p₊:(0, 1), k₊:([-0.9166666666666667, -0.75, -0.583333333333 └─19: Newton onestep - residual norm: 1.2488e-03, Δa norm: 1.4208e-02, compu └─20: Newton onestep - residual norm: 5.3362e-06, Δa norm: 3.5955e-04, com └─21: 📌 Newton onestep - residual norm: 2.5164e-09, Δa norm: 1.9453e-06 Each step has its embedded shape and strain distribution: 18: Refinement - More control points for more accuracy 31

Details of the shape computation The output of auto_allsteps(dom) was like this: 15: Initial state - domain: [-1, 1]×[0.2, 0.3] └─16: Newton onestep - residual norm: 4.7580e-02, Δa norm: 2.2337e-02, computation └─17: Newton onestep - residual norm: 2.7644e-03, Δa norm: 1.3057e-02, computati └─18: Refinement - p₊:(0, 1), k₊:([-0.9166666666666667, -0.75, -0.583333333333 └─19: Newton onestep - residual norm: 1.2488e-03, Δa norm: 1.4208e-02, compu └─20: Newton onestep - residual norm: 5.3362e-06, Δa norm: 3.5955e-04, com └─21: 📌 Newton onestep - residual norm: 2.5164e-09, Δa norm: 1.9453e-06 Each step has its embedded shape and strain distribution: 19: Newton onestep - Relax the shape to decrease the strain energy 32

Details of the shape computation The output of auto_allsteps(dom) was like this: 15: Initial state - domain: [-1, 1]×[0.2, 0.3] └─16: Newton onestep - residual norm: 4.7580e-02, Δa norm: 2.2337e-02, computation └─17: Newton onestep - residual norm: 2.7644e-03, Δa norm: 1.3057e-02, computati └─18: Refinement - p₊:(0, 1), k₊:([-0.9166666666666667, -0.75, -0.583333333333 └─19: Newton onestep - residual norm: 1.2488e-03, Δa norm: 1.4208e-02, compu └─20: Newton onestep - residual norm: 5.3362e-06, Δa norm: 3.5955e-04, com └─21: 📌 Newton onestep - residual norm: 2.5164e-09, Δa norm: 1.9453e-06 Each step has its embedded shape and strain distribution: 20: Newton onestep - Relax the shape to decrease the strain energy 33

Details of the shape computation The output of auto_allsteps(dom) was like this: 15: Initial state - domain: [-1, 1]×[0.2, 0.3] └─16: Newton onestep - residual norm: 4.7580e-02, Δa norm: 2.2337e-02, computation └─17: Newton onestep - residual norm: 2.7644e-03, Δa norm: 1.3057e-02, computati └─18: Refinement - p₊:(0, 1), k₊:([-0.9166666666666667, -0.75, -0.583333333333 └─19: Newton onestep - residual norm: 1.2488e-03, Δa norm: 1.4208e-02, compu └─20: Newton onestep - residual norm: 5.3362e-06, Δa norm: 3.5955e-04, com └─21: 📌 Newton onestep - residual norm: 2.5164e-09, Δa norm: 1.9453e-06 Each step has its embedded shape and strain distribution: 21: Newton onestep - Relax the shape to decrease the strain energy 34

Initial shape determination - ➀ How can we calculate this initial state of the Newton-Raphson method? 15: Initial state - Approximate the embedded shape (somehow) Next theorem is useful to construct the initial shape 35

Initial shape determination - ➁ Theorem (Embedding approximation) Let C[0] be the center curve of M[0] , κ[0] be its geodesic curvature, B be the breadth from center curve of M[0] . Similarly, let C[t] be the center curve of M[t] , κ[t] be its planer curvature. If the breadth B is sufficiently small, then the following approximation is satisfied. g[t] |C ≈ g[0] |C κ[t] ≈ κ[0] (Isometric on the center curve) (The curvature is conserved) 36

Strip breadth trade-offs - ➀ Paraboloid case: devide the surface into 10 strips • Large breadth: ⟨0⟩ • Increases strain E11 in paper medium* • Causes non-smooth surface model • Small breadth: • Complicates assembly process Identifying optimal breadth is essential for assembling smooth surface. ⟨0⟩ ⟨0⟩ * E11 is a coefficient of the Green's strain tensor E = Eij θ ⟨0⟩i ⊗ θ ⟨0⟩j where {θ ⟨0⟩i } is a dual orthonormal frame of the reference state M[0] . 37

Strip breadth trade-offs - ➁ Theorem (Strain approximation) The stress state is approximately one-dimensional, and the principal ⟨0⟩ strain E11 can be approximated by ( ) 1 1 ⟨0⟩ E11 ≈ K[0] B 2 r2 − 2 3 if the breadth B is sufficiently small. Where K[0] is a Gaussian curvature of the surface, r is a normalized (|r| ≤ 1) breadth-directional coordinate. show_strain(dom) # predict strain, should be < 1% 38

Contents 1. Overview 2. Our formulation framework 3. More topics on the weaving theory 4. Results 39

Results • Paraboloid • Hyperbolic paraboloid • Helicatenoid • Stereographic projection 40

Paraboloid z = x2 + y 2 41

Hyperbolic paraboloid z = x2 − y 2 • The surface is doubly ruled. • Pulling the surface reminds us a straight line on the surface. 42

Helicatenoid Catenoid   cosh(v) cos(u)   p[0] (u, v) =  cosh(v) sin(u)  v Helicoid   sinh(v) cos(u)   p[0] (u, v) =  sinh(v) sin(u)  u 43

Helicatenoid It deforms to each other! 44

Stereographic projection   2u 2 + v2 + 1 u     2v   p[0] (u, v) =  2  2 u + v + 1  2  2 u +v −1 u2 + v 2 + 1 45

Thank you! The following packages helped ElasticSurfaceEmbedding.jl a lot! • FastGaussQuadrature.jl • ForwardDiff.jl • Luxor.jl • IntervalSets.jl • BasicBSpline.jl • Documenter.jl Special Gifts! Weaving kits for: • Paraboloid • Hyperbolic paraboloid • Stereographic projection 46

Backup slides • Why weaving? • Swapping states • Coordinate system and notations • References 47

Why weaving? - ➀ There are several ways to create surfaces: • Solid models (e.g. 3D printing) a • Combining cross-sections b • Combining paper strips • Combining strips in one-direction c • Weaving strips in hexiagonal shape d • (Weaving strips in checker pattern) a Twitter: Twitter: c Twitter: d Twitter: b @henryseg @onomahedron @jmitani @alisonmartin57 48

• https://twitter.com/henryseg
• https://twitter.com/onomahedron
• https://twitter.com/jmitani
• https://twitter.com/alisonmartin57

Why weaving? - ➁ Benefits of weaving strips in checker pattern: • Easy to split the surface along with its coordinates • Enough strength • Deformable • Totally thin 49

Swapping states The embedding Φ is characterized by minimizing the strain enery. • Physically correct way (left): • Φ : M[0] → S • Flat material deforms into curved strip • Computationally easy way (right): • Φ : M[0] → E2 • Curved strip deforms into flat material • Tension / Compression are swapped • Strain energies of Φ are equivalent 50

Coordinate system and notations 51

References • Yuto Horikawa and Ryuichi Tarumi. “Weaving paper strips for designing of general curved surface with geometrical elasticity”. In: arXiv:2211.06372 [cs] (2022). • Arthur Ogawa. “Helicatenoid”. In: Mathematica Journal 2.2 (1992). • Javier Bonet and Richard D Wood. “Nonlinear Continuum Mechanics for Finite Element Analysis”. In: (2008). • Shigeyuki Morita. Geometry of Differential Forms. Translations of Mathematical Monographs v. 201. Providence, R.I: American Mathematical Society, 2001. ISBN : 978-0-8218-1045-3. 52