# 2 Neumann Boundary Conditions

Hello

how do I add 2 Neumann boundary conditions?
My idea was:

u1 = ufl.VersuchsFunktion(V)
v1 = ufl.TestFunktion(V)
dObs = Measure("ds", domain=mesh, subdomain_data=ft, subdomain_id=obstacle2_marker)
dObt = Measure("ds", domain=mesh, subdomain_data=ft, subdomain_id=obstacle1_marker)
L = ufl.inner(E1,v1) * dObt + ufl.inner(E2, v1) * dObs
u1 = fem.Function(V)
problem1 = fem.petsc.LinearProblem(a, L, bcwall, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
p = problem1.solve()

but it does not work like this.

1 Like

But there are not 2 Neumann boundary conditions executed at the same time, are they?

I believe all 4 boundary conditions are executed, but I cannot yet understand how.

Is there no easier way on my example?

What your doing in your suggested example should work (up to some some spelling).

If you can provide a minimal working example that does not do what you would expect people can help you.

What do you mean by â€śup to some some spellingâ€ť?

I am trying to make a minimal working example.

These are German? names for the test and trial functon, and are thus not in the ufl namespace.

1 Like

Hi
I am trying to solve a problem with different materials.
The code i am using is:

from __future__ import print_function

from dolfin import *
from fenics import *
# from boxfield import *
import numpy as np
import matplotlib.pyplot as plt

def solver_bcs(kappa, f, boundary_conditions, Nx, Ny,
degree=1,
subdomains=[],
linear_solver='Krylov',
abs_tol=1E-5,
rel_tol=1E-3,
max_iter=1000):
"""
Solve -div(kappa*grad(u) = f on (0, 1) x (0, 1) with 2*Nx*Ny Lagrange
elements of specified degree and u = u_D on the boundary. This version
of the solver uses a specified combination of Dirichlet, Neumann, and
Robin boundary conditions.

"""

# Create mesh and define function space
mesh = UnitSquareMesh(Nx, Ny)
V = FunctionSpace(mesh, 'P', degree)

# Check if we have subdomains
if subdomains:
if not isinstance(kappa, (list, tuple, np.ndarray)):
raise TypeError(
'kappa must be array if we have sudomains, not %s'
% type(kappa))
materials = CellFunction('size_t', mesh)
materials.set_all(0)
for m, subdomain in enumerate(subdomains[1:], 1):
subdomain.mark(materials, m)

kappa_values = kappa
V0 = FunctionSpace(mesh, 'DG', 0)
kappa  = Function(V0)
help = np.asarray(materials.array(), dtype=np.int32)
kappa.vector()[:] = np.choose(help, kappa_values)
else:
if not isinstance(kappa, (Expression, Constant)):
raise TypeError(
'kappa is type %s, must be Expression or Constant'
% type(kappa))

# Define boundary subdomains
tol = 1e-14

class BoundaryX0(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0, tol)

class BoundaryX1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 1, tol)

class BoundaryY0(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0, tol)

class BoundaryY1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 1, tol)

# Mark boundaries
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
boundary_markers.set_all(9999)
bx0 = BoundaryX0()
bx1 = BoundaryX1()
by0 = BoundaryY0()
by1 = BoundaryY1()
bx0.mark(boundary_markers, 0)
bx1.mark(boundary_markers, 1)
by0.mark(boundary_markers, 2)
by1.mark(boundary_markers, 3)

# Redefine boundary integration measure
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

# Collect Dirichlet conditions
bcs = []
for i in boundary_conditions:
if 'Dirichlet' in boundary_conditions[i]:
bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'],
boundary_markers, i)
bcs.append(bc)

debug1 = False
if debug1:

# Print all vertices that belong to the boundary parts
for x in mesh.coordinates():
if bx0.inside(x, True): print('%s is on x = 0' % x)
if bx1.inside(x, True): print('%s is on x = 1' % x)
if by0.inside(x, True): print('%s is on y = 0' % x)
if by1.inside(x, True): print('%s is on y = 1' % x)

# Print the Dirichlet conditions
print('Number of Dirichlet conditions:', len(bcs))
if V.ufl_element().degree() == 1:  # P1 elements
d2v = dof_to_vertex_map(V)
coor = mesh.coordinates()
for i, bc in enumerate(bcs):
print('Dirichlet condition %d' % i)
boundary_values = bc.get_boundary_values()
for dof in boundary_values:
print('   dof %2d: u = %g' % (dof, boundary_values[dof]))
if V.ufl_element().degree() == 1:
print('    at point %s' %
(str(tuple(coor[d2v[dof]].tolist()))))

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)

