I had some indications that in one of my dolfinx-scripts the Neumann BC is not fulfilled, so I barely modified one of the tutorials from @dokken to create a MWE and was surprised to find the same indication there as well: The Neumann BC does not seem to be fulfilled.
MWE (see Setting multiple Dirichlet, Neumann, and Robin conditions — FEniCSx tutorial):
import numpy as np
import pyvista
from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction,
div, dot, dx, grad, inner, lhs, rhs)
from dolfinx.io import XDMFFile
from dolfinx.plot import create_vtk_mesh
mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
u_ex = lambda x: 1 + x[0]**2 + 2*x[1]**2
x = SpatialCoordinate(mesh)
# Define physical parameters and boundary condtions
s = u_ex(x)
f = -div(grad(u_ex(x)))
n = FacetNormal(mesh)
g = -dot(n, grad(u_ex(x)))
# g = 0.0
kappa = Constant(mesh, ScalarType(1))
r = Constant(mesh, ScalarType(1000))
# Define function space and standard part of variational form
V = FunctionSpace(mesh, ("CG", 1))
u, v = TrialFunction(V), TestFunction(V)
F = kappa * inner(grad(u), grad(v)) * dx - inner(f, v) * dx
boundaries = [(1, lambda x: np.isclose(x[0], 0)),
(2, lambda x: np.isclose(x[0], 1)),
(3, lambda x: np.isclose(x[1], 0)),
(4, lambda x: np.isclose(x[1], 1))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, 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 = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
with XDMFFile(mesh.comm, "output/facet_tags.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_meshtags(facet_tag)
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
class BoundaryCondition():
def __init__(self, type, marker, values):
self._type = type
if type == "Dirichlet":
u_D = Function(V)
u_D.interpolate(values)
facets = facet_tag.find(marker)
dofs = locate_dofs_topological(V, fdim, facets)
self._bc = dirichletbc(u_D, dofs)
elif type == "Neumann":
self._bc = inner(values, v) * ds(marker)
elif type == "Robin":
self._bc = values[0] * inner(u-values[1], v)* ds(marker)
else:
raise TypeError("Unknown boundary condition: {0:s}".format(type))
@property
def bc(self):
return self._bc
@property
def type(self):
return self._type
# Define the Dirichlet condition
boundary_conditions = [BoundaryCondition("Dirichlet", 1, u_ex),
BoundaryCondition("Dirichlet", 2, u_ex),
BoundaryCondition("Robin", 3, (r, s)),
BoundaryCondition("Neumann", 4, g)]
bcs = []
for condition in boundary_conditions:
if condition.type == "Dirichlet":
bcs.append(condition.bc)
else:
F += condition.bc
# Solve linear variational problem
a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
solution_trace_norm = assemble_scalar(form((dot(n, grad(uh))-g)**2*ds(4)))**0.5
print(f"Neumann BC trace: {solution_trace_norm}")
The output I receive should be close to zero, but is:
Neumann BC trace: 7.800058502092114
Even if I set the Neumann BC to g=0, the output is still too high with
Neumann BC trace: 1.366936098344996
I am using the version fenics-dolfinx 0.6.0 via conda-forge on Ubuntu via WSL2.
Is this expected behaviour?
Did I do any mistakes evaluating (last line of the code) the Neumann BC?