Hello, I have two questions

1. How is the mesh made with Fenics ? If for example I take a sheet of paper and I draw a square, how will he proceed to create the mesh with horizontal and/or vertical lines? (like in Eshelby problem)

2. I want to plot the mesh with a color code to locate the inclusion and the matrix in the Eshelby Problem. With the code below I get the image that I show again below. But the inclusion is not colored.

``````from dolfinx.plot import create_vtk_mesh
pyvista.set_jupyter_backend("none")
grid = pyvista.UnstructuredGrid(*create_vtk_mesh(domain))
plotter = pyvista.Plotter()
plotter.show_bounds(grid='front', location='outer', all_edges=True)
plotter.view_xy()
plotter.show()
``````

What I want to achieve is something like this: I don’t get this result.

The source code for mesh generation is here:

The code you have supplied is not complete, in the sense that one cannot reproduce the plots, as you haven’t supplied the code that defines your `domain` variable.

Note that you can send in your vertex and cells in any order you would like, see: Mesh creation in serial and parallel — FEniCSx Documentation

Here is a complete code :

``````import dolfinx # FEM in python
import matplotlib.pyplot as plt
import ufl # variational formulations
import numpy as np
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import gmsh # Mesh generation
import pyvista
import extract

# Geometry
R_i = 1.0 # Radius of the inclusion
R_e = 6.9  # Radius of the matrix (whole domain)

# Material
E_m = 0.8 # Young's modulus in matrix
nu_m = 0.35 # Poisson's ratio in matrix
E_i = 11.0 # Young's modulus of inclusion
nu_i = 0.3 # Poisson's ratio in inclusion

#CREATE THE MESH WITH GMSH

mesh_size = R_i/5
mesh_order = 1

mesh_comm = MPI.COMM_WORLD
model_rank = 0
gmsh.initialize()
facet_names = {"inner_boundary": 1, "outer_boundary": 2}
cell_names = {"inclusion": 1, "matrix": 2}
model = gmsh.model()
model.setCurrent("Disk")
gdim = 2 # geometric dimension of the mesh
inner_disk = gmsh.model.occ.addDisk(0, 0, 0, R_i, aspect_ratio * R_i)
outer_disk = gmsh.model.occ.addDisk(0, 0, 0, R_e, R_e)
whole_domain = gmsh.model.occ.fragment([(gdim, outer_disk)], [(gdim, inner_disk)])
gmsh.model.occ.synchronize()

# Add physical tag for bulk
inner_domain = whole_domain
outer_domain = whole_domain
model.setPhysicalName(gdim, inner_domain, "Inclusion")
model.setPhysicalName(gdim, outer_domain, "Matrix")

# Add physical tag for boundaries
lines = gmsh.model.getEntities(dim=1)
inner_boundary = lines
outer_boundary = lines
gmsh.option.setNumber("Mesh.CharacteristicLengthMin",mesh_size)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax",mesh_size)
model.mesh.generate(gdim)
gmsh.option.setNumber("General.Terminal", 1)
model.mesh.setOrder(mesh_order)
gmsh.option.setNumber("General.Terminal", 0)

# Import the mesh in dolfinx
from dolfinx.io import gmshio
domain, cell_tags, facet_tags = gmshio.model_to_mesh(model, mesh_comm, model_rank, gdim=gdim)
domain.name = "composite"
cell_tags.name = f"{domain.name}_cells"
facet_tags.name = f"{domain.name}_facets"
gmsh.finalize()

cell_names["matrix"]
cell_names["inclusion"]

from dolfinx.plot import create_vtk_mesh
pyvista.set_jupyter_backend("none")

grid = pyvista.UnstructuredGrid(*create_vtk_mesh(domain))
plotter = pyvista.Plotter()
plotter.show_bounds(grid='front', location='outer', all_edges=True)
plotter.view_xy()
plotter.show()
``````

I would strongly consider looking at the tutorials:
https://jsdokken.com/dolfinx-tutorial/chapter3/em.html

As it shows how to plot the cell markers

``````

from dolfinx.plot import create_vtk_mesh
pyvista.set_jupyter_backend("static")
pyvista.start_xvfb()
grid = pyvista.UnstructuredGrid(*create_vtk_mesh(domain))
plotter = pyvista.Plotter()
num_local_cells = domain.topology.index_map(domain.topology.dim).size_local
grid.cell_data["Marker"] = cell_tags.values[cell_tags.indices<num_local_cells]
grid.set_active_scalars("Marker")