Locate_dofs_topological behave abnormally

dolfinx
1

Dear all,

I used dolfinx to solve the Poisson equation for the problem with two electric poles, which are set to 100 V and -100 V respectively. The rest of the boundaries are Neumann ones. However I can not get the correct solution. My code is displayed as following

import pygmsh
import numpy as np
import meshio
from mpi4py import MPI
from ufl import  dx, grad, inner, TestFunction, TrialFunction
from dolfinx.fem import Function, FunctionSpace, dirichletbc, Constant
from dolfinx.fem import locate_dofs_topological
from dolfinx.io import XDMFFile
from dolfinx.fem.petsc import LinearProblem
from petsc4py.PETSc import ScalarType

def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh
with pygmsh.geo.Geometry() as geom:
    msh = meshio.read("rfq_pole.msh")

    # Create and save one file for the mesh, and one file for the facets
    line_mesh = create_mesh(msh, "line", prune_z=True)
    triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
    meshio.write("mesh.xdmf", triangle_mesh)
    meshio.write("mt.xdmf", line_mesh)


with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf_infile:
    mesh = xdmf_infile.read_mesh(name="Grid")
    mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)
    mvc_subdomain = xdmf_infile.read_meshtags(mesh, "Grid")

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
with XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf_infile:
    mvc_boundaries = xdmf_infile.read_meshtags(mesh, "Grid")

V = FunctionSpace(mesh, ('CG', 2))
u = TrialFunction(V)
v = TestFunction(V)

pos_pole = Function(V)
with pos_pole.vector.localForm() as bc_local:
    bc_local.set(100)

neg_pole = Function(V)
with neg_pole.vector.localForm() as bc_local:
    bc_local.set(-100)

tdim = mesh.topology.dim
fdim = tdim - 1
print(mvc_boundaries.values)
pos_pole_facets = np.where(mvc_boundaries.values == 17)[0]
pos_pole_dofs = locate_dofs_topological(V, fdim, pos_pole_facets)
pos_bc = dirichletbc(pos_pole, pos_pole_dofs)
neg_pole_facets = np.where(mvc_boundaries.values == 18)[0]
neg_pole_dofs = locate_dofs_topological(V, fdim, neg_pole_facets)
neg_bc = dirichletbc(neg_pole, neg_pole_dofs)

a = inner(grad(u), grad(v))*dx
L = inner(Constant(mesh, ScalarType(0)), v)*dx 
problem = LinearProblem(a, L, bcs=[pos_bc, neg_bc])
uh = problem.solve()

with XDMFFile(MPI.COMM_WORLD, "pos.xdmf", 'w',
              encoding=XDMFFile.Encoding.HDF5) as file:
    file.write_mesh(mesh)
    file.write_function(uh)

Here is my geo file:

// Gmsh project created on Fri May 06 10:02:11 2022
SetFactory("OpenCASCADE");
//+
Point(3) = {0, 0.4, 0, 0.05};
//+
Point(4) = {0.375, 0.775, 0, 0.05};
//+
Circle(1) = {3, 1, 4};
//+
Line(2) = {4, 5};
//+
Line(3) = {5, 6};
//+
Line(4) = {6, 3};
//+
Point(7) = {0.6, 0, 0, 0.05};
//+
Point(9) = {0.975, 0.375, 0, 0.05};
//+
Circle(5) = {7, 8, 9};
//+
Line(6) = {9, 10};
//+
Line(7) = {10, 11};
//+
Line(8) = {11, 7};
//+
Point(12) = {0, 0, 0, 0.05};
//+
Line(9) = {3, 12};
//+
Line(10) = {12, 7};
//+
Point(13) = {0, 0.775, 0, 1.0};
//+
Circle(11) = {3, 13, 4};
//+
Point(14) = {0.975, 0, 0, 1.0};
//+
Circle(12) = {7, 14, 9};
//+
Line(13) = {5, 10};
//+
Point(15) = {0.375, 0.375, 0, 1.0};
//+
Point(16) = {0.375, 2.375, 0, 0.15};
//+
Point(17) = {2.375, 0.375, 0, 0.15};
//+
Line(13) = {16, 4};
//+
Line(14) = {9, 17};
//+
Circle(15) = {16, 15, 17};
//+
Physical Curve("neumann", 16) = {9, 10, 15};
//+
Physical Curve("positive_pole", 17) = {11, 13};
//+
Physical Curve("negative_pole", 18) = {12, 14};
//+
Curve Loop(1) = {13, -11, 9, 10, 12, 14, -15};
//+
Plane Surface(1) = {1};
//+
Physical Surface("vac", 19) = {1};

I have already tested the code with unit square mesh, and it gave the correct result. Therefore, the error must occur at the stage of importing mesh from gmsh. I really appreciate your help.

2

The issue is that your code finding the facets to set the DirichletBC on does not make sense:

Here you are assuming that mvc_boundaries contains all facet indices of the mesh, which is not the case.
You get the expected result if you use:

pos_pole_facets = mvc_boundaries.indices[mvc_boundaries.values == 17]
print(pos_pole_facets)
pos_pole_dofs = locate_dofs_topological(V, fdim, pos_pole_facets)
pos_bc = dirichletbc(pos_pole, pos_pole_dofs)
neg_pole_facets = mvc_boundaries.indices[mvc_boundaries.values == 18]
3

Dear dokken,
It works and thank you very much.