# How to set 'boundary condition' inside the boundary? **URL:** \ **Category:** General\ **Created:** [May 24, 2023, 5:38am UTC](https://fenicsproject.discourse.group/t/how-to-set-boundary-condition-inside-the-boundary/11253 "2023-05-24T05:38:40Z")\ **Posts on this page:** 12\ **Page:** 1
**Author:** ![psj1866](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/psj1866/32/4987_2.png) [@psj1866](https://fenicsproject.discourse.group/u/psj1866)\ **Post date:** [May 24, 2023, 5:38am UTC](https://fenicsproject.discourse.group/t/how-to-set-boundary-condition-inside-the-boundary/11253/1 "2023-05-24T05:38:40Z")
Hello, I am solving the coupling equation for displacement(u) and (p). The current topic is “Cryer’s problem” with such boundary conditions : ![Cryers-problem-of-a-clay-sphere](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/0/0b61ae1b6ba19a41a2d976263642470fd2f9d74c.png) - 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 : ```auto # 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) ``` ```auto 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 : ```auto 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, dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace) 3. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function, dofs: numpy.ndarray[numpy.int32]) 4. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function, 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, dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace) 3. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function, dofs: numpy.ndarray[numpy.int32]) 4. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function, 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)], ``` Could you help me to implement this setting? ALSO, If you can suggest better function for defining center compared to this code : ```auto 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. ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [May 24, 2023, 6:44am UTC](https://fenicsproject.discourse.group/t/how-to-set-boundary-condition-inside-the-boundary/11253/2 "2023-05-24T06:44:43Z")
> [@psj1866](#): > > ```auto > 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)) > > ``` Here you are mixing approaches, how about using ```python 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)) ``` ---
**Author:** ![psj1866](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/psj1866/32/4987_2.png) [@psj1866](https://fenicsproject.discourse.group/u/psj1866)\ **Post date:** [May 24, 2023, 7:06am UTC](https://fenicsproject.discourse.group/t/how-to-set-boundary-condition-inside-the-boundary/11253/3 "2023-05-24T07:06:17Z")
Dear Dokken, Thanks for your help, it works! Also, could you comment on this one? ```auto def center(x): return x[0]**2 + x[1]**2 < (0.1*R)**2 ``` My mesh is like this : ![Screenshot from 2023-05-24 16-05-21](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/4/4dcdc98d0eac6453b921cf87276735717c857d41.png) Always thanks. ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [May 24, 2023, 7:13am UTC](https://fenicsproject.discourse.group/t/how-to-set-boundary-condition-inside-the-boundary/11253/4 "2023-05-24T07:13:21Z")
You should really mark the center of your mesh during mesh creation, rather than doing it using geometrical information after mesh creation. ---
**Author:** ![psj1866](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/psj1866/32/4987_2.png) [@psj1866](https://fenicsproject.discourse.group/u/psj1866)\ **Post date:** [May 24, 2023, 7:21am UTC](https://fenicsproject.discourse.group/t/how-to-set-boundary-condition-inside-the-boundary/11253/5 "2023-05-24T07:21:25Z")
Thanks for your reply. I think I understand what you mean. Here is the code for mesh generation (by following the “Electromagnetics example” demo code) ```auto 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() whole_domain = gmsh.model.addPhysicalGroup(gdim, [background], 1) gmsh.model.occ.synchronize() # Generate mesh gmsh.model.mesh.field.add("Distance", 1) gmsh.model.mesh.field.add("Threshold", 2) 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.add("Min", 5) 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? ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [May 24, 2023, 7:57am UTC](https://fenicsproject.discourse.group/t/how-to-set-boundary-condition-inside-the-boundary/11253/6 "2023-05-24T07:57:02Z")
Do you actually want to have a vertex in the center of the mesh? If so I would try to explicitly enforce that. ---
**Author:** ![psj1866](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/psj1866/32/4987_2.png) [@psj1866](https://fenicsproject.discourse.group/u/psj1866)\ **Post date:** [May 24, 2023, 9:03am UTC](https://fenicsproject.discourse.group/t/how-to-set-boundary-condition-inside-the-boundary/11253/7 "2023-05-24T09:03:55Z")
Thanks for your reply. How can I find some materials to implement such settings? ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [May 24, 2023, 9:32am UTC](https://fenicsproject.discourse.group/t/how-to-set-boundary-condition-inside-the-boundary/11253/8 "2023-05-24T09:32:56Z")
[http://onelab.info/pipermail/gmsh/2019/013265.html](http://onelab.info/pipermail/gmsh/2019/013265.html) ---
**Author:** ![dparsons](https://avatars.discourse-cdn.com/v4/letter/d/5e9695/32.png) [@dparsons](https://fenicsproject.discourse.group/u/dparsons)\ **Post date:** [May 25, 2023, 8:21am UTC](https://fenicsproject.discourse.group/t/how-to-set-boundary-condition-inside-the-boundary/11253/9 "2023-05-25T08:21:41Z")
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. ---
**Author:** ![psj1866](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/psj1866/32/4987_2.png) [@psj1866](https://fenicsproject.discourse.group/u/psj1866)\ **Post date:** [May 27, 2023, 12:11pm UTC](https://fenicsproject.discourse.group/t/how-to-set-boundary-condition-inside-the-boundary/11253/10 "2023-05-27T12:11:51Z")
Dear Dokken, Thanks to your reply, I applied the marker at the center as follows : ```auto # 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 center = gmsh.model.occ.addPoint(0, 0, 0) gmsh.model.occ.synchronize() background_domain = gmsh.model.addPhysicalGroup(gdim, [background], 1) center_domain = gmsh.model.addPhysicalGroup(0, [center], 1) # Generate mesh gmsh.model.mesh.field.add("Distance", 1) gmsh.model.mesh.field.add("Threshold", 2) 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.add("Min", 5) 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 : ```auto center_domain = gmsh.model.addPhysicalGroup(0, [center], 1) ``` How should I fix my code? Always thanks. ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [May 27, 2023, 1:41pm UTC](https://fenicsproject.discourse.group/t/how-to-set-boundary-condition-inside-the-boundary/11253/11 "2023-05-27T13:41:16Z")
The link > [@dokken](#): > > [http://onelab.info/pipermail/gmsh/2019/013265.html](http://onelab.info/pipermail/gmsh/2019/013265.html) I suggested didnt suggest adding a physical tag, but using the embed command, as follows: ```python 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 center = gmsh.model.occ.addPoint(0, 0, 0) gmsh.model.occ.synchronize() gmsh.model.mesh.embed(0, [center], 2, background) background_domain = gmsh.model.addPhysicalGroup(gdim, [background], 1) # Generate mesh gmsh.model.mesh.field.add("Distance", 1) gmsh.model.mesh.field.add("Threshold", 2) 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.add("Min", 5) 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))) ``` Please read through the links I post before asking a question. ---
**Author:** ![psj1866](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/psj1866/32/4987_2.png) [@psj1866](https://fenicsproject.discourse.group/u/psj1866)\ **Post date:** [May 27, 2023, 2:15pm UTC](https://fenicsproject.discourse.group/t/how-to-set-boundary-condition-inside-the-boundary/11253/12 "2023-05-27T14:15:21Z")
Dear Dokken, Sorry for my poor understanding of your suggested document. I didn’t come up with using locate\_entities, but only concentrate on imposing a marker in the center 😅 Thanks for spending your valuable time on my question!