How to impose Dirichlet boundary condition at a point

General
1

Hello I want to assign Dirichlet condition in the center of my mesh. But the solution doesn’t seem to agree with the Dirichlet condition assigned.

The system has 2 unknowns both in second order. 2 bc were weakly imposed. The other two are (u_x, u_y) = (0, 0) in the center of my mesh. The correct solution should be (u_x, u_y) = (x, y). But the actual solution it returns me seems to be this correct solution shifted by some points. And the Dirichlet condition seems to be assigned somewhere else instead of the center. From my understanding, the strong imposition should at least ensure the value of the center. Can you help me to see what is wrong? Here is an old but related post.

I used

# bc at center
def center(x):
    return np.logical_and(np.isclose(x[0], 0), np.isclose(x[1], 0))
boundary_facets = dolfinx.mesh.locate_entities(mesh, 0, center)
u_center = dolfinx.fem.Function(V)
u_center.x.array[:] = 0

bc_center = dolfinx.fem.dirichletbc(u_center, boundary_facets)

to locate center point and apply center bc. I think in the solution at least this point should [0, 0] but it’s not.

The complete code

import basix.ufl
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import dolfinx.nls.petsc
import dolfinx.log
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import ufl
import multiphenicsx.fem
import multiphenicsx.fem.petsc
import viskex
from petsc4py.PETSc import ScalarType
import matplotlib.pyplot as plt

# hexogonal mesh
a = 6.
lcar = 1.2

gmsh.initialize()
gmsh.model.add("mesh")

p0 = gmsh.model.geo.addPoint(a, 0, 0, lcar)
p1 = gmsh.model.geo.addPoint(a/2, np.sqrt(3.)*a/2, 0, lcar)
p2 = gmsh.model.geo.addPoint(-a/2, np.sqrt(3)*a/2, 0, lcar)
p3 = gmsh.model.geo.addPoint(-a, 0, 0, lcar)
p4 = gmsh.model.geo.addPoint(-a/2, -np.sqrt(3)*a/2, 0, lcar)
p5 = gmsh.model.geo.addPoint(a/2, -np.sqrt(3)*a/2, 0, lcar)

l0 = gmsh.model.geo.addLine(p0, p1)
l1 = gmsh.model.geo.addLine(p1, p2)
l2 = gmsh.model.geo.addLine(p2, p3)
l3 = gmsh.model.geo.addLine(p3, p4)
l4 = gmsh.model.geo.addLine(p4, p5)
l5 = gmsh.model.geo.addLine(p5, p0)

line_loop = gmsh.model.geo.addCurveLoop([l0, l1, l2, l3, l4, l5])
domain = gmsh.model.geo.addPlaneSurface([line_loop])

p_center = gmsh.model.geo.addPoint(0, 0, 0, lcar)

gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [l0, l1, l2, l3, l4, l5], 1)
gmsh.model.addPhysicalGroup(2, [domain], 0)
gmsh.model.addPhysicalGroup(0, [p_center], 3)

gmsh.model.mesh.embed(0, [p_center], 2, domain)

gmsh.model.mesh.generate()

mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(
    gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2
)
gmsh.finalize()

boundaries_full = boundaries.indices[boundaries.values == 1]

bb_tree = dolfinx.geometry.bb_tree(mesh, mesh.topology.dim)

V = dolfinx.fem.VectorFunctionSpace(mesh, ("Lagrange", 2))

vxvy = ufl.TestFunction(V)
(vx, vy) = vxvy[0], vxvy[1]
duxuy = ufl.TrialFunction(V)
uxuy = dolfinx.fem.Function(V)
(ux, uy) = uxuy[0], uxuy[1]

Y = 1.
nu_p = 0.3
tau = 0

beta = 2.
R_0 = 10

gbar_ij = ufl.as_matrix([[(ux.dx(0))**2 + (uy.dx(0))**2, (ux.dx(0)) * (ux.dx(1)) + (uy.dx(0)) * (uy.dx(1))],
                         [(ux.dx(0)) * (ux.dx(1)) + (uy.dx(0)) * (uy.dx(1)), (ux.dx(1))**2 + (uy.dx(1))**2]])

