Boundary conditions mixed elements

Linear Algebra
1

I am trying to impose boundary conditions in a special way on a MixedElement space.

Consider : Lq=Bf. Imposing L.zeroRows(dofs) and B.zeroRows(dofs,0) should be equivalent to enforcing q[dofs]=0. Indeed PETSc garantees L[dofs,dofs]=1.

A toy problem follows. I’m trying to solve on unit 2D square \Omega using mixed element space (u,v) :
\begin{cases}\Delta u=1\text{ with }u_{\partial\Omega}=0\\\Delta v=0\text{ with }v_{\partial\Omega}=0\end{cases}

Notice that v=0 is a trivial solution to the second equation and the problem is decoupled. I manage to get the first part working on its own, not as part of an element space.

Component wise Dirichlet BC is taken from this tutorial.

import ufl
import numpy as np
import dolfinx as dfx
from petsc4py import PETSc as pet
from dolfinx.io import XDMFFile
from mpi4py.MPI import COMM_WORLD

mesh=dfx.UnitSquareMesh(COMM_WORLD,100,100)

FE_scalar=ufl.FiniteElement("Lagrange",mesh.ufl_cell(),1)
V=dfx.FunctionSpace(mesh,FE_scalar*FE_scalar)
U, _ = V.sub(0).collapse(collapsed_dofs=True)

u,v=ufl.TrialFunction(V),ufl.TestFunction(V)
w,z=ufl.TrialFunction(U),ufl.TestFunction(U)
v1,_=ufl.split(v)

form=ufl.inner(ufl.grad(u),ufl.grad(v))*ufl.dx
B_form=ufl.inner(w,v1)*ufl.dx

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

dofs=np.empty(0,dtype=np.int32)
for i in range(2):
	U_i =V.sub(i)
	U_ic=U_i.collapse()
	# Compute DoFs
	np.union1d(dofs,dfx.fem.locate_dofs_geometrical((U_i, U_ic), boundary))

L = dfx.fem.assemble_matrix(form);   L.assemble(); L.zeroRows(dofs)
B = dfx.fem.assemble_matrix(B_form); B.assemble(); B.zeroRows(dofs,0)

KSP = pet.KSP().create()
KSP.setOperators(L)
KSP.setFromOptions()

x,y=B.createVecs()
f=dfx.Function(V)
x.set(1)
B.mult(x,y)
KSP.solve(y,f.vector)

f1,f2=f.split()
with XDMFFile(COMM_WORLD, "sanity_check_bcs_1.xdmf","w") as xdmf:
	xdmf.write_mesh(mesh)
	xdmf.write_function(f1)
with XDMFFile(COMM_WORLD, "sanity_check_bcs_2.xdmf","w") as xdmf:
	xdmf.write_mesh(mesh)
	xdmf.write_function(f2)

This code works and gives f2=0, but produces a f1 that does not satisfy f1=0 on the boundary.

f1
f12252×1194 87.3 KB

Did I miss anything ? I know there are other ways to impose boundary conditions in FENICsx but I would prefer to use B here.

1 Like
2

Consider essentially the same code without a MixedElement :

import ufl
import numpy as np
import dolfinx as dfx
from petsc4py import PETSc as pet
from dolfinx.io import XDMFFile
from mpi4py.MPI import COMM_WORLD

mesh=dfx.UnitSquareMesh(COMM_WORLD,100,100)

FE_scalar=ufl.FiniteElement("Lagrange",mesh.ufl_cell(),1)
V=dfx.FunctionSpace(mesh,FE_scalar)

u,v=ufl.TrialFunction(V),ufl.TestFunction(V)

form=ufl.inner(ufl.grad(u),ufl.grad(v))*ufl.dx
B_form=ufl.inner(u,v)*ufl.dx

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

dofs=dfx.fem.locate_dofs_geometrical(V, boundary)

L = dfx.fem.assemble_matrix(form);   L.assemble(); L.zeroRows(dofs)
B = dfx.fem.assemble_matrix(B_form); B.assemble(); B.zeroRows(dofs,0)

KSP = pet.KSP().create()
KSP.setOperators(L)
KSP.setFromOptions()

x,y=B.createVecs()
f=dfx.Function(V)
x.set(1)
B.mult(x,y)
KSP.solve(y,f.vector)

with XDMFFile(COMM_WORLD, "sanity_check_bcs.xdmf","w") as xdmf:
	xdmf.write_mesh(mesh)
	xdmf.write_function(f)

This one works with no boundary problem :

f
f1121×1194 82.9 KB

3

I guess I’m going to answer my own question. There are two problems with the original code :

  1. Add dofs= in front of np.union1d
  2. Only take the first part of the dfx.fem.locate_dofs_geometrical