Identifying subdomains from level set function on fenicsx

dolfinx
1

Hello,

I’m trying to understand the best way to identify subdomains defined from a level set function on dolfinx. Briefly speaking, I have a reference domain D and the domain of interest Omega is defined as x in D such that phi(x) < 0.

When an analytical expression of phi is know, i’m using “locate_entities” from dolfinx.mesh (“strategy1” in the example bellow). However, assuming that the expression of phi is not known, the only approach I was able to figure out was using ufl conditional (“strategy2” in the example). Although error with conditional is smaller, the behavior with mesh refinement seems a bit strange.

Is there any problems with this use of “conditional”? I wonder if there is a better approach when no analytical expression of phi is know.

I appreciate any help!

import numpy as np
from dolfinx import fem
from dolfinx.mesh import (locate_entities, create_rectangle, CellType, meshtags)
from mpi4py import MPI
from ufl import Measure, conditional, ge,dx
from matplotlib import pyplot as plt
import dolfinx.cpp as _cpp

def tag_subdomains(msh, subdomains): # Identifies and marks subdomains accoring to locator function
    cell_indices, cell_markers = [], [] # List for facet indices and respective markers
    fdim = msh.topology.dim
    for (marker, locator) in subdomains:
        cells = locate_entities(msh, fdim, locator)
        cell_indices.append(cells)
        cell_markers.append(np.full_like(cells, marker))
    cell_indices = np.hstack(cell_indices).astype(np.int32)
    cell_markers = np.hstack(cell_markers).astype(np.int32)
    sorted_cells = np.argsort(cell_indices)
    cell_tag = meshtags(msh, fdim, cell_indices[sorted_cells], cell_markers[sorted_cells]) 
    return cell_tag

tol = 1e-8
# simple example with circle inside a square
def phi(x):
    return np.sqrt((x[0] - 0.5)**2 + (x[1]-0.5)**2) - radius
    
subdomains = [(1, lambda x: phi(x) >= tol), 
                (0, lambda x: phi(x) <= tol)]

N_list = np.array([5,10,20,40,80,160,320])
err = []; err2 = []
for N in N_list:
    comm = MPI.COMM_WORLD
    msh = create_rectangle(comm=MPI.COMM_WORLD,
                        points=((0.0, 0.0), (1.0, 1.0)), n=(N, N),
                        cell_type=CellType.triangle, diagonal = _cpp.mesh.DiagonalType.crossed)

    # Strategy 1: analytical expression
    cell_tag = tag_subdomains(msh, subdomains)
    dx = Measure("dx", domain=msh, subdomain_data=cell_tag)
    one = fem.Constant(msh, 1.0)
    area = 1 - np.pi*radius**2
    err.append(np.abs(area - fem.assemble_scalar(fem.form(one*dx(1)))))

    # Strategy 2: if no analytical expression is known
    dx = Measure("dx", domain=msh)
    V = fem.FunctionSpace(msh, ("CG", 1))
    phi_ = fem.Function(V)
    phi_.interpolate(phi)
    err2.append(np.abs(area - fem.assemble_scalar(fem.form(conditional(ge(phi_, tol), 1, 0)*dx))))

plt.plot(1/N_list,err, label = 'strategy 1'); plt.plot(1/N_list,err2, label = 'strategy 2'); plt.grid()
plt.yscale("log"); plt.xscale("log") 
plt.xlabel('h'); plt.ylabel('error'); plt.legend()

error

2

First of all, you should look at the cell markers, as they clearly show that your approach might have some faults (adding this in your loop):

    with XDMFFile(msh.comm, f"tag_{N}.xdmf", "w") as xdmf:
        xdmf.write_mesh(msh)
        xdmf.write_meshtags(cell_tag)

coarse_tags
coarse_tags1531×789 26.2 KB

As you can see from the figure, you have a layer where some vertices satisfies one criteria, while others satisfies the other, thus no tagged cell.

