# Defining subdomains from gmsh on CG elements

**URL:** <https://fenicsproject.discourse.group/t/defining-subdomains-from-gmsh-on-cg-elements/8629>\
**Category:** dolfinx\
**Created:** [June 27, 2022, 10:08am UTC](https://fenicsproject.discourse.group/t/defining-subdomains-from-gmsh-on-cg-elements/8629 "2022-06-27T10:08:53Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![fpradelli](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/fpradelli/32/3291_2.png) [@fpradelli](https://fenicsproject.discourse.group/u/fpradelli)\
**Post date:** [June 27, 2022, 10:08am UTC](https://fenicsproject.discourse.group/t/defining-subdomains-from-gmsh-on-cg-elements/8629/1 "2022-06-27T10:08:53Z")

</div>

Hi everyone,  
This great [FEniCSx tutorial](https://jorgensd.github.io/dolfinx-tutorial/chapter3/subdomains.html#subdomains-defined-from-external-mesh-data) 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:

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

```

to:

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

```

I get an error:

```auto
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:

```auto
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

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [June 27, 2022, 12:48pm UTC](https://fenicsproject.discourse.group/t/defining-subdomains-from-gmsh-on-cg-elements/8629/2 "2022-06-27T12:48:09Z")

</div>

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.