Hello everyone,

I’m using the example from the dolfin-adjoint documentation about topology optimization of the heat conduction problem.

I would like to add a **custom inequality constraint** to it, but I did not succeed up to now …

The constrained value is a scalar resulting from multiple operations such as the result of a variational form, integral over the domain, etc…

Is there a way to implement such a constraint, with its gradient/jacobian computed with the tools for automatic derivative provided by dolfin-adjoint?

*My first tentative and first error are the following:*

I tried to use the function `UFLInequalityConstraint()`

, but it returns an error stating:

`ValueError: Form must depend on control`

. It seems like the dependencies on the previous operations are not accounted for. Aren’t they recorded on the `Tape`

?

*My second tentative and second error are the following:*

I tried to create my own `InequalityConstraint`

class (named `CustomConstraint`

) as in the example but without success either … I get the following error:

`AttributeError: 'numpy.ndarray' object has no attribute 'block_variable'`

Here is a minimal working example. To have the first or second error appear, please set `flag_showError1`

or `flag_showError2`

to `True`

individually:

```
from fenics import *
from fenics_adjoint import *
import matplotlib.pyplot as plt
try:
from pyadjoint import ipopt # noqa: F401
except ImportError:
print("""This example depends on IPOPT and Python ipopt bindings. \
When compiling IPOPT, make sure to link against HSL, as it \
is a necessity for practical problems.""")
raise
# turn off redundant output in parallel
parameters["std_out_all_processes"] = False
flag_showError1 = False
flag_showError2 = False
Vlim = Constant(0.4) # volume bound on the control
Qlim = 1e-4
p = Constant(5) # power used in the solid isotropic material with penalisation (SIMP) rule, to encourage the control solution to attain either 0 or 1
eps = Constant(1.0e-3) # epsilon used in the solid isotropic material
alpha = Constant(1.0e-8) # regularisation coefficient in functional
def k(a):
return eps + (1 - eps) * a ** p
n = 250
mesh = UnitSquareMesh(n, n)
A = FunctionSpace(mesh, "CG", 1) # function space for control
P = FunctionSpace(mesh, "CG", 1) # function space for solution
class WestNorth(SubDomain):
"""The top and left boundary of the unitsquare, used to enforce the Dirichlet boundary condition."""
def inside(self, x, on_boundary):
return (x[0] == 0.0 or x[1] == 1.0) and on_boundary
bc = [DirichletBC(P, 0.0, WestNorth())] # the Dirichlet BC; the Neumann BC will be implemented implicitly by dropping the surface integral after integration by parts
f = interpolate(Constant(1.0e-2), P) # the volume source term for the PDE
def forward(a):
"""Solve the forward problem for a given material distribution a(x)."""
T = Function(P)
v = TestFunction(P)
F = inner(grad(v), k(a) * grad(T)) * dx - f * v * dx
solve(F == 0, T, bc, solver_parameters={"newton_solver": {"absolute_tolerance": 1.0e-7, "maximum_iterations": 20}})
return T
#-------------------------------------------------
class CustomConstraint(InequalityConstraint):
"""A class that enforces the volume constraint g(a) = V - a*dx >= 0."""
def __init__(self, Qlim):
self.Qlim = float(Qlim)
self.tmpvec = Function(A)
def function(self, m):
from pyadjoint.reduced_functional_numpy import set_local
set_local(self.tmpvec, m)
T = forward(self.tmpvec)
Q = assemble(inner(grad(T),k(self.tmpvec)*grad(T)) * dx)
return [-self.Qlim + Q]
def jacobian(self, m):
from pyadjoint.reduced_functional_numpy import set_local
set_local(self.tmpvec, m)
T = forward(self.tmpvec)
Q = assemble(inner(grad(T),k(self.tmpvec)*grad(T)) * dx)
dQdm = compute_gradient(Q, m)
return [dQdm]
def output_workspace(self):
return [0.0]
def length(self):
"""Return the number of components in the constraint vector (here, one)."""
return 1
#-------------------------------------------------
if __name__ == "__main__":
a = interpolate(Vlim, A) # initial guess.
T = forward(a) # solve the forward problem once.
J = assemble(f * T * dx + alpha * inner(grad(a), grad(a)) * dx)
m = Control(a)
Jhat = ReducedFunctional(J, m)
lb = 0.0; ub = 1.0
volume_constraint = UFLInequalityConstraint((Vlim - a)*dx, m)
if flag_showError1:
losses_constraint = UFLInequalityConstraint((-Qlim + inner(grad(forward(T)),k(T)*grad(forward(T))))*dx, m)
problem = MinimizationProblem(Jhat, bounds=(lb, ub), constraints=[volume_constraint,losses_constraint])
elif flag_showError2:
problem = MinimizationProblem(Jhat, bounds=(lb, ub), constraints=[volume_constraint,CustomConstraint(Qlim)])
else:
problem = MinimizationProblem(Jhat, bounds=(lb, ub), constraints=volume_constraint)
solver = IPOPTSolver(problem, parameters={"acceptable_tol": 1.0e-3, "maximum_iterations": 100})
a_opt = solver.solve()
T_opt = forward(a_opt)
pl, ax = plt.subplots(); fig = plt.gcf(); fig.set_size_inches(16, 4)
plt.subplot(1, 2, 1); p = plot(a_opt,title='Topology',mode='color',vmin=0,vmax=1.0); p.set_cmap("Greys"); cbar = plt.colorbar(p);
plt.subplot(1, 2, 2); p = plot(T_opt,title='Temperature',mode='color'); p.set_cmap("coolwarm"); cbar = plt.colorbar(p);
```

Is there a way to fix my errors or a more efficient was to implement such a constraint?

Thank you very much for your time!