So I am trying to port this notebook Topology Optmisation to Fenicsx.
The code I managed to put together so far is below, for your convenience I place the old, commented out Fenics lines as well as the Fenicsx “translation”.
######### DOLFIN ##########
#%matplotlib notebook
#from dolfin import *
#import matplotlib.pyplot as plt
#import numpy as np
########## DOLFIN ##########
########## DOLFINX ##############
%matplotlib notebook
from dolfinx import *
import matplotlib.pyplot as plt
import numpy as np
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from mpi4py import MPI
from dolfinx.fem import FunctionSpace
from dolfinx.fem import VectorFunctionSpace
from dolfinx import fem
from dolfinx import mesh
########## DOLFINX #############
# Algorithmic parameters
niternp = 20 # number of non-penalized iterations
niter = 80 # total number of iterations
pmax = 4 # maximum SIMP exponent
exponent_update_frequency = 4 # minimum number of steps between exponent update
tol_mass = 1e-4 # tolerance on mass when finding Lagrange multiplier
thetamin = 0.001 # minimum density modeling void
########## DOLFIN ##########
# Problem parameters
#thetamoy = 0.4 # target average material density
#E = Constant(1)
#nu = Constant(0.3)
#lamda = E*nu/(1+nu)/(1-2*nu)
#mu = E/(2*(1+nu))
#f = Constant((0, -1)) # vertical downwards force
########## DOLFIN ##########
########## DOLFINX #############
# Problem parameters
thetamoy = 0.4 # target average material density
E = fem.Constant(domain, ScalarType(1))
nu = fem.Constant(domain, ScalarType(0.3))
lamda = E*nu/(1+nu)/(1-2*nu)
mu = E/(2*(1+nu))
f = fem.Constant(domain, ScalarType([0,-1])) # vertical downwards force
########## DOLFINX #############
########## DOLFIN ##########
# Mesh
#mesh = RectangleMesh(Point(-2, 0), Point(2, 1), 50, 30, "crossed")
# Boundaries
#def left(x, on_boundary):
# return near(x[0], -2) and on_boundary
#def load(x, on_boundary):
# return near(x[0], 2) and near(x[1], 0.5, 0.05)
#facets = MeshFunction("size_t", mesh, 1)
#AutoSubDomain(load).mark(facets, 1)
#ds = Measure("ds", subdomain_data=facets)
########## DOLFIN ##########
########## DOLFINX #############
L = 5
domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0], [L, 1, 1]], [20, 5, 5], mesh.CellType.hexahedron)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
########## DOLFINX #############
# Boundaries
def left(x, on_boundary):
return near(x[0], -2) and on_boundary
def load(x, on_boundary):
return near(x[0], 2) and near(x[1], 0.5, 0.05)
facets = MeshFunction("size_t", mesh, 1)
AutoSubDomain(load).mark(facets, 1)
ds = Measure("ds", subdomain_data=facets)
########## DOLFIN ##########
# Function space for density field
#V0 = FunctionSpace(mesh, "DG", 0)
# Function space for displacement
#V2 = VectorFunctionSpace(mesh, "CG", 2)
# Fixed boundary condtions
#bc = DirichletBC(V2, Constant((0, 0)), left)
########## DOLFIN ##########
########## DOLFINX #############
# Function space for density field
V0 = FunctionSpace(domain, ("DG", 0))
# Function space for displacement
V2 = VectorFunctionSpace(domain, ("CG", 1))
# Fixed boundary condtions
#bc = fem.Constant(domain, ScalarType([0,0]), left)
########## DOLFINX #############
########## DOLFIN ##########
p = Constant(1) # SIMP penalty exponent
exponent_counter = 0 # exponent update counter
lagrange = Constant(1) # Lagrange multiplier for volume constraint
thetaold = Function(V0, name="Density")
thetaold.interpolate(Constant(thetamoy))
coeff = thetaold**p
theta = Function(V0)
volume = assemble(Constant(1.)*dx(domain=mesh))
avg_density_0 = assemble(thetaold*dx)/volume # initial average density
avg_density = 0.
########## DOLFIN ##########
#p = Constant(1) # SIMP penalty exponent
#exponent_counter = 0 # exponent update counter
#lagrange = Constant(1) # Lagrange multiplier for volume constraint
#thetaold = Function(V0, name="Density")
#thetaold.interpolate(Constant(thetamoy))
#coeff = thetaold**p
#theta = Function(V0)
########## DOLFIN ##########
########## DOLFINX #############
p = fem.Constant(domain, ScalarType(1)) # SIMP penalty exponent
exponent_counter = 0 # exponent update counter
lagrange = fem.Constant(domain, ScalarType(1))# Lagrange multiplier for volume constraint
thetaold = fem.Function(V0, name="Density")
thetaold.interpolate(fem.Constant(domain, ScalarType(thetamoy)))
coeff = thetaold**p
theta = Function(V0)
########## DOLFINX #############
The problem lies in
thetaold.interpolate(fem.Constant(domain, ScalarType(thetamoy)))
I get an error message
Couldn't map 'c_5' to a float, returning ufl object without evaluation.
like in this post, couldnt-map-float-returning-object-without-evaluation/
I tried the suggested solution
thetaold.x.set(thetamoy)
but I am not so sure it is correct, what does for example .x stand for?
The next step where I have issues is the Fenics lines
volume = assemble(Constant(1.)*dx(domain=mesh))
avg_density_0 = assemble(thetaold*dx)/volume # initial average density
avg_density = 0.
I tried to rewrite in Fenicsx as
################# DOLFINX#############
volume = fem.assemble(fem.Constant(domain, ScalarType(1))*dx(domain=domain))
avg_density_0 = fem.assemble(thetaold*dx)/volume # initial average density
avg_density = 0.
################# DOLFINX#############
but I get the error message
name 'dx' is not defined
and I am at a loss. In the legacy code ds = Measure("ds", subdomain_data=facets) is defined, but not dx
Also in the Fenicsx Poisson tutorial,
Implementation — FEniCSx tutorialFenicsx Poisson tutorial
dx is not defined.
Any hint would be most appreciated, I am really keen to port this example to Fenicsx, thanks