Assign subdomain id to the measure

mesh
1

Hello everyone,

I have the following problem. I have a mesh of a brain slice and within this slice there is another region that is supposed to mark a tumor region. Since I am building a model to describe tumor progression in the brain, I am feeling my way forward step by step. What I want to do now is to set up a differential measure for the tumor region itself. In doing so, I’m having trouble assigning the subdomain ID to this measure. Here is the relevant part of my code so far:

with XDMFFile(MPI.COMM_WORLD, "Domain.xdmf", "r") as xdmf:
    msh = xdmf.read_mesh(name="Grid")
    healthy_unhealthy = xdmf.read_meshtags(msh, name="Grid")

with XDMFFile(MPI.COMM_WORLD, "Domain_boundary.xdmf", "r") as xdmf:
    boundary_msh = xdmf.read_mesh(name="Grid")
    boundary = xdmf.read_meshtags(boundary_msh, name="Grid")

In the snippet above, healthy_unhealthy contains the two areas of interest. The boundary_msh contains the boundaries.
Within the given brain slice, I have 6 physical groups describing the subdomains, as seen below.

healthy = healthy_unhealthy.find(1)
unhealthy = healthy_unhealthy.find(2)
outer = boundary.find(3)
inner1 = boundary.find(4)
inner2 = boundary.find(5)
tumor_bound = boundary.find(6)

Now i want to define the differential measure for the tumor region as

dx_unhealthy = ufl.Measure("dx", domain=msh, subdomain_data=facet_tag, subdomain_id=?)

for the weak formulation

F = (dt * d * inner(grad(u), grad(v)) + dt * inner(inner(v, u), (K_u - u))) * dx_unhealthy + inner(u, v) * dx_unhealthy - inner(u0, v) * dx_unhealthy

The question I have now is how can I assign the subdomain for the tumor region (unhealthy = healthy_unhealthy.find(2)) to my measure?
I found some similar problems, but after a few attempts I decided to post in the forum. Please let me know if I make any major mistakes or forget important things.

Any help would be great, thanks in advance.
Julian

2

Should be as simple as

dx = ufl.Measure("dx", domain=msh, subdomain_data=celltags)
dx_unhealthy = dx(2)

where, from your snippets, celltags seems to be the variable healthy_unhealthy.

3

If the cells that are healthy and or unhealthy changes over time, i would use a DG-0 function as a marker (1 in unhealthy cells, 0 in healthy cells). Once:
Mixed dimension assembly in C++ and form independence for subdomains by jorgensd · Pull Request #3262 · FEniCS/dolfinx · GitHub Is merged there will be another approach that will be slightly more efficient.

4

Thanks, this seems to work. Now the occuring error is the known one for the solver.

Traceback (most recent call last):
  File "/home/jz97x/Mesh_25_06/Reaction_Diffusion_2D_V1_mesh_submesh.py", line 167, in <module>
    r = solver.solve(u)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/nls/petsc.py", line 47, in solve
    n, converged = super().solve(u.x.petsc_vec)
RuntimeError: Failed to successfully call PETSc function 'KSPSolve'. PETSc error code is: 73, Object is in wrong state

I set up the solver as following:

problem = NonlinearProblem(F, u, bcs=[bc])
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual" #incremental
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2

file = XDMFFile(MPI.COMM_WORLD, "output_V1_submesh.xdmf", "w") 
file.write_mesh(msh)

u0.x.array[:] = u.x.array

while t < T:
    t += dt 
    r = solver.solve(u)
    print(f"Step {int(t / dt)}: num iterations: {r[0]}")
    u0.x.array[:] = u.x.array 
    file.write_function(u, t)
file.close()

with boundary conditions as

facets = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),marker=lambda x: np.logical_or(np.isclose(x[0], 0.0),np.isclose(x[0], 2.0)))
dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)
bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)

Do you have experience with this error?

5

This usually means that you have not used u in your variational form:

Please see

or provide a minimal reproducible example

6

I´m defining

u.interpolate(initial_cond)

with a function for the initial condition after defining

u = fem.Function(V)

Without the interpolate(initial_cond) I´m getting the mesh as an output, but with no solution of the equation. This is the error, because I´m redefining the u. Is there a way that I can use u.interpolate but without redefining u? I have tried to introduce a new variable.

7

u.interpolate is not redefining u. As you have only provided code snippets and not a reproducible example, I cannot Give you any further guidance.

8

Sorry, I wanted to keep it as short and clear as possible, here is the full code :slight_smile: If required, I can also provide the mesh.

