Problems about setting boundary conditions with two different ways

dolfinx
1

Hi all,

I am using trying to set Dirichlet boundary conditions in a 2D domain, for example, I wanna fix the displacement u_x and u_y on the bottom. Here I tried two different ways to implement this step:

import numpy as np

import ufl
from dolfinx import fem, io, mesh, plot
from ufl import ds, dx, grad, inner
from dolfinx.io import XDMFFile

from dolfinx.fem import FunctionSpace, Function, Constant
from ufl import SpatialCoordinate, FiniteElement, MixedElement, TestFunctions, split

from mpi4py import MPI
from petsc4py.PETSc import ScalarType

Lx, Ly = 1.0, 1.0
nx, ny = 2, 2
# define the domain with quadrilateral mesh
msh = mesh.create_rectangle(comm=MPI.COMM_WORLD, 
                            points=((0.0, 0.0), (Lx, Ly)), 
                            n=(nx, ny),
                            cell_type=mesh.CellType.quadrilateral)

vector_element = ufl.VectorElement("CG", msh.ufl_cell(), 1) # vector displacement 
scalar_element = ufl.FiniteElement("CG", msh.ufl_cell(), 1) # scalar pressure

element = ufl.MixedElement(vector_element, scalar_element)
mixed_space = fem.FunctionSpace(msh, element)
M = mixed_space

V_displacement = M.sub(0)  # displacement
P_pressure = M.sub(1)  # pressure

ndim = msh.geometry.dim
fdim = ndim -1

# msh.topology.create_connectivity(fdim, msh.topology.dim)

def bottom(x):
    return np.isclose(x[1], 0, 1e-6)

def top(x):
    return np.isclose(x[1], Ly, 1e-6)

The first way is:

V, submap = V_displacement.collapse()

dofs_bottom = fem.locate_dofs_geometrical(V, bottom)
u_xy_zero = np.array((0,)*ndim, dtype=ScalarType)
BC_bottom_xy= fem.dirichletbc(u_xy_zero, dofs_bottom, V)

which could be run in dolfinx

The second way fails to interpolate the value on the bottom, could someone tell me why?

u1_bc = fem.Function(V)
u1 = lambda x: np.zeros_like(x, dtype=ScalarType)
u1_bc.interpolate(u1)
2

If the domain is 3D, then u1 could be successfully interpolated to the function u1_bc.

3

The input x to this function is always (3, number of points). Thus for a vector function in 2D this should be
u1 = lambda x: np.zeros((2, x.shape[1]), dtype=ScalarType)
You could also encode it with
u1 = lambda x: np.zeros((msg.geometry.dim, x.shape[1]), dtype=ScalarType)

4

Thanks dokken, could I also write in the following way?

u1_bc = fem.Function(V)
u1 = lambda x: np.array([[0], [0]], dtype =ScalarType)
u1_bc.interpolate(u1)

so that I can define any value for u_x and u_y on the boundary.

5

No.
The function you send i to interpolate should be vectorized, ie the input is (3, num_points), and output of shape (2, num_points) for a vector space in 2D

7

got you, then how about this one to set u_y as 1 on the top in 2D:

u1_bc = fem.Function(V)
def set_uy_top(x):
    uy_top = np.zeros((2, x.shape[1]))
    uy_top[1] = 1
    return ScalarType(uy_top)
u1_bc.interpolate(set_uy_top)

and the applied value is:

image
image896×518 18.5 KB

8

I dont understand your question here. You have interpolated (0,1) into your space.

9

Sorry, the following is the complete code for impose the bcs:

Lx, Ly  = 1.0, 1.0
nx, ny  = 4, 4
# define the domain with quadrilateral mesh
msh = mesh.create_rectangle(comm=MPI.COMM_WORLD, 
                            points=((0.0, 0.0), (Lx, Ly)), 
                            n=(nx, ny),
                            cell_type=mesh.CellType.quadrilateral)

vector_element = ufl.VectorElement("CG", msh.ufl_cell(), 1) # vector displacement 
scalar_element = ufl.FiniteElement("CG", msh.ufl_cell(), 1) # scalar pressure

element = ufl.MixedElement(vector_element, scalar_element)
mixed_space = fem.FunctionSpace(msh, element)
M = mixed_space

V_displacement = M.sub(0)  # displacement
P_pressure = M.sub(1)  # pressue

ndim = msh.geometry.dim
fdim = ndim -1

def bottom(x):
    return np.isclose(x[1], 0, 1e-6)

def top(x):
    return np.isclose(x[1], Ly, 1e-6)

V, submap = V_displacement.collapse()

dofs_bottom = fem.locate_dofs_geometrical(V, bottom)
u_xy_zero = np.array((0,)*ndim, dtype=ScalarType)
BC_bottom_xy= fem.dirichletbc(u_xy_zero, dofs_bottom,V)

u1_bc = fem.Function(V)
def set_uy_top(x):
    uy_top = np.zeros((ndim,x.shape[1]))
    uy_top[1] = 1
    return ScalarType(uy_top)
u1_bc.interpolate(set_uy_top)

dofs_top = fem.locate_dofs_geometrical((M.sub(0), V), top)
BC_top_y = fem.dirichletbc(u1_bc, dofs_top, M.sub(0).sub(1))

bcs = [BC_bottom_xy, BC_top_y]