Function on dirichletBC in a mixed functionspace

dolfinx
1

Hi all,

Building upon other posts: splitting mixed function spaces and DirichletBC depending on both time and another function.

I would like to apply a function to a dirichletbc within a mixed function space. The following mwe runs in dolfinx v0.7.0, however resulting in 0 everywhere.

from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.io import XDMFFile
import dolfinx.fem as fem
import dolfinx.fem.petsc
from dolfinx.nls.petsc import NewtonSolver
import ufl
import dolfinx.mesh
import numpy as np
from dolfinx import log

my_mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 20, 20)
fdim = my_mesh.topology.dim - 1

elements = ufl.FiniteElement("CG", my_mesh.ufl_cell(), 1)
mixed_element = ufl.MixedElement([elements, elements])
V = fem.FunctionSpace(my_mesh, mixed_element)
u = fem.Function(V)
u_n = fem.Function(V)

c_1, c_2 = ufl.split(u)
v_c_1, v_c_2 = ufl.TestFunctions(V)

C1, C1_to_V = V.sub(0).collapse()
C2, C2_to_V = V.sub(1).collapse()

indices_L = dolfinx.mesh.locate_entities_boundary(
    my_mesh, fdim, lambda x: np.isclose(x[0], 0)
)
indices_R = dolfinx.mesh.locate_entities_boundary(
    my_mesh, fdim, lambda x: np.isclose(x[0], 1)
)
left_dofs_c1 = fem.locate_dofs_topological(V.sub(0), fdim, indices_L)
right_dofs_c1 = fem.locate_dofs_topological(V.sub(0), fdim, indices_R)

x = ufl.SpatialCoordinate(my_mesh)
my_bc_value = 100 * x[1]
value_function = fem.Function(C1)
bc_expr = fem.Expression(my_bc_value, C1.element.interpolation_points())
value_function.interpolate(bc_expr)

bcs = [
    fem.dirichletbc(value=value_function, dofs=left_dofs_c1),
    fem.dirichletbc(
        fem.Constant(my_mesh, PETSc.ScalarType(0)),
        right_dofs_c1,
        V.sub(0),
    ),
]

dx = ufl.Measure("dx", domain=my_mesh)
F = ufl.dot(ufl.grad(c_1), ufl.grad(v_c_1)) * dx
F += ufl.dot(ufl.grad(c_2), ufl.grad(v_c_2)) * dx

problem = dolfinx.fem.petsc.NonlinearProblem(F, u, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
log.set_log_level(log.LogLevel.INFO)

output_c1_xdmf = XDMFFile(MPI.COMM_WORLD, "fenics_MWE.xdmf", "w")
output_c1_xdmf.write_mesh(my_mesh)

solver.solve(u)
output_c1_xdmf.write_function(u.sub(0))

Is there a better way to do this?

Thanks in advance,
James

2

Actually, think I found a solution already, by defining the boundary dofs using a tuple of the sub space and the collapsed function space:

fem.locate_dofs_topological((V.sub(0), C1), fdim, indices_L)

and adding the function space sub space to the dirichletbc:

fem.dirichletbc(value_function, left_dofs_c1, V.sub(0))

This seemed to fix the issue.

Heres the full code for anyone interested:

from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.io import XDMFFile
import dolfinx.fem as fem
import dolfinx.fem.petsc
from dolfinx.nls.petsc import NewtonSolver
import ufl
import dolfinx.mesh
import numpy as np
from dolfinx import log

my_mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 20, 20)
fdim = my_mesh.topology.dim - 1

elements = ufl.FiniteElement("CG", my_mesh.ufl_cell(), 1)
mixed_element = ufl.MixedElement([elements, elements])
V = fem.functionspace(my_mesh, mixed_element)
u = fem.Function(V)
u_n = fem.Function(V)

c_1, c_2 = ufl.split(u)
v_c_1, v_c_2 = ufl.TestFunctions(V)

C1, C1_to_V = V.sub(0).collapse()
C2, C2_to_V = V.sub(1).collapse()

indices_L = dolfinx.mesh.locate_entities_boundary(
    my_mesh, fdim, lambda x: np.isclose(x[0], 0)
)
indices_R = dolfinx.mesh.locate_entities_boundary(
    my_mesh, fdim, lambda x: np.isclose(x[0], 1)
)
left_dofs_c1 = fem.locate_dofs_topological((V.sub(0), C1), fdim, indices_L)
right_dofs_c1 = fem.locate_dofs_topological(V.sub(0), fdim, indices_R)

x = ufl.SpatialCoordinate(my_mesh)
my_bc_value = 100 * x[1]
value_function = fem.Function(C1)
bc_expr = fem.Expression(my_bc_value, C1.element.interpolation_points())
value_function.interpolate(bc_expr)

bcs = [
    fem.dirichletbc(value_function, left_dofs_c1, V.sub(0)),
    fem.dirichletbc(
        fem.Constant(my_mesh, PETSc.ScalarType(0)),
        right_dofs_c1,
        V.sub(0),
    ),
]

dx = ufl.Measure("dx", domain=my_mesh)
F = ufl.dot(ufl.grad(c_1), ufl.grad(v_c_1)) * dx
F += ufl.dot(ufl.grad(c_2), ufl.grad(v_c_2)) * dx

problem = dolfinx.fem.petsc.NonlinearProblem(F, u, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
log.set_log_level(log.LogLevel.INFO)

output_c1_xdmf = XDMFFile(MPI.COMM_WORLD, "fenics_MWE.xdmf", "w")
output_c1_xdmf.write_mesh(my_mesh)

solver.solve(u)
output_c1_xdmf.write_function(u.sub(0))
output_c1_xdmf.close()
1 Like