Boundary conditions about Mixed space

General
1

I have a question on how to apply Dirichlet condition in Dolfinx for a mixed finte element scheme.The problem is as follows:

The function space I defined is

Ve=ufl.VectorElement("CG",domain.ufl_cell(),2)
Pe=ufl.FiniteElement("CG",domain.ufl_cell(),1)
We=ufl.MixedElement([Ve,Pe])
W=fem.FunctionSpace(domain,We)
w=ufl.TrialFunction(W)
u,p=ufl.split(w)

where domain is a 2d unit square.Now I want to apply Dirichlet condition to u such that u = 0 on x=0,x=1 and y=0 while u=1 on y=1 .I implement this in this way:

def no_slip(x):
    return np.logical_or(np.logical_or(np.isclose(x[0],0),np.isclose(x[0],1)),
                         np.isclose(x[1],0))
noslip_v=np.zeros(domain.geometry.dim,dtype=PETSc.ScalarType)
facets=mesh.locate_entities_boundary(domain,1,no_slip)

bc_noslip=fem.dirichletbc(noslip_v,fem.locate_dofs_topological(W.sub(0),1,facets),W.sub(0))

and it sends the error like

Invoked with: array([0., 0.]), array([ 0,  1,  2,  3,  4,  5,  6,  7, 10, 11, 22, 23, 26, 27, 29, 30, 31,
       32, 36, 37, 40, 41, 52, 53, 56, 57], dtype=int32), FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), VectorElement(FiniteElement('Lagrange', triangle, 2), dim=2))
RuntimeError: Creating a DirichletBC using a Constant is not supported when the Constant size is not equal to the block size of the constrained (sub-)space. Use a fem::Function to create the fem::DirichletBC.

I think noslip_v is dim =2 and W.sub(0) is also dim =2, so why do they not match?
And maybe one way to solve this is to use a Function used in dirichlebc() and how?But I wonder since dirichletbc() can accept a Constant in its first parameter, is there a way in which I can still use np.zeros as above in dirichletbc() without errors?(Constant in dirichletbc cannot work in MixedElement space?)

2

Please make the example reproducible. Currently you have only posted snippets of code, which doesn’t make it easy for people to reproduce the problem.

An example of how to make your code reproducible is:

import dolfinx
from mpi4py import MPI
import ufl
from petsc4py import PETSc
import numpy as np
domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
Ve = ufl.VectorElement("CG", domain.ufl_cell(), 2)
Pe = ufl.FiniteElement("CG", domain.ufl_cell(), 1)
We = ufl.MixedElement([Ve, Pe])
W = dolfinx.fem.FunctionSpace(domain, We)


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


noslip_v = np.zeros(domain.geometry.dim, dtype=PETSc.ScalarType)
facets = dolfinx.mesh.locate_entities_boundary(domain, 1, no_slip)
bc_noslip = dolfinx.fem.dirichletbc(
    noslip_v, dolfinx.fem.locate_dofs_topological(W.sub(0), 1, facets), W.sub(0))

The issue here is that the block size of W.sub(0) is 1. You can check this by looking at W.sub(0).dofmap.bs. This means that it has unrolled all it’s dofs with respect to the dimension of the vector space.

To create a bc for your case, I would:

from IPython import embed
import dolfinx
from mpi4py import MPI
import ufl
from petsc4py import PETSc
import numpy as np
domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
Ve = ufl.VectorElement("CG", domain.ufl_cell(), 2)
Pe = ufl.FiniteElement("CG", domain.ufl_cell(), 1)
We = ufl.MixedElement([Ve, Pe])
W = dolfinx.fem.FunctionSpace(domain, We)


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


# Set zero BC
noslip_v = np.zeros(domain.geometry.dim, dtype=PETSc.ScalarType)
V, _ = W.sub(0).collapse()
u_bc = dolfinx.fem.Function(V)
u_bc.x.set(0)
facets = dolfinx.mesh.locate_entities_boundary(domain, 1, no_slip)
bc_noslip = dolfinx.fem.dirichletbc(
    u_bc, dolfinx.fem.locate_dofs_topological((W.sub(0), V), 1, facets), W.sub(0))


w = dolfinx.fem.Function(W)
# Fill with non-zero data to illustrate input
w.x.set(2)

dolfinx.fem.petsc.set_bc(w.vector, [bc_noslip])

with dolfinx.io.VTXWriter(domain.comm, "output.bp", [w.sub(0).collapse()]) as vtx:
    vtx.write(0)
3

Thank you for reply and sorry for the incomplete codes.
If I understand right, when using MixedElement(), a constant ,for exmaple, like here noslip_v, is just a constant of block 2 but W has two spaces with each one being 1 block.Instead, if I just use V=VectorElement(…), and supply noslip here to dirichletbc() as follows

noslip_v = np.zeros(domain.geometry.dim, dtype=PETSc.ScalarType)
facets = dolfinx.mesh.locate_entities_boundary(domain, 1, no_slip)
bc_noslip = dolfinx.fem.dirichletbc(
    noslip_v, dolfinx.fem.locate_dofs_topological(V ,1, facets), V)

then it will be okay since noslip_v and V have the same block?

4

I’m not sure what you are trying to do in this code.
V is the collapsed sub-space of W. It has no relation to W, and cannot be used in a matrix whose unknown is from W. That is why I created the boundary condition above on W.sub(0).

5

Thanks.I just want to compare different ways of implementing dirichletbc in mixed finite element method.

By this I mean if at beginning I use

Ve=ufl.VectorElement("CG",domain.ufl_cell(),2)
V=fem.FunctionSpace(domain,Ve)

then I can supply noslip_v=np.zeros(domain.geometry.dim,dtype=PETSc.ScalarType) directly to dirichletbc like dirichletbc(noslip_v,...,V) since this time nothing is unmatched.

6

Yes, if you work on the vector element outside the mixed element context, you can use constants to assign boundary conditions. See for instance: dolfinx/test_bcs.py at main · FEniCS/dolfinx · GitHub

7

Thanks a lot.I am much more clear about it.