Dirichlet boundary conditions not working with DG0 vector function space in dolfinx

dolfinx
1

Hi everyone,

I’m a beginner working with dolfinx, and I’m currently trying to solve a coupled ODE system using a vector-valued DG0 (discontinuous Galerkin) space. I apply Dirichlet boundary conditions at the left end of a 1D interval (x = 0), but they don’t seem to take effect when I use DG0 elements. However, when I switch to a CG1 (continuous Galerkin) space, the same boundary conditions work as expected.

The ODE system I’m solving is:

u_1' = -(\alpha + \beta_1) u_1 + \alpha u_2 + \beta_1 + S_1 \sin(P_1 x),
u_2' = \alpha u_1 - (\alpha + \beta_2) u_2 + S_2 \sin(P_2 x),
u_1(0) = 0.75, \quad u_2(0) = 0.25.

Below is a simplified version of my code. I’m applying fixed values at x = 0 using dirichletbc on V.sub(0) and V.sub(1). This works for CG1 but not for DG0 — the solution at the boundary doesn’t match the expected values.

Could someone kindly explain why this is happening? Is this behavior expected when using DG elements in dolfinx?

Any explanation or advice would be greatly appreciated. Thank you!

Here it’s my code:

import numpy as np
import matplotlib.pyplot as plt
import dolfinx
from mpi4py import MPI
from petsc4py import PETSc
from ufl import *
import basix
from dolfinx.fem import petsc

# Create mesh
T = 800.0
n = 80000
mesh = dolfinx.mesh.create_interval(MPI.COMM_WORLD, n, [0.0, T])

# DG0 vector-valued function space
vector_element = basix.ufl.element("DG", "interval", 0, shape=(2, 1))
V_vec = dolfinx.fem.functionspace(mesh, vector_element)

# Parameters
alpha = 0.1
beta1 = 0.01
beta2 = 0.01
S1 = 0.005
P1 = pi / 180
S2 = 0.3
P2 = pi

# Boundary condition at x = 0
facets = dolfinx.mesh.locate_entities_boundary(mesh, dim=0, marker=lambda x: np.isclose(x[0], 0.0))
dofs_1 = dolfinx.fem.locate_dofs_topological(V=V_vec, entity_dim=0, entities=facets)
dofs_2 = dolfinx.fem.locate_dofs_topological(V=V_vec, entity_dim=0, entities=facets)

bc_u1 = dolfinx.fem.dirichletbc(value=0.75, dofs=dofs_1, V=V_vec.sub(0))
bc_u2 = dolfinx.fem.dirichletbc(value=0.25, dofs=dofs_2, V=V_vec.sub(1))

# Define variational problem
u = TrialFunction(V_vec)
v = TestFunction(V_vec)
x = SpatialCoordinate(mesh)
A = as_matrix([[alpha + beta1, -alpha], [-alpha, alpha + beta2]])
f = as_matrix([[beta1 + S1 * sin(P1 * x[0])], [S2 * sin(P2 * x[0])]])

a = inner(dot(A, u), v) * dx + inner(jump(u), v('+')) * dS
L = inner(f, v) * dx

# Solve
problem = petsc.LinearProblem(a, L, bcs=[bc_u1, bc_u2])
uh = problem.solve()

# Extract solution
x_coords = mesh.geometry.x[:, 0]
u1 = uh.x.array[::2]
u2 = uh.x.array[1::2]

print(u1[0], u2[0])  # Values at left boundary

# Plot
plt.figure(figsize=(8, 5))
plt.plot(x_coords[:-1], u1, "r-", label="u1(x)", linewidth=2)
plt.plot(x_coords[:-1], u2, "b--", label="u2(x)", linewidth=2)
plt.xlabel("x")
plt.ylabel("Value")
plt.title("Components of Vector Function")
plt.grid(True)
plt.legend()
plt.show()
2

This is down to the fundamental definition of the discontinuous basis that you’re employing. In the case of "DG" of degree 0 this aligns with a discontinuous Lagrange basis of degree 0. This basis function consists only of piecewise constant functions defined in each cell. As such, defining your boundary conditions with:

dofs_1 = dolfinx.fem.locate_dofs_topological(V=V_vec, entity_dim=0, entities=facets)

and crucially the component entity_dim=0 which corresponds to vertices, will not yield any degree of freedom. This is because the DG0 basis degrees of freedom are associated with the cells of the mesh.

To resolve your issue you could either redesign your DoF location component to find cell degrees of freedom. Or, more naturally, impose Dirichlet boundary data weakly as is typical for a discontinuous Galerkin formation.

3

Thank you very much for your detailed and insightful explanation.
Your response really helped me understand the root cause of the issue and clarified the role of cell-based degrees of freedom in DG0 spaces. As a beginner, I truly appreciate the time and effort you took to guide me—it’s been extremely helpful! I’ll look into implementing the boundary conditions weakly as you suggested.

