Specific DOF of a Function as control variable

dolfin-adjoint
3

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 *
from fenics_adjoint import *

def fix_dofs(func,fixed_indexes,fixed_values):             
    vec = func.vector()
    vec[fixed_indexes] = fixed_values
    return Function(func.function_space(), vec)
  • Overloading the function with the following module fix_dofs_overloaded.py:
from fenics import *
from fenics_adjoint import *
from pyadjoint import Block
from pyadjoint.overloaded_function import overload_function

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
        self.add_dependency(func)

        #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'

    def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx, prepared=None):
        adj_input = adj_inputs[0] 
        x = inputs[0].vector() 

        output = adj_input
        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[0],self.fixed_indexes,self.fixed_values) 

fix_dofs = overload_function(fix_dofs, FixedDofsBlock)
  • The main script / MWE (Minimal Working Example):
from fenics import *
from fenics_adjoint import *
import numpy as np   
from numpy.random import rand

from fix_dofs_overloaded import fix_dofs

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, 'CG', 1)
uh = project(Expression("sin(pi*x[0])*cos(pi*x[1])",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)[0]
free_indexes = np.where(dof_V_x < 0.5)[0]

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

J = assemble(inner(grad(u), grad(u)) * dx) # Some Functional 

m = Control(uh)
Jhat = ReducedFunctional(J, m)

h = Function(V)
h.vector()[:] = rand(h.vector().local_size())
taylor_test(Jhat, uh, h)