Problem with functional - dolfin-adjoint - FEniCS Project

相关文章推荐
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
 newbyJulia:
\int_{\Omega) |\nabla \phi|^2 dx = - \frac{1}{|\Gamma_4} * I_2 * \int_{\Omega_4} \phi ds
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
# 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,
        ((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):        
        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)
        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)
        J = assemble(0.5 * dolfin.dot(dolfin.grad(u_h), dolfin.grad(v))*dx - (leitfaehigkeit*(1/area_gamma1)*((-sol[3,1]))*v*ds(4)))
        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.dolfin_adjoint.solve()
       # 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")
        a = dot(grad(phi), grad(v))*dx
        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)
            #--------PLOT Gradienten
            plot(-grad(u_h))
            #------------------------
            #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?)
            elec_field = dolfin.project(-grad(u_h))
            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:

    Read before posting: How do I get my question answered? 
    The following is a list of well constructed questions adhering to the guidelines above: 
Getting cell dofs for vector valued function space in dolfinx
Problem with 3D nedelec element 
What is the “default” linear solver/preconditioner? 
Setting Dirichlet boundary condition/constraint on 3D edges 
              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 dolfin_adjoint 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)
a = dot(grad(phi), grad(v))*dx
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)
J = assemble(0.5 * dolfin.dot(dolfin.grad(u_h1), dolfin.grad(v))*dx) - float(leitfaehigkeit*(1/area_gamma1)*(-z[3,1]))*assemble(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()
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
f = interpolate(Expression("x[0]+x[1]", degree=1), W)
J = assemble(0.5 * dolfin.dot(dolfin.grad(u_h1), dolfin.grad(v))*dx) - float(leitfaehigkeit*(1/area_gamma1)*(-z[3,1]))*assemble(v*ds(4))
control = Control(f)
rf = ReducedFunctional(J, control)
You are also mixing dolfin commands and dolfin_adjoint commands, i.e.
 newbyJulia:
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)
this solve will not be overloaded, as you are explicitly using the dolfin.solve function.
              Hello Mr. Dokken,
thank you for all the good comments about my problem.

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 dolfin_adjoint 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,
        ((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):        
        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)
        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ß
        a = dot(grad(phi), grad(v))*dx
        # 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
            J = assemble(0.5 * dot(grad(phi), grad(v))*dx - (leitfaehigkeit*(1/area_gamma1)*(-sol[3,i]))*v*dss(4))
            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!
              As your example is far from minimal, it is quite hard for me to help you.

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:
 newbyJulia:
(leitfaehigkeit*(1/area_gamma1)*(-sol[3,i]))*v*dss(4)
to be your control. You cannot have an integral as a control, as it depends on several variables.
-\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.
 newbyJulia:
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.
∫Ωγ∇ϕ∇vdx=1|Γ4|I2(t)∫Γ4v(x)dx.