How to generate multiple subdomains based on each cell

General
1

Dear all,
I want to have multiple subdomains based on each cell. For example, I have 100 cells in my mesh, then I need to have 100 subdomains, in order to do the integration around each cell. I know this can be done by a For loop. However, when there are thousands of cells, it will take a lot of time to do the for loop at every time step. The code attached below is my thought to do that, but it fails.

  1. Can anyone help me with this problem?
  2. Or if you have a better idea to do the integration around each point, please let me know?
  3. Is there a way to find the index of points during “interpolate”?

Thanks a lot for your time and help!

import numpy as np
from mpi4py import MPI
import ufl
from dolfinx import fem, io, mesh, plot, common, geometry
from ufl import ds, dx, grad, inner

# Create mesh and define function space
msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                            points=((0.0, 0.0), (1, 1)), n=(10, 10),
                            cell_type=mesh.CellType.quadrilateral)

V = fem.FunctionSpace(msh, ("CG", 1))
Q = fem.FunctionSpace(msh, ("DG", 0))

P = fem.Function(Q)
T = fem.Function(V)
T.interpolate(lambda x: x[0] + 2 * x[1])

R = 0.02

def f_(x):
    a = np.array(x[0])
    b = np.array(x[1])

    subcells = mesh.locate_entities(msh, msh.topology.dim, lambda y: np.sqrt((y[0]-a)**2 + (y[1]-b)**2) <= R)
    subdm = mesh.meshtags(msh, msh.topology.dim, subcells, np.array([point_index], dtype = np.int32)) 
    #here the point_index should be the index of each point and the length should be equal to subcells
    dx_subdm = ufl.Measure("dx", subdomain_data = subdm)    
    
    return fem.assemble_scalar(fem.form(T * dx_subdm(point_index)))

P.interpolate(f_)
print(P.x.array)
2

See for instance: Is there a more efficient way to many integrals over subdomains? - #2 by dokken

3

Hello Dokken,

Thanks a lot. That one really helps.

One more question:

If the function I want to apply to the subdomain f is like this 1-r/R, do you have any good idea to compile this function?

R is a fixed number that is the border of the subdomain. r is the distance of each point in the subdomain to the midpoint (a cell).

Finally, I want to apply fem.form(f * dx)

4

You can compute the midpoint of any cell with compute_midpoints.
Something along the lines of

x = ufl.SpatialCoordinate(mesh)
tdim = mesh.topology.dim
num_cells_local = mesh.topology.index_map(tdim).size_local
midpoints = dolfinx.mesh.compute_midpoints(mesh, tdim, np.arange(num_cells_local, dtype=np.int32))
midpoint = dolfinx.fem.Constant(mesh, np.zeros(2, dtype=PETSc.ScalarType))
r = ufl.sqrt((x[0]-midpoint[0])**2 + (x[1]-midpoint[1])**2)
R = dolfinx.fem.Constant(mesh, some_fixed_value)
f = fem.form((1- r/R)*ufl.dx)
for cell in range(num_cells_local):
    midpoint_i = midpoints[cell]
    midpoint.value[:] = midpoint_i
    print(dolfinx.fem.assemble_scalar(f))
5

Hello Dokken,

Thanks a lot for your time and help. I think that’s what needed and I will try that.