# Collect Neumann integrals
integrals_N = []
for i in boundary_conditions:
if 'Neumann' in boundary_conditions[i]:
if boundary_conditions[i]['Neumann'] != 0:
g = boundary_conditions[i]['Neumann']
integrals_N.append(g*v*ds(i))

# Collect Robin integrals
integrals_R_a = []
integrals_R_L = []
for i in boundary_conditions:
if 'Robin' in boundary_conditions[i]:
r, s = boundary_conditions[i]['Robin']
integrals_R_a.append(r*u*v*ds(i))
integrals_R_L.append(r*s*v*ds(i))

# Simpler Robin integrals
integrals_R = []
for i in boundary_conditions:
if 'Robin' in boundary_conditions[i]:
r, s = boundary_conditions[i]['Robin']
integrals_R.append(r*(u - s)*v*ds(i))

# Sum integrals to define variational problem
L = f*v*dx - sum(integrals_N) + sum(integrals_R_L)

# Simpler variational problem
sum(integrals_R) - f*v*dx + sum(integrals_N)
a, L = lhs(F), rhs(F)

# Set linear solver parameters
if linear_solver == 'Krylov':
prm.linear_solver = 'gmres'
prm.preconditioner = 'ilu'
prm.krylov_solver.absolute_tolerance = abs_tol
prm.krylov_solver.relative_tolerance = rel_tol
prm.krylov_solver.maximum_iterations = max_iter
else:
prm.linear_solver = 'lu'

# Compute solution
u = Function(V)
solve(a == L, u, bcs, solver_parameters=prm)

return u

def demo_bcs():
"Compute and plot solution using a combination of boundary conditions"

# Define manufactured solution in sympy and derive f, g, etc.
import sympy as sym
x, y = sym.symbols('x[0], x[1]')            # needed by UFL
u = 1 + x**2 + 2*y**2                       # exact solution
u_e = u                                     # exact solution
u_00 = u.subs(x, 0)                         # restrict to x = 0
u_01 = u.subs(x, 1)                         # restrict to x = 1
f = -sym.diff(u, x, 2) - sym.diff(u, y, 2)  # -Laplace(u)
f = sym.simplify(f)                         # simplify f
g = -sym.diff(u, y).subs(y, 1)              # compute g = -du/dn
r = 1000                                    # Robin data, arbitrary
s = u                                       # Robin data, u = s

# Collect variables
variables = [u_e, u_00, u_01, f, g, r, s]

# Turn into C/C++ code strings
variables = [sym.printing.ccode(var) for var in variables]

# Turn into FEniCS Expressions
variables = [Expression(var, degree=2) for var in variables]

# Extract variables
u_e, u_00, u_01, f, g, r, s = variables

# Define boundary conditions
boundary_conditions = {0: {'Dirichlet': u_00},   # x = 0
1: {'Dirichlet': u_01},   # x = 1
2: {'Robin':     (r, s)}, # y = 0
3: {'Neumann':   g}}      # y = 1

# Set linear solver parameters
prm = LinearVariationalSolver.default_parameters()
if linear_solver == 'Krylov':
prm.linear_solver = 'gmres'
prm.preconditioner = 'ilu'
prm.krylov_solver.absolute_tolerance = abs_tol
prm.krylov_solver.relative_tolerance = rel_tol
prm.krylov_solver.maximum_iterations = max_iter
else:
prm.linear_solver = 'lu'

# Compute solution
kappa = Constant(1)
Nx = Ny = 8
u = solver_bcs(kappa, f, boundary_conditions, Nx, Ny,
degree=1, linear_solver='direct')

# Compute maximum error at vertices
mesh = u.function_space().mesh()
vertex_values_u_e = u_e.compute_vertex_values(mesh)
vertex_values_u = u.compute_vertex_values(mesh)
error_max = np.max(np.abs(vertex_values_u_e -
vertex_values_u))
print('error_max =', error_max)

# Save and plot solution
vtkfile = File('poisson_extended/solution_bcs.pvd')
vtkfile << u
plot(u)

demo_bcs()


However, I always get this Error:

NameError                                 Traceback (most recent call last)
Cell In[93], line 1
----> 1 demo_bcs()

Cell In[90], line 37, in demo_bcs()
35     # Set linear solver parameters
36 prm = LinearVariationalSolver.default_parameters()
---> 37 if linear_solver == 'Krylov':
38     prm.linear_solver = 'gmres'
39     prm.preconditioner = 'ilu'