4

Hi Nate,

Thank you again for your previous explanation — it was really helpful!

After adjusting my weak form as you suggested, I was able to impose the Dirichlet boundary condition using a large penalty term. However, I noticed that changing the boundary value only affects the initial value (i.e., at the left boundary), but does not seem to influence the rest of the domain. I’m wondering if this might be due to a mistake or limitation in my current weak form?

I’m still new to DG methods, so I really appreciate your patience and guidance. Apologies for the repeated questions, and thank you once again for your kind support!

Here it’s my week form:

\int AUdx +\int ([U],v^{+} )dS+1e8*\int (u,v)ds = \int fvdx+1e8*\int (g,v)ds
g=(u_1(0),u_2(0)) 
import numpy as np
import matplotlib.pyplot as plt
import dolfinx
from mpi4py import MPI
from petsc4py import PETSc
from ufl import *
import basix
from dolfinx.fem import petsc

# Create a 1D mesh over the interval [0, T] with n cells
T = 800.0
n = 80000
mesh = dolfinx.mesh.create_interval(MPI.COMM_WORLD, n, [0.0, T])

# Define a discontinuous Galerkin (DG0) vector-valued function space (2 components)
vector_element = basix.ufl.element("DG", "interval", 0, shape=(2, 1))
V_vec = dolfinx.fem.functionspace(mesh, vector_element)

# Define physical and forcing parameters
alpha = 0.1
beta1 = 0.01
beta2 = 0.01
S1 = 0.005
P1 = pi / 180
S2 = 0.3
P2 = pi

# Mark left boundary (x = 0) with id 1 for boundary integration
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, 0, lambda x: np.isclose(x[0], 0))
boundary_domains = [(1, boundary_facets)]
ds = Measure("ds", domain=mesh, subdomain_data=boundary_domains)

# Define trial and test functions
u = TrialFunction(V_vec)
v = TestFunction(V_vec)
x = SpatialCoordinate(mesh)

# Prescribed Dirichlet boundary value at x = 0
n1 = as_matrix([[0.75], [0.25]])

# (Unused) alternative value, could be for right side
n2 = as_matrix([[0.25], [0.75]])

# Define system matrix A and source term f(x)
A = as_matrix([[alpha + beta1, -alpha], [-alpha, alpha + beta2]])
f = as_matrix([[beta1 + S1 * sin(P1 * x[0])], [S2 * sin(P2 * x[0])]])

# Variational form:
# - Volume term: inner product of A * u and v over domain
# - Jump term: weak coupling between elements
# - Strong enforcement of boundary condition at x=0 via large penalty on ds(1)
a = inner(dot(A, u), v) * dx + inner(jump(u), v('+')) * dS + 1e8 * inner(u, v) * ds(1)
L = inner(f, v) * dx + 1e8 * inner(n1, v) * ds(1)

# Solve the linear system
problem = petsc.LinearProblem(a, L, bcs=[], petsc_options={"ksp_type": "cg"})
uh = problem.solve()

# Extract degrees of freedom: u1 and u2 components
x_coords = mesh.geometry.x[:, 0]
u1 = uh.x.array[::2]
u2 = uh.x.array[1::2]

# Print values at the left boundary
print(u1[0], u2[0])

# Plot the solution components
plt.figure(figsize=(8, 5))
plt.plot(x_coords[:-1], u1, "r-", label="u1(x)", linewidth=2)
plt.plot(x_coords[:-1], u2, "b--", label="u2(x)", linewidth=2)
plt.xlabel("x")
plt.ylabel("Value")
plt.title("Components of Vector Function")
plt.grid(True)
plt.legend()
plt.show()

Best regards,
huanwei

5

Unfortunately I don’t have time to do this in detail for you, but I think you may need to audit your discontinuous Galerkin formulation. It doesn’t seem quite right. It’s crucial that the numerical flux of u can propagate through the domain correctly (the “flow of information”).

Some examples which may help:

6

You can also consider

for a dg advection diffusion problem with a divergence conforming velocity field.

7

Thank you, Nate and Dokken, for providing the valuable information and for the time you’ve invested—it’s greatly appreciated!

8

The weak form is modified to

a = inner(dot(A, u), v) * dx + inner(u('+'), jump(v)) * dS + 1e8*inner(u,v)*ds(1)
L = inner(f, v) * dx + 1e8*inner(n1,v)*ds(1)

where the last term enforces a strong boundary condition.
For visualization purposes, values on the right boundary are omitted, as no boundary condition is imposed there.

9

Just be careful. As far as I’m aware, u("+") is not guaranteed to be evaluated in the neighbouring cell of the facet you expect. And naturally this formulation won’t extend to higher dimensions.

10

To ensure the correct side being used, see for instance: fenics-in-the-wild/src/emi_primal_single.py at dokken/emi-models · scientificcomputing/fenics-in-the-wild · GitHub