Hello Everyone,

In my displacement controlled code below, I tried to store the calculated traction and displacement in each iteration to create the Force-Displacement diagram, but I could not find a way to do this.

If you have any suggestions, it would be very helpful.

Thanks for your help in advance!

```
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
# Create mesh and define function space
mesh = Mesh("FinerTriangularMesh05.xml")
V = VectorFunctionSpace(mesh, "Lagrange", 2)
def boundary_bot(x, on_boundary):
tol = 1E-7
return on_boundary and near(x[1], 0, tol)
bc_bot = DirichletBC(V, [0,0,0], boundary_bot)
def boundary_top(x, on_boundary):
tol = 1E-14
return on_boundary and near(x[1], 9, tol)
bc_top = DirichletBC(V, [0.0,-1.8,0], boundary_top)
bcs = [bc_bot, bc_top]
# Define functions
du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement from previous iteration
T = Constant((0.0, 0.0, 0.0)) # Traction force on the boundary
boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
AutoSubDomain(boundary_top).mark(boundary_subdomains, 1)
dss = ds(subdomain_data=boundary_subdomains)
# Kinematics
d = u.geometric_dimension()
I = Identity(d) # Identity tensor
F = variable(I + grad(u)) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
# Invariants of deformation tensors
Ic = tr(C)
J = det(F)
# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
TS = TensorFunctionSpace(mesh, "Lagrange", 2)
# Total potential energy
Pi = psi*dx
# Compute first variation of Pi (directional derivative about u in the direction of v)
dPi = derivative(Pi, u, v)
# Compute Jacobian of F
J = derivative(dPi, u, du)
# Solve variational problem
solve(dPi == 0, u, bcs, J=J,
form_compiler_parameters=ffc_options)
storedu=[]
storedu.append(u)
print(storedu)
#Piola-Kirchhoff Stress
def S(u):
pi = psi
stress = diff(psi, F)
return stress
#Traction
Stress = project(S(u), TS)
traction = Function(V)
e2 = Constant((0, 1, 0))
traction = project(dot(Stress,e2), V)
storedtraction = assemble(traction[1]*dss) #top surface traction
Traclist.append(storedtraction)
print(Traclist)
print(storedtraction)
# Save solution in VTK format
file = File("displacement.pvd");
file << u;
file = File("stress.pvd");
file << Stress;
file = File("traction.pvd");
file << traction;
# # Plot and hold solution
# plot(u, mode = "displacement", interactive = True)
print("Done")
# plot(storedtraction,storedu)
# plt.show()
[quote="HaticeGokcenGuner, post:1, topic:7351, full:true"]
Hello Everyone,
In my displacement controlled code below, I tried to store the calculated traction and displacement in each iteration to create the Force-Displacement diagram, but I could not find a way to do this.
If you have any suggestions, it would be very helpful.
Thanks for your help in advance!
---
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
# Create mesh and define function space
mesh = Mesh("FinerTriangularMesh05.xml")
V = VectorFunctionSpace(mesh, "Lagrange", 2)
def boundary_bot(x, on_boundary):
tol = 1E-7
return on_boundary and near(x[1], 0, tol)
bc_bot = DirichletBC(V, [0,0,0], boundary_bot)
def boundary_top(x, on_boundary):
tol = 1E-14
return on_boundary and near(x[1], 9, tol)
bc_top = DirichletBC(V, [0.0,-1.8,0], boundary_top)
bcs = [bc_bot, bc_top]
# Define functions
du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement from previous iteration
T = Constant((0.0, 0.0, 0.0)) # Traction force on the boundary
boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
AutoSubDomain(boundary_top).mark(boundary_subdomains, 1)
dss = ds(subdomain_data=boundary_subdomains)
# Kinematics
d = u.geometric_dimension()
I = Identity(d) # Identity tensor
F = variable(I + grad(u)) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
# Invariants of deformation tensors
Ic = tr(C)
J = det(F)
# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
TS = TensorFunctionSpace(mesh, "Lagrange", 2)
# Total potential energy
Pi = psi*dx
# Compute first variation of Pi (directional derivative about u in the direction of v)
dPi = derivative(Pi, u, v)
# Compute Jacobian of F
J = derivative(dPi, u, du)
# Solve variational problem
solve(dPi == 0, u, bcs, J=J,
form_compiler_parameters=ffc_options)
storedu=[]
storedu.append(u)
print(storedu)
#Piola-Kirchhoff Stress
def S(u):
pi = psi
stress = diff(psi, F)
return stress
#Traction
Stress = project(S(u), TS)
traction = Function(V)
e2 = Constant((0, 1, 0))
traction = project(dot(Stress,e2), V)
storedtraction = assemble(traction[1]*dss) #top surface traction
Traclist.append(storedtraction)
# Save solution in VTK format
file = File("displacement.pvd");
file << u;
file = File("stress.pvd");
file << Stress;
file = File("traction.pvd");
file << traction;
print("Done")
# plot(storedtraction,storedu)
# plt.show()
```