Yarn Modelization


I am a begginner on Fenics and for my studies, I am trying to do some yarn modelization. To help me do it I use a model presented in “Real time simulation of inextensible surgical thread using a kirchhoff rod model with force output for haptic feedback applications” by Zhujiang Wang.

Here what I already did but it is clearly not functional:

from dolfin import *
from fenics import *
from ufl import Jacobian, diag
from plotting import plot
import numpy as np
import matplotlib.pyplot as plt

#Beam Caractéristique
L = 1.
M = Constant(0.01)
k_b = Constant(0.)
k_t = Constant(0.02)

N_L = 20
mesh = IntervalMesh(N_L, 0, L)

TotalTime = 10
TimeStep = Constant(0.001)
NumStep = TotalTime/TimeStep

#Finite Elements
Ue = VectorElement("CG", mesh.ufl_cell(), 1, dim=3)
Space = FunctionSpace(mesh, Ue*Ue)

u_ = TestFunction(Space)
du = TrialFunction(Space)
(w_, theta_) = split(u_)
(dw, dtheta) = split(du)

def both_ends(x, on_boundary):
    return on_boundary

def left_end(x, on_boundary):
    return near(x[0], 0) and on_boundary

def right_end (x, on_boundary):
    return near (x[0], L) and on_boundary

#Boundary Conditions
bc = [DirichletBC(Space.sub(1), Constant(0.), both_ends),
      DirichletBC(Space.sub(2), Constant(0.), both_ends),
      DirichletBC(Space.sub(0), Constant(0.), left_end)]

for i in range(N_L):
    e = x(i+1) - x(i)
    e /= sqrt(dot(x(i+1) - x(i), x(i+1) - x(i)))

#Matérial Frame {T, M1, M2}
ez = as_vector([0, 0, 1])

for i in range(N_L):
    M1 = cross(e(i), ez)
    M1 /= sqrt(dot(M1(i), M1(i)))
    M2 = cross(e(i), M1(i))
    M2 /= sqrt(dot(M2(i), M2(i)))

#Bishop Frame {T, U, V}
#définition de U
#définition de V

TotalTorsion = np.arcos(dot(M1(N_L), U(N_L))) - np.arcos(dot(M1(1), U(1)))

li = L / N_L

for i in range(N_L):
    psi = dot(x(i+1) - x(i), x(i+1) - x(i)) - pow(li, 2)
limoy = sqrt(dot(x(i) - x(i-1), x(i) - x(i-1))) + sqrt(dot(x(i+1) - x(i), x(i+1) - x(i)))

for i in range(2, N_L):
    Lmoy += limoy(i)
for i in range(2, N_L):
    Energy += k_b/limoy(i)*(2/(1 + dot(e(i-1), e(i))) - 1) + k_t*pow(TotalTorsion, 2)/Lmoy(i)

Phi = theta_[1, 2]
Beta = cos(Phi)
BetaPoint = NoIdea

for i in range(2, N_L):
    PowerDiss += k_b*pow(BetaPoint, 2)/2

My first question is about IntervalMesh. Isn’t it supposed to create some vertices x ? When I try to run my code, I have an error saying it’s not defined. Is there anything already made to have the straight edges of the mesh ?

Is someone know how to create a Bishop Frame ?

How can I do some time derivative ?

first, I think you are mixing numpy functions with UFL operators.
I suggest you start with very simple models first to understand properly how FEniCS works. Time-dependent problems are treated here for instance:
Regarding beam models, you can find demos here:
Also, definition of local frames is discussed here: