Cannot apply Neumann BC to 1D Fenicsx problem

dolfinx
1

Greetings.

I am trying to solve the following 1D problem in Fenicsx. I have been successful in the past with the 2019 version but not yet successful with the 2023 version.

The problem is u_xx + x+1 = 0 with BCs: u’(0) = 0, u(1) = 0. This is what I have so far. It is the RHS term that is throwing up errors. Thank you for your attention.

import numpy as np

import ufl
from dolfinx import fem, io, mesh, plot
from ufl import ds, dx, grad, inner, Measure

import dolfinx.fem.petsc

from mpi4py import MPI
from petsc4py.PETSc import ScalarType

# 1. Define domain and mesh
domain = mesh.create_interval(comm=MPI.COMM_WORLD, points=(0.0, 1.0), nx=5)

# 2. Define function space (for shape functions)
V = fem.FunctionSpace(domain, ("Lagrange", 1))

# 3. Define boundary and boundary conditions

u1 = ScalarType(0)

def boundary_D(x):
    return np.isclose(x[0], 1)

dof1 = fem.locate_dofs_geometrical(V, boundary_D)
bc1 = fem.dirichletbc(value=ScalarType(u1), dofs=dof1, V=V)

up0 = 0

# 4. Define u, w function space and solve for uh
u = ufl.TrialFunction(V)
w = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(domain)
f = x[0]+1

#ds = Measure("ds", domain=mesh, subdomain_data=dof0)

LHS = inner(grad(w), grad(u)) * dx
RHS = inner(w, f) * dx - up0 * w * ds

problem = fem.petsc.LinearProblem(LHS, RHS, bcs=[bc1], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

This is the error I get:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[3], line 39
     36 #ds = Measure("ds", domain=mesh, subdomain_data=dof0)
     38 LHS = inner(grad(w), grad(u)) * dx
---> 39 RHS = inner(w, f) * dx - up0 * w * ds
     41 problem = fem.petsc.LinearProblem(LHS, RHS, bcs=[bc1], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
     42 uh = problem.solve()

File ~/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ufl/measure.py:386, in Measure.__rmul__(self, integrand)
    384     domain, = domains
    385 elif len(domains) == 0:
--> 386     raise ValueError("This integral is missing an integration domain.")
    387 else:
    388     raise ValueError("Multiple domains found, making the choice of integration domain ambiguous.")

ValueError: This integral is missing an integration domain.

It seems this has to do with how up0 * w * ds is being parsed by dolfinx. But I am not sure how to fix it.

2

I was able to fix the issue. Defining up0 as the following helped.

up0 = fem.Constant(domain, ScalarType(0))