Unable to apply dirichlet BC in mixed formulation

dolfinx
1

Hi,
Attached is the minimum working example.
Solution in uh.x.array only contains inf values.
Perhaps, my understanding about mixed formulation needs some guidance.
Please enlighten.

Thanks,

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


norm_lst = []
sigma_lst = []


def boundary_D(x):
    log1 = np.logical_or(np.isclose(x[0], 0), np.isclose(x[0], 1))
    log2 = np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1))
    return np.logical_or(log1, log2)


for i in [8, 16, 32, 64, 128]:
    domain = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                                   points=[(-1, -1), (1, 1)],
                                   n=(i, i),
                                   cell_type=mesh.CellType.triangle)

    BDM = ufl.FiniteElement(family="BDM", cell=domain.ufl_cell(), degree=1)
    DG = ufl.FiniteElement(family="DG", cell=domain.ufl_cell(), degree=0)
    W = fem.FunctionSpace(mesh=domain, element=ufl.MixedElement([BDM, DG]))

    (sigma, u) = ufl.TrialFunctions(function_space=W)
    (tau, v) = ufl.TestFunctions(function_space=W)

    f = fem.Constant(domain=domain, c=ScalarType(1.))

    a = (-ufl.dot(sigma, tau) + ufl.div(tau)*u - ufl.div(sigma)*v) * ufl.dx
    L = -f * v * ufl.dx

# ---------------------------------------------------------------------------
    V0, submap = W.sub(0).collapse()

    u_zero = fem.Function(V=V0)
    u_zero.x.array[:] = 0.0

    # determine boundary DOFs
    boundary_facets = mesh.locate_entities_boundary(
        domain, domain.topology.dim - 1, boundary_D)
    boundary_dofs = fem.locate_dofs_topological(
        (W.sub(0), V0), domain.topology.dim - 1, boundary_facets)
    bc = fem.dirichletbc(value=u_zero, dofs=boundary_dofs, V=W)
# ---------------------------------------------------------------------------
    problem = fem.petsc.LinearProblem(a=a, L=L,
                                      bcs=[bc],
                                      petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    uh = problem.solve()
2

This problem looks like a similar formulation to Add demo of mixed Poisson problem by jhale · Pull Request #2502 · FEniCS/dolfinx · GitHub
Have you considered having a look at it?

2 Likes
3

Thanks @dokken for you quick help. Appreciated!