# Problem with functional

Hello everybody,

i have the following problem with my functional where the inhomogeneous Neumann condition is coupled in the weak formulation.

\int_{\Omega} |\nabla \phi|^2 dx = - \frac{1}{|\Gamma_4} * I_2 * \int_{\Omega_4} \phi ds

I_2 was computed through the Rosenbrock-Wanner method and it is a numpy.ndarray.
When i am trying to optimize the problem with

J = assemble(0.5 * dolfin.dot(dolfin.grad(u_h), dolfin.grad(v))*dx - (leitfaehigkeit*(1/area_gamma1)*((-sol[3,1]))*v*ds(4)))

control = Control(f)

rf = ReducedFunctional(J, control)

problem = MoolaOptimizationProblem(rf)
f_moola = moola.DolfinPrimalVector(f)
solver2 = moola.NewtonCG(problem, f_moola, options={'gtol': 1e-9,
'maxiter': 20,
'display': 3,
'ncg_hesstol': 0})

sol2 = solver2.solve()


I get the following error in the line solver2.solve()

AttributeError: 'dolfin.cpp.la.Matrix' object has no attribute 'block_variable'


Can somebody please explain me what is wrong?

Greetings

Please write this with Latex formatting, i.e. \int_\Omega (done by $encapsulation). Also, please add a minimal code example reproducing hte error, as this makes it alot easier for others to help you, and makes the problem more useful for others that might face the same problem Thank you dokken. I will do it. import numpy as np import matplotlib.pylab as plt import math from numpy.linalg import solve, norm from dolfin import * from dolfin_adjoint import * from mshr import * import pdb import moola # Form compiler options parameters["form_compiler"]["cpp_optimize"] = True parameters["form_compiler"]["optimize"] = True class ode(object): # Konstruktor def __init__(self, l_1:int, l_2:int, k:float, r_2:float, c_1:float, u_omega:float): self.l_1 = l_1 self.l_2 = l_2 self.k = k self.r_2 = r_2 self.c_1 = c_1 self.u_omega = u_omega def rhs(self,q1,p1,q2,p2,l_12): f = [p1, (-self.l_2/((self.l_1*self.l_2-l_12**2)*self.c_1))*q1+(l_12/((self.l_1*self.l_2-l_12**2)*self.c_1))*q2, p2, ((l_12*self.r_2)/((self.l_1*self.l_2-l_12**2)))*p1-((self.l_1*self.r_2)/(self.l_1*self.l_2-l_12**2))*p2] return f def semiImpRK(self,A,q1_i,p1_i,q2_i,p2_i,i,h): # k1 l_12 = self.k*np.sqrt(self.l_1*self.l_2) b1 = prob.rhs(q1_i,p1_i,q2_i,p2_i,l_12) k1 = np.linalg.solve(A,b1) # k2 b2 = prob.rhs(q1_i+0.5*h*k1[0],p1_i+0.5*h*k1[1],q2_i+0.5*h*k1[2],p2_i+0.5*h*k1[3],l_12) k2 = np.linalg.solve(A,b2) # yip1 yip1 = np.zeros(4) yip1[0] = q1_i + h*k2[0] yip1[1] = p1_i + h*k2[1] yip1[2] = q2_i + h*k2[2] yip1[3] = p2_i + h*k2[3] print("Ergebnis des " + str(i) + "ten Durchlaufs des Algorithmus\n") print(str(yip1)) return yip1[0], yip1[1], yip1[2], yip1[3] def problem(self): l_12 = self.k*np.sqrt(self.l_1*self.l_2) jac = np.array([[0, 1, 0, 0], [-self.l_2/((self.l_1*self.l_2-l_12**2)*self.c_1), 0, l_12/((self.l_1*self.l_2-l_12**2)*self.c_1), 0], [0, 0, 0, 1], [0, (l_12*self.r_2)/(self.l_1*self.l_2-l_12**2), 0, -(self.l_1*self.r_2)/(self.l_1*self.l_2-l_12**2)]]) t0 = 0.0 T = 100 # Maximale SchweiĂzeit 15ms N = 1000 h = T/float(N) # Die Jacobimatrix muss nur einmal berechnet werden, deshalb liegt diese Berechnung auĂerhalb der Schleife n = 4 a = 1.0/(2.0+math.sqrt(2.0)) I = np.identity(n) A = I - a*h*jac # Anfangswerte z = np.zeros((4,N)) z[0,0] = self.c_1 * self.u_omega z[1,0] = 0.0 z[2,0] = 0.0 z[3,0] = 0.0 t = np.zeros(N+1) sol = np.zeros((4,N+1)) # t = 0 sol[0,0], sol[0,1], sol[0,2], sol[0,3] = prob.semiImpRK(A,z[0,0],z[1,0],z[2,0],z[3,0],0,h) # t > 0 for i in range(1,N+1): q1_i, p1_i, q2_i, p2_i = sol[0,i-1], sol[1,i-1], sol[2,i-1], sol[3,i-1] sol[0,i], sol[1,i], sol[2,i], sol[3,i] = prob.semiImpRK(A,q1_i,p1_i,q2_i,p2_i,i,h) t[i] = t[i-1] + h # Plot------------------------------------------------- fig1 = plt.figure(1) fig1, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, sharex=True) ax1.plot(t, sol[0,:]) ax1.set_ylabel('A*s') ax1.set_title('Q_1') ax1.grid(True) ax2.plot(t, sol[1,:]) ax2.set_ylabel('A') ax2.set_title('I_1') ax2.grid(True) ax3.plot(t, sol[2,:]) ax3.set_ylabel('A*s') ax3.set_title('Q_2') ax3.grid(True) ax4.plot(t, sol[3,:]) ax4.set_ylabel('A') ax4.set_xlabel('Time$t\$ [s]')
ax4.set_title('I_2')
ax4.grid(True)
plt.show()
plt.savefig("LĂ¶sung_odeSol.png")

