# Computing area of triangles gives different results in parallel

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
[&](const auto a)
{ return selectedDofs.contains(a); }))
{
// Get the coordinates of the three points
std::valarray<PetscScalar> p0 = {xx[3 * conn->links(i) + 0],
std::valarray<PetscScalar> p1 = {xx[3 * conn->links(i) + 0],
std::valarray<PetscScalar> p2 = {xx[3 * conn->links(i) + 0],

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.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 `i`th row is the geometry indices of the the `i`th `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,
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 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? 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!