Hello,
I have a target polynomial curve and I want to achieve that curve in my modeling. I shared my working code below.
I changed my strain energy density function by integrating my polynomial, aiming to match my final result with the desired curve. Also, I’ve shared both the target graph and the one my model produced. You can observe that my model deviates from the target, particularly towards the end. I couldn’t find the error that caused this. I would be grateful if you could help me.
import gmsh
import math
import sys
from dolfin import *
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline
from scipy import interpolate
import meshio
import numpy as np
t = 1
h = 1
max_disp = -0.25
last_point_of_spline = 0.25
gmsh.initialize()
commit = ("versuch")
N = 10
a, b, c, d = 0, 0, 0, 0
lc = 1 # mesh
def createGeometryAndMesh(t, lc):
# Clear all models and create a new one
gmsh.initialize()
gmsh.clear()
gmsh.model.add("t3")
# in units mm
# cube points:
P1 = gmsh.model.occ.addPoint(0, 0, 0, lc)
P2 = gmsh.model.occ.addPoint(t, 0, 0, lc)
P3 = gmsh.model.occ.addPoint(0, h, 0, lc)
P4 = gmsh.model.occ.addPoint(t, h, 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)
gmsh.model.occ.synchronize()
pg1 = gmsh.model.addPhysicalGroup(2, [Surface1])
gmsh.model.setPhysicalName(2, pg1, "1")
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) + "geoforstl" + ".geo_unrolled")
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
gmsh.write("t3.opt");
gmsh.finalize()
createGeometryAndMesh(t, lc)
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)
# Set log level
# set_log_level(LogLevel.WARNING)
# 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) + "/MSH.pvd") << mesh
File("./Result--" + str(commit) + "MSH2.pvd") << mf
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], h, 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)
# E = (F.T + F - 2*I)
E = (C - I) / 2
return Ic, J, F, C, I, E
def poly(x_values, y_values):
# Fit a cubic polynomial to the given data points
coeffs = np.polyfit(x_values, y_values, len(x_values) - 1)
print(coeffs)
# print(coeffs.convert().coe f)
if len(coeffs) == 4:
a, b, c, d = coeffs
if len(coeffs) == 3:
a = 0
b = coeffs[0]
c = coeffs[1]
d = coeffs[2]
if len(coeffs) == 2:
a = 0
b = 0
c = coeffs[0]
d = coeffs[1]
return a, b, c, d
x_values = [0, 0.1, 0.15, last_point_of_spline]
y_values = [0, 0.19, 0.19, 0.6]
x_values = [0, 0.08, 0.22, last_point_of_spline]
y_values = [0, 0.35, 0.09, 0.1]
# Solve variational problem
timestep = np.linspace(0, max_disp, N, endpoint=True)
stepu = []
force = []
normal = Constant((0, 1, 0))
for i in timestep:
u.assign(u0)
disp.assign(Constant([0.0, i, 0.0]))
print("u", u((0, h, 0))[1])
print("i", i)
stepu.append(i)
Ic, J, F, C, I, E = kinematics(u)
# Define the custom strain energy density function
a, b, c, d = poly(x_values, y_values)
psi = (a / 4) * (E[1, 1]) ** 4 - (b / 3) * (E[1, 1]) ** 3 + (c / 2) * (E[1, 1]) ** 2 - d * E[
1, 1] # + 1* E[0,0] **2
stress = diff(psi, F)
# 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)
# Define the nonlinear variational problem with Jacobian
solve(dPi == 0, u, bcs,
form_compiler_parameters=ffc_options)
u0.assign(u)
E_ = Function(V_tensor, name='E')
E_.assign(project(E, V_tensor))
stress_calculation = Function(V_tensor, name='stress')
stress_calculation.assign(project(stress, V_tensor))
storeforce = assemble((dot(stress, normal)[1]) * dss(1)) # Traction in the y-direction
force.append(storeforce)
K = -1
dff = [x * K for x in stepu]
dff2 = [y * K for y in force]
plt.plot(dff, dff2, 'o-', label='Model Simulation')
x_spline = np.linspace(0, last_point_of_spline, N * 5, endpoint=True)
y_spline = a * x_spline ** 3 + b * x_spline ** 2 + c * x_spline + d
plt.plot(x_values, y_values, 'o', label='Data Points')
plt.plot(x_spline, y_spline, label='Target Curve')
plt.legend()
plt.xlabel("X")
plt.ylabel("Y")