g_ij = ufl.Identity(2)

gbarij = ufl.as_matrix([[((ux.dx(1))**2 + (uy.dx(1))**2) / (((uy.dx(1)) * (ux.dx(0))-(ux.dx(1)) * (uy.dx(0)))**2), - ((ux.dx(1)) * (ux.dx(0)) + (uy.dx(1)) * uy.dx(0)) / (((uy.dx(1)) * (ux.dx(0)) - (ux.dx(1)) * (uy.dx(0)))**2)],
                       [- ((ux.dx(1)) * (ux.dx(0)) + (uy.dx(1)) * (uy.dx(0))) / (((uy.dx(1)) * (ux.dx(0)) - (ux.dx(1)) * (uy.dx(1)))**2), ((ux.dx(0))**2 + (uy.dx(0))**2) / (((uy.dx(1)) * (ux.dx(0)) - (ux.dx(1)) * (uy.dx(0)))**2)]])

gij = ufl.Identity(2)

g_ij = ufl.as_tensor([[dolfinx.fem.Constant(mesh, ScalarType(1)), dolfinx.fem.Constant(mesh, ScalarType(0))],
                      [dolfinx.fem.Constant(mesh, ScalarType(0)), dolfinx.fem.Constant(mesh, ScalarType(1))]])
gij = ufl.as_tensor([[dolfinx.fem.Constant(mesh, ScalarType(1)), dolfinx.fem.Constant(mesh, ScalarType(0))],
                      [dolfinx.fem.Constant(mesh, ScalarType(0)), dolfinx.fem.Constant(mesh, ScalarType(1))]])

i, j, k, l = ufl.indices(4)

Gammai_jk = ufl.as_tensor((1./2.) * gij[i, j] * (g_ij[j, k].dx(l) + g_ij[j, l].dx(k) - g_ij[k, l].dx(j)), (i,k,l))
# Gammai_jk = ufl.as_tensor(0 , (i,j,k))
Gammabari_jk = ufl.as_tensor((1./2.) * gbarij[i, j] * (gbar_ij[j, k].dx(l) + gbar_ij[j, l].dx(k) - gbar_ij[k, l].dx(j)), (i, k, l))

Aijkl = ufl.as_tensor((Y / (1 - nu_p**2)) * (nu_p * gij[i, j] * gij[k, l] + (1 - nu_p) * gij[i, k] * gij[j, l]), (i,j,k,l))

u_ij = ufl.as_tensor((1./2.) * (g_ij[i, j] - gbar_ij[i, j]), (i, j))

sigmaij = ufl.as_tensor(Aijkl[i,j,k,l] * u_ij[k,l], (i,j))

F = - sigmaij[i, j] * (vxvy[j].dx(i)) * ufl.dx + (Gammabari_jk[i,j,k] * sigmaij[j,k]) * vxvy[i] * ufl.dx

# bc at center
def center(x):
    return np.logical_and(np.isclose(x[0], 0), np.isclose(x[1], 0))
boundary_facets = dolfinx.mesh.locate_entities(mesh, 0, center)
u_center = dolfinx.fem.Function(V)
u_center.x.array[:] = 0

bc_center = dolfinx.fem.dirichletbc(u_center, boundary_facets)

problem = dolfinx.fem.petsc.NonlinearProblem(F, uxuy, bcs=[bc_center])

