Hello 
I wrote an article to ask a question about a problem that occurred in the process of converting an existing ‘A 55-line topology optimization’ dolfin code to a dolfinx code.
The existing code is structured as follows, so it is being translated without destroying the structure as much as possible.
from dolfin import *
import numpy as np, sklearn.metrics.pairwise as sp
# A 55 LINE TOPOLOGY OPTIMIZATION CODE ---------------------------------
def main(nelx, nely, volfrac, penal, rmin):
sigma = lambda _u: 2.0 * mu * sym(grad(_u)) + lmbda * tr(sym(grad(_u))) * Identity(len(_u))
psi = lambda _u: lmbda / 2 * (tr(sym(grad(_u))) ** 2) + mu * tr(sym(grad(_u)) * sym(grad(_u)))
xdemf = XDMFFile("output/density.xdmf")
mu, lmbda = Constant(0.4), Constant(0.6)
# PREPARE FINITE ELEMENT ANALYSIS ---------------------------------
mesh = RectangleMesh(Point(0, 0), Point(nelx, nely), nelx, nely, "right/left")
U = VectorFunctionSpace(mesh, "P", 1)
D = FunctionSpace(mesh, "DG", 0)
u, v = TrialFunction(U), TestFunction(U)
u_sol, density_old, density = Function(U), Function(D), Function(D, name="density")
density.vector()[:] = volfrac
# DEFINE SUPPORT ---------------------------------
support = CompiledSubDomain("near(x[0], 0.0, tol) && on_boundary", tol=1e-14)
bcs = [DirichletBC(U,Constant((0.0, 0.0)), support)]
# DEFINE LOAD ---------------------------------
load_marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
CompiledSubDomain("x[0]==l && x[1]<=2", l=nelx).mark(load_marker, 1)
ds = Measure("ds")(subdomain_data=load_marker)
F = dot(v, Constant((0.0, -1.0))) * ds(1)
# SET UP THE VARIATIONAL PROBLEM AND SOLVER ---------------------------------
K = inner(density**penal*sigma(u), grad(v)) * dx
problem = LinearVariationalProblem(K, F, u_sol, bcs=bcs)
solver = LinearVariationalSolver(problem)
# PREPARE DISTANCE MATRICES FOR FILTER ---------------------------------
midpoint = [cell.midpoint().array()[:] for cell in cells(mesh)]
distance_mat = rmin - sp.euclidean_distances(midpoint, midpoint)
distance_mat[distance_mat < 0] = 0
distance_sum = distance_mat.sum(1)
# START ITERATION ---------------------------------
loop, change = 0, 1
while change > 0.01 and loop < 2000:
loop = loop + 1
density_old.assign(density)
# FE-ANALYSIS ---------------------------------
solver.solve()
# OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS ---------------------------------
objective = project(density ** penal * psi(u_sol), D).vector()[:]
sensitivity = -penal * (density.vector()[:]) ** (penal-1) * project(psi(u_sol), D).vector()[:]
# FILTERING/MODIFICATION OF SENSITIVITIES ---------------------------------
sensitivity = np.divide(distance_mat @ np.multiply(density.vector()[:], sensitivity), np.multiply(density.vector()[:], distance_sum))
# DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD ---------------------------------
l1, l2, move = 0, 100000, 0.2
while l2 - l1 > 1e-4:
l_mid = 0.5 * (l2 + l1)
density_new = np.maximum(0.001, np.maximum(density.vector()[:] - move, np.minimum(1.0, np.minimum(density.vector()[:] + move, density.vector()[:] * np.sqrt(-sensitivity/l_mid)))))
l1, l2 = (l_mid, l2) if sum(density_new) - volfrac*mesh.num_cells() > 0 else (l1, l_mid)
density.vector()[:] = density_new
# PRINT RESULTS ---------------------------------
change = norm(density.vector() - density_old.vector(), norm_type="linf", mesh=mesh)
print("it.: {0} , obj.: {1:.3f} Vol.: {2:.3f}, ch.: {3:.3f}".format(loop, sum(objective), sum(density.vector()[:]) / mesh.num_cells(), change))
xdemf.write(density, loop)
I read and applied the dolfinx implementation as much as possible, but I couldn’t figure out how to proceed with the translation in # DEFINE LOAD.
So far the code I’ve tried to translate is:
import ufl
import numpy as np, sklearn.metrics.pairwise as sp
from dolfinx import fem, mesh, io
from ufl import dx, sym, grad, tr, Identity, Constant
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
# A 55 LINE TOPOLOGY OPTIMIZATION CODE ---------------------------------
def main(nelx, nely, volfrac, penal, rmin):
sigma = lambda _u: 2.0 * mu * sym(grad(_u)) + lmd * tr(sym(grad(_u))) * Identity(len(_u))
psi = lambda _u: lmd / 2 * (tr(sym(grad(_u))) ** 2) + mu * tr(sym(grad(_u)) * sym(grad(_u)))
# xdemf = io.XDMFFile("output/density.xdmf") **# How to convert to XDMF has not been translated yet.**
# mu, lmbda = Constant(0.4), Constant(0.6)
mu, lmd = ScalarType(0.4), ScalarType(0.6)
# PREPARE FINITE ELEMENT ANALYSIS ---------------------------------
# mesh = RectangleMesh(Point(0, 0), Point(nelx, nely), nelx, nely, "right/left")
msh = mesh.create_rectangle(MPI.COMM_WORLD, ((0.0, 0.0), (nelx, nely)), (nelx, nely))
U = fem.VectorFunctionSpace(msh, ("CG", 1))
D = fem.FunctionSpace(msh, ("DG", 0))
u, v = ufl.TrialFunction(U), ufl.TestFunction(U)
u_sol, density_old, density = fem.Function(U), fem.Function(D), fem.Function(D, name="density")
density.vector()[:] = volfrac
# DEFINE SUPPORT ---------------------------------
# support = CompiledSubDomain("near(x[0], 0.0, tol) && on_boundary", tol=1e-14)
# bcs = [DirichletBC(U,Constant((0.0, 0.0)), support)]
def left(x):
return np.isclose(x[0], 0.0)
u_zero = np.array((0,) * msh.geometry.dim, dtype=ScalarType)
bc_l = fem.dirichletbc(u_zero, fem.locate_dofs_geometrical(U, left), U)
bcs = [bc_l]
# DEFINE LOAD ---------------------------------
# load_marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
# MeshFunction: mesh entities with integer markers
tdm = msh.topology.dim
fdm = tdm - 1
msh.topology.create_connectivity(fdm, tdm)
load_facets = mesh.exterior_facet_indices(msh.topology)
load_marker = *# From this point onwards, the code could not be converted.*
CompiledSubDomain("x[0]==l && x[1]<=2", l=nelx).mark(load_marker, 1)
ds = ufl.Measure("ds", subdomain_data=load_marker)
f = ufl.dot(v, Constant((0.0, -1.0))) * ds(1)
# SET UP THE VARIATIONAL PROBLEM AND SOLVER ---------------------------------
k = ufl.inner(density ** penal * sigma(u), grad(v)) * dx
problem = fem.petsc.LinearProblem(k, f, u_sol, bcs)
solver = problem.solver
# PREPARE DISTANCE MATRICES FOR FILTER ---------------------------------
midpoint = [cell.midpoint().array()[:] for cell in cells(mesh)]
distance_mat = rmin - sp.euclidean_distances(midpoint, midpoint)
distance_mat[distance_mat < 0] = 0
distance_sum = distance_mat.sum(1)
# START ITERATION ---------------------------------
loop, change = 0, 1
while change > 0.01 and loop < 2000:
loop += 1
density_old.assign(density)
# FE-ANALYSIS ---------------------------------
solver.solve()
# OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS ---------------------------------
objective = project(density ** penal * psi(u_sol), D).vector()[:]
sensitivity = -penal * (density.vector()[:]) ** (penal-1) * project(psi(u_sol), D).vector()[:]
# FILTERING/MODIFICATION OF SENSITIVITIES ---------------------------------
sensitivity = np.divide(distance_mat @ np.multiply(density.vector()[:], sensitivity), np.multiply(density.vector()[:], distance_sum))
# DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD ---------------------------------
l1, l2, move = 0, 1e9, 0.2
while (l2-l1)/(l1+l2) > 1e-3:
l_mid = 0.5*(l2+l1)
density_new = np.maximum(0.001, np.maximum(density.vector()[:] - move, np.minimum(1.0, np.minimum(density.vector()[:] + move, density.vector()[:] * np.sqrt(-sensitivity/l_mid)))))
l1, l2 = (l_mid, l2) if sum(density_new) - volfrac*mesh.num_cells() > 0 else (l1, l_mid)
density.vector()[:] = density_new
# PRINT RESULTS ---------------------------------
change = norm(density.vector() - density_old.vector(), norm_type="linf", mesh=mesh)
print("it.: {0} , obj.: {1:.3f} Vol.: {2:.3f}, ch.: {3:.3f}".format(loop, sum(objective), sum(density.vector()[:]) / mesh.num_cells(), change))
# xdemf.write(density, loop) **# This part is also commented out as it is an XDMF part.**
I read the implementation(Implementation — FEniCSx tutorial), but the following error occurred in the currently installed dolfinx environment.
“AttributeError: module ‘dolfinx.mesh’ has no attribute ‘exterior_facet_indices’”
Hence,
May I ask together if there is a way to solve the problem and if there are errors in the translation so far?
(Among the codes I have modified, the ones treated with ‘#’ are related to the existing code or description.)
