I have a “tricky” mixed problem that I can’t (easily) solve using PETSc/Tao/etc. I’m minimizing an energy functional, \Pi(\mathbf{q}), \mathbf{q}:=[\mathbf{u | a | b | c}] subject to some dirichlet boundary conditions on the displacement field \mathbf{u}, and additional inequality and bound constraints on the extra degrees of freedom a,b,c:

min \Pi(\mathbf{q}), \mathbf{q}:=[\mathbf{u | a | b | c}] \\ s.t. \\ u = u_D \forall x \in \Gamma_D \\a_i >0 \\ b_i >0 \\c_i >0 \\ c_i < \pi,\\ a_i >b_i \\

I got ok results with a DIY primal-dual interior point method (with dynamic relaxation), but I need something more robust.

I was wondering if there is an easy way to hack the IPOPT interface provided by dolfin adjoint. It seems like it should work, but I am not sure how to pass in my boundary conditions and bounds. Inspecting the `pyipopt`

interface to IPOPT https://github.com/live-clones/dolfin-adjoint/blob/master/dolfin_adjoint/optimization/ipopt_solver.py, it seems like the bounds correspond to `lb`

and `ub`

. I’m not sure what `nconstraints`

is (does it include bounds and constraints or just constraints?)…

For the modified hyperelasticity demo, how can I impose the Dirichlet boundary conditions /interface to IPOPT? I appreciate any tips. Thanks

```
from dolfin import *
from dolfin_adjoint import *
mesh = UnitCubeMesh(24, 16, 16)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
# Mark boundary subdomians
left = CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)
# Define Dirichlet boundary (x = 0 or x = 1)
c = Expression(("0.0", "0.0", "0.0"), degree=1)
r = Expression(("scale*0.0",
"scale*(y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta) - x[1])",
"scale*(z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta) - x[2])"),
scale = 0.5, y0 = 0.5, z0 = 0.5, theta = pi/3, degree=1)
bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, r, right)
bcs = [bcl, bcr]
du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement from previous iteration
B = Constant((0.0, -0.5, 0.0)) # Body force per unit volume
T = Constant((0.1, 0.0, 0.0)) # Traction force on the boundary
d = u.geometric_dimension()
I = Identity(d) # Identity tensor
F = I + grad(u) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
Ic = tr(C)
J = det(F)
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds
F = derivative(Pi, u, v)
J = derivative(F, u, du)
# Solve using newton
solve(F == 0, u, bcs, J=J) #works
file = File("newton.pvd");
file << u
# solve using ipopt
Jhat = ReducedFunctional(assemble(Pi), Control(u))
problem = MinimizationProblem(Jhat,
bounds=None, # default to +- inf
constraints=None # bcs go here???
)
solver = IPOPTSolver(problem, parameters={"maximum_iterations": 50,
"acceptable_tol": 1.0e-8})
u = solver.solve() # obviously does not work
file = File("ipopt.pvd");
file << u
```