FEniCS / FEniCS/dolfinx

More than one mesh tag (per integral type) for a single form

Open
#3,663 0 comments 0 reactions 0 assignees View on GitHub
enhancement
Dominant language
C++
Stars
1.2k
Forks
261
Avg merge
1d 15h
Merged PRs (30d)
65

Description

### Describe new/missing feature

In some problems it may be handy to define forms using different mesh tags (whose entities may intersect).
For instance, imagine that you want to build a form in which you assemble the same operator everywhere, but you have to include some extra terms for a few cells/facets (see the example below).
First possibility would be to define a single mesh tags object with two tags: for cells/facets with and without extra terms.
However, including additional terms to your form, the number of different cell/facet groups grows exponentially.

The most natural solution would be define several mesh tags and use them all in the same form (see the code below). Unfortunately, computing quantities with two (or more) `MeshTags` and the same `IntegralType` is not ideal as they cannot be combined in a single DOLFINx form.
Two (or more) forms need to be generated, which involves calling the assembler more than once and generating several (matrices, vectors) that need to be summed.

I am looking for possible alternatives:

1. Define mesh tags that contain more than one value per cell. This way a single `MeshTags` object would be needed and therefore 1 single form is created.
Just looking at the code in `MeshTags.h`, I can’t see any particular reason why this shouldn’t work (apart from the built-in tests for checking for duplicates). I am afraid that `MeshTags` is used in a lot of places and extesions and breaking this “no duplicate entities” principle might blow up somewhere else.
As mentioned in this [slack discussion](https://fenicsproject.slack.com/archives/C08BPUCUAT0/p1741779201232879), another possibility would be associate every entity in the `MeshTags` to a `std::pair` of `T`

2. Allow more than one `MeshTags` to be used in the same DOLFINx form. At the end of the pipeline, the (C++) `Form` object simply stores (for each integral type) a container mapping tag ids to the associated entities (these ids and entities could initially come from different `MeshTags` objects). So again, apart from the [Python code](https://github.com/FEniCS/dolfinx/blob/4627a17522e73020dc7982e1a030da5297605f09/python/dolfinx/fem/forms.py#L327C1-L329C70) that checks that only one `MeshTags` object is present in the form, I don’t see any particular reason why this shouldn’t work.
In fact, I introduced some changes in [`forms.py`](https://github.com/FEniCS/dolfinx/blob/4627a17522e73020dc7982e1a030da5297605f09/python/dolfinx/fem/forms.py) and everything worked smoothly.

### Suggested user interface

```python3
from mpi4py import MPI
import ufl
from dolfinx import fem, mesh

domain = mesh.create_unit_square(MPI.COMM_WORLD, 2, 1, mesh.CellType.quadrilateral)

x = ufl.SpatialCoordinate(domain)
term_0 = 1 + x[0]
term_1 = 1 + x[1]
dx = ufl.dx(domain=domain)

# Current version
tags_0 = mesh.meshtags(domain, 2, np.array([0, 1]), np.array([0, 1]))
dx_0 = dx(subdomain_data=tags_0)

ufl_form_0 = term_0 * (dx_0(0) + dx_0(1)) + term_1 * dx_0(1)
form_0 = fem.form(ufl_form_0)

# Proposed version
tags_1 = mesh.meshtags(domain, 2, np.array([0, 1]), 1)
tags_2 = mesh.meshtags(domain, 2, np.array([1]), 2)

dx_1 = dx(subdomain_data=tags_1)
dx_2 = dx(subdomain_data=tags_2)

ufl_form_1 = term_0 * dx_1(1) + term_1 * dx_2(2)
form_1 = fem.form(ufl_form_1) # Error
```

Contributor guide

Open the contributing guide

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.