NameError: name 'linear_solver' is not defined


You never send linear_solver as an argument to demo_bcs(), and thus the variable does not exist.

Thank you for your answer, I am actually running the solver_bcs and demo_bcs functions in ft10_possion_extended and The code i am using is shown in the code below:

from __future__ import print_function

from fenics import *
from boxfield import *
import numpy as np

def solver_bcs(kappa, f, boundary_conditions, Nx, Ny,
degree=1,
subdomains=[],
linear_solver='Krylov',
abs_tol=1E-5,
rel_tol=1E-3,
max_iter=1000):
"""
Solve -div(kappa*grad(u) = f on (0, 1) x (0, 1) with 2*Nx*Ny Lagrange
elements of specified degree and u = u_D on the boundary. This version
of the solver uses a specified combination of Dirichlet, Neumann, and
Robin boundary conditions.

"""

# Create mesh and define function space
mesh = UnitSquareMesh(Nx, Ny)
V = FunctionSpace(mesh, 'P', degree)

# Check if we have subdomains
if subdomains:
if not isinstance(kappa, (list, tuple, np.ndarray)):
raise TypeError(
'kappa must be array if we have sudomains, not %s'
% type(kappa))
materials = CellFunction('size_t', mesh)
materials.set_all(0)
for m, subdomain in enumerate(subdomains[1:], 1):
subdomain.mark(materials, m)

kappa_values = kappa
V0 = FunctionSpace(mesh, 'DG', 0)
kappa  = Function(V0)
help = np.asarray(materials.array(), dtype=np.int32)
kappa.vector()[:] = np.choose(help, kappa_values)
else:
if not isinstance(kappa, (Expression, Constant)):
raise TypeError(
'kappa is type %s, must be Expression or Constant'
% type(kappa))

# Define boundary subdomains
tol = 1e-14

class BoundaryX0(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0, tol)

class BoundaryX1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 1, tol)

class BoundaryY0(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0, tol)

class BoundaryY1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 1, tol)

# Mark boundaries
# boundary_markers = FacetFunction('size_t', mesh)
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)######################äż®ć”ą
boundary_markers.set_all(9999)
bx0 = BoundaryX0()
bx1 = BoundaryX1()
by0 = BoundaryY0()
by1 = BoundaryY1()
bx0.mark(boundary_markers, 0)
bx1.mark(boundary_markers, 1)
by0.mark(boundary_markers, 2)
by1.mark(boundary_markers, 3)

# Redefine boundary integration measure
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

# Collect Dirichlet conditions
bcs = []
for i in boundary_conditions:
if 'Dirichlet' in boundary_conditions[i]:
bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'],
boundary_markers, i)
bcs.append(bc)

debug1 = False
if debug1:

# Print all vertices that belong to the boundary parts
for x in mesh.coordinates():
if bx0.inside(x, True): print('%s is on x = 0' % x)
if bx1.inside(x, True): print('%s is on x = 1' % x)
if by0.inside(x, True): print('%s is on y = 0' % x)
if by1.inside(x, True): print('%s is on y = 1' % x)

# Print the Dirichlet conditions
print('Number of Dirichlet conditions:', len(bcs))
if V.ufl_element().degree() == 1:  # P1 elements
d2v = dof_to_vertex_map(V)
coor = mesh.coordinates()
for i, bc in enumerate(bcs):
print('Dirichlet condition %d' % i)
boundary_values = bc.get_boundary_values()
for dof in boundary_values:
print('   dof %2d: u = %g' % (dof, boundary_values[dof]))
if V.ufl_element().degree() == 1:
print('    at point %s' %
(str(tuple(coor[d2v[dof]].tolist()))))

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)

# Collect Neumann integrals
integrals_N = []
for i in boundary_conditions:
if 'Neumann' in boundary_conditions[i]:
if boundary_conditions[i]['Neumann'] != 0:
g = boundary_conditions[i]['Neumann']
integrals_N.append(g*v*ds(i))

# Collect Robin integrals
integrals_R_a = []
integrals_R_L = []
for i in boundary_conditions:
if 'Robin' in boundary_conditions[i]:
r, s = boundary_conditions[i]['Robin']
integrals_R_a.append(r*u*v*ds(i))
integrals_R_L.append(r*s*v*ds(i))

