How to set 'boundary condition' inside the boundary?

[psj1866]

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.

[dokken]

@psj1866

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

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))

[psj1866]

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.

[dokken]

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

[psj1866]

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)

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?

[dokken]

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

[psj1866]

Thanks for your reply.
How can I find some materials to implement such settings?

[dokken]

http://onelab.info/pipermail/gmsh/2019/013265.html

[dparsons]

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.

[psj1866]

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
    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 :

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

How should I fix my code?
Always thanks.

[dokken]

The link

@dokken

http://onelab.info/pipermail/gmsh/2019/013265.html

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
    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.

[psj1866]

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!