Hi
I am marking all interior facets with the ID of the cell where the facet belongs. The logic I am using when looping over all facets is:
- If facet is on boundary, skip it
- If the two neighboring cells of each facet have the same subdomain ID, I mark the facet with the ID
- If they differ, I leave it as is
Here’s the piece of Cpp code string (I am writing for pybind11):
void mark_1d_facets_2d_mesh(const Mesh & mesh,
MeshFunction<std::size_t> & cell_subdomains,
MeshFunction<std::size_t> & facet_subdomains) {
const std::size_t dim = mesh.topology().dim();
mesh.init(dim - 1);
mesh.init(dim - 1, dim);
for (FacetIterator facet(mesh); !facet.end(); ++facet) {
if (facet->exterior())
continue;
auto i1 = facet->entities(dim)[0];
auto i2 = facet->entities(dim)[1];
auto id1 = cell_subdomains[i1];
auto id2 = cell_subdomains[i2];
if (id1 == id2)
facet_subdomains[facet->index()] = id1;
}
}
Here is how I construct my mesh and test it
from dolfin import *
mesh = UnitSquareMesh(6, 6)
subs = MeshFunction('size_t', mesh, 2)
facets = MeshFunction('size_t', mesh, 1)
CompiledSubDomain('x[0]<=0.5').mark(subs, 11)
CompiledSubDomain('x[0]>=0.5').mark(subs, 22)
CompiledSubDomain('on_boundary').mark(facets, 3)
CompiledSubDomain('!on_boundary && near(x[0], 0.5)').mark(facets, 33)
mark_1d_facets_2d_mesh(mesh, subs, facets)
print(Process RANK: facets: Counter(facets.array()), flush=True)
Theoretically, there will be 45 facets with ID==11, 45 facets with ID==22, 6 facets with ID==33, and 24 facets with ID==3, and I get the correct solution in serial:
Process 0: facets: Counter({11: 45, 22: 45, 3: 24, 33: 6})
However, when I use MPI with 2 cores, what I got was:
Process 0: facets: Counter({22: 29, 11: 16, 3: 13, 33: 3, 0: 1})
Process 1: facets: Counter({11: 31, 22: 18, 3: 11, 33: 3, 0: 1})
Well, apparently, this is not correct answer for the interior facets: 11 counts as 47 and 22 also counts as 47. And there are 0, that are not supposed to be
Thanks in advance
Victor