# How to set 'boundary condition' inside the boundary?

Hello,
I am solving the coupling equation for displacement(u) and (p).
The current topic is “Cryer’s problem” with such boundary conditions :

• u(r=0) = 0
• p(r=R) = 0
where R is the radius of the sphere.

I could deal with pressure at the boundary, as it is laterally boundary condition.
However, for the center, I am not sure how to implement this ‘boundary condition’, u(r=0)=0

What I have tried is as follows :

# Define function spaces
u_elem     = VectorElement('CG', domain.ufl_cell(), 2) # displacement(2D)
p_elem     = FiniteElement('CG', domain.ufl_cell(), 1) # pore water pressure(scalar)
V    = FunctionSpace(domain, MixedElement([u_elem, p_elem])) #Define vector function space for [u,p]

#Set space for u and p from V
U = FunctionSpace(domain, u_elem)
P = FunctionSpace(domain, p_elem)

def center(x):
return x[0]**2 + x[1]**2 < (0.1*R)**2

cells_center = locate_entities(domain, dim, center)

#Dirichlet BC for u : center
U, _ = V.sub(0).collapse()
u_bc = Function(U)
u_bc.x.set(0)

#center : u = 0
boundary_dofs_u_center = locate_dofs_topological((V.sub(0), U), dim, cells_center)
BC_u_center = dirichletbc(np.full_like(boundary_dofs_u_center, 0, dtype=ScalarType)[0], boundary_dofs_u_center, V.sub(0))


which followed the demo code of ‘’ Defining subdomains for different materials".

The error message is as follows :

TypeError                                 Traceback (most recent call last)
File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/bcs.py:125, in DirichletBCMetaClass.__init__(self, value, dofs, V)
124 try:
--> 125     super().__init__(_value, dofs, V)  # type: ignore
126 except TypeError:

TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
1. dolfinx.cpp.fem.DirichletBC_float64(g: numpy.ndarray[numpy.float64], dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace)
2. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Constant<double>, dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace)
3. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double>, dofs: numpy.ndarray[numpy.int32])
4. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double>, dofs: List[numpy.ndarray[numpy.int32][2]], V: dolfinx::fem::FunctionSpace)

Invoked with: array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0.]), [array([1796, 1797, 1888, 1889, 1901, 1902, 2041, 2042, 2043, 2044, 2045,
2046, 2133, 2134, 2137, 2138, 2140, 2141], dtype=int32), array([1270, 1271, 1260, 1261, 1268, 1269, 1460, 1461, 1462, 1463, 1456,
1457, 1482, 1483, 1478, 1479, 1480, 1481], dtype=int32)], FunctionSpace(Mesh(VectorElement(Basix element (P, triangle, 1, equispaced, unset, False), 2), 0), VectorElement(FiniteElement('Lagrange', triangle, 2), dim=2))

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
Cell In[17], line 20
16 #center : u = 0
17 #boundary_dofs_u_center = locate_dofs_topological((V.sub(0), U), bdim, boundary_facets_center)
18 #BC_u_center = dirichletbc(u_bc, boundary_dofs_u_center, V.sub(0))
19 boundary_dofs_u_center = locate_dofs_topological((V.sub(0), U), dim, cells_center)
---> 20 BC_u_center = dirichletbc(np.full_like(boundary_dofs_u_center, 0, dtype=ScalarType)[0], boundary_dofs_u_center, V.sub(0))
22 #Dirichlet BC for p : surface
23 P, _ = V.sub(1).collapse()

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/bcs.py:183, in dirichletbc(value, dofs, V)
180     raise AttributeError("Boundary condition value must have a dtype attribute.")
182 formcls = type("DirichletBC", (DirichletBCMetaClass, bctype), {})
--> 183 return formcls(value, dofs, V)

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/bcs.py:127, in DirichletBCMetaClass.__init__(self, value, dofs, V)
125         super().__init__(_value, dofs, V)  # type: ignore
126     except TypeError:
--> 127         super().__init__(_value, dofs, V._cpp_object)  # type: ignore
128 else:
129     super().__init__(_value, dofs)

TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
1. dolfinx.cpp.fem.DirichletBC_float64(g: numpy.ndarray[numpy.float64], dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace)
2. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Constant<double>, dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace)
3. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double>, dofs: numpy.ndarray[numpy.int32])
4. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double>, dofs: List[numpy.ndarray[numpy.int32][2]], V: dolfinx::fem::FunctionSpace)

Invoked with: array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0.]), [array([1796, 1797, 1888, 1889, 1901, 1902, 2041, 2042, 2043, 2044, 2045,
2046, 2133, 2134, 2137, 2138, 2140, 2141], dtype=int32), array([1270, 1271, 1260, 1261, 1268, 1269, 1460, 1461, 1462, 1463, 1456,
1457, 1482, 1483, 1478, 1479, 1480, 1481], dtype=int32)], <dolfinx.cpp.fem.FunctionSpace object at 0x7f68ce2abdb0>


Could you help me to implement this setting?
ALSO, If you can suggest better function for defining center compared to this code :

def center(x):
return x[0]**2 + x[1]**2 < (0.1*R)**2


Please let me know. It seems that np.isclose might not work depending on mesh cell size.

Thanks.

Here you are mixing approaches, how about using

boundary_dofs_u_center = locate_dofs_topological((V.sub(0), U), dim, cells_center)
BC_u_center = dirichletbc(u_bc, boundary_dofs_u_center, V.sub(0))


