Defining subdomains from gmsh on CG elements

dolfinx
1

Hi everyone,
This great FEniCSx tutorial clearly shows how to define a dolfinx function from gmsh data.

However, playing around with the code I found that this solution is dependent on the choice of the element. If I change:

Q = dolfinx.fem.FunctionSpace(mesh, ("DG", 0))

to:

Q = dolfinx.fem.FunctionSpace(mesh, ("CG", 1))

I get an error:

IndexError: index 108 is out of bounds for axis 0 with size 103

Is there any way to define a CG function from mesh data?

Here you can find a reduced version of the tutorial code, modified to reproduce the issue:

import gmsh
import meshio
import dolfinx
from mpi4py import MPI
import numpy as np
from petsc4py import PETSc

# get comm and rank
comm = MPI.COMM_WORLD
proc = comm.Get_rank()

gmsh.initialize()
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)
    gmsh.model.mesh.generate(2)
    gmsh.write("mesh.msh")
gmsh.finalize()


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)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh


if proc == 0:
    # Read in mesh
    msh = meshio.read("mesh.msh")

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

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
    ct = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)

# define function
Q = dolfinx.fem.FunctionSpace(mesh, ("CG", 1))  # I changed this
kappa = dolfinx.fem.Function(Q)
bottom_cells = ct.find(bottom_marker)
kappa.x.array[bottom_cells] = np.full_like(bottom_cells, 1, dtype=PETSc.ScalarType)
top_cells = ct.find(top_marker)
kappa.x.array[top_cells] = np.full_like(top_cells, 0.1, dtype=PETSc.ScalarType)

Thank you

2

A CG-1 function based on cell data is not well defined (a vertex can be part of multiple cells with different values in each cell).

If you still want to do this, you would have to:

  1. Read in the data as DG-0
  2. Interpolate (not well defined) or project (well defined) the DG-0 data into a CG-1 space.