return sol

class pde(ode):

def solverIntern(self,a,v,dx,ds,leitfaehigkeit,area_gamma1,sol,f,u_h,bc2,i):

L = f*v*dx + (leitfaehigkeit*(1/area_gamma1)*(-sol[3,i]))*v*ds(4)
dolfin.solve(a == L, u_h, bc2)

return u_h

def achsenSetzen(self):
# Das Achsen-Objekt des Diagramms in einer Variablen ablegen:
ax = plt.gca()
# Die obere und rechte Achse unsichtbar machen:
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
# Die linke Diagrammachse auf den Bezugspunkt '0' der x-Achse legen:
ax.spines['left'].set_position(('data',0))
# Die untere Diagrammachse auf den Bezugspunkt '0' der y-Achse legen:
ax.spines['bottom'].set_position(('data',0))
# Ausrichtung der Achsen-Beschriftung festlegen:
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')

return ax

def energyfunctional(self,phi,v,ds,leitfaehigkeit,area_gamma1,mesh,u_h,f,sol):

W = FunctionSpace(mesh, "DG", 0)
f = interpolate(Expression("x[0]+x[1]", degree=1), W)

alpha = Constant(1e-6)
print("J = " + str(J))
control = Control(f)

rf = ReducedFunctional(J, control)

problem = MoolaOptimizationProblem(rf)
f_moola = moola.DolfinPrimalVector(f)
pdb.set_trace()
solver2 = moola.NewtonCG(problem, f_moola, options={'gtol': 1e-9,
'maxiter': 20,
'display': 3,
'ncg_hesstol': 0})

# sol2 = solver2.solve()
f_opt = sol['control'].data

plot(f_opt, title="f_opt")
plt.savefig("f_opt.png")

# Define the expressions of the analytical solution
f_analytic = Expression("1/(1+alpha*4*pow(pi, 4))*w", w=w, alpha=alpha, degree=3)
u_analytic = Expression("1/(2*pow(pi, 2))*f", f=f_analytic, degree=3)
#We can then compute the errors between numerical and analytical solutions.

f.assign(f_opt)
solve(F == 0, u_h, bc2)
control_error = errornorm(f_analytic, f_opt)
state_error = errornorm(u_analytic, u_h)
print("h(min):           %e." % mesh.hmin())
print("Error in state:   %e." % state_error)
print("Error in control: %e." % control_error)

def __init__(self,l_1,l_2,k,r_2,c_1,u_omega):
super().__init__(l_1,l_2,k,r_2,c_1,u_omega)

def solver(self,f,degree,sol):

fig3 = plt.figure(3)
# Gebiet und Mesh erstellen------------------------------------------------
domain_vertices = [Point(2.25, 1.75),
Point(1.0, 3.0),
Point(-1.0, 1.0),
Point(4.0, -4.0),
Point(5.0, -3.0),
Point(2.0, 0.0),
Point(0.5, 0.0)]

