Computing area of triangles gives different results in parallel

dolfinx
1

Hello!
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
        selectedDofs.insert(i);
    }
}

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(),
                    conn->links(i).end(),
                    [&](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?

2

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)
print(geometry_vertices)

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

3

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,
        mesh->topology()->dim() - 1,
        [&](const auto &x)
        {
            std::vector<std::int8_t> marked;
            marked.reserve(x.extent(1));
            for (const auto &j :
                    std::ranges::views::iota(static_cast<size_t>(0),
                                            x.extent(1)))
            {
                marked.push_back(criterion(x(0, j), x(1, j), x(2, j)));
            }

            return marked;
        });

contactFaces.erase(
        std::remove_if(
                contactFaces.begin(),
                contactFaces.end(),
                [n = mesh->topology()
                                ->index_map(mesh->topology()->dim() - 1)
                                ->size_local()](const auto i)
                { return i >= n; }),
        contactFaces.end());

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:

4

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:

5

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!