What is the more succint way to get the "cell" and "facet" tags?

mesh
1

Hello,

Is there is a more compact way of getting “facet” and “cell” tags (topology.dim - 1 and topology.dim, respectively)? I did this on DOLFINx 0.8.0. Thanks.

[Edit 1]
I just want to be able to access all the facet and cell tags in the domain without creating a boundary condition or function (e.g. to be used with locate_entities).

from dolfinx import mesh
from mpi4py import MPI
import numpy as np

# Mesh
msh = mesh.create_unit_square(MPI.COMM_WORLD, 5, 5)
gdim = msh.topology.dim

# Facets?
# https://fenicsproject.discourse.group/t/
#   transfer-meshtags-to-submesh-in-dolfinx/8952/7
fdim = msh.topology.dim - 1
msh.topology.create_connectivity(fdim, gdim)
facet_map = msh.topology.index_map(fdim)
facet_total = facet_map.size_local + facet_map.num_ghosts
facet_tags = np.zeros(facet_total, dtype=np.int32)

# Cells?
cell_map = msh.topology.index_map(gdim)
cell_total = cell_map.size_local + cell_map.num_ghosts
cell_tags = np.zeros(cell_total, dtype=np.int32)

# Left boundary (just to test)
left_facets = mesh.locate_entities_boundary(
    msh, fdim, lambda x: np.isclose(x[0], 0))
#  Mark boundary
left_tags = mesh.meshtags(msh, fdim, left_facets,
                          np.full_like(left_facets, np.int32(1)))
left_tags.indices == left_facets
2

I’m not sure what you’re asking for here. You can generate all of the facet indices with, e.g.,

facet_idxs = np.arange(facet_total, dtype=np.int32)
1 Like
3

As @nate said, it is not clear what you are asking for here. Do you only want to mark a subset of facets with a single value?
If so one can simplify your code quite alot.

1 Like
4

Dear @nate and @dokken,

Summary

If I have this, how do I get all the facet and cell tags? (are there any?). If not how do I create them?

  msh = mesh.create_unit_square(MPI.COMM_WORLD, 5, 5)

Details

When creating a mesh with an external program (e.g. Gmsh), one usually ends with (1) a mesh, (2) cell tags and (3) facet tags:

::: [from_gmsh_object_to_mesh_in_dolfinx.py]

mesh_gmsh, cell_tags, facet_tags = dolfinx.io.gmshio.model_to_mesh(
      gmsh.model, mesh_comm, model_rank, gdim=2)

:::

(a similar result can be achieved with XDMFFile.read_meshtags and XDMFFile.read_mesh).

Recently, I was debugging some code, and tried to create just a square with msh = mesh.create_unit_square.

msh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 5, 5)

The rest of the code required a dolfinx.mesh.MeshTags object (what [from_gmsh_object_to_mesh_in_dolfinx.py] shows as facet_tags). In particular

ufl.Measure("ds", domain=msh, subdomain_data=facet_tags)

Thus, I needed the tags, but I didn’t know how to get them. The documentation of mesh.MeshTags states that one should not do

tmp = mesh.MeshTags(msh)

(or so I understood :P; doing that results in something weird–I tried anyway :DD ). I knew that one can add tags to the facets by means of functions (as shown above), but I was wondering if the tags are already there; if they can be accessed. Without reading the documentation, I assumed that dolfinx.mesh.exterior_facet_indices would only give the facets around the domain, not all the tags (including the ones inside).

My question is just how do I get the tags of all the facets (without a def foo or lambda). That’s why I think it’s a simpleton’s question, I just don’t know how q: .

If you don’t mind, now I would also be interested in reading about this:

Thank you very much!

5

I’m still not sure I understand your question…

Mesh tags don’t belong to the mesh. Mesh tags are just a way to associate mesh entity indices of a given dimension with given values.

The in-built mesh generation, e.g., dolfinx.mesh.create_unit_square, dolfinx.mesh.create_unit_cube, etc, are just meshes. They don’t have any mesh tag information to share because why would they? It’s up to the user to decide how they want to annotate the mesh with information specific to their problem.

The Gmsh interface provided by dolfinx returns the entity tags as created in the mesh generation process. This is particularly useful as one typically works with non trivial domain descriptions and using geometric location of entities, as you describe above, is not feasible.

1 Like
6

Dear @nate. I think that you got the question very well:

I wanted tags which don’t exist. So, if I want to use

I have to create them first, most likely, with something like this

from dolfinx import mesh
from mpi4py import MPI
import numpy as np

# Mesh
msh = mesh.create_unit_square(MPI.COMM_WORLD, 5, 5)
gdim = msh.topology.dim

# Facets
fdim = msh.topology.dim - 1
msh.topology.create_connectivity(fdim, gdim)
facet_map = msh.topology.index_map(fdim)
facet_total = facet_map.size_local + facet_map.num_ghosts
# Indices of all facets
all_facets_idx = np.arange(facet_total, dtype=np.int32)
# # The above is equivalent to, but possibly better (due to MPI) than
# all_facets_idx = mesh.locate_entities(msh, fdim, lambda x: np.full(x.shape[1], True))
facet_tags = mesh.meshtags(msh, fdim, all_facets_idx, values=5)

Now, I can access the tags and see that they have a value of 5 (and change them if I need):

print(facet_tags.values)

[5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
 5 5 5 5 5 5 5 5 5 5 5]

I just didn’t know. It was through the process of asking this question that I came to the answer :- ) . (I will leave the shorter solution of creating a simple array with a value of 5 to the interested reader). Thanks! (Sorry for the confusion. I hope I didn’t take too much of your time).

(Should I mark my own answer as the solution?)