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 ?