Recreating the Poisson subdomain problem in Dolfin - with external mesh

General
1

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.

2

See for instance: How to define different materials / importet 3D geometry (gmsh) - #2 by dokken

3

@dokken I was just reading this in the dolfin docs, https://fenicsproject.org/olddocs/dolfin/latest/python/demos/subdomains-poisson/documentation.html?highlight=subdomains . Is this still a valid way to do things? I can see how this maps to the example that you gave, with the MeshFunction and Subdomain class. I was just about to try and run this code and see if it works :). Thanks so much for your help and for providing a link to this example. It helps to have multiple examples.

4

The approach is still valid in legacy dolfin.

1 Like