domain = Polygon(domain_vertices)
mesh = generate_mesh(domain,20)

# Plot und Speicherung des nicht konvexen Gebiets---------------------------------
plt.plot([0.5, 2.25, 1.0, -1.0, 4.0, 5.0, 2.0, 0.5], [0.0, 1.75, 3.0, 1.0, -4.0, -3.0, 0.0, 0.0])

# Das Achsen-Objekt des Diagramms in einer Variablen ablegen:
ax = prob_pde.achsenSetzen()
plt.savefig("Domain.png")

# Plot und Speicherung der Triangulierung------------------------------------------
plt.figure(4)
p = plot(mesh, title = 'Finite Element Mesh')
ax = prob_pde.achsenSetzen()
plt.savefig("Mesh.png")

# Randbedingungen-----------------------------------------------------------
class Left(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (between(x[1], (-0.999, 0.999)) and between(x[0], (3.999, -3.999)))

class Right(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and ((between(x[1], (2.25, 1.75)) and  between(x[0], (1.001, 2.999))) and (between(x[1], (0.5, 0.0)) and between(x[0], (2.249, 1.749))) and (between(x[1], (1.999, 0.0)) and between(x[0], (0.501, 0.0))) and (between(x[1], (4.999, -2.999)) and between(x[0], (2.0, 0.0))))

class Bottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (between(x[1], (5.0, -3.0)) and between(x[0], (4.0, -4.0)))

class Top(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (between(x[1], (1.0, 3.0)) and between(x[0], (-1.0, 1.0)))

# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()

# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)

boundaries.set_all(4)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

# Defining the function space for electric field
V_phi = FunctionSpace(mesh, "CG", degree)
#V_phi = VectorFunctionSpace(mesh, "Lagrange", degree)
# pdb.set_trace()

# Define Dirichlet boundary condition at top boundary
bc2 = DirichletBC(V_phi, 0.0, boundaries, 2)
#Take a look at the solution:
dx = Measure('dx', mesh)
ds=Measure('ds', domain=mesh, subdomain_data=boundaries)
area_gamma1 = dolfin.assemble(1*ds(4))
leitfaehigkeit = 1.43e-7

# Permeability (should be 6.3e-3)
mu = 1e-5

# Defining the trial function
phi = TrialFunction(V_phi)

# Defining the test function
v = TestFunction(V_phi)

u_h1 = Function(V_phi, name = "Displacement")

for i in range(0,10):
u_h = prob_pde.solverIntern(a,v,dx,ds,leitfaehigkeit,area_gamma1,sol,f,u_h1,bc2,i)
p = plot(u_h,title='Finite element solution at t_' + str(i))
ax = prob_pde.achsenSetzen()
plot(mesh)
cbar = plt.colorbar(p)
cbar.set_label('elektrisches Potential [V]', rotation=270)
#------------------------
#plt.savefig('FiniteESol' + str(i)+ '.png')
plt.savefig("FiniteESol_gmres" + str(i)+ '.png')
plt.clf()

# Find the electric field by taking the negative of the
# gradient of the electric potential. Project this onto
# the function space for evaluation. (L_2 projection?)
plt.figure(5)
p1 = plot(elec_field)
cbar = plt.colorbar(p1)
cbar.set_label('elektrisches Potential [V]', rotation=270)
plt.show()
plt.savefig('solution electric_field' + str(i) + '.png')
plt.clf()

prob_pde.energyfunctional(phi,v,ds,leitfaehigkeit,area_gamma1,mesh,u_h,f,sol)

def problem2(self, sol):

f = Constant(0.0)
u = prob_pde.solver(f,1,sol)

print("Test ENDE")

if __name__ == "__main__":

# Initialisierung einer Instanz der Klasse
prob = ode(l_1=5.0,
l_2=1000.0,
k=0.5,
r_2=20.0,
c_1=8.0,
u_omega=160)

prob_pde = pde(l_1=5.0,
l_2=1000.0,
k=0.5,
r_2=20.0,
c_1=8.0,
u_omega=160)

sol = prob.problem()
prob_pde.problem2(sol)


Maybe it isnât that bad if you can see everything.

Please reduce your code by removing everything not needed to reproduce the error. This includes plotting etc.

Note that you should try to remove a variable and see if the error is still there, and keep removing variables and reducing complexity of the code until it is a smaller example, such as:

Hello dokken,

here is a minimal working example. I just replaced the Array z with some random numbers.

import numpy as np
import matplotlib.pylab as plt
import math
from numpy.linalg import solve, norm
from dolfin import *
from mshr import *
import pdb
import moola

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True

z = np.random.rand(4,1000)

domain_vertices = [Point(2.25, 1.75),
Point(1.0, 3.0),
Point(-1.0, 1.0),
Point(4.0, -4.0),
Point(5.0, -3.0),
Point(2.0, 0.0),
Point(0.5, 0.0)]

domain = Polygon(domain_vertices)
mesh = generate_mesh(domain,20)

class Left(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (between(x[1], (-0.999, 0.999)) and between(x[0], (3.999, -3.999)))

class Right(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and ((between(x[1], (2.25, 1.75)) and  between(x[0], (1.001, 2.999))) and (between(x[1], (0.5, 0.0)) and between(x[0], (2.249, 1.749))) and (between(x[1], (1.999, 0.0)) and between(x[0], (0.501, 0.0))) and (between(x[1], (4.999, -2.999)) and between(x[0], (2.0, 0.0))))

class Bottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (between(x[1], (5.0, -3.0)) and between(x[0], (4.0, -4.0)))

class Top(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (between(x[1], (1.0, 3.0)) and between(x[0], (-1.0, 1.0)))

left = Left()
top = Top()
right = Right()
bottom = Bottom()

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)

boundaries.set_all(4)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

degree = 2
V_phi = FunctionSpace(mesh, "CG", degree)
bc2 = DirichletBC(V_phi, 0.0, boundaries, 2)
dx = Measure('dx', mesh)
ds=Measure('ds', domain=mesh, subdomain_data=boundaries)
area_gamma1 = dolfin.assemble(1*ds(4))
leitfaehigkeit = 1.43e-7

phi = TrialFunction(V_phi)
v = TestFunction(V_phi)

u_h1 = Function(V_phi, name = "Displacement")

f = Constant(0)
for i in range(1,10):
L = f*v*dx + (leitfaehigkeit*(1/area_gamma1)*(-z[3,i]))*v*ds(4)
dolfin.solve(a == L, u_h1, bc2)

x = dolfin.SpatialCoordinate(mesh)
W = FunctionSpace(mesh, "DG", 0)
f = interpolate(Expression("x[0]+x[1]", degree=1), W)
control = Control(f)

rf = ReducedFunctional(J, control)

problem = MoolaOptimizationProblem(rf)
f_moola = moola.DolfinPrimalVector(f)
solver2 = moola.NewtonCG(problem, f_moola, options={'gtol': 1e-9,
'maxiter': 20,
'display': 3,
'ncg_hesstol': 0})

sol2 = solver2.solve()
f_opt = sol['control'].data


I am getting the following error:

 func_value = self.scale * self.functional.block_variable.checkpoint
AttributeError: 'dolfin.cpp.la.Vector' object has no attribute 'block_variable'


Greetings

So I do not understand why f is your control, as f is not part of your functional:

You are also mixing dolfin commands and dolfin_adjoint commands, i.e.

this solve will not be overloaded, as you are explicitly using the dolfin.solve function.

Hello Mr. Dokken,

Iâve attached a new MWE to ask you to look over it again please.

Of course you were right about the f and the solve. Now i am doing the following:

import numpy as np
import matplotlib.pylab as plt
import math
from numpy.linalg import solve, norm
from dolfin import *
from mshr import *
import pdb
import moola
from ufl import nabla_div
import dijitso
import pkgconfig
import scipy.io
import sympy as sy

# Form compiler options
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True

class ode(object):

# Konstruktor
def __init__(self, l_1:int, l_2:int, k:float, r_2:float, c_1:float, u_omega:float):
self.l_1 = l_1
self.l_2 = l_2
self.k = k
self.r_2 = r_2
self.c_1 = c_1
self.u_omega = u_omega

def rhs(self,q1,p1,q2,p2,l_12):
f = [p1,
(-self.l_2/((self.l_1*self.l_2-l_12**2)*self.c_1))*q1+(l_12/((self.l_1*self.l_2-l_12**2)*self.c_1))*q2,
p2,
((l_12*self.r_2)/((self.l_1*self.l_2-l_12**2)))*p1-((self.l_1*self.r_2)/(self.l_1*self.l_2-l_12**2))*p2]
return f

def semiImpRK(self,A,q1_i,p1_i,q2_i,p2_i,i,h):
# k1
l_12 = self.k*np.sqrt(self.l_1*self.l_2)
b1 = prob.rhs(q1_i,p1_i,q2_i,p2_i,l_12)
k1 = np.linalg.solve(A,b1)

# k2
b2 = prob.rhs(q1_i+0.5*h*k1[0],p1_i+0.5*h*k1[1],q2_i+0.5*h*k1[2],p2_i+0.5*h*k1[3],l_12)
k2 = np.linalg.solve(A,b2)

# yip1
yip1 = np.zeros(4)
yip1[0] = q1_i + h*k2[0]
yip1[1] = p1_i + h*k2[1]
yip1[2] = q2_i + h*k2[2]
yip1[3] = p2_i + h*k2[3]

return yip1[0], yip1[1], yip1[2], yip1[3]

def problem(self):

l_12 = self.k*np.sqrt(self.l_1*self.l_2)

jac = np.array([[0, 1, 0, 0],
[-self.l_2/((self.l_1*self.l_2-l_12**2)*self.c_1), 0, l_12/((self.l_1*self.l_2-l_12**2)*self.c_1), 0],
[0, 0, 0, 1],
[0, (l_12*self.r_2)/(self.l_1*self.l_2-l_12**2), 0, -(self.l_1*self.r_2)/(self.l_1*self.l_2-l_12**2)]])

t0 = 0.0
T = 100 # Maximale SchweiĂzeit 15ms
N = 1000
h = T/float(N)

# Die Jacobimatrix muss nur einmal berechnet werden, deshalb liegt diese Berechnung auĂerhalb der Schleife
n = 4
a = 1.0/(2.0+math.sqrt(2.0))
I = np.identity(n)
A = I - a*h*jac

# Anfangswerte
z = np.zeros((4,N))
z[0,0] = self.c_1 * self.u_omega
z[1,0] = 0.0
z[2,0] = 0.0
z[3,0] = 0.0

t = np.zeros(N+1)
sol = np.zeros((4,N+1))

# t = 0
sol[0,0], sol[0,1], sol[0,2], sol[0,3] = prob.semiImpRK(A,z[0,0],z[1,0],z[2,0],z[3,0],0,h)
# t > 0
for i in range(1,N+1):
q1_i, p1_i, q2_i, p2_i = sol[0,i-1], sol[1,i-1], sol[2,i-1], sol[3,i-1]
sol[0,i], sol[1,i], sol[2,i], sol[3,i] = prob.semiImpRK(A,q1_i,p1_i,q2_i,p2_i,i,h)
t[i] = t[i-1] + h

return sol

class pde(ode):

def __init__(self,l_1,l_2,k,r_2,c_1,u_omega):
super().__init__(l_1,l_2,k,r_2,c_1,u_omega)

def solver(self,f,degree,sol):

mesh = UnitSquareMesh(10, 10)

class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)

class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0)

class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)

class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 1.0)

# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()

# Initialize mesh function for boundary domains
boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)

