Check ownership of Facet

Question

Having

f = Facet(mesh, 42)

can I determine whether f belongs also to some other process? Is this information stored somewhere in mesh partitioning data or do I need to perform some MPI communication by myself?

Answer 1

Use the MeshEntity function

std::size_t num_global_entities(std::size_t dim) const

From Python, try:

dim = mesh.topology().dim()
mesh.init(dim-1, dim)
f = Facet(mesh, 42)
num_incident = f.num_global_entities(dim)
num_local = f.num_entities(dim)

If num_incident == 2, then two cells (globally) share the facet. On the other hand if num_local == 1, then there is only one incident cell on this partition. Therefore f is shared iff num_incident != num_local.

Comments

edited to give a full answer

Answer 2

In C++ you can do:

const std::size_t dim = mesh.topology().dim() - 1;
std::map<unsigned int, std::set<unsigned int> >& shared_facets =  mesh.topology().shared_entities(dim);

to get a map of shared facets. Not sure about Python.

Comments

In python mesh.topology().shared_entities(dim) returns useless

<Swig Object of type 'std::map< unsigned int,std::set< unsigned int > > *' at address>

with recent master.

---

Useless? What are you talking about? Simply a matter of:

from dolfin import *

mesh = UnitSquareMesh(4, 4)
code = """
#include "dolfin.h"
namespace dolfin{

  std::vector<unsigned int> get_keys(std::map<unsigned int,std::set<unsigned int > >& themap )
  {
    std::vector<unsigned int> keys;
    for (std::map<unsigned int, std::set<unsigned int> >::iterator it=themap.begin(); it!= themap.end(); it++)
       keys.push_back(it->first);
    return keys;
  }  

  std::vector<unsigned int> get_vals(std::map<unsigned int,std::set<unsigned int > >& themap , unsigned int i)
  {
    std::vector<unsigned int> vals;
    for (std::set<unsigned int>::iterator it=themap[i].begin(); it!= themap[i].end(); it++)
      vals.push_back(*it);
    return vals;
  }  
}
"""
cm = compile_extension_module(code)
m = mesh.topology().shared_entities(0)
keys = cm.get_keys(m)
for k in keys:
    print "Vertex ", k, " shared with = ", cm.get_vals(m, k)[:]

Or perhaps Johan making a map would be a better idea:-)

---

OK, I take it back - it's definitely useful. But typemap would be convenient -> filed as issue 73