# Simpler Robin integrals
integrals_R = []
for i in boundary_conditions:
if 'Robin' in boundary_conditions[i]:
r, s = boundary_conditions[i]['Robin']
integrals_R.append(r*(u - s)*v*ds(i))

# Sum integrals to define variational problem
L = f*v*dx - sum(integrals_N) + sum(integrals_R_L)

# Simpler variational problem
sum(integrals_R) - f*v*dx + sum(integrals_N)
a, L = lhs(F), rhs(F)

# Set linear solver parameters
if linear_solver == 'Krylov':
prm.linear_solver = 'gmres'
prm.preconditioner = 'ilu'
prm.krylov_solver.absolute_tolerance = abs_tol
prm.krylov_solver.relative_tolerance = rel_tol
prm.krylov_solver.maximum_iterations = max_iter
else:
prm.linear_solver = 'lu'

# Compute solution
u = Function(V)
solve(a == L, u, bcs, solver_parameters=prm)

return u

def demo_bcs():
"Compute and plot solution using a combination of boundary conditions"

# Define manufactured solution in sympy and derive f, g, etc.
import sympy as sym
x, y = sym.symbols('x[0], x[1]')            # needed by UFL
u = 1 + x**2 + 2*y**2                       # exact solution
u_e = u                                     # exact solution
u_00 = u.subs(x, 0)                         # restrict to x = 0
u_01 = u.subs(x, 1)                         # restrict to x = 1
f = -sym.diff(u, x, 2) - sym.diff(u, y, 2)  # -Laplace(u)
f = sym.simplify(f)                         # simplify f
g = -sym.diff(u, y).subs(y, 1)              # compute g = -du/dn
r = 1000                                    # Robin data, arbitrary
s = u                                       # Robin data, u = s

# Collect variables
variables = [u_e, u_00, u_01, f, g, r, s]

# Turn into C/C++ code strings
variables = [sym.printing.ccode(var) for var in variables]

# Turn into FEniCS Expressions
variables = [Expression(var, degree=2) for var in variables]

# Extract variables
u_e, u_00, u_01, f, g, r, s = variables

# Define boundary conditions
boundary_conditions = {0: {'Dirichlet': u_00},   # x = 0
1: {'Dirichlet': u_01},   # x = 1
2: {'Robin':     (r, s)}, # y = 0
3: {'Neumann':   g}}      # y = 1

# Compute solution
kappa = Constant(1)
Nx = Ny = 8
u = solver_bcs(kappa, f, boundary_conditions, Nx, Ny,
degree=1, linear_solver='direct')

# Compute maximum error at vertices
mesh = u.function_space().mesh()
vertex_values_u_e = u_e.compute_vertex_values(mesh)
vertex_values_u = u.compute_vertex_values(mesh)
error_max = np.max(np.abs(vertex_values_u_e -
vertex_values_u))
print('error_max =', error_max)

# Save and plot solution
vtkfile = File('poisson_extended/solution_bcs.pvd')
vtkfile << u
plot(u)

demo_bcs()


The first problem I am having is shown below:

ModuleNotFoundError                       Traceback (most recent call last)
Cell In[52], line 4
1 from __future__ import print_function
3 from fenics import *
----> 4 from boxfield import *
5 import numpy as np
7 def solver_bcs(kappa, f, boundary_conditions, Nx, Ny,
8                degree=1,
9                subdomains=[],
(...)
12                rel_tol=1E-3,
13                max_iter=1000):

ModuleNotFoundError: No module named 'boxfield'

NameError                                 Traceback (most recent call last)
Cell In[51], line 1
----> 1 demo_bcs()

Cell In[50], line 38, in demo_bcs()
36 kappa = Constant(1)
37 Nx = Ny = 8
---> 38 u = solver_bcs(kappa, f, boundary_conditions, Nx, Ny,
39                degree=1, linear_solver='direct')
41 # Compute maximum error at vertices
42 mesh = u.function_space().mesh()

Cell In[49], line 162, in solver_bcs(kappa, f, boundary_conditions, Nx, Ny, degree, subdomains, linear_solver, abs_tol, rel_tol, max_iter)
160     prm.krylov_solver.maximum_iterations = max_iter
161 else:
--> 162     prm.linear_solver = 'lu'
164 # Compute solution
165 u = Function(V)

NameError: name 'prm' is not defined


You need to install the boxfield package. Please note that boxfield and legacy FEniCS is deprecated.

Similar error as the first one, you have not defined prm inside solver_bcs. Please look carefully through your code before posting, as the error message literally tells you what is wrong.