# Yarn Modelization

Hello,

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.
http://dx.doi.org/10.1016/j.ijsolstr.2017.02.017

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)

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

#Temp
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)

#Function
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)]
"""

#Edges
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 ?

Hi,
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:
https://fenicsproject.org/pub/tutorial/html/._ftut1006.html
Regarding beam models, you can find demos here: