# Specific DOF of a Function as control variable

Hello,

Using dolfin_adjoint, is it possible to have, as control variables, some specific DOF of a function?
I know that we can have all of the DOF as control variables but what if I want to keep some of them fixed and only use the free ones as control variables?

Here is a minimal example showing what I would like to do:

``````from dolfin import *
import numpy as np

mesh = UnitSquareMesh(50,50)
V = FunctionSpace(mesh, "CG", 2)

n = V.dim()
d = mesh.geometry().dim()

# Recover the coordinates of the dof
dof_V_coordinates = V.tabulate_dof_coordinates()
dof_V_coordinates.resize((n, d))
dof_V_x = dof_V_coordinates[:, 0]
dof_V_y = dof_V_coordinates[:, 1]

# Some function as initialization
initial = Expression("sin(pi*x)*cos(pi*x)",degree=3)
u = interpolate(initial, V)

# find some indexes of the free and fixed dof
fixed_indexes = np.where(dof_V_x >= 0.5)
free_indexes = np.where(dof_V_x < 0.5)

fixed_value = 0.4
u.vector()[fixed_indexes] =  fixed_value

# Working example:
m = Control(u) #Control variables (all the dofs of the Function u)
Jhat = ReducedFunctional(J, m)
dJdm = Jhat.derivative()

# Non-working example:
m_1 = Control(u.vector()[free_indexes]) # Control variables (some dofs of the Function u)
Jhat_1 = ReducedFunctional(J, m_1)
dJdm_1 = Jhat_1.derivative()
``````

I get: AttributeError: ‘numpy.ndarray’ object has no attribute ‘block_variable’ at the line where m_1 is declared. Does it mean that I cannot achieve this without using Custom functions? And if yes what are your suggestions?
To have it working, just comment the last three lines that are just there to better explain what I am trying to achieve.

The simplest solution is to modify the gradient of the full problem, aka:

``````m = Control(u) #Control variables (all the dofs of the Function u)
Jhat = ReducedFunctional(J, m)
dJdm = Jhat.derivative()
dJdm.vector()[fixed_indexes] = 0
``````

However, this implies that you cannot use the built-in optimization methods.

To do this properly, you need to make an intermediate variable, which is the one that will be used in your expression. I have added a minimal working example of this strategy below:

``````from dolfin import *
import numpy as np

mesh = UnitSquareMesh(4,4)
V = FunctionSpace(mesh, "CG", 2)

# # Some function as initialization
initial = Expression("sin(pi*x)*cos(pi*x)",degree=3)
uh = interpolate(initial, V)
uh.rename("u","")
fixed_value = 0.4

class fixed_indices(SubDomain):
def inside(self, x, on_boundary):
return x >= 0.5
mf = MeshFunction("size_t", mesh, 2, 0)
fixed_domain = fixed_indices()
fixed_domain.mark(mf, 1)
dx_c = Measure("dx", domain=mesh, subdomain_data=mf)
u, v = TrialFunction(V), TestFunction(V)
a = inner(u,v)*dx
l = inner(uh, v)*dx_c(0) + inner(Constant(fixed_value), v)*dx_c(1)

u = Function(V)
solve(a==l, u)

m = Control(uh) #Control variables
Jhat = ReducedFunctional(J, m)
dJdm = Jhat.derivative()
print(dJdm.vector().get_local())
``````

which yields

``````[ 0.2934843   0.60295222  0.07689859 -0.42947131 -0.1527415   0.90639199
0.19458717  1.4761873  -0.52695545 -0.75811268  0.50790983 -0.03113176
-0.99850619 -3.74725313 -0.39626439 -0.34956218  0.80677139 -0.0541857
-0.38557608 -0.34352392  0.2649522   0.09599802 -1.75726921 -2.9235389
-1.60314787  0.          0.          0.         -0.34876289 -0.13449154
0.2601133   0.20486509 -0.92068632  0.33012212  0.0787908   0.31704804
-1.19407037 -0.84793381  0.          0.          0.          0.
0.          0.          0.         -0.63982149  0.56348174  0.37797151
-0.02842044  1.38792627  0.9383983   0.19644893  0.          0.
0.          0.          0.          0.          0.          0.
-0.17723889  1.08410744  1.55911393  0.          0.          0.
0.          0.          0.          0.          0.          0.
0.          0.          0.          0.          0.          0.           0.          0.]
``````
1 Like

Thank you very much for your answer. I indeed wanted to use the built-in optimization methods.
I understand the method you proposed with the intermediate variable.
However, I was hoping for something less costly than solving a PDE each time I want to fix some DOFs. I plan on integrating it into some optimization based on a large number of degrees of freedom.

After quite some time, I think that I found a “more efficient” way to implement it with a custom function as introduced in dolfin_adjoint documentation.
Here is how I proceeded:

• Custom function in a file/module named fix_dofs.py :
``````from fenics import *

def fix_dofs(func,fixed_indexes,fixed_values):
vec = func.vector()
vec[fixed_indexes] = fixed_values
return Function(func.function_space(), vec)
``````
``````from fenics import *

from fix_dofs import fix_dofs

backend_fix_dofs = fix_dofs

class FixedDofsBlock(Block):
def __init__(self, func, fixed_indexes, fixed_values, **kwargs):
super(FixedDofsBlock, self).__init__()
self.kwargs = kwargs

#variables that do not appear in the dependencies but are still used for the computation
self.fixed_indexes = fixed_indexes
self.fixed_values = fixed_values

def __str__(self):
return 'FixedDofsBlock'

x = inputs.vector()

output[self.fixed_indexes] = 0.0 # the derivative of the fixed dofs is null
return output

def recompute_component(self, inputs, block_variable, idx, prepared):
return backend_fix_dofs(inputs,self.fixed_indexes,self.fixed_values)

``````
• The main script / MWE (Minimal Working Example):
``````from fenics import *
import numpy as np
from numpy.random import rand

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, 'CG', 1)
uh = project(Expression("sin(pi*x)*cos(pi*x)",degree=3), V)

# Recover the coordinates of the dof
dof_V_coordinates = V.tabulate_dof_coordinates()
dof_V_x = dof_V_coordinates[:, 0]
dof_V_y = dof_V_coordinates[:, 1]
# find some indexes of the free and fixed dof
fixed_indexes = np.where(dof_V_x >= 0.5)
free_indexes = np.where(dof_V_x < 0.5)

fixed_value= 0.0
u = fix_dofs(uh,fixed_indexes,fixed_value)