Define bilinear form with subdomains

General
1

Hi all,
I am attempting to solve a problem with an expression which is only non-zero in certain subdomain. To be specific, I need to add a term in the bilinear form : some_function * u * v * dx(0) where dx(0) is the ufl.Measure for the corresponding subdomain. However, it only works if I define the entire bilinear form using the defined measure dx(index) without ufl.dx, i.e. on top of

# works
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx(0) + ufl.inner(ufl.grad(u), ufl.grad(v)) * dx(1) + x[0] * ufl.inner(u, v) * dx(0)
# error
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx + x[0] * ufl.inner(u, v) * dx(0)

Is there a better approach? This is inconvenient as the term inner(grad(u), grad(v)) is the same in all subdomains. Below is the full script to reproduce the issue.
Thanks a lot !!

import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import mesh, fem
from dolfinx.fem import FunctionSpace
from petsc4py.PETSc import ScalarType

domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = FunctionSpace(domain, ('CG', 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(domain)

Q = FunctionSpace(domain, ("DG", 0))
Omega = {
    0: lambda x: x[1] <= 0.5,
    1: lambda x: x[1] >= 0.5,
}

# mark subdomains
tdim = domain.topology.dim
domain_indices, domain_markers = [], []
for tag, omega in Omega.items():
    entities = mesh.locate_entities(domain, tdim, omega)
    domain_indices.append(entities)
    domain_markers.append(np.full_like(entities, tag))

domain_indices = np.hstack(domain_indices).astype(np.int32)
domain_markers = np.hstack(domain_markers).astype(np.int32)
sorted_domain = np.argsort(domain_indices)
subdomain_tag = mesh.meshtags(domain, tdim, domain_indices[sorted_domain], domain_markers[sorted_domain])

dx = ufl.Measure('dx', domain=domain, subdomain_data=subdomain_tag)

a = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx(0) + ufl.inner(ufl.grad(u), ufl.grad(v)) * dx(1) 
# this line will cause an error when solving the problem
# a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx 
a += x[0] * ufl.inner(u, v) * dx(0)
L = fem.Constant(domain, ScalarType(1)) * v * ufl.dx

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

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

What about using:

a = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx 
a += x[0] * ufl.inner(u, v) * dx(0)
1 Like