Hallo, I want to calculate the traction and strain curve for a beam with a nonlinear spline to model the strain energy. Basically I want to have a specific spline response in my final plot. I could not find how to do that, I attached the code I am working on. Could you please help me? Thanks.
I am doing a huge mistake probably but here is my code:
import gmsh
import sys
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline
# 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}
un = "trial"
b = 1
lc = 1e-1 * 10
def createGeometryAndMesh(b):
# Clear all models and create a new one
gmsh.initialize()
gmsh.clear()
gmsh.model.add("t3")
P1 = gmsh.model.occ.addPoint(-b , 0, 0, lc)
P2 = gmsh.model.occ.addPoint(b, scale *0, 0, lc)
P3 = gmsh.model.occ.addPoint(-b, b , 0, lc)
P4 = gmsh.model.occ.addPoint(b , b , 0, lc)
L1 = gmsh.model.occ.addLine(P1, P2)
L2 = gmsh.model.occ.addLine(P2, P4)
L3 = gmsh.model.occ.addLine(P4, P3)
L4 = gmsh.model.occ.addLine(P3, P1)
Curve1 = gmsh.model.occ.addCurveLoop([L1, L2, L3, L4])
Surface1 = gmsh.model.occ.addPlaneSurface([Curve1], 1)
# Synchronize the CAD model
gmsh.model.occ.synchronize()
# Create physical groups
pg1 = gmsh.model.addPhysicalGroup(2, [Surface1])
gmsh.model.setPhysicalName(2, pg1, "beam")
# Set mesh size at all points
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc)
gmsh.model.occ.synchronize()
gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
gmsh.model.mesh.generate()
def checkForEvent():
action = gmsh.onelab.getString("ONELAB/Action")
if len(action) and action[0] == "check":
gmsh.onelab.setString("ONELAB/Action", [""])
createGeometryAndMesh()
gmsh.graphics.draw()
return True
# Refine the mesh
gmsh.model.mesh.refine()
gmsh.write("t3.opt");
gmsh.finalize()
createGeometryAndMesh(b)
# Importing mesh from gmsh and defining surface and boundary markers
msh = meshio.read("./Result--" + str(un)+ "geo.msh")
for key in msh.cell_data_dict["gmsh:physical"].keys():
if key == "triangle":
triangle_data = msh.cell_data_dict["gmsh:physical"][key]
for cell in msh.cells:
if cell.type == "tetra":
tetra_cells = cell.data
elif cell.type == "triangle":
triangle_cells = cell.data
triangle_mesh = meshio.Mesh(points=msh.points,
cells=[("triangle", triangle_cells)],
cell_data={"name_to_read": [triangle_data]})
meshio.write("./Result--" + str(un)+ "mesh.xdmf",
triangle_mesh)
meshio.write("mesh.xdmf", triangle_mesh)
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
print(mesh.topology().dim() - 1)
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1, 0)
class Boundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0.0)
b_c = Boundary()
b_c.mark(boundary_markers, 3)
File("./Result--" + str(un)+ "MSH3.pvd") << boundary_markers
# Compiling subdomains
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx = Measure("dx", domain=mesh, subdomain_data=mf, subdomain_id=1)
V = VectorFunctionSpace(mesh, "CG", 2)
W_tensor = TensorFunctionSpace(mesh, "CG", 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-7
return on_boundary and near(x[1], b, tol)
disp = Constant([0.0, 0.0, 0.0])
bc_top = DirichletBC(V, disp, boundary_top)
bcs = [bc_bot, bc_top]
du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement from previous iteration
u0 = Function(V) # Initial displacement
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)
def kinematics(u):
# Kinematics
d = u.geometric_dimension() # Space 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)
return Ic, J, F, C
from scipy import interpolate
# Strain data points and energy values
strain_values = [0.0, 0.6, 0.19, 0.25]
energy_values = [0.0, 0.03, 0.005, 0.047]
if any(strain_values[i] >= strain_values[i+1] for i in range(len(strain_values)-1)):
# If there are repeated or non-unique values, preprocess the data
unique_strain_values, unique_indices = np.unique(strain_values, return_index=True)
unique_energy_values = [energy_values[i] for i in unique_indices]
strain_values = list(unique_strain_values)
energy_values = list(unique_energy_values)
# Calculate the spline for strain-energy interpolation
spline = CubicSpline(strain_values, energy_values)
# Solve variational problem
store_u = np.linspace(0, -0.4, 50 , endpoint=True)
normal = Constant((0, 1, 0))
# Define the stress tensor function space
V_tensor = TensorFunctionSpace(mesh, "CG", 2)
force = []
for i in store_u:
disp.assign(Constant([0.0, i, 0.0]))
Ic, J, F, C = kinematics(u)
# Define the custom strain energy density function
def psi(F, u, i):
strain = i
energy = Constant(spline(strain))
stress = diff(energy, F)
return energy, stress
psi, stress = psi(F, u, i)
# 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 second variation of Pi (directional derivative about du in the direction of v)
ddPi= derivative(dPi, u, du)
# Define the nonlinear variational problem with Jacobian
problem = NonlinearVariationalProblem(dPi u, bcs, J=ddPi)
solver = NonlinearVariationalSolver(problem)
solver.solve()
# Create a new Function with the same FunctionSpace as stress to store the stress tensor
stress_calculation = Function(W_tensor, name='Stress')
stress_calculation.assign(stress) # Assign the stress tensor to stress_calculation directly
# Compute traction (force) on the boundary
traction = dot(stress, normal)
force.append(assemble((traction[1]) * dss(1)) )