solver = dolfinx.nls.petsc.NewtonSolver(mpi4py.MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"
solver.rtol = 1e-10
solver.max_it = 100
solver.report = True

ksp = solver.krylov_solver
opts = petsc4py.PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

uxuy.interpolate(lambda x: [(x[0]) * 2, (x[1])])

dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
n, converged = solver.solve(uxuy)
assert (converged)
print(f"Number of interations: {n:d}")

viskex.dolfinx.plot_vector_field(uxuy, "uxuy") # plot solution

points_center = np.zeros((3, 1))
cell_candidates_center = dolfinx.geometry.compute_collisions_points(bb_tree, points_center.T)
colliding_cells_center = dolfinx.geometry.compute_colliding_cells(mesh, cell_candidates_center, points_center.T)
cells_center = []
points_on_proc_center = []
for i, point in enumerate(points_center.T):
    if len(colliding_cells_center.links(i)) > 0:
        points_on_proc_center.append(point)
        cells_center.append(colliding_cells_center.links(i)[0])

uxuy.eval(points_on_proc_center, cells_center)
2

Boundary conditions get applied based on the DOF number, not the geometric entity ID.

You may want to try something like

dof_center = dolfinx.fem.locate_dofs_geometrical(V, center)

bc_center = dolfinx.fem.dirichletbc(u_center, dof_center)

although I haven’t tested it with your example.

1 Like
3

as @francesco-ballarin says you need to send in the degrees of freedom to dirichletbc. Irregardless of this, if you are looking for a single point (and looking in something than 1D), a facet does have more than two vertices, and both vertices of a facet cannot satisfy this criterion

4

Francesco and dokken thank you very much! This is exactly the issue.

But as I type print(dolfinx.mesh.locate_entities(mesh, 0, center)) and print(dolfinx.fem.locate_dofs_geometrical(V, center)), they both return me a tag: the first one returns array([54], dtype=int32) and the second array([183], dtype=int32). The input for the first one is mesh, so the output should be about this geometrical mesh. That should be the geometrical point. The second input is the function space. So the output should be related to the function thing. That’s the dofs. But how the function dolfinx.fem.dirichletbc() differentiate one from the other given they are in the same type. Is it in the tag [#], the smaller number represents the geometrical point, the larger number represents the dofs?

Here’s the revised code for this session if you are interested:

# bc at center
def center(x):
    return np.logical_and(np.isclose(x[0], 0), np.isclose(x[1], 0))
dof_center = dolfinx.fem.locate_dofs_geometrical(V, center)
u_center = dolfinx.fem.Function(V)
u_center.x.array[:] = 0

bc_center = dolfinx.fem.dirichletbc(u_center, dof_center)
5

But how the function dolfinx.fem.dirichletbc() differentiate one from the other given they are in the same type

It doesn’t. help(dolfinx.fem.dirichletbc) clearly states:

Help on function dirichletbc in module dolfinx.fem.bcs:

dirichletbc(value: 'typing.Union[Function, Constant, np.ndarray]', dofs: 'numpy.typing.NDArray[np.int32]', V: 'typing.Optional[dolfinx.fem.FunctionSpace]' = None) -> 'DirichletBC'
    Create a representation of Dirichlet boundary condition which
    is imposed on a linear system.
    
    Args:
        value: Lifted boundary values function. It must have a ``dtype``
        property.
        dofs: Local indices of degrees of freedom in function space to
            which boundary condition applies. Expects array of size
            (number of dofs, 2) if function space of the problem, ``V``,
            is passed. Otherwise assumes function space of the problem
            is the same of function space of boundary values function.
        V: Function space of a problem to which boundary conditions are applied.
    
    Returns:
        A representation of the boundary condition for modifying linear
        systems.

If you pass in an entity index rather than a DOF index, it will apply the BC to the incorrect physical location (the one which has a DOF associated to that index), as you were observing in the original post.

1 Like
6

I see. Thank you very much!

8

One more note is that, in this implementation, it is essential to embed a point exactly at where you want to put the boundary condition. Otherwise the degree of freedom is a null array. It is unclear why it’s the case though.

9

If you don’t have a mesh point at the origin (= the point you would like to find), you can loosen up the tolerances in np.isclose(...), and use a tolerance comparable to the mesh size. This may end up finding more than one point, and then you should decide which one to use (the first one? the one closest to the origin? a random one? up to you)

2 Likes