Hi, I’m trying to get a specific plot with the block analysis below. I have made my stress and strain energy definitions according to the spline I have defined, but still the graph I obtained does not match the graph I am trying to obtain. I couldn’t find what I’m missing, I would be very happy if you could help. thanks
Here, how I defined stress:
def psi(u, i, Ic, J, F, C, I):
def strain_energy_density(strain):
a, b, c, d = calculate_spline(x_values, y_values)
energy = a * strain **3 + b * strain **2 + c * strain + d
stress = (3 * a * strain ** 2 + b* 2 * strain +c) * grad(u)
return energy, stress
strain = 0.5 * inner(grad(u), grad(u))
energy, stress = strain_energy_density(strain)
return energy, stress
and this is the full code:
import gmsh
import numpy as np
from dolfin import *
import matplotlib.pyplot as plt
import meshio
import itertools
from scipy.interpolate import CubicSpline
from scipy import interpolate
b = 1
gmsh.initialize()
commit = ("zwei")
N = 50
lc = 1 # mesh
def createGeometryAndMesh(b,lc):
# Clear all models and create a new one
gmsh.initialize()
gmsh.clear()
gmsh.model.add("t3")
P1 = gmsh.model.occ.addPoint(0 , 0, 0, lc)
P2 = gmsh.model.occ.addPoint(b, 0, 0, lc)
P3 = gmsh.model.occ.addPoint(0, 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, "main_geometry")
# 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()
gmsh.write("./Result--" + str(commit)+ "geo.msh")
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, lc)
# Importing mesh from gmsh and defining surface and boundary markers
msh = meshio.read("./Result--" + str(commit)+ "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(commit)+ "mesh.xdmf",
triangle_mesh)
meshio.write("mesh.xdmf", triangle_mesh)
# 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()
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)
File("./Result--" + str(commit)+ "MSH.pvd") << mesh
File("./Result--"+ str(commit) + "MSH2.pvd") << mf
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(commit)+ "MSH3.pvd") << boundary_markers
# Compiling subdomains
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx_main = Measure("dx", domain=mesh, subdomain_data=mf, subdomain_id=1)
V = VectorFunctionSpace(mesh, "CG", 2) # Lagrange
# Define the stress tensor function space
V_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], scale * b, tol)
disp = Constant([0.0, 0.0, 0.0])
bc_top = DirichletBC(V, disp, 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
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, I
# Elasticity parameters
E1, nu1 = 1e-20, 0.0 # void
mu1, lmbda1 = Constant(E1 / (2 * (1 + nu1))), Constant(
E1 * nu1 / ((1 + nu1) * (1 - 2 * nu1)))
E2, nu2 = 30.5, 0.4 # main
mu2, lmbda2 = Constant(E2 / (2 * (1 + nu2))), Constant(
E2 * nu2 / ((1 + nu2) * (1 - 2 * nu2)))
def calculate_spline(x_values, y_values):
# Fit a cubic polynomial (degree 3) to the given data points
coeffs = np.polyfit(x_values, y_values, 3)
a, b, c, d = coeffs
return a, b, c, d
# Define the spline values
x_values = [0, 0.06, 0.19, 0.25]
y_values = [0, 0.3, 0.05, 0.07]
# Solve variational problem
store_u = np.linspace(0, -0.25, N, endpoint=True)
traction = []0))
fileu = open("./Result--" + str(commit)+ "/displacement.txt", "w")
filet = open("./Result--" + str(commit)+ "/traction" + ".txt", "w")
recu = []
strain_list = []
for i in store_u:
disp.assign(Constant([0.0, i, 0.0]))
u.assign(u0)
Ic, J, F, C, I = kinematics(u)
# Define the custom strain energy density function
def psi(u, i, Ic, J, F, C, I):
def strain_energy_density(strain):
a, b, c, d = calculate_spline(x_values, y_values)
energy = a * strain **3 + b * strain **2 + c * strain + d
stress = (3 * a * strain ** 2 + b* 2 * strain +c) * grad(u)
return energy, stress
strain = 0.5 * inner(grad(u), grad(u))
energy, stress = strain_energy_density(strain)
return energy, stress
psi, stress = psi(u, i, Ic, J, F, C, 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()
u0.assign(u)
# Create a new Function with the appropriate FunctionSpace to store the stress tensor
stress_calculation = Function(V_tensor, name='Stress')
stress_calculation.assign(project(stress, V_tensor))
# Compute traction (force) on the boundary
traction = dot(stress_calculation, normal)
traction_save = assemble((traction[1]) * dss(1))
traction.append(traction_save)
recu.append(i)
np.savetxt("./Result--" + str(commit)+ "/traction" +".txt", Traclist)
np.savetxt("./Result--" + str(commit)+ "/displacement.txt", recu)
# Plot the force-displacement curve
df = np.linspace(0, -0.25, N)
K = -1
dff2 = [y * K for y in Traction]
# Calculate spline coefficients
a, b, c, d = calculate_spline(x_values, y_values)
# Generate x values for the spline plot
x_spline = np.linspace(0, 0.25, N)
# Calculate corresponding y values using the cubic polynomial
y_spline = a * x_spline**3 + b * x_spline**2 + c * x_spline +d
y_spline2 = 3 * a * x_spline ** 2 + b* 2 * x_spline +c
# Plot the data points and the spline
plt.plot(x_values, y_values, 'o', label='Data Points')
plt.plot(x_spline, y_spline, label='Cubic Spline')
plt.plot(x_spline, dff2, label='Model' )
plt.legend()
plt.grid(True)
plt.xlabel("Displacement")
plt.ylabel("Force")
plt.savefig('Graph.png')