You can alter this with the following code:

import numpy as np
from dolfinx import fem
from dolfinx.mesh import (locate_entities, create_rectangle, CellType, meshtags)
from mpi4py import MPI
from ufl import Measure, conditional, ge,dx
from matplotlib import pyplot as plt
import dolfinx.cpp as _cpp
from dolfinx.io import XDMFFile
radius = 0.25
def tag_subdomains(msh, levelset): # Identifies and marks subdomains accoring to locator function
    dim = msh.topology.dim
    num_cells = msh.topology.index_map(dim).size_local
    cells = np.arange(num_cells, dtype=np.int32)
    markers = np.zeros_like(cells)
    markers[locate_entities(msh, dim, levelset)] = 1
    cell_tag = meshtags(msh, dim, cells, markers)
    return cell_tag

tol = 1e-8
# simple example with circle inside a square
def phi(x):
    return np.sqrt((x[0] - 0.5)**2 + (x[1]-0.5)**2) - radius
    
subdomains = lambda x: phi(x) >= -tol

N_list = np.array([5,10,20,40,80,160,320])
err = []; err2 = []
for N in N_list:
    comm = MPI.COMM_WORLD
    msh = create_rectangle(comm=MPI.COMM_WORLD,
                        points=((0.0, 0.0), (1.0, 1.0)), n=(N, N),
                        cell_type=CellType.triangle, diagonal = _cpp.mesh.DiagonalType.crossed)

    # Strategy 1: analytical expression
    cell_tag = tag_subdomains(msh, subdomains)
    dx = Measure("dx", domain=msh, subdomain_data=cell_tag)
    one = fem.Constant(msh, 1.0)
    area = 1 - np.pi*radius**2
    err.append(np.abs(area - fem.assemble_scalar(fem.form(one*dx(1)))))

    # Strategy 2: if no analytical expression is known
    dx = Measure("dx", domain=msh)
    V = fem.FunctionSpace(msh, ("CG", 1))
    phi_ = fem.Function(V)
    phi_.interpolate(phi)
    err2.append(np.abs(area - fem.assemble_scalar(fem.form(conditional(ge(phi_, tol), 1, 0)*dx))))
    with XDMFFile(msh.comm, f"tag_{N}.xdmf", "w") as xdmf:
        xdmf.write_mesh(msh)
        xdmf.write_meshtags(cell_tag)
plt.plot(1/N_list,err, label = 'strategy 1'); plt.plot(1/N_list,err2, label = 'strategy 2'); plt.grid()
plt.yscale("log"); plt.xscale("log") 
plt.xlabel('h'); plt.ylabel('error'); plt.legend()
plt.savefig("rates.png")

giving

tags2
tags21531×789 22 KB

What do you mean by this? You are trying to fit your area integration with a non-matching mesh, using the whole cell, giving the non-smooth boundary as shown in the figures.

If the analytical level set is known, I would mesh the boundary.

Conditionals simply evaluate phi(x) at every quadrature point, making it the most accurate representation of an Nth order polynomial on the given interval.
This means that you can change the quadrature degree of the integral to improve the error estimates (up to a point):

err = []; err2 = []
qd = 3
for N in N_list:
    comm = MPI.COMM_WORLD
    msh = create_rectangle(comm=MPI.COMM_WORLD,
                        points=((0.0, 0.0), (1.0, 1.0)), n=(N, N),
                        cell_type=CellType.triangle, diagonal = _cpp.mesh.DiagonalType.crossed)

    # Strategy 1: analytical expression
    cell_tag = tag_subdomains(msh, subdomains)
    dx = Measure("dx", domain=msh, subdomain_data=cell_tag)
    one = fem.Constant(msh, 1.0)
    area = 1 - np.pi*radius**2
    err.append(np.abs(area - fem.assemble_scalar(fem.form(one*dx(1)))))

    # Strategy 2: if no analytical expression is known
    dx = Measure("dx", domain=msh, metadata={"quadrature_degree":qd})
    V = fem.FunctionSpace(msh, ("CG", 1))
    phi_ = fem.Function(V)
    phi_.interpolate(phi)
    err2.append(np.abs(area - fem.assemble_scalar(fem.form(conditional(ge(phi_, tol), 1, 0)*dx))))
    with XDMFFile(msh.comm, f"tag_{N}.xdmf", "w") as xdmf:
        xdmf.write_mesh(msh)