boundary_subdomains.set_all(4)
left.mark(boundary_subdomains, 1)
top.mark(boundary_subdomains, 2)
right.mark(boundary_subdomains, 3)
bottom.mark(boundary_subdomains, 4)

dss=Measure('ds', domain=mesh, subdomain_data=boundary_subdomains)
area_gamma1 = dolfin.assemble(1*dss(4))
# specific electric conductivity of carbon steel at 20 C
leitfaehigkeit = 1.43e-7

V_phi = FunctionSpace(mesh, "CG", degree)

# Dirichlet boundary condition for the equation of potential
c = Constant(0.0)
bc2 = DirichletBC(V_phi, c, boundary_subdomains, 2)
dx = Measure('dx', mesh)

# Test and trial functions for equation of potential
phi = TrialFunction(V_phi)
v = TestFunction(V_phi)

f = Constant(0.0)

# time parameters
T = 15            # final time should be 15ms
num_steps = 100     # number of time steps
dt = T / num_steps # time step size
t = 0.0

# bilinearform for the equation of potentialĂ

# repleasment functions for the equation of potential in the discrete space
u = Function(V_phi)
u_n = Function(V_phi)

for i in range(num_steps):

# equation of potential not time dependet-------------------------
L = f*v*dx + (leitfaehigkeit*(1/area_gamma1)*(-sol[3,i]))*v*dss(4)
# Compute solution
solve(a==L, u, bc2)

