Assigning function values using MeshFunction in parallel

Legacy DOLFIN
1

Hi,

I have a script for assigning DG function values based on a MeshFunction which works in serial.

Is there a parallel safe way to achieve the same behaviour? Currently the last two lines are not compatible with MPI since they operate on .vector().

Thanks,
Kiri

import numpy as np
from dolfin import *

# Dummy subdomain function
class Omega0(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[0] <= 0.5 else False


mesh = UnitSquareMesh(10, 10, diagonal='crossed')
mf = MeshFunction("size_t", mesh, mesh.topology().dim())
U = FunctionSpace(mesh, 'DG', 0)
u = Function(U)

# Values to be assigned based on subdomain values
vals = [100, -100]

left = Omega0()
left.mark(mf, 1)

##### Is there a parallel safe equivalent of the following? #######
help_func = np.asarray(mf.array(), dtype=np.int32)
u.vector()[:] = np.choose(help_func, vals)
3

To me it is unclear what you mean is not parallel safe with this?
The only thing you are missing (i believe), is a
u.vector().apply("insert") to update potential ghost cells, depending on the ghost mode.

Please note that your problem, as it is currently written is not showing any errors when executed in parallel.

4

Thanks for your reply @Dokken.

When running my script using MPI I was getting an error on np.choose which I assumed was due to calling .array() in parallel. After your message I have been debugging further and it seems that the values in the mf.array() are behaving weirdly in parallel.

I my use case I have an imported mesh and MeshFunction where the above script runs fine in serial, e.g. m_f \in [0,1] where m_f is mf.array(), so I get np.min(mf.array()) = 0 and np.max(mf.array()) = 1 in serial.

In parallel however, I get np.min(mf.array()) = 0 and np.min(mf.array()) = 18446744073709551615. I assume the latter is some sort of overflow error?

Do you know what could be causing this and how I can avoid this? I can share a mesh and MeshFunction if you think that would help in the debugging process.

Best,
Kiri

5

That is because the default value in MeshFunction is -1, which is equivalent to a huge number.
Have you tried with
mf = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
to initialize all entries of the meshfunction?

6

Thanks for your quick reply!

That seems like it could potentially be an issue. However, I am loading in a MeshFunction from file so I would have assumed there would be no need to set default values since all values should be pre-defined from the .xdmf file?

If my assumption is incorrect, could you please point me to the syntax for assigning default values when loading a MeshFunction from file?

My current script for loading my mesh and MeshFunction:

# Load mesh
mesh = Mesh()
with XDMFFile(comm, mesh_path) as infile:
            infile.read(mesh)

# Load mesh function
mvc1 = MeshValueCollection("size_t", mesh, gdim)
    with XDMFFile(meshfunc_filepath) as infile:
        infile.read(mvc1, "name_to_read")

vols = MeshFunction("size_t", mesh, mvc1) #  MeshFunction("size_t", mesh, mvc1, value=0) does not work

For additional context, I am generating my MeshFunction by converting a .gmsh file using the function shown below (this is done in serial). Could this be the issue?

def convert_msh_to_xdmf(gmsh_mesh_filename: str, gdim: int, scaling: float = 1):
    """Converts gmsh mesh file (.msh) to .xdmf format for use in FEniCS"""
    mesh = meshio.read(gmsh_mesh_filename)

    assert gdim == 3

    facet = "triangle"
    element = "tetra"
    
    element_data = []
    facet_data = []

    for cell in mesh.cells:
        if cell.type == element:
            element_data.append(cell.data)

        elif cell.type == facet:
            facet_data.append(cell.data)

    element_data = np.vstack(element_data)

    facet_data = np.vstack(facet_data)

    points = scaling * mesh.points

    for key in mesh.cell_data_dict["gmsh:physical"].keys():
        if key == element:
            elements = mesh.cell_data_dict["gmsh:physical"][key]

        elif key == facet:
            facets = mesh.cell_data_dict["gmsh:physical"][key]

    element_mesh = meshio.Mesh(
        points=points,
        cells=[(element, element_data)],
        cell_data={"name_to_read": [elements]},
    )

    facet_mesh = meshio.Mesh(
        points=points, cells=[(facet, facet_data)], cell_data={"name_to_read": [facets]}
    )

    name = gmsh_mesh_filename.split(".")[0]
    meshio.write(f"{name}.xdmf", element_mesh)
    meshio.write(f"{name}_facets.xdmf", facet_mesh)

    return len(points), len(elements)

Thanks for your help @dokken!

7

a meshvalue collection does not contain all sub entities marked (think interior facets of your domain).
What you can do, is to indentify what values in mf.array() that are out of your range or negative, and replace them with appropriate values (0).

8

That’s interesting, thank you! The only thing I am not clear on is how to choose what to replace the out-of-range value(s) with, or does this not matter?

E.g. if m_f \in [0,1] but there is a -1 in mf.array(), how know whether this should be assigned 0 or 1?

9

I would choose 0. It would indicate indices that are not marked in your input mesh.

10

Hi @dokken, I have been tinkering with this but have still not managed to get the assigning to work correctly when using MPI.

My current issue is that for some reason len(mf.array()) != len(u.vector()[:]) resulting in Unable to set local values of PETSc vector. error on the u.vector()[:] = np.choose(help_func, vals) line.

Given mf is defined on a cell-wise manner and u is in DG0 I would have assumed these should have identical lengths on each rank.

Do you have any suggestions for things I could try to debug this further?

Many thanks,
K

11

You could ignore the ghost entries of the mesh function, and only set data to the first M values (len(u.vector()[:])), then use vector().apply("insert") afterwards to update ghost values.

1 Like
12

Thanks @dokken - does this imply the entries in mf.array()[M:], where M = len(u.vector()[:]), are always guaranteed to be ghost values?

13

Yes. it does.
as one has seen before there is a one to one correspondence between cell indexing and a DG-0 space, where ghosts are located at the end on each process.

14

Amazing - thanks for clarifying!