Dear Dokken,
Thanks for your help, it works!
Also, could you comment on this one?

def center(x):
return x[0]**2 + x[1]**2 < (0.1*R)**2


My mesh is like this :

Always thanks.

You should really mark the center of your mesh during mesh creation, rather than doing it using geometrical information after mesh creation.

I think I understand what you mean.
Here is the code for mesh generation (by following the “Electromagnetics example” demo code)

rank = MPI.COMM_WORLD.rank

gmsh.initialize()
R = 1
gdim = 2  # Geometric dimension of the mesh
model_rank = 0
mesh_comm = MPI.COMM_WORLD
if mesh_comm.rank == model_rank:

# Define geometry for background
background = gmsh.model.occ.addDisk(0, 0, 0, R, R)
gmsh.model.occ.synchronize()

gmsh.model.occ.synchronize()

# Generate mesh

gmsh.model.mesh.field.setNumber(2, "IField", 1)
#gmsh.model.mesh.field.setNumber(2, "LcMin", 0.1)
gmsh.model.mesh.field.setNumber(2, "LcMax", 0.1)
#gmsh.model.mesh.field.setNumber(2, "DistMin", 1)
#gmsh.model.mesh.field.setNumber(2, "DistMax", 1)

gmsh.model.mesh.field.setNumbers(5, "FieldsList", [2])
gmsh.model.mesh.field.setAsBackgroundMesh(5)

gmsh.option.setNumber("Mesh.Algorithm", 7)
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.optimize("Netgen")

domain, ct, _ = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim)
gmsh.finalize()


Here I set the tag for the circle as 1.
Can I use this tag 1 or should define ‘center’ again?

Do you actually want to have a vertex in the center of the mesh? If so I would try to explicitly enforce that.

How can I find some materials to implement such settings?

1 Like

For a different perspective, note that you’ve got spherical symmetry. You could just use a 1D domain along the coordinate r, you don’t need to run a 2D or 3D calculation. In that case the boundary condition at r=0 can be handled naturally.

1 Like

Dear Dokken,
Thanks to your reply, I applied the marker at the center as follows :

# https://github.com/FEniCS/dolfinx/blob/main/python/demo/demo_gmsh.py
rank = MPI.COMM_WORLD.rank

gmsh.initialize()
R = 1
gdim = 2  # Geometric dimension of the mesh
model_rank = 0
mesh_comm = MPI.COMM_WORLD
if mesh_comm.rank == model_rank:

# Define geometry for background
background = gmsh.model.occ.addDisk(0, 0, 0, R, R)
gmsh.model.occ.synchronize()

# Define geometry for center
gmsh.model.occ.synchronize()

# Generate mesh

gmsh.model.mesh.field.setNumber(2, "IField", 1)
#gmsh.model.mesh.field.setNumber(2, "LcMin", 0.1)
gmsh.model.mesh.field.setNumber(2, "LcMax", 0.1)
#gmsh.model.mesh.field.setNumber(2, "DistMin", 1)
#gmsh.model.mesh.field.setNumber(2, "DistMax", 1)

gmsh.model.mesh.field.setNumbers(5, "FieldsList", [2])
gmsh.model.mesh.field.setAsBackgroundMesh(5)

gmsh.option.setNumber("Mesh.Algorithm", 7)
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.optimize("Netgen")

domain, ct, _ = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim)
gmsh.finalize()

with XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_meshtags(ct)

x = SpatialCoordinate(domain)

dim = domain.topology.dim
bdim = dim - 1


However, If I run this code, the Python kernel is broken, even without any error message. It seems that the code for the marking center would have some problems :

center_domain = gmsh.model.addPhysicalGroup(0, [center], 1)


How should I fix my code?
Always thanks.

I suggested didnt suggest adding a physical tag, but using the embed command, as follows:

import gmsh
from mpi4py import MPI
from dolfinx.io.gmshio import model_to_mesh
from dolfinx.io import XDMFFile
# https://github.com/FEniCS/dolfinx/blob/main/python/demo/demo_gmsh.py
rank = MPI.COMM_WORLD.rank

gmsh.initialize()
R = 1
gdim = 2  # Geometric dimension of the mesh
model_rank = 0
mesh_comm = MPI.COMM_WORLD
if mesh_comm.rank == model_rank:

# Define geometry for background
background = gmsh.model.occ.addDisk(0, 0, 0, R, R)

gmsh.model.occ.synchronize()

# Define geometry for center

gmsh.model.occ.synchronize()
gmsh.model.mesh.embed(0, [center], 2, background)

# Generate mesh

gmsh.model.mesh.field.setNumber(2, "IField", 1)
gmsh.model.mesh.field.setNumber(2, "LcMin", 0.01)
gmsh.model.mesh.field.setNumber(2, "LcMax", 0.1)
gmsh.model.mesh.field.setNumber(2, "DistMin", 0.2)
gmsh.model.mesh.field.setNumber(2, "DistMax", 0.7)

gmsh.model.mesh.field.setNumbers(5, "FieldsList", [2])
gmsh.model.mesh.field.setAsBackgroundMesh(5)

gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.optimize("Netgen")

domain, ct, _ = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim)
gmsh.finalize()

with XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_meshtags(ct)
import numpy as np
from dolfinx.mesh import locate_entities
print(locate_entities(domain, 0, lambda x: np.isclose(x[0],0) & np.isclose(x[1], 0)))