# Update previous solution
u_n.assign(u)

t+=dt

f_con = (leitfaehigkeit*(1/area_gamma1)*(-sol[3,i]))*v*dss(4)
control = Control(f_con)

def problem2(self, sol):

f = Constant(0.0)
degree = 1
u = prob_pde.solver(f,degree,sol)

if __name__ == "__main__":

# Initialisierung einer Instanz der Klasse
prob = ode(l_1=5.0,
l_2=1000.0,
k=0.5,
r_2=20.0,
c_1=8.0,
u_omega=160)

prob_pde = pde(l_1=5.0,
l_2=1000.0,
k=0.5,
r_2=20.0,
c_1=8.0,
u_omega=160)

sol = prob.problem()
prob_pde.problem2(sol)


I am very sorry that I did not remove the ODE class from the MWE, but I think that this is important for my problem.

The error is

AttributeError: 'Form' object has no attribute 'block_variable'.


I want my control to be the voltage drop over the boundary.
May you tell me if there is any possibility for this?

Thank you!

You cannot use a whole ufl.Form as a control variable. You need to choose a variable (or create a variable your functional is dependent on, that is your Control).

1 Like

Do you mean something like the coordinates of my boundary?
But what if i do not know how this boundary is going to look like after my process?

Please formulate the mathematical problem you would like to solve, with functional and PDEs (in strong form).
as of your problem above, you are saying that you would like the integral:

