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)
a = ufl.inner(ufl.grad(u1), ufl.grad(v1)) * ufl.dx 
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.

See for instance Setting multiple Dirichlet, Neumann, and Robin conditions — FEniCSx tutorial

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
    a = kappa*dot(grad(u), grad(v))*dx + sum(integrals_R_a)
    L = f*v*dx - sum(integrals_N) + sum(integrals_R_L)

    # Simpler variational problem
    F = kappa*dot(grad(u), grad(v))*dx + \
        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
    a = kappa*dot(grad(u), grad(v))*dx + sum(integrals_R_a)
    L = f*v*dx - sum(integrals_N) + sum(integrals_R_L)

    # Simpler variational problem
    F = kappa*dot(grad(u), grad(v))*dx + \
        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.