Computing area of triangles gives different results in parallel

I ran into a strange issue with a code like the following, that tries to compute the area of the triangles that satisfy a certain condition.

// The list of point coordinates in the mesh
auto xx = mesh->geometry().x();

// The dofs that satisfy a certain criterion
std::set<std::size_t> selectedDofs;

// Loop on all the points; could be using mesh->geometry().index_map().local_size() too
for (int i = 0; i < xx.size() / 3; ++i)
    // Shortcuts for a point's coordinates
    auto &x = xx[3 * i + 0];
    auto &y = xx[3 * i + 1];
    auto &z = xx[3 * i + 2];

    // If this point satisfies a certain criterion
    if (criterion(x, y, z) == true)
        // Add it to the list of selected dofs

double area = 0;

// Loop on all the boundary facets given by the indices of a boundary MeshTag
for (const auto i : mt_boundaries.indices())
    // If all the dofs connected to this facet were selected
    if (std::all_of(conn->links(i).begin(),
                    [&](const auto a)
                    { return selectedDofs.contains(a); }))
        // Get the coordinates of the three points
        std::valarray<PetscScalar> p0 = {xx[3 * conn->links(i)[0] + 0],
                                         xx[3 * conn->links(i)[0] + 1],
                                         xx[3 * conn->links(i)[0] + 2]};
        std::valarray<PetscScalar> p1 = {xx[3 * conn->links(i)[1] + 0],
                                         xx[3 * conn->links(i)[1] + 1],
                                         xx[3 * conn->links(i)[1] + 2]};
        std::valarray<PetscScalar> p2 = {xx[3 * conn->links(i)[2] + 0],
                                         xx[3 * conn->links(i)[2] + 1],
                                         xx[3 * conn->links(i)[2] + 2]};

        contactArea += computeArea(p1, p2, p3);

With this code I get the correct area when running in serial but it counts fewer triangles [and thus yields a smaller area] when I run it in parallel. The triangles that are selected form a contiguous region and I suspect this region somehow crosses inter-process boundaries and some triangles are lost, somehow.

Can anyone see anything obviously wrong?

You are making some very strong assumptions on the connection between the mesh topology and the mesh geometry.

Is there as specific reason for not using locate_entities_boundary? You could even do this on the mesh vertices.

Here is a minimal computation in Python illustrating one way of doing what you want:

import dolfinx
from mpi4py import MPI
import numpy as np
import ufl
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 6, 6, 6)

def criterion(x):
    return x[0] <= 0.5 + 1e-13

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

facets = dolfinx.mesh.locate_entities_boundary(
    mesh, mesh.topology.dim-1, criterion)
num_local_facets = mesh.topology.index_map(mesh.topology.dim-1).size_local
local_facets = facets[facets < num_local_facets]  # Count only once

ft = dolfinx.mesh.meshtags(mesh, mesh.topology.dim-1,
                           local_facets, np.ones_like(local_facets))
area = dolfinx.fem.form(
    1*ufl.ds(domain=mesh, subdomain_data=ft, subdomain_id=1))
local_area = dolfinx.fem.assemble_scalar(area)

print(mesh.comm.rank, mesh.comm.allreduce(local_area, MPI.SUM), local_area)

However, if you insist on getting the mesh vertices, you could use entities_to_geometry

which would give you the vertices of each facet directly, so you can send them into xx to get each of the points and use your own algorithm.

geometry_vertices = dolfinx.cpp.mesh.entities_to_geometry(
    mesh, mesh.topology.dim-1, local_facets, False)

where the ith row is the geometry indices of the the ith local_facets.

Hello Jørgen! Thanks for the quick and very helpful reply. I tried the method you suggested, with something like

auto contactFaces = mesh::locate_entities_boundary(
        mesh->topology()->dim() - 1,
        [&](const auto &x)
            std::vector<std::int8_t> marked;
            for (const auto &j :
                marked.push_back(criterion(x(0, j), x(1, j), x(2, j)));

            return marked;

                [n = mesh->topology()
                                ->index_map(mesh->topology()->dim() - 1)
                                ->size_local()](const auto i)
                { return i >= n; }),

std::vector<int> values(contactFaces.size(), 1);

auto mt_contact =
        mesh::MeshTags<int>(mesh->topology(), 2, contactFaces, values);

// Then integrate surface area...

but I’m getting the same exact number of triangles as with my previous method and the same total surface area :thinking:

If you could make a minimal reproducible example, I could help you further.
Running the code I supplied in the previous post, gives a consistent result on one or more processes, so there is either something wrong with the mesh, or there is something I’m not grasping in the C++ code.

For the sake of simple debugging, code you move it to Python?:smiley:

I found the solution. In my MWE above I did not show that I was moving the mesh nodes by writing into mesh->geometry().x(). I disabled the mesh modification and it started to give me good numbers, so I tracked it down to me not moving all the nodes.
Quite confusing explanation, I know, but at least it works now. Thanks for the help!