Hello,
My code below is doing optimization using Ipopt but it is sending me back my initial estimate, not optimizing. I simply want to implement the following problem.
Also, how can I add filter values and sensitivity on Ipopt?
I couldn’t find what I did wrong, and Ipopt gave me output after a certain iteration:
“The current filter has 0 entries.
minimal step size ALPHA_MIN = 0.000000E+00”
I couldn’t find my problem, I will be glad if you can help. Thank you so much!
Here is my code,
l1 = 1e-1
l2 = 1e-9
normal = Constant((0, 1))
E1, nu1 = 26.0, 0.25 # base material
E2, nu2 = 1e-10, 0.0 # void
mu1, lmbda1 = Constant(E1 / (2 * (1 + nu1))), Constant(E1 * nu1 / ((1 + nu1) * (1 - 2 * nu1)))
mu2, lmbda2 = Constant(E2 / (2 * (1 + nu2))), Constant(E2 * nu2 / ((1 + nu2) * (1 - 2 * nu2)))
max_disp = -0.14
f_desired = -0.019
p = 10
eps = Constant(1e-6)
c1, c2 = 1, 1e+1
chi = project(Expression("x[1]", name='Control', degree=3), W)
def H(a):
return 0.5 * (1 + a / (1e-10 + (a **2 ) ** (1/2)))
# Create mesh and define function space
n = 32
mesh = UnitSquareMesh(n, n, "crossed")
W = FunctionSpace(mesh, "CG", 1)
V_vector = VectorFunctionSpace(mesh, "CG", 1)
W_tensor = TensorFunctionSpace(mesh, "CG", 1)
force_FS = VectorFunctionSpace(mesh, "Real", 0)
df = TestFunction(force_FS)
disp = Constant([0.0, max_disp])
bcs, dss = bcs_(mesh, V_vector, disp)
u0 = Function(V_vector)
u = Function(V_vector)
Ic, J_, v, F = kinematics(u)
# Stored strain energy density (compressible neo-Hookean model)
Psi = (eps - (1 - eps) * ((H((chi - 0.5), l2)) ** 3)) * psi(mu1, Ic, J_, lmbda1)
# 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)
step_size = np.linspace(0, max_disp, 50)
for j in step_size:
disp = Constant([0.0, j])
bcs, dss = bcs_(mesh, V_vector, disp)
solve(dPi== 0, u, bcs,
form_compiler_parameters=ffc_options)
# Piola-Kirchhoff Stress
stress = diff(Psi, F)
stress_calculation = (project(stress, W_tensor))
file2 << stress_calculation
t = dot(stress, normal)
t = assemble(-t[1] * dss(1))
max_error = ((t - f_desired) ** p) ** (1 / p)
J = c1 * max_error + assemble((inner(grad(chi), grad(chi)) * dx))
print("max_error", max_error)
print("J", J)
control = Control(chi)
volume_constraint = UFLInequalityConstraint(
H((chi - 0.5), l2) / 1 * dx, control)
rf = ReducedFunctional(J, control, eval_cb_post=eval_cb)
problem = MinimizationProblem(rf, bounds=(0.0, 1.0), constraints=volume_constraint)
parameters = {"acceptable_tol": 1.0e-150, "maximum_iterations": 100}
solver = IPOPTSolver(problem, parameters=parameters)
f_opt = solver.solve()