Interpolation issue for topology optimisation

dolfinx
1

Good morning,

I am trying to convert the FEniCS code contained in this article to FEniCSx, but i have a really hard time with these lines (in the full code in the paper’s annex, lines 42&49):

sensitivity = -penal * (density.vector()[:]) ** (penal - 1) * project(psi(u_sol), D).vector()[:]
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)))))

Here the u_sol is definined on CG1 space and density on DG0, thus needing projection.
My tries always end up with an interpolation error resulting in infinite recursion or type incompatibilities.
Does someone know how efficiently wrap these together?
Thanks in advance.

2

Without having a minimal code in DOLFINx reproducing your error message, it is hard to give you any guidance. If you add your converted code (as far as you have gotten) and supply the error message, people might be more inclined to help you.

3

Good afternoon,
I tried this solution but the compliance doesn’t change at all with every iteration. There must be an issue somewhere:

#IMPORTS
import numpy as np
from dolfinx import fem, cpp, io, mesh
from mpi4py import MPI
from petsc4py import PETSc
import ufl

#PARAMETERS
NITER_TOTAL=300
TOL=1e4
THETAMIN=0.001
NOMFICHIER="opti_test.xdmf"
M=0.2
ETA=0.5
PENAL=2
LMIN_INIT, LMAX_INIT = 0, 1e5
DENS_CIBLE=0.3
WEIGHT=1
LAMBDA=0.6
MU=0.4
L, W, H=0.303, 0.19, 0.064
EMP=0.06
AXE=0.025
BOTTOM_MARKER=1
AXES_MARKER=2
NX, NY, NZ = 30, 19, 6

#MESH CREATION
points=[np.array([0, 0, 0]), np.array([L, W, H])]
domain=mesh.create_box( 
    MPI.COMM_WORLD, 
    points,
    [NX,NY,NZ],
    cell_type=cpp.mesh.CellType.hexahedron,
    ghost_mode=mesh.GhostMode.shared_facet
    )
dim=domain.topology.dim

#BOUNDARY MARKERS
def axes(x):
    return np.logical_and(
               np.logical_or(np.sqrt( (x[0]-EMP)**2 + (x[2]-H/2)**2 ) <=AXE,
                             np.sqrt( (x[0]-L+EMP)**2 + (x[2]-H/2)**2 ) <=AXE),
               np.logical_or(np.isclose(x[1], 0, atol=1e-6), 
                             np.isclose(x[1], W, atol=1e-6))
           )
def bottom(x):
    return np.isclose(x[2], 0, atol=1e-6)

bottom_facets=mesh.locate_entities(domain, dim-1, bottom)
axes_facets=mesh.locate_entities(domain, dim-1, axes)
facet_indices=[bottom_facets,axes_facets]
facet_indices=np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers=[np.full(len(bottom_facets), BOTTOM_MARKER), np.full(len(axes_facets), AXES_MARKER)]
facet_markers=np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets=np.argsort(facet_indices)
facet_tag=mesh.meshtags(
    domain, 
    dim-1, 
    facet_indices[sorted_facets],
    facet_markers[sorted_facets]
    )

#LOADS AND BC
U=fem.VectorFunctionSpace(domain, ("CG", 1))
u=fem.Function(U, name="Déplacement")

bcdofs=fem.locate_dofs_geometrical(U, axes)
u_0=np.array((0,)*dim, dtype=PETSc.ScalarType)
bc=fem.dirichletbc(u_0, bcdofs, U)

ds= ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)
bottom_surf=fem.assemble_scalar(fem.form(fem.Constant(domain, 1.)*ds(1)))

def eps(v):
    return ufl.sym(ufl.grad(v))
def sigma(v):
    return (LAMBDA*ufl.nabla_div(v)*ufl.Identity(dim)+2*MU*eps(v))
def dE(u, v):
    return ufl.inner(sigma(u), eps(v))/2


#DENSITY SET UP
T=fem.FunctionSpace(domain, ("DG", 0))

density_old=fem.Function(T)
density=fem.Function(T, name="Density")
def homogenous(x):
    return np.full(len(x[0]), DENS_CIBLE)
density.interpolate(homogenous)
VOL=fem.assemble_scalar(fem.form(fem.Constant(domain, 1.)*ufl.dx))
frac_init=fem.assemble_scalar(fem.form(density_old*ufl.dx))/VOL

#VARIATONAL PROBLEM
u_=ufl.TestFunction(U)
du=ufl.TrialFunction(U)
a=ufl.inner(density**PENAL*sigma(u_), eps(du))*ufl.dx
b=ufl.dot((WEIGHT/bottom_surf)*ufl.FacetNormal(domain), u_)*ds(1)

solver=PETSc.KSP().create(domain.comm)
OPTIONS=PETSc.Options()
OPTIONS["ksp_type"]="cg"
OPTIONS["ksp_atol"]=1.0e-10
OPTIONS["pc_type"]="ilu"
solver.setFromOptions()
compliance_history=[]
compliance_old=1e30

file = io.XDMFFile(domain.comm, NOMFICHIER, 'w')
file.write_mesh(domain)

#OPTIMISATION LOOP
for i in range(NITER_TOTAL):

    density_old.interpolate(density)

    #SOLVE
    A=fem.petsc.assemble_matrix(fem.form(a), bcs=[bc])
    A.assemble()

    B=fem.petsc.assemble_vector(fem.form(b))
    fem.petsc.apply_lifting(B, [fem.form(a)], bcs=[[bc]])
    B.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
    fem.petsc.set_bc(B, [bc])

    solver.setOperators(A)
    iter_solve=solver.solve(B, u.vector)
    u.x.scatter_forward()

    file.write_function(density_old, i)
    file.write_function(u, i)

    #OPTIMISATION CRITERIA
    sensitivity=-PENAL*density**(PENAL-1)*dE(u,u)

    #LAGRANGE MULTIPLIER FINDER (DICHOTOMY)
    lmin, lmax=LMIN_INIT, LMAX_INIT
    while lmax-lmin>TOL:
        lmoy=(lmin+lmax)/2
        densitynew=fem.Function(T)
        densitynew.interpolate(fem.Expression((-sensitivity/lmoy)**ETA, T.element.interpolation_points))
        criterion=np.maximum(THETAMIN,
                             np.maximum(density_old.vector-M,
                                        np.minimum(1,
                                                   np.minimum(density_old.vector+M,
                                                               density.vector
                                                             )
                                                   )
                                       )         
                            )
        densitynew.vector.setArray(criterion)
        frac=fem.assemble_scalar(fem.form(densitynew*ufl.dx))/VOL
        if  frac-frac_init>0:
            lmin=lmoy
        else:
            lmax=lmoy
    density.vector.setArray(densitynew.vector.getArray())

    compliance=fem.assemble_scalar(fem.form(ufl.action(b,u)))
    compliance_history.append(compliance)
    print("Iteration {}: compliance=".format(i), compliance)

file.close()

Please excuse le if the code isn’t the most minimal, just don’t want to break the code already working.
Thanks again!