Applying constant vector for dirichletBC in dolfinx

dolfinx
1

Hi,

In old FEniCS, we could apply DirichletBC using constant vector as folllows.

from dolfin import *

mesh = UnitSquareMesh(20, 20)
# Build function space
P2 = VectorElement("CG", mesh.ufl_cell(), 2)
P1 = FiniteElement("CG", mesh.ufl_cell(), 1)
TH = P2 * P1
VQ = FunctionSpace(mesh, TH)
bc = DirichletBC(VQ.sub(0), Expression(("1", "0"), degree=0), "std::abs(x[1])>1-1e-12 && on_boundary")

I would like to do the same thing in dolfinx, but I could only find an example of using either scalar or function in the tutorial.

How could I apply constant vector that has different value (x-component is 1 and y-component is 0) using dolfinx?

1 Like
2

Not sure if this is completely correct, but it may give you some useful hints on how to solve the problem:

import numpy as np
from dolfinx import mesh
from dolfinx.fem import Constant, Function
from dolfinx.fem import FunctionSpace, dirichletbc
from dolfinx.fem locate_dofs_topological
from dolfinx.mesh import locate_entities_boundary
from mpi4py import MPI
from petsc4py import PETSc
from ufl import FiniteElement, MixedElement, VectorElement

domain = mesh.create_unit_square(
    MPI.COMM_WORLD, 20, 20, mesh.CellType.quadrilateral)

P2 = VectorElement("CG", domain.ufl_cell(), 2)
P2_space = FunctionSpace(domain, P2)

P1 = FiniteElement("CG", domain.ufl_cell(), 1)
TH = FunctionSpace(domain, MixedElement([P2, P1]))

def on_boundary(x):
    return np.full(x.shape[1], True)

def bc_function(x):
    return np.logical_and(on_boundary, np.abs(x[1]) > 1-1e-12)

boundary_facets = locate_entities_boundary(
    domain, domain.topology.dim-1, bc_function)

boundary_dofs = locate_dofs_topological(
    P2_space, domain.topology.dim-1, boundary_facets)

f = Constant(domain, PETSc.ScalarType((0.0, 1.0)))
bc = dirichletbc(f, boundary_dofs, P2_space)
1 Like
3

What you are doing here is actually not applying a Constant in the meaning you are thinking. The expression has to be interpolated into the space (CG-2). So even if you are sending in a constant value, it computes this as a function under the hood.

Other than that, it is worth noting that there are only certain sub-spaces that support constant valued bcs.
Under the hood of DOLFINx, there is only a minor memory increase (size of a vector) when using functions instead of a constant. It does not affect the actual assembly time.
To use a function you could:

V, _ = VQ.sub(0).collapse(0)
u_bc = dolfinx.fem.Function(V)
def c(x):
    vals = np.zeros((2, x.shape[1]))
    vals[0] = 1
    return vals
u_bc.interpolate(c)
2 Likes
4

You can also use Expression in DOLFINx for this purpose:

u_bc = dolfinx.fem.Function(V)
bc_val = dolfinx.fem.Constant(mesh, PETSc.ScalarType([1,0]))
bc_expr = dolfinx.fem.Expression(bc_val, V.element.interpolation_points())
u_bc.interpolate(bc_expr)
3 Likes
5

Hi,

Using the suggested approach indeed helps defining the dirichlet boundary with a Constant (bc_val) but it does not seem to update u_bc when I redefine bc_val.value. If this is the expected behavior then is there a way to accomplish what I want?

Below is an MWE evaluated on 0.8.0 in which I try to do it the way I thought is right except it does not work:

import dolfinx, gmsh, ufl, mpi4py, basix, numpy as np

gdim = 2
mesh = dolfinx.mesh.create_unit_square(mpi4py.MPI.COMM_WORLD, 10, 10)
mesh.topology.create_connectivity(gdim-1, gdim)

P1 = basix.ufl.element("Lagrange", mesh.basix_cell(), 1)
P2 = basix.ufl.element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))
ME = basix.ufl.mixed_element([P2, P1]) # mixed space
W = dolfinx.fem.functionspace(mesh, ME)
V2, _ = W.sub(0).collapse()
V1, _ = W.sub(1).collapse()

# apply dirichlet boundary at the inlet
inlet_marker = lambda x: np.isclose(x[0], 0)
inlet_value = dolfinx.fem.Constant(mesh, np.array((1.0, 0), dtype=dolfinx.default_scalar_type))
inlet_dofs = dolfinx.fem.locate_dofs_geometrical((W.sub(0), V2), inlet_marker)
dbc = [dolfinx.fem.dirichletbc(inlet_value, inlet_dofs)]
6

I do not understand what is not working in this post, please add the remaining code that illustrates your issue.