to be your control. You cannot have an integral as a control, as it depends on several variables.

Oh wow. Thank you.

I do not mean that you should implement it. I would simply like to see that mathematical problem you are trying to solve.

Ok.

-\gamma\Delta\phi = 0 \text{ in } \Omega\\ \phi = 0 \text{ on } \Gamma_1\\ -\gamma \frac{\partial \phi}{\partial n} = 0 \text{ on } \Gamma_2 \text{ and } \Gamma_3\\ -\gamma \frac{\partial \phi}{\partial n} = -\frac{1}{|\Gamma_4|}I_2(t) \text{ on } \Gamma_4

whereby $$\phi$$ is the potential, n is the outer normal and $$I_2(t)$$ is a current.
Defining V as

V:=\{v \in H^1|v=0 \text { on } \Gamma_1\}

and multiplying with the test function v and doing the integration yields for the weak formulation

\int_{\Omega} \gamma \nabla \phi \nabla v dx = \frac{1}{|\Gamma_4|}I_2(t)\int_{\Gamma_4} v(x) dx.

Now we are applying a little trick by saying v should be $$\phi$$ which yields now

\int_{\Omega} \gamma |\nabla \phi|^2 dx = \frac{1}{|\Gamma_4|}I_2(t)\int_{\Gamma_4} \phi(x) dx.

This would give an optimal control problem with the energy-functional like

\min_{\phi \in V} J(\phi) := \min_{\phi \in V}(\frac{1}{2}\int_{\Omega} \gamma |\nabla \phi |^2 dx - \frac{1}{|\Gamma_4|}I_2(t)\int_{\Gamma_4} \phi(x) dx)

under the constrains of the potential equation.

I attached a picture of my area and hope that this is understandable (I really canât draw).

The marked area is to be melted. I was told to only look at the equation of potential and not at the PDEs for the stress and the heat. In an isotropic material J denotes the density of the electric current, E = - \nabla \phi and \phi is the electric potential. With Ohmâs law we have

J = - \gamma(u)\nabla \phi

where u is the temperature.

For the optimal shape of the bulk i should maybe look at the map to the x-axes.

So with this PDE, what is it that you want to use as a control parameter? The current optimal control problem only states \phi as a Control variable, i.e. a classical energy minimization problem.

Yes, this is true.

Indeed i want to minimize the energy, which depends on

\phi.

I think this would be in my code

control = Control(u)


But i do not want to minimize the given energy from the capacitor and at the same time i want to have an optimal shape of my bulk.

This is not the intention of dolfin-adjoint, i would just solve it with the weak form:

Thank you Mr. Dokken.