Bug in `locate_subdomain_entities`

[RemDelaporteMathurin]

https://github.com/festim-dev/FESTIM/blob/013937a860d72252436118dabcbf75b80ebc5cce/src/festim/subdomain/volume_subdomain.py#L32-L43

This will always return all the cells of the mesh, not the subdomain cells!!

We should use meshtags.find()

[RemDelaporteMathurin]

I guess this works with the mixed domain approach since each subdomain has its own submesh.... this will not work with HydrogenTransportMaterial with several materials

[RemDelaporteMathurin]

Actually running this with HydrogenTransportMaterial with several materials crashes:

from mpi4py import MPI

import dolfinx
import dolfinx.fem.petsc
import numpy as np

import festim as F

def generate_mesh(n=20):

    def bottom_boundary(x):
        return np.isclose(x[1], 0.0)

    def top_boundary(x):
        return np.isclose(x[1], 1.0)

    def half(x):
        return x[1] <= 0.5 + 1e-14

    mesh = dolfinx.mesh.create_unit_square(
        MPI.COMM_WORLD, n, n, dolfinx.mesh.CellType.triangle
    )

    # Split domain in half and set an interface tag of 5
    gdim = mesh.geometry.dim
    tdim = mesh.topology.dim
    fdim = tdim - 1
    top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary)
    bottom_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, bottom_boundary)
    num_facets_local = (
        mesh.topology.index_map(fdim).size_local
        + mesh.topology.index_map(fdim).num_ghosts
    )
    facets = np.arange(num_facets_local, dtype=np.int32)
    values = np.full_like(facets, 0, dtype=np.int32)
    values[top_facets] = 1
    values[bottom_facets] = 2

    bottom_cells = dolfinx.mesh.locate_entities(mesh, tdim, half)
    num_cells_local = (
        mesh.topology.index_map(tdim).size_local
        + mesh.topology.index_map(tdim).num_ghosts
    )
    cells = np.full(num_cells_local, 4, dtype=np.int32)
    cells[bottom_cells] = 3
    ct = dolfinx.mesh.meshtags(
        mesh, tdim, np.arange(num_cells_local, dtype=np.int32), cells
    )
    all_b_facets = dolfinx.mesh.compute_incident_entities(
        mesh.topology, ct.find(3), tdim, fdim
    )
    all_t_facets = dolfinx.mesh.compute_incident_entities(
        mesh.topology, ct.find(4), tdim, fdim
    )
    interface = np.intersect1d(all_b_facets, all_t_facets)
    values[interface] = 5

    mt = dolfinx.mesh.meshtags(mesh, mesh.topology.dim - 1, facets, values)
    return mesh, mt, ct

mesh, mt, ct = generate_mesh(10)

my_model = F.HydrogenTransportProblemDiscontinuousChangeVar()
my_model = F.HydrogenTransportProblem()
my_model.mesh = F.Mesh(mesh)
my_model.volume_meshtags = ct
my_model.facet_meshtags = mt

material_bottom = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0)
material_top = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0)

top_domain = F.VolumeSubdomain(4, material=material_top)
bottom_domain = F.VolumeSubdomain(3, material=material_bottom)

top_surface = F.SurfaceSubdomain(id=1)
bottom_surface = F.SurfaceSubdomain(id=2)
my_model.subdomains = [
    bottom_domain,
    top_domain,
    top_surface,
    bottom_surface,
]

H = F.SpeciesChangeVar("H", mobile=True)

my_model.species = [H]

my_model.boundary_conditions = [
    F.FixedConcentrationBC(top_surface, value=1, species=H),
    F.FixedConcentrationBC(bottom_surface, value=0, species=H),
]

my_model.temperature = 500.0

my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)
max_top = F.MaximumVolume(field=H, volume=top_domain)
max_bottom = F.MaximumVolume(field=H, volume=bottom_domain)
my_model.exports = [max_top, max_bottom]

my_model.initialise()
my_model.run()

assert max_top.value != max_bottom.value

Produces:

Traceback (most recent call last):
  File "/home/remidm/FESTIM/mwe_bug.py", line 103, in <module>
    my_model.run()
  File "/home/remidm/FESTIM/src/festim/problem.py", line 151, in run
    self.post_processing()
  File "/home/remidm/FESTIM/src/festim/hydrogen_transport_problem.py", line 774, in post_processing
    export.compute()
  File "/home/remidm/FESTIM/src/festim/exports/maximum_volume.py", line 30, in compute
    self.value = np.max(self.field.solution.x.array[indices])
                        ~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^
IndexError: index 121 is out of bounds for axis 0 with size 121