Dirichletbc type error

I get the following type error when trying to define dirichletbc.

The code is a modified version of this example Setting multiple Dirichlet, Neumann, and Robin conditions

The mathematics:

Equation and problem definition

For a domain \Omega \subset \mathbb{R}^n with boundary \partial \Omega = \Gamma_{D} \cup \Gamma_{N} \cup \Gamma_{M}, the Poisson equation with
particular boundary conditions reads:

\nabla^{2}u =0 \quad {\rm in} \ \Omega, \\ u = 1 \quad {\rm on} \ \Gamma_{D}, \\ \nabla u \cdot n = 0 \quad {\rm on} \ \Gamma_{N}, \\ \nabla u \cdot n = h \quad {\rm on} \ \Gamma_{M}
a(u, v) := \int_{\Omega} \nabla u \cdot \nabla v \, {\rm d} x + \int_{\Gamma_{M}} r u v \, {\rm d} x \\ L(v) := \int_{\Omega} 0 v \, {\rm d} x + \int_{\Gamma_{M}} h v \, {\rm d} s + \int_{\Gamma_{N}} 0 v \, {\rm d} s.
  • \Omega = [0,2] \times [0,1] (a rectangle)
  • \Gamma_{D} = \{(x, 0) \subset \partial \Omega\}
  • \Gamma_{M} = \{(2, y) \cup(x, 1) \subset \partial \Omega\}
  • \Gamma_{N} = \{(0, y) \subset \partial \Omega\}
  • h = 0.1*(u-0) = r*(u-0)

The full code:

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, create_rectangle, CellType
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_rectangle(comm=MPI.COMM_WORLD,
                            points=((0.0, 0.0), (2.0, 1.0)), n=(32, 16),
                            cell_type=CellType.triangle)

u_ex =Constant(mesh, ScalarType(1)) #dirichlet
x = SpatialCoordinate(mesh)

# Define physical parameters and boundary condtions
s = Constant(mesh, ScalarType(0))
f = Constant(mesh, ScalarType(0))
n = FacetNormal(mesh)
g = Constant(mesh, ScalarType(0))
kappa = Constant(mesh, ScalarType(1))
r = Constant(mesh, ScalarType(0.1))

# 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[1], 0)),
              (2, lambda x: np.logical_or(np.isclose(x[0], 2),np.isclose(x[1], 1))),
              (3, lambda x: np.isclose(x[0], 0))]

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])

ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)

# use a simple class to define bc

class BoundaryCondition():
    def __init__(self, tipo, marker, values):
        self._type = tipo
        
        if tipo == "Dirichlet":
            facets = facet_tag.find(marker)
            dofs = locate_dofs_topological(V, fdim, facets)
            
            self._bc = dirichletbc(value=Constant(mesh, ScalarType(1)), dofs=dofs)
            
        elif tipo == "Neumann":
                self._bc = inner(values, v) * ds(marker)
        elif tipo == "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("Robin", 2, (r, s)),
                       BoundaryCondition("Neumann", 3, g)]

Thank you!

There is no function space input here as a third argument. See for instance: Component-wise Dirichlet BC — FEniCSx tutorial
or just the error message you posted in the screenshot above. The 1. and second supported input arguments has a V as the third input parameter.

Thank you so much for the previous solution.

I’m still quite confused by the functioning of the following example, a modification of the tutorial Stokes equations with Taylor-Hood elements in which I’m trying to impose different condition (parabolic inlet, no slip, pressure equal 0)

The following code do not raise error:

P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V, Q = FunctionSpace(mesh, P2), FunctionSpace(mesh, P1)

# parabolic inlet
u_in = Function(V)
u_in.interpolate(velocity_in)
where=ft.indices[ft.values==3]

dofs = locate_dofs_topological(V, fdim, where)
bc0 = dirichletbc(u_in, dofs)

# noslip
noslip = Constant(mesh,(0.0,0.0))
where=np.append(ft.indices[ft.values==1], ft.indices[ft.values==2])
dofs = locate_dofs_topological(V, fdim, where)
bc1 = dirichletbc(noslip, dofs, V)

# pressure=0 out
pres=Constant(mesh,(0.0))
where=ft.indices[ft.values==4]
dofs = locate_dofs_topological(Q, fdim, where)
bc2 = dirichletbc(pres, dofs, Q)

bcs = [bc0,bc1,bc2]

While if I change this line I get a error, even if it’s the same syntax:

bc1 = dirichletbc(u_in, dofs, V)

I can’t see why

When u_in is a function in the function space you are constraining, you dont need to supply the function space (as there is an identity map).
If you apply BCs to a subspace, or by using a constant, you have to tell the bc constructor what space you are constraining to

1 Like

Oh right, thank you!

I’m now having the following problem with assemble_matrix:
" RuntimeError: Diagonal sub-block (1, 1) cannot be ‘None’ and have DirichletBC applied. Consider assembling a zero block. "

previous working code, plus this:

# Define variational problem
(u, p) = TrialFunction(V), TrialFunction(Q)
(v, q) = TestFunction(V), TestFunction(Q)
f = Constant(mesh, (ScalarType(0), ScalarType(0)))

a = form([[inner(grad(u), grad(v)) * dx, inner(p, div(v)) * dx],
          [inner(div(u), q) * dx, None]])
L = form([inner(f, v) * dx, inner(Constant(mesh, ScalarType(0)), q) * dx])

I then run the nested_iterative_solver(), or block_iterative_solver(), or block_direct_solver() as written in Stokes equations using Taylor-Hood elements and get the error written above.

What am I doing wrong?

As the error states, you have no block relating the pressure test and trial functions, but you have a boundary condition on pressure. Thus you would have to supply a block of zeros instead of supplying None as the last argument in a.