# Unknown ufl objective type FiniteElement

**URL:** <https://fenicsproject.discourse.group/t/unknown-ufl-objective-type-finiteelement/12660>\
**Category:** General\
**Created:** [November 6, 2023, 11:44am UTC](https://fenicsproject.discourse.group/t/unknown-ufl-objective-type-finiteelement/12660 "2023-11-06T11:44:08Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![Mohammed\_Hameed](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mohammed_hameed/32/5604_2.png) [@Mohammed\_Hameed](https://fenicsproject.discourse.group/u/Mohammed_Hameed)\
**Post date:** [November 6, 2023, 11:44am UTC](https://fenicsproject.discourse.group/t/unknown-ufl-objective-type-finiteelement/12660/1 "2023-11-06T11:44:08Z")

</div>

Hello,  
The error is: when I run the code in python, this error appears as you can see in the attached screenshot.

 ![fe](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/e/e4b4ca0a1a5b81aef08c2e08617fe7be530e4b65.png)

This is the whole code:  
‘’'Python  
“”"  
This code is made available under a BSD 3-clause license. See LICENSE for further information.

This script is a minimal working example of how to perform gradient-based topology optimisation  
based on a homogenisation interpolation and GCMMA with FEniCS. This exaple minimises the temporal  
average variance of the spatial average temperature at the heat source boundary by placing phase  
change material (PCM) and highly thermally conductive material (HCM).

To run this you need

- legacy FEniCS (sudo ap), to install in VSC in Linux use following commands:  
sudo apt-get install software-properties-common  
sudo add-apt-repository ppa:fenics-packages/fenics  
sudo apt-get update  
sudo apt-get install fenics
- dolfin adjoint ([dolfin-adjoint about](https://www.dolfin-adjoint.org/en/latest/)), to install in VSC in Linux use following commands:  
pip install git+https://github.com/dolfin-adjoint/pyadjoint.git@2019.1.0  
git clone [Bitbucket](https://bitbucket.org/dolfin-adjoint/pyadjoint.git)
- MMA/GCMMA code (the MMA.py file) from Arjee Deetman ([GitHub - arjendeetman/GCMMA-MMA-Python: Python implementation of the Method of Moving Asymptotes (MMA)](https://github.com/arjendeetman/GCMMA-MMA-Python)), download from github

If you do not have a linux computer available a virtual machine with Linux can be used.  
I have had good experiences with WMware and Oracle VM Virtualbox.

By Mark Christensen  
28-07-23

Disclaimer:  
The author does not guarantee that the code is free from errors.  
Furthermore, the author shall not be liable in any event caused by the use of the program.  
“”"

# #####################################################

# LIBRARIES

# #####################################################

from fenics import \*  
from fenics\_adjoint import \* # Note: It is important that fenics\_adjoint is called after fenics  
import matplotlib.pyplot as plt  
import numpy as np  
import statistics as stat

from mma import gcmmasub,subsolv,kktcheck,asymp,concheck,raaupdate # Import MMA functions

set\_log\_level(LogLevel.ERROR) # Only print message if error message

# #####################################################

# PREPROCESSING

# #####################################################

# ----- PHYSICAL MODEL DATA -----

# Material constants

khcm = 10 # Thermal conductivity og HCM [W/(m_K)]  
cphcm = 1 # Specific heat capacity of HCM [J/(K_kg)]  
rhohcm = 1 # Mass density density of HCM [kg/m³]

kpcm = 1e-3_khcm # Thermal conductivity og PCM [W/(m_K)]  
cppcm = 1 # Specific heat capacity of PCM [J/(K\*kg)]  
rhopcm = 1 # Mass density density of PCM [kg/m³]  
Tmelt = 0.5 # Melting temperature of PCM [K]  
dTmelt = 0.5 # Melting temperture range of PCM [K]  
Lheat = 10 # Latent heat of fusion [J/kg]

# Thermal data

qheat = 1 # Heat production rate in electronics [W]  
w = 1 # Ocsilation frequency of heat source [Hz]  
hconv = 5 # Het atransfer coefficient [W/(m²\*K)]  
Tinf = 0 # Temperature of surroundings [K]  
Tinitial = 0 # Initial temperature [K]

# Temporal data

tfin = 20 # Final time [s]  
num\_steps = 500 # Number of time steps  
dt = tfin/num\_steps # Time step size [s]

# Finite element data

nx = ny = 100 # Number of element in x and y direction  
lx = ly = 1 # Dimension of the 2D heatsink [m]  
lz = 1 # Out of plane thickness of the 2D heatsink [m]  
degree\_phys = 1 # Element order for physical problem

# ----- OPTIMIZATION ELEMENT DATA -----

# Loop

niter = 300 # Max number of optimisation iterations  
tol\_opt = 1E-3 # Convergence tolerance, based on the absolute error

# Problem parameters

volfrac = 0.3 # Maximum allowable volume fraction used in the optimisation

# Filter data

r = 0.01 # Length scale for the Helmholzt PDE filterer

# ----- DEFINING DOMAIN AND FUNCTION SPACES -----

# Create mesh

mesh = RectangleMesh(Point(-lx/2, -ly/2), Point(lx/2, ly/2), nx, ny, “crossed”)

# Function space

V0 = FunctionSpace(mesh, “DG”, 0) # Discontinous function space  
Vphys = FunctionSpace(mesh, “CG”, degree\_phys) # Continous function space

# ----- BOUNDARY CONDITIONS -----

# Define the Heat source and convection boundaries were the Neumann BC are a aplied.

class HeatSourceBoundary(SubDomain):  
def inside(self, x, on\_boundary):  
tol\_bc = 1E-14  
return on\_boundary and near(x[1], -ly/2, tol\_bc) and x[0]+lx/4\>=-tol\_bc and x[0]-lx/4\<=tol\_bc  
HS\_boundary = HeatSourceBoundary()

# Applying Neuman boundary condition to simulate the heat loss due to convection

class ConvectionBoundary(SubDomain):  
def inside(self, x, on\_boundary):  
tol\_bc = 1E-14  
return on\_boundary and near(x[1], ly/2, tol\_bc)  
C\_boundary = ConvectionBoundary()

# Marking the facets (boundaries) that coresspond to HeatSourceBoundary and ConvectionBoundary

boundaries = MeshFunction(‘size\_t’, mesh, mesh.topology().dim()-1) # store all facets in array  
HS\_boundary.mark(boundaries, 1) # Mark all facest that are on the HS\_boundary with 1  
C\_boundary.mark(boundaries, 2) # Mark all facest that are on the c\_boundary with 2

# Denfine the line intrgrals

dsHS = ds(subdomain\_id=1, subdomain\_data=boundaries)  
dsC = ds(subdomain\_id=2, subdomain\_data=boundaries)

# Define volume and area for later use

voltot = assemble(Constant(lz)\*dx(domain=mesh)) # Compute total volume Of design domain  
AHS = assemble(Constant(lz)\*dsHS(domain=mesh)) # Compute crossectionl area of heat source

# #####################################################

# FUNCTIONS

# #####################################################

# Apparent heat capacity method

def cppcmmod(T\_,cppcm\_, Lheat\_):  
“”"  
Incoorporates the latent heat of fusion into the heat capacity with  
the aparent heat capacity method with heaviside step functions.  
This function assumes no change in spefic heat capacity due to phase change.

```
Inputs: T_ is the temperture field, cppcm_ is the heat capacity without 
phase change, Lheat_ is the latent heat of fusion.
Outputs: Field of the modified heat capacity
"""
k_H = 25 # Steepness of Heviside step function 
return cppcm_ + Lheat_ / (dTmelt) * \
    (
        1 / (1 + exp(-2 * k_H * (T_ - (Tmelt - dTmelt/2))))
        - (1 / (1 + exp(-2 * k_H * (T_ - (Tmelt + dTmelt/2)))))
    )

```

def PDEfilter(u\_, r\_, mesh\_):  
“”"  
Applies ering to a field by solving af PDEfilter type PDE  
based on the work of B.S. Lazarov and O. Sigmund DOI: 10.1002/nme.3072

```
Inputs: u_ is the unfiltered field, r_ is the characteristic length
Outputs: filtered field u_tilde
"""
Vold = u_.function_space() # Saving function space from u_
VCG = FunctionSpace(mesh_, "CG", 1) # Defining continous function space

# Defining variational problem
u_tilde = TrialFunction(VCG)
vf = TestFunction(VCG)
# HELMHOLTZ-TYPE PDE (u_ is projected to a continouos function space)
F = u_tilde*vf*dx + r_**2*dot(grad(u_tilde), grad(vf))*dx - project(u_, VCG)*vf*dx(mesh_)
# a: unknowns, L: knowns
a, L = lhs(F), rhs(F)

u_temp = Function(VCG)
u_tilde = Function(Vold)

# Compute solution and project onto inital function space
solve(a == L, u_temp)
u_tilde.assign(project(u_temp, Vold))

return u_tilde

```

def forward(rho\_tilde\_):  
“”"  
Solves the forward problem (solving PDE) using homogenisation

```
Inputs: rho_tilde is the filtered material density variable field
Outputs: Variables used for computing the objective function
"""
# ----- DEFINING PROBLEM -----
# Material interpolations
def rhocppcmphys(rho_,T_):
    return (rhohcm*cphcm)*(rho_)+(1-rho_)*(rhopcm*cppcmmod(T_,cppcm, Lheat))

def kphys(rho_):
    a = 1-(1-rho_)**0.5
    return 1/(a/khcm + (1-a)/(kpcm*(1-a)+khcm*a))

# Defining the flucturating heat source
qelec = Expression("qheat/(AHS)*(1+sin(2*DOLFIN_PI*w*t))",
                    qheat=qheat, AHS=AHS, lz=lz, w=w, t=0, degree=0)

# Define initial values
T_ini = Constant(Tinitial)
T_n = interpolate(T_ini, Vphys)

# Defining the trial and test functions
T = TrialFunction(Vphys)
v = TestFunction(Vphys)

# Weak form rho_k=0 -> PCM, rho_k = 1 -> HCM
F = rhocppcmphys(rho_tilde_,T_n)*(T-T_n)/dt*v*dx \
    + kphys(rho_tilde_)*dot(grad(T),grad(v))*dx \
    - qelec*v*dsHS \
    + hconv*(T-Tinf)*v*dsC
# a: unknowns, L: knowns
a, L = lhs(F), rhs(F)

# ------ SOLVING THE PROBLEM -----

# Initialise the field for the Temperature
T = Function(Vphys) 

# Set initial condistions
T_n.assign(interpolate(T_ini, Vphys))
T.assign(interpolate(T_ini, Vphys))
t = 0
qelec.t = t

# Save the average temperture at the heat source
T_elec_array = [assemble(T*dsHS)/AHS]   

# ------ Solving Physics -----
for timestep in range(num_steps):
    t += float(dt) # Update time
    qelec.t = t # Update heat source
    solve(a == L, T) # Solve problem for current time step with lieaner solver

    T_n.assign(T) # Update previous temperature
    T_elec_array.append(assemble(T*dsHS)/AHS)

# Compute the variance of the average temperture at the heat source over time
T_time_mean = stat.mean(T_elec_array) # Average over time
T_var = []

for timesteps in range(len(T_elec_array)):
    T_var.append((T_elec_array[timesteps] - T_time_mean)**2)

return T_var

```

def Optimization(rho\_,rho\_tilde\_, rf\_f0\_, rf\_h\_,f0\_history\_,Mnd\_history\_):  
“”"  
Performs gradient based topology optimazation using the Methode of moving assymptotes  
based on a MMA and GCMMA script by Krister Svanberg ([https://people.kth.se/~krille/mmagcmma.pdf](https://people.kth.se/~krille/mmagcmma.pdf))  
which is translated to Python Arjeen Deetman ([GitHub - arjendeetman/GCMMA-MMA-Python: Python implementation of the Method of Moving Asymptotes (MMA)](https://github.com/arjendeetman/GCMMA-MMA-Python))

```
Inputs: rho is the density variable field, rho_tilde is the filtered density variable field, 
rf_f0 is the reduced functional for f0, rf_h is the reduced functional for h, f0_history is 
a empty list for storing the history of f0, and Mnd_history is a empty list for storing the
history of h.
Outputs: rho and rho_tilde at the final design and f0_history, Mnd_history
"""
with pyadjoint.stop_annotating() as _: # Stop saving opterations with dolfin-adjoint
    
    # ------ INITIALASING LOOP -----
    # Preallocate rho and rho_tilde for each iteration
    rho_k = rho_
    rho_tilde_k = rho_tilde_

    # Initialize GCMMA parameters
    m = 1 # Number of constraints
    n = rho.vector()[:].size # Number of design variables
    xval = rho.vector()[:].reshape(-1, 1) # Initial design variables
    epsimin = 0.0000001
    eeen = np.ones((n,1))
    eeem = np.ones((m,1))
    zeron = np.zeros((n,1))
    zerom = np.zeros((m,1))
    xold1 = xval.copy()
    xold2 = xval.copy()
    xmin = 0*eeen.copy()
    xmax = 1*eeen.copy()
    low = xmin.copy()
    upp = xmax.copy()
    c = 1000*eeem
    d = eeem.copy()
    a0 = 1
    a = zerom.copy()
    raa0 = 0.01
    raa = 0.01*eeem
    raa0eps = 0.000001
    raaeps = 0.000001*eeem
    outeriter = 0	

    # Initilize preliminary design 
    f0_k = rf_f0_(rho_k)
    h_k = rf_h_(rho_k)
    Mnd = assemble(4*rho_tilde_k*(1-rho_tilde_k)*dx)/(lx*ly)*100
    # Calculate Sensitivities
    df0drho = Function(V0, name="Object function Sensitivity")
    dhdrho = Function(V0, name="Constraint function Sensitivity")
    df0drho.assign(rf_f0_.derivative())                     
    dhdrho.assign(rf_h_.derivative())

    # Scaling functions and storing them in matrices for the MMA solver
    xval = rho_k.vector()[:].reshape(-1, 1) # Saving design variable in xval
    scale = np.abs(f0_k/10)
    f0val = 1 + f0_k/scale # Scaling f0 at xval
    df0dx = df0drho.vector()[:].reshape(-1, 1)/scale # nx1 matrix of sensitivities of f0 at xval
    fval = h_k # Array of constraint functions at xval
    dfdx = dhdrho.vector()[:].reshape(m, n)  

    # Initialize the outer iteration and convergence counter
    outit = 0
    cc = 0 

    # ----- OPTIMIZATION LOOP -----
    while (outit < niter):  
        outit += 1
        outeriter += 1

        if outit > 1: # update prev values for convergence check
            f0_kprev = f0_k
            Mndprev = Mnd
        
        # The parameters low, upp, raa0 and raa are calculated:
        low,upp,raa0,raa= \
            asymp(outeriter,n,xval,xold1,xold2,xmin,xmax,low,upp,raa0,raa,raa0eps,raaeps,df0dx,dfdx)
        # The MMA subproblem is solved at the point xval:
        xmma,ymma,zmma,lam,xsi,eta,mu,zet,s,f0app,fapp= \
            gcmmasub(m,n,iter,epsimin,xval,xmin,xmax,low,upp,raa0,raa,f0val,df0dx,fval,dfdx,a0,a,c,d)

        # Evaluate the functionals with the new deisgn
        rho_k.vector()[:] = xmma.flatten() # Overwrite rho with the new rho from the MMA solver
        f0_knew = rf_f0_(rho_k)
        h_knew = rf_h_(rho_k)
        f0valnew = np.array([1 + f0_knew/scale]) # Scaling f0 at xval
        fvalnew = np.array([h_knew])

        # It is checked if the approximations are conservative:
        conserv = concheck(m,epsimin,f0app,f0valnew,fapp,fvalnew)

        # While the approximations are non-conservative (conserv=0), repeated inner iterations are made:
        innerit = 0
        if conserv == 0:
            while conserv == 0 and innerit < 2:
                innerit += 1
                # New values on the parameters raa0 and raa are calculated:
                raa0,raa = raaupdate(xmma,xval,xmin,xmax,low,upp,f0valnew,fvalnew,f0app,fapp,raa0, \
                    raa,raa0eps,raaeps,epsimin)
                # The GCMMA subproblem is solved with these new raa0 and raa:
                xmma,ymma,zmma,lam,xsi,eta,mu,zet,s,f0app,fapp = gcmmasub(m,n,iter,epsimin,xval,xmin, \
                    xmax,low,upp,raa0,raa,f0val,df0dx,fval,dfdx,a0,a,c,d)
                
                # Evaluate the functionals with the new deisgn
                rho_k.vector()[:] = xmma.flatten() # Overwrite rho with the new rho from the MMA solver
                f0_knew = rf_f0_(rho_k)
                h_knew = rf_h_(rho_k)
                f0valnew = np.array([1 + f0_knew/scale]) # Scaling f0 at xval
                fvalnew = np.array([h_knew])
                # It is checked if the approximations have become conservative:
                conserv = concheck(m,epsimin,f0app,f0valnew,fapp,fvalnew)

        # Updat vectors:
        xold2 = xold1.copy()
        xold1 = xval.copy()
        xval = xmma.copy()

        # Re-calculate function values and gradients of the objective and constraints functions
        rho_k.vector()[:] = xmma.flatten() # Overwrite rho with the new rho from the MMA solver
        rho_tilde_k.assign(PDEfilter(rho_k, r, mesh))
        f0_k = rf_f0_(rho_k)
        h_k = rf_h_(rho_k)
        df0drho.assign(rf_f0_.derivative())                     
        dhdrho.assign(rf_h_.derivative())

        # Scaling functions and storing them in matrices for the MMA solver
        f0val = 1 + f0_k/scale # Scaling f0 at xval
        df0dx = df0drho.vector()[:].reshape(-1, 1)/scale # nx1 matrix of sensitivities of f0 at xval
        fval = h_k # Array of constraint functions at xval
        dfdx = dhdrho.vector()[:].reshape(m, n)  

        # Store f0
        f0_history_.append(f0_k) # Save f0 in f0 history

        Mnd = assemble(4*rho_tilde_k*(1-rho_tilde_k)*dx)/(lx*ly)*100
        Mnd_history_.append(Mnd)

        # Chect for convergence
        if outit > 1:
             # Check if converged
            f0ichange = abs((f0_k-f0_kprev)/(f0_k))
            Mndchange = abs((Mnd-Mndprev)/(Mnd) )

            if f0ichange <= tol_opt and Mndchange <= tol_opt:
                cc += 1
                if cc > 3:
                    print("Tolerance Reached!!")
                    break 
            else:
                cc = 0
            print("Iteration " + str(outeriter) + "." + str(innerit) + ": f0 = " + str(f'{f0_k:.3f}') + ", Mnd = " + str(f'{Mnd:.2f}')+", f0ichange = " + str("{:.2e}".format(f0ichange))+", Mndchange = " + str("{:.2e}".format(Mndchange)))               

        
        # ----- PLOTS While Solving -----               
        plt.figure(1) # plot of current design
        plt.clf()
        plt.colorbar(plot(rho_tilde_k, cmap="Greys", vmin=0, vmax=1))
        plt.xlim([-lx/2, lx/2])
        plt.ylim([-ly/2, ly/2])
        plt.title("rho_tilde at iter " + str(outit)+": f0 = "+ str(f'{f0_history_[-1]:.2f}')+ ", Mnd = " + str(f'{Mnd_history[-1]:.2f}'))
        plt.pause(0.05)

        plt.figure(2) # plot of f0 and Mnd history
        plt.clf()
        plt.plot(np.arange(1, len(f0_history_)+1)-1, f0_history_, color="red")
        ax = plt.gca()
        plt.xlabel("Number of iterations")
        plt.ylabel("Objective", color="red")
        plt.tick_params(axis="y", which="both", labelcolor="red")
        plt.title("Convergence history")
        ax2 = ax.twinx()
        ax2.plot(np.arange(1, len(Mnd_history_)+1)-1, Mnd_history_, color="blue")
        ax2.set_ylabel("Measure of non-discreteness", color="blue")
        ax2.tick_params(axis="y", which="both", labelcolor="blue")
        plt.pause(0.05)

return rho_k,rho_tilde_k, f0_history_, Mnd_history_

```

if **name** == ‘ **main** ’:  
# ----- INITIALIZE DENSITIES -----  
rho = interpolate(Constant(volfrac),V0) # Define initial design  
rho\_tilde = PDEfilter(rho, r, mesh) # Apply filter

```
T_var = forward(rho_tilde) # Solve forward problem (physical model)

# Objective and Constraint function
h = assemble(rho_tilde*lz*dx)/(volfrac*voltot) - 1
f0 = sum(T_var)/len(T_var)

# Reduced functionals
control = Control(rho)  
rf_f0 = ReducedFunctional(f0, control)  
rf_h = ReducedFunctional(h, control)  

# Initiate f0 and Mnd historry lists
f0_history = []
Mnd_history = []

rho,rho_tilde, f0_history, Mnd_history = Optimization(rho, rho_tilde, rf_f0, rf_h, (f0_history), (Mnd_history))

# ----- PLOTS -----  
plt.figure(1) #The filal design
plt.clf()
plt.colorbar(plot(rho_tilde, cmap="Greys", vmin=0, vmax=1))
plt.xlim([-lx/2, lx/2])
plt.ylim([-ly/2, ly/2])
plt.title("rho_tilde at iter " + str(iter)+": f0 = "+ str(f'{f0_history[-1]:.2f}')+ ", Mnd = "+ str(f'{Mnd_history[-1]:.2f}'))

plt.figure(2) #The f0 and Mnd history
plt.clf()
plt.plot(np.arange(1, len(f0_history)+1)-1, f0_history, color="red")
ax = plt.gca()
plt.xlabel("Number of iterations")
plt.ylabel("Objective", color="red")
plt.tick_params(axis="y", which="both", labelcolor="red")
plt.title("Convergence history")
ax2 = ax.twinx()
ax2.plot(np.arange(1, len(Mnd_history)+1)-1, Mnd_history, color="blue")
ax2.set_ylabel("Measure of non-discreteness", color="blue")
ax2.tick_params(axis="y", which="both", labelcolor="blue")
plt.show()

```

‘’’

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [November 13, 2023, 6:41pm UTC](https://fenicsproject.discourse.group/t/unknown-ufl-objective-type-finiteelement/12660/2 "2023-11-13T18:41:15Z")

</div>

You have a mixture of system wide packages installed with `apt` (in `/usr/lib/petsc`), and manually installed packages in `/home/yourname/.local`. Make a backup of `/home/yourname/.local` just in case there are important files there, then delete `/home/yourname/.local` and try again.