![rates_20|640x480](upload://frHmlLrvjPD258GoC9vyCjjfACN.png)
![rates_10|640x480](upload://nNYbQgVah1W057QjKonIY0lDawN.png)
![rates_3|640x480](upload://8l8H83WJg2cg5rUTJKULAHJnto4.png)

        xdmf.write_meshtags(cell_tag)
plt.plot(1/N_list,err, label = 'strategy 1'); plt.plot(1/N_list,err2, label = f'strategy 2 (degree: {qd})'); plt.grid()
plt.yscale("log"); plt.xscale("log") 
plt.xlabel('h'); plt.ylabel('error'); plt.legend()
plt.savefig(f"rates_{qd}.png")

Using different qd yields:
rates_3
rates_10
rates_20

2 Likes
3

Thank you very much!

4

Hello! May I ask another question? If I should make another post, just let me know.

I’m using ufl conditional to identify subdomains from a level set function, as described in the above discussion. However, I’m facing a memory issue (Error: error code 55). The thing is: i’m working with a shape optimization algorithm, so that i have to run several iterations. For each iteration I have to use the level set function to identify the subdomains and solve an EDP on the subdomains

Since the level set function/subdomains change, I’m initializing the LinearProblem for the EDP inside the loop. I belive that is what is causing the memory issue. However, I cannot figure out a solution. I’ve read other posts suggesting the use of fem.Constants, for example, but it doesn’t seem to apply with the use of conditional for mapping the subdomains. Do you have any sugestions?

Here is a simple example that reproduces my problem. In my computer the error appears after ~380 iterations. (In real examples i don’t have the analytical expression of the level set function, so that mapping the boundary is not an option)

import numpy as np
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
                         form, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_rectangle, CellType, locate_entities
from dolfinx.plot import create_vtk_mesh
from ufl import (Measure, SpatialCoordinate, TestFunction, TrialFunction,
                 dx, grad, inner, le, ge, conditional)
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import dolfinx.cpp as _cpp

def phi(x):
    return np.sqrt((x[0] - 0.5)**2 + (x[1]-0.5)**2) - radius

radius = 0.25; tol =  1e-8
msh = create_rectangle(comm=MPI.COMM_WORLD,
                    points=((0.0, 0.0), (1.0, 1.0)), n=(200, 100),
                    cell_type=CellType.triangle, diagonal = _cpp.mesh.DiagonalType.crossed)
comm = MPI.COMM_WORLD
V = FunctionSpace(msh, ("CG", 1))
dx = Measure("dx", domain=msh)

kappa = 0.1
eps = 10E-3

er = []; i = 0
r_list = np.array([0.1, 0.2, 0.3, 0.4])
for radius in np.linspace(0.1, 0.4, 10000):
    
    print(i)

    # Identifying the two-phases/subdomains from level-set
    phi_ = Function(V)
    phi_.interpolate(phi)
    interior = conditional(le(phi_, tol), 1, 0)
    exterior = conditional(ge(phi_, tol), 1, 0)

    # Solving two-phase problem
    u, v = TrialFunction(V), TestFunction(V)
    a = interior * inner(kappa*grad(u), grad(v)) * dx \
        + eps * exterior * inner(kappa*grad(u), grad(v)) * dx
    x = SpatialCoordinate(msh)
    L =  interior * Constant(msh, ScalarType(1)) * v * dx \
        + eps * exterior * Constant(msh, ScalarType(1)) * v * dx

    dofs = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
    bcs = [dirichletbc(ScalarType(1), dofs, V)]

    problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    uh = problem.solve()

    i+=1
5

It is simple to rewrite your problem to not having linear-problem inside the loop:

import dolfinx.cpp as _cpp
import numpy as np
from dolfinx.fem import (Constant, Function, FunctionSpace, dirichletbc,
                         locate_dofs_geometrical)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import CellType, create_rectangle
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import (Measure, SpatialCoordinate, TestFunction, TrialFunction,
                 conditional, dx, ge, grad, inner, le)


def phi(x):
    return np.sqrt((x[0] - 0.5)**2 + (x[1]-0.5)**2) - radius

radius = 0.25; tol =  1e-8
msh = create_rectangle(comm=MPI.COMM_WORLD,
                    points=((0.0, 0.0), (1.0, 1.0)), n=(200, 100),
                    cell_type=CellType.triangle, diagonal = _cpp.mesh.DiagonalType.crossed)
comm = MPI.COMM_WORLD
V = FunctionSpace(msh, ("CG", 1))
dx = Measure("dx", domain=msh)

kappa = 0.1
eps = 10E-3

er = []; i = 0
r_list = np.array([0.1, 0.2, 0.3, 0.4])

# Identifying the two-phases/subdomains from level-set
phi_ = Function(V)
interior = conditional(le(phi_, tol), 1, 0)
exterior = conditional(ge(phi_, tol), 1, 0)

# Solving two-phase problem
u, v = TrialFunction(V), TestFunction(V)
a = interior * inner(kappa*grad(u), grad(v)) * dx \
    + eps * exterior * inner(kappa*grad(u), grad(v)) * dx
x = SpatialCoordinate(msh)
L =  interior * Constant(msh, ScalarType(1)) * v * dx \
    + eps * exterior * Constant(msh, ScalarType(1)) * v * dx

dofs = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
bcs = [dirichletbc(ScalarType(1), dofs, V)]

problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})


for radius in np.linspace(0.1, 0.4, 10000):
    phi_.interpolate(phi)
    print(i)

    uh = problem.solve()
    print(uh.x.array)
    i+=1
1 Like
6

Hello @dokken, thank you for your quick reply! The results obtained with the two approaches seem to be different, however (i’ve checked by saving uh.x.array[:] in both cases).

The left-hand side of the system also depends on the level set function. Therefore, wouldn’t it be necessary to redefine the solver?

7

Calling solve re-assembles the matrix and vector (i.e. with the updated levelset phi_, which was updated with interpolation inside your loop.

To highlight this, see the following codes, that both produce the same result

import dolfinx.cpp as _cpp
import numpy as np
from dolfinx.fem import (Constant, Function, FunctionSpace, dirichletbc,
                         locate_dofs_geometrical)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import CellType, create_rectangle
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import (Measure, SpatialCoordinate, TestFunction, TrialFunction,
                 conditional, dx, ge, grad, inner, le)


def phi(x):
    return np.sqrt((x[0] - 0.5)**2 + (x[1]-0.5)**2) - radius

radius = 0.25; tol =  1e-8
msh = create_rectangle(comm=MPI.COMM_WORLD,
                    points=((0.0, 0.0), (1.0, 1.0)), n=(200, 100),
                    cell_type=CellType.triangle, diagonal = _cpp.mesh.DiagonalType.crossed)
comm = MPI.COMM_WORLD
V = FunctionSpace(msh, ("CG", 1))
dx = Measure("dx", domain=msh)

kappa = 0.1
eps = 10E-3

er = []; i = 0
r_list = np.array([0.1, 0.2, 0.3, 0.4])

# Identifying the two-phases/subdomains from level-set
phi_ = Function(V)
interior = conditional(le(phi_, tol), 1, 0)
exterior = conditional(ge(phi_, tol), 1, 0)

# Solving two-phase problem
u, v = TrialFunction(V), TestFunction(V)
a = interior * inner(kappa*grad(u), grad(v)) * dx \
    + eps * exterior * inner(kappa*grad(u), grad(v)) * dx
x = SpatialCoordinate(msh)
L =  interior * Constant(msh, ScalarType(1)) * v * dx \
    + eps * exterior * Constant(msh, ScalarType(1)) * v * dx

dofs = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
bcs = [dirichletbc(ScalarType(1), dofs, V)]

problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})


