Hi all,
I want my force to reach the values I want with the applied displacement. I determine these desired points in the polynomial. So I minimize my function to find the polynomial coefficients that will allow it to pass my target points in the force-displacement graph.
However I am having an error like:
vector = vector._cpp_object
^^^^^^^^^^^^^^^^^^
AttributeError: ‘float’ object has no attribute ‘_cpp_object’
The error starts with minimization.
I have divided the code into pieces and left my minimal example below. Even though I optimize only the first constant in the example, I will later ask it to optimize all constants and many polynomials will be included in my geometry.
I can provide you more information if needed.
Thank you!
Geometry-Boundary-Function Definitions (No error here)
from dolfin import *
from dolfin_adjoint import *
import gmsh
import numpy as np
import matplotlib.pyplot as plt
import meshio
import moola
a = 0.1
t = 1
h = 1
x_dir = 1.5
y_dir = h + a
scale = 1
N = 10
lc = 1/2
max_disp = -0.3
last_point_of_spline = 0.25
Traclist = []
Traclist2 = []
normal = Constant((0, 1, 0))
recu = []
J = 0
multiply = 10
x_valuesA = [0, 0.1, 0.225, last_point_of_spline]
y_valuesA = [0, multiply*0.32, multiply*0.1, multiply*0.15]
desired_response_x = [0, 0.1, 0.225, last_point_of_spline]
desired_response_y = [0, multiply*0.36, multiply*0.02, multiply*0.07]
def createGeometryAndMesh(a, t, lc):
# Clear all models and create a new one
gmsh.initialize()
gmsh.clear()
gmsh.model.add("t3")
def geopoints(dx, dy, surf_point):
P1 = gmsh.model.occ.addPoint(0 + dx, dy, 0, lc)
P2 = gmsh.model.occ.addPoint(t + dx, dy, 0, lc)
P3 = gmsh.model.occ.addPoint(dx, h + dy, 0, lc)
P4 = gmsh.model.occ.addPoint(t + dx, h + dy, 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], 11 + surf_point)
return Surface1
# Entity creation operations
SurfaceA = geopoints(0, 0, 2000)
gmsh.model.occ.synchronize()
# One group for each of the unit squares
pg1 = gmsh.model.addPhysicalGroup(2, [SurfaceA])
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) + "geoforstl" + ".geo_unrolled")
gmsh.write("./Result--" + str(commit) + "geo.msh")
# colors:
gmsh.option.setNumber("Geometry.PointNumbers", 1)
gmsh.option.setColor("Geometry.Points", 255, 165, 0)
gmsh.option.setColor("General.Text", 255, 255, 255)
gmsh.option.setColor("Mesh.Points", 255, 0, 0)
rr, gg, bb, aaa = gmsh.option.getColor("Geometry.Points")
gmsh.option.setColor("Geometry.Surfaces", rr, gg, bb, aaa)
gmsh.onelab.set("""[
{"type":"number",
"name":"Parameters/Twisting angle",
"values":[90],
"min":0,
"max":120,
"step":1}]""")
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(a, t, 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) + "/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_A = Measure("dx", domain=mesh, subdomain_data=mf, subdomain_id=1)
V = VectorFunctionSpace(mesh, "CG", 2) # Lagrange
W_tensor = TensorFunctionSpace(mesh, "CG", 2)
V_tensor = TensorFunctionSpace(mesh, "CG", 2)
W = FunctionSpace(mesh, "R", 0)
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, max_disp, 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)
Calculation starts here
Error is also here
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 = (C-I)/2
return Ic, J, F, C, I, E
def poly(x_values, y_values):
# Fit a cubic polynomial (degree 3) to the given data points
coeffs = np.polyfit(x_values, y_values, len(x_values) -1)
print(coeffs)
# print(coeffs.convert().coef)
if len(coeffs) == 4:
c1, c2, c3, c4= coeffs
if len(coeffs) == 3:
c1 = 0
c2 = coeffs[0]
c3 = coeffs[1]
c4 = coeffs[2]
if len(coeffs) == 2:
c1 = 0
c2 = 0
c3 = coeffs[0]
c4 = coeffs[1]
return c1, c2, c3, c4
for i in np.linspace(0, max_disp, N):
u.assign(u0)
disp.assign(Constant([0.0, i, 0.0]))
print("i", i)
Ic, J, F, C, I, E = kinematics(u)
# Define the custom strain energy density function
def psi(x_values, y_values, num):
aA, bA, cA, dA= poly(x_values, y_values)
aA = interpolate(Expression(("a_constant"), a_constant=aA, degree=1), W) # W = FunctionSpace(mesh, "R", 0)
Ic, J, F, C, I, E = kinematics(u)
psi = (aA((0,h,0))/ 4) * (E[1, 1]) ** 4 - (bA / 3) * (E[1, 1]) ** 3 + (cA / 2) * (E[1, 1]) ** 2 - dA * E[1, 1] + 1E+4 * E[0,0] **2 + 1E+4 * E[0,1] **2
stress = diff(psi, F)
return psi, stress, aA, bA, cA, dA
psiA, stressA, aA, bA, cA, dA = psi(x_valuesA, y_valuesA, 1)
a_d, b_d, c_d, d_d = poly(desired_response_x, desired_response_y)
target_y = (a_d) * (i) ** 3 + (b_d) * (i) ** 2 + (c_d ) * (i) + d_d
psi = psiA
PiA = psiA * dx_A
Pi = PiA
dPiA = derivative(PiA, u, v)
dPi = derivative(Pi, u, v)
stress = stressA
solve(dPi == 0, u, bcs,
form_compiler_parameters=ffc_options)
print("u", u((0, h, 0))[1])
traction = dot(stress, normal)
stored_traction = assemble((traction[1]) * dss(1)) # Traction in the y-direction
Traclist.append(stored_traction)
"""
ERROR STARTS HERE
"""
J = (target_y - stored_traction) ** 2
# aA, bA, cA, dA = Function(W), Function(W), Function(W),Function(W)
# aA = interpolate(Expression(('a_constant'), a_constant=aA, degree=1), W)
# bA = interpolate(Expression(('b_constant'), b_constant=bA, degree=1), W)
# cA = interpolate(Expression(('c_constant'), c_constant=cA, degree=1), W)
# dA = interpolate(Expression(('d_constant'), d_constant=dA, degree=1), W)
control = Control(aA) # , Control(bA), Control(cA), Control(dA)
rf = ReducedFunctional(J, control)
problem = MoolaOptimizationProblem(rf)
n_moola = moola.DolfinPrimalVector(aA)
# d_moola = moola.DolfinPrimalVector(bA)
# k_moola = moola.DolfinPrimalVector(cA)
# l_moola = moola.DolfinPrimalVector(dA)
m_moola = moola.DolfinPrimalVectorSet([n_moola]) # , d_moola, k_moola, l_moola
solver = moola.BFGS(problem, m_moola, options={'jtol': 0,
'gtol': 1e-9,
'Hinit': "default",
'maxiter': 10,
'mem_lim': 10})
sol = solver.solve()
m_opt = sol['control'].data
print(m_opt)
