A question arose during the process of converting existing 55-line code of topology optimization to dolfinx

Hello :slight_smile:

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.)

I guess you are using dolfinx v0.4.1.
The tutorial is using the development version (0.4.2.0).

If you want to look at tutorials for the 0.4.1 version, see the notebooks in: GitHub - jorgensd/dolfinx-tutorial at v0.4.0

I’m currently not at a computer, so I can’t say much about your code. However, I would suggest you add in some spaces to make the code more readable, as it is currently quite dense.

1 Like

Thank you! I will edit the contents.

P.S) In addition, the Linux (Ubuntu 22.04 LTS) workstation was constructed and the dolfinx environment was installed in our Lab.
As in this case, If the version does not match, a version update seems to be necessary. I will also ask if I need to reinstall FEniCSx after uninstalling it, or if there is an update command in Ubuntu environment. :slight_smile:

to dokken

Sorry for the late additional question. :disappointed_relieved:

Reading the tutorial you recommended, ‘MeshFunction’ was replaced with ‘MeshTags’ in the dolfinx library. It is difficult to find a tutorial that uses MeshTags. How should I use it?
Previously, ‘MeshFunction’ was used in dolfin to display mesh entities with integer markers in the subdomain part obtained with load_marker.

P.S) Could you possibly get an answer to the question I asked above?

As I said previously, there are plenty of examples using MeshTags in the 0.4.0 version of the tutorial, see for instance: dolfinx-tutorial/subdomains.ipynb at v0.4.0 · jorgensd/dolfinx-tutorial · GitHub
or
dolfinx-tutorial/robin_neumann_dirichlet.ipynb at v0.4.0 · jorgensd/dolfinx-tutorial · GitHub

As I said in the previous post, exterior_facet_indices are only available on main, and you need to use the appropriate syntax from 0.4.0, you can see the changes in the diff of: Remove duplicate functionality for finding exterior facets (`mesh:: compute_boundary_facets` removed) by garth-wells · Pull Request #2203 · FEniCS/dolfinx · GitHub