for radius in np.linspace(0.1, 0.4, 4):
    phi_.interpolate(phi)

    uh = problem.solve()
    print(i, radius, np.linalg.norm(uh.x.array))
    i+=1

produces the same result as

import numpy as np
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
                         form, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_rectangle, CellType, locate_entities
from dolfinx.plot import create_vtk_mesh
from ufl import (Measure, SpatialCoordinate, TestFunction, TrialFunction,
                 dx, grad, inner, le, ge, conditional)
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import dolfinx.cpp as _cpp

def phi(x):
    return np.sqrt((x[0] - 0.5)**2 + (x[1]-0.5)**2) - radius

radius = 0.25; tol =  1e-8
msh = create_rectangle(comm=MPI.COMM_WORLD,
                    points=((0.0, 0.0), (1.0, 1.0)), n=(200, 100),
                    cell_type=CellType.triangle, diagonal = _cpp.mesh.DiagonalType.crossed)
comm = MPI.COMM_WORLD
V = FunctionSpace(msh, ("CG", 1))
dx = Measure("dx", domain=msh)

kappa = 0.1
eps = 10E-3

er = []; i = 0
for radius in np.linspace(0.1, 0.4, 4):
    
    # Identifying the two-phases/subdomains from level-set
    phi_ = Function(V)
    phi_.interpolate(phi)
    interior = conditional(le(phi_, tol), 1, 0)
    exterior = conditional(ge(phi_, tol), 1, 0)

    # Solving two-phase problem
    u, v = TrialFunction(V), TestFunction(V)
    a = interior * inner(kappa*grad(u), grad(v)) * dx \
        + eps * exterior * inner(kappa*grad(u), grad(v)) * dx
    x = SpatialCoordinate(msh)
    L =  interior * Constant(msh, ScalarType(1)) * v * dx \
        + eps * exterior * Constant(msh, ScalarType(1)) * v * dx

    dofs = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
    bcs = [dirichletbc(ScalarType(1), dofs, V)]

    problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    uh = problem.solve()
    print(i, radius, np.linalg.norm(uh.x.array))
    i+=1
2 Likes
8

Hello, dokken. is it possible to change these codes to parallel version ？ When it is runned by mpirun, the indexes returned by locate_entities will exceed the size of marker array.

Traceback (most recent call last):
  File "/mnt/g/workdir/MKM/NO_mkm/test_a.py", line 35, in <module>
Traceback (most recent call last):
  File "/mnt/g/workdir/MKM/NO_mkm/test_a.py", line 35, in <module>
    cell_tag = tag_subdomains(msh, subdomains)
    cell_tag = tag_subdomains(msh, subdomains)
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/mnt/g/workdir/MKM/NO_mkm/test_a.py", line 15, in tag_subdomains
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/mnt/g/workdir/MKM/NO_mkm/test_a.py", line 15, in tag_subdomains
    markers[locate_entities(msh, dim, levelset)] = 1
    markers[locate_entities(msh, dim, levelset)] = 1
    ~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IndexError: index 53 is out of bounds for axis 0 with size 50
    ~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IndexError: index 51 is out of bounds for axis 0 with size 50

and the write_meshtags function needs an argument “x”, what is it ?

#this will cause error in the latest fenicsx
xdmf.write_meshtags(cell_tag)