from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import numpy as np
import ufl
from dolfinx import fem, io, mesh, plot, default_real_type, default_scalar_type, log
from dolfinx.fem.petsc import LinearProblem, NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import XDMFFile, gmshio
from dolfinx.mesh import CellType, create_unit_square, meshtags
from ufl import dx, grad, inner
import pyvista as pv
import pyvistaqt as pvqt
from dolfinx import *

#       ∂ u 
#       ─── - div(d ∇ u) = u * (K - u)
#       ∂ t

dt = 0.02
t = 0.0
T = 0.1

def custom_function(x):
    #print('Test')
    #print(x[0],x[1])
    dummy = np.zeros(x[0].shape)
    for index, (X, Y) in enumerate(zip(x[0],x[1])):
        if X <= 0.1 and X >= 0.08 and Y >= -0.07 and Y <= -0.05:  
            dummy[index] = 1 + X**10 + Y**5
        else:
            dummy[index] = 10 + (X**5)/10 + (X**4)/5 + (X**3)
    return dummy

def initial_cond(x):
    # print('Test')
    # print(x[0],x[1])
    dummy = np.zeros(x[0].shape)
    for index, (X, Y) in enumerate(zip(x[0],x[1])): #replace with tumor region
        #if X <= 0.106 and X >= 0.098 and Y >= -0.102 and Y <= -0.095:
        if X <= 0.1 and X >= 0.08 and Y >= -0.07 and Y <= -0.05:
        #if X <= 0.075 and X >= 0.073 and Y >= -0.07 and Y <= -0.05:    
            dummy[index] = 10 + X**2 + Y**2
        else:
            dummy[index] = 0
    return dummy

with XDMFFile(MPI.COMM_WORLD, "Domain.xdmf", "r") as xdmf:
    msh = xdmf.read_mesh(name="Grid")
    healthy_unhealthy = xdmf.read_meshtags(msh, name="Grid")

with XDMFFile(MPI.COMM_WORLD, "Domain_boundary.xdmf", "r") as xdmf:
    boundary_msh = xdmf.read_mesh(name="Grid")
    boundary = xdmf.read_meshtags(boundary_msh, name="Grid")
        
healthy = healthy_unhealthy.find(1)
unhealthy = healthy_unhealthy.find(2)
outer = boundary.find(3)
inner1 = boundary.find(4)
inner2 = boundary.find(5)
tumor_bound = boundary.find(6)

V = fem.functionspace(msh, ("Lagrange", 1))

celltags = healthy_unhealthy
dx = ufl.Measure("dx", domain=msh, subdomain_data=celltags)
dx_unhealthy = dx(2)

#print("Healthy subdomain:", healthy_unhealthy.find(1))
#print("Unhealthy subdomain:", healthy_unhealthy.find(2))

facets = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),marker=lambda x: np.logical_or(np.isclose(x[0], 0.0),np.isclose(x[0], 2.0)))
dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)
bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)

x = ufl.SpatialCoordinate(msh)

K_u = fem.Function(V)
K_u.interpolate(custom_function)

u = fem.Function(V) 
u0 = fem.Function(V) 
v = ufl.TestFunction(V)

u.interpolate(initial_cond) 

d = 0.04

F = (dt * d * inner(grad(u), grad(v)) + dt * inner(inner(v, u), (K_u - u))) * dx_unhealthy + inner(u, v) * dx_unhealthy - inner(u0, v) * dx_unhealthy

problem = NonlinearProblem(F, u, bcs=[bc])
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual" #incremental
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2

file = XDMFFile(MPI.COMM_WORLD, "output_V1_submesh.xdmf", "w") 
file.write_mesh(msh)

u0.x.array[:] = u.x.array 

while t < T:
    t += dt 
    r = solver.solve(u)
    print(f"Step {int(t / dt)}: num iterations: {r[0]}")
    u0.x.array[:] = u.x.array 
    file.write_function(u, t)
file.close()

# u plotting
import pyvista
cells, types, x = plot.vtk_mesh(V)
grid = pyvista.UnstructuredGrid(cells, types, x)
print(u.x.array.real)
grid.point_data["u"] = u.x.array.real
grid.set_active_scalars("u")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
warped = grid.warp_by_scalar()
plotter.add_mesh(warped)
plotter.show()
9

It seems like you are trying to restrict the problem to only the “unhealthy” tissue.
To do so you should use the sub mesh creation by Joe Dean in DOLFINx: GitHub - jpdean/mixed_domain_demos
or multiphenicsx: https://multiphenics.github.io/ by Francesco Ballarin.