I am trying to recreate the demo problem from @dokken 's tutorial for subdomains. That tutorial uses an external mesh and solves the poisson equation with variable coefficients. However, I was trying to develop this demo in FEniCs dolfin
instead of dolfinx
, so that I can go back and forth between the old and new version.
I checked the existing posts, and I don’t think this is a duplicate. Sorry if I am wrong about that.
I have read through the very comprehensive post.Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio - #3 by dokken , a number of times. and so I have a sense of how to create the mesh and import that mesh into Fenics. But I am lost when it comes to using this imported mesh to set the boundary conditions and function spaces in Fenics. In dolfinx
there is a function to find()
a value on a mesh, but that function does not exist in dolfin
. I was hoping someone could help me complete this example. I really just need help with connecting the boundary conditions and the function value for the kappa
function. My MWE is a bit long, but hopefully the code is very generic and so most people are familiar with it. The most frequent comment I saw on the existing posts was a request to include the mesh file, or code that generated the error, so I wanted to be complete.
I created the mesh using the same code as in the dolfinx
tutorial.
import gmsh
import numpy as np
import pyvista
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
pyvista.start_xvfb()
gmsh.initialize()
proc = MPI.COMM_WORLD.rank
top_marker = 2
bottom_marker = 1
left_marker = 1
if proc == 0:
# We create one rectangle for each subdomain
gmsh.model.occ.addRectangle(0, 0, 0, 1, 0.5, tag=1)
gmsh.model.occ.addRectangle(0, 0.5, 0, 1, 0.5, tag=2)
# We fuse the two rectangles and keep the interface between them
gmsh.model.occ.fragment([(2,1)],[(2,2)])
gmsh.model.occ.synchronize()
# Mark the top (2) and bottom (1) rectangle
top, bottom = None, None
for surface in gmsh.model.getEntities(dim=2):
com = gmsh.model.occ.getCenterOfMass(surface[0], surface[1])
if np.allclose(com, [0.5,0.25, 0]):
bottom = surface[1]
else:
top = surface[1]
gmsh.model.addPhysicalGroup(2, [bottom], bottom_marker)
gmsh.model.addPhysicalGroup(2, [top], top_marker)
# Tag the left boundary
left = []
for line in gmsh.model.getEntities(dim=1):
com = gmsh.model.occ.getCenterOfMass(line[0], line[1])
if np.isclose(com[0], 0):
left.append(line[1])
gmsh.model.addPhysicalGroup(1, left, left_marker)
gmsh.model.mesh.generate(2)
gmsh.write("mesh.msh")
gmsh.finalize()
Then I imported this mesh into dolfin
using the code from the previously posted thread and the tutorial. This part works.
import gmsh
import pygmsh
import meshio
import numpy as np
from dolfin import *
# This function works with meshio 5.4.3 and later.
def create_mesh(mesh: meshio.Mesh, cell_type: str, data_name: str = "name_to_read",
prune_z: bool = False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:, :2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
data_name: [cell_data]})
return out_mesh
msh = meshio.read('mesh.msh')
triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)
line_mesh = create_mesh(msh, "line", data_name="Boundary", prune_z=True)
meshio.write("mf.xdmf", line_mesh)
# import the meshes
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mf.xdmf") as infile:
infile.read(mvc, "Boundary")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
#File("dolfin_mesh.pvd").write(mesh)
#File("dolfin_cellfunc.pvd").write(mf)
#File("dolfin_boundary.pvd").write(boundary_markers)
The problem is I don’t know how to connect the pieces after this point. The code below uses the find()
function which does not exist in dolfin
according to the error message posted below. Here was my attempt.
top_marker = 2
bottom_marker = 1
left_marker = 1
Q = FunctionSpace(mesh, "DG", 0)
kappa = Function(Q)
bottom_cells = mesh.find(bottom_marker) # <----SOURCE OF THE PROBLEM
kappa.x.array[bottom_cells] = np.full_like(bottom_cells, 1, dtype=ScalarType)
top_cells = mesh.find(top_marker)
kappa.x.array[top_cells] = np.full_like(top_cells, 0.1, dtype=ScalarType)
V = FunctionSpace(mesh, "CG", 1)
u_bc = Function(V)
left_facets = ft.find(left_marker)
left_dofs = locate_dofs_topological(V, mesh.topology.dim-1, left_facets)
bcs = [dirichletbc(ScalarType(1), left_dofs, V)]
The error message I get is.
AttributeError: 'dolfin.cpp.mesh.Mesh' object has no attribute 'find'
I just need some help to find the right way to connect the imported mesh to the dirichlet boundary conditions, and the kappa
function. Any help would be appreciated.