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.
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!