1D implementation of div-conforming (BDM) elements

[drew-parsons]

As I understand it, the div-conforming elements like BDM are suitable for representing flux, especially under conditions where a Dirichlet flux (controlled flux) is desired.

There is a BDMCF example in the dolfinx demos for a Mixed formulation for the Poisson equation, demo_mixed-poisson.py. I raised a question about adapting it for Dirichlet flux (controlling the normal component of the flux at the boundary) on the FEniCS discourse site at https://fenicsproject.discourse.group/t/setting-scalar-function-for-bdm-boundary-condition-with-controlled-flux-mixed-poisson/13508

Those demos used a square 2D domain. The demo can be trivially extended to 3D (a cubic domain).

But in some cases the symmetry of a system of interest (a cubic domain, for instance, with periodicity or null BCs in the lateral dimensions) means only a 1D solution is needed. It would be convenient if the same code base could easily switch between systems of different dimensions, using the same class of elements for consistency.

But the div-conforming elements, or at least BDMCF and RT as implemented in basix, do not support for 1D domains. If you try it, you get an error

ValueError: Unknown element family: BDMCF with cell type interval

I'm curious to know if it's a fundamental contradiction to try to set div-conforming elements in 1D. Or does the definition reduce to a trivial case in 1D that has simply not yet been implemented in basix? I don't understand the internal structure of the div-conforming elements well enough to say which it is.

If it's the latter, should the 1D case be added to basix for the div-conforming elements?

The use-case would be allowing a solver to handle different dimensionalities with the same code and same family of elements.
A simple demo is

# div-conforming elements with controlled boundary flux

# cf. https://fenicsproject.discourse.group/t/setting-scalar-function-for-bdm-boundary-condition-with-controlled-flux-mixed-poisson/13508
# adapted from dolfinx demo Mixed formulation for the Poisson equation (demo_mixed-poisson.py)
# https://docs.fenicsproject.org/dolfinx/v0.9.0/python/demos/demo_mixed-poisson.html

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np

from basix.ufl import element, mixed_element
from dolfinx import fem, io, mesh, default_scalar_type
from dolfinx.fem import Function
from dolfinx.fem.petsc import LinearProblem
from ufl import (CellDiameter, FacetNormal, Measure, TestFunctions, TrialFunctions,
                 div, inner)

###########################################
# Choose the dimensions of the test system.

# dim=2 and 3 successfully generate a solution over a square or cube
# dim=1 fails since BDMCF is not implemented for 1D

dim=1
#dim=2
#dim=3

###########################################

if dim==1:
    # simple 1D interval mesh with BDMCF elements gets
    # ValueError: Unknown element family: BDMCF with cell type interval
    domain = mesh.create_unit_interval(MPI.COMM_WORLD, 50)

    # alternative 1D interval mesh with mesh.CellType.interval gets
    # TypeError: Cannot interpret 'CellType.interval' as a data type
    # Is this error expected?
    #domain = mesh.create_unit_interval(MPI.COMM_WORLD, 50, mesh.CellType.interval)
elif dim==3:
    domain = mesh.create_unit_cube(MPI.COMM_WORLD, 20, 20, 20, mesh.CellType.hexahedron)
else:
    # 2d system
    domain = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32, mesh.CellType.quadrilateral)

right_boundary = lambda x: np.isclose(x[0], 1.0)
left_boundary = lambda x: np.isclose(x[0], 0.0)
right_bc_i = 1
left_bc_i = 2

boundaries = [(right_bc_i, right_boundary),
              (left_bc_i, left_boundary)]

facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = mesh.locate_entities(domain, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = Measure("ds", domain=domain, subdomain_data=facet_tag)


k = 1
Q_el = element("BDMCF", domain.basix_cell(), k)
P_el = element("DG", domain.basix_cell(), k - 1)
V_el = mixed_element([Q_el, P_el])
V = fem.functionspace(domain, V_el)

(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)

f = 0

dx = Measure("dx", domain)
a = inner(sigma, tau) * dx + inner(u, div(tau)) * dx + inner(div(sigma), v) * dx
L = -inner(f, v) * dx

Q, _ = V.sub(0).collapse()

# The boundary condition of the flux is
# constant outflow out of both the left and the right boundaries
# (+1 in the x-direction on the right, -1 in the x-direction on the left)
boundary_flux_value = 1

# impose BC weakly
bcs = []

# weak condition
n = FacetNormal(V.mesh)

# Nitsche penalty parameter
h = CellDiameter(domain)
alpha = fem.Constant(domain, default_scalar_type(10))
Nitsche_penalty = alpha/h

# weak Dirichlet boundary condition as scalar normal projection
boundary_flux = fem.Constant(domain, default_scalar_type(boundary_flux_value))
a += (Nitsche_penalty*inner(sigma, n) * inner(tau,n)  * ds(right_bc_i)
      + Nitsche_penalty*inner(sigma, n) * inner(tau,n)  * ds(left_bc_i))
L += (Nitsche_penalty*boundary_flux * inner(tau,n)  * ds(right_bc_i)
      + Nitsche_penalty*boundary_flux * inner(tau,n)  * ds(left_bc_i))
    
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu",
                                                       "pc_factor_mat_solver_type": "mumps"})
try:
    w_h = problem.solve()
except PETSc.Error as e:  # type: ignore
    if e.ierr == 92:
        print("The required PETSc solver/preconditioner is not available. Exiting.")
        print(e)
        exit(0)
    else:
        raise e

sigma_h, u_h = w_h.split()

with io.XDMFFile(domain.comm, f"out_mixed_poisson_weak_dim{dim}/u.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_function(u_h)

# interpolate vector solution to CG elements for visualisation
Vl = fem.functionspace(domain, ("Lagrange", k, (domain.geometry.dim,)))
sigmal_h = Function(Vl)
sigmal_h.interpolate(sigma_h)
with io.VTXWriter(domain.comm, f"out_mixed_poisson_weak_dim{dim}/sigmal.bp", sigmal_h) as file:
    file.write(0)

A successful 2D square or 3D cube solution is obtained if dim is set to 2 or 3. But dim=1 gets the ValueError

[jpdean]

In 1D, $H(\mathrm{div})=H^1$, so Lagrange elements can be used

[drew-parsons]

I suspected that might be the case. The definition itself is simple then.

To clarify the use-case, the aim would be to allow the code to be as simple as possible without extra conditional statements. i.e. setting the element to "BDMCF" in just one line, without having to add an alternative special case for domain.topology.dim==1