Constructing subdomains from knowing the indices of the cells in the mesh

Question

In the examples I've seen so far that deals with constructing subdomains out of a mesh, the subdomains could always be defined fairly simple like, for example:

class Omega0(SubDomain):
   def inside(self, x, on_boundary):
      return True if x[1] <= 0.5 else False

However, in my case, I just can't define inside in this way. My mesh is slightly more complex and I cringe at the thought of having to check the point against all of the subdomain's tetrahedron cell (by this method for example).

My mesh is constructed out of the geometry from another program, so I have the indices for the cells of the various subdomains. Assuming DOLFIN doesn't rearrange the cells and scramble the indices when I load my mesh, there should be some simple way of just marking these cells by their index in the loaded Mesh (or MeshFunction?) to their corresponding subdomain. However I can't seem to figure out how to do this.

As an example, say if have 5 cells in my mesh, and indices 0-2 is in subdomain A and indices 3-4 is in subdomain B.

Answer 1

I would worry about the indices matching up as a reliable assumption, particularly in parallel. Also, some DOFs which are on the edge may not be in the mesh for a high order element?

Perhaps you could keep a sorted list of facets along each boundary? That way you can simply check if the current coordinate lies along any of the facets. Pad in DOLFIN_EPS for floating math, and maybe cache the result for each coordinate if the checks are slow.

Just an idea :)

Overall, though, I think it would be nice to specify the boundary sub domains in a mesh directly as is done in CGNS files.

---

Another option (from the comments) is to load a mesh for the sub domain and check if it contains the DOF. Something like:

    bool foundPoint = false;
    CellIterator cell(subDomainMesh);
    for(; !cell.end(); ++cell)
    {
        if (cell->contains(point))
        {
            foundPoint = true;
        }
    }

Comments

As an aside, I find this issue a bit more encumbering when doing mesh refinement along complex domains.

---

Do you suppose I could generate mesh files for each individual subdomain and somehow combine them in DOLFIN as well as marking the subdomains? Because I could easily generate such subdomain-specific .xml files. The problem then is how to solve any interconnectivity between the meshes.

---

I think you could - just instantiate a new mesh for the subdomain and do a check if it contains the spatial coordinate of the DOF. You could do this by iterating through the subdomain's cells and calling cell.contains(point). This is probably the way to go? :)

---

I solved it with an extended SubDomain class:

# Class for defining subdomains
class DolfinSubDomain(SubDomain):
   def __init__(self, file):
      SubDomain.__init__(self)
      self.mesh = Mesh(file)
   def inside(self, point, on_boundary):
      for c in cells(self.mesh):
         if (c.contains(Point(point[0], point[1], point[2]))):
            return True
      return False

And initiated my subdomains with:

mesh = Mesh('complete_mesh.xml')
subdomains = MeshFunction("size_t", mesh, 3)
subdomain0 = DolfinSubDomain('subdomain_0.xml')
subdomain1 = DolfinSubDomain('subdomain_1.xml')
subdomain0.mark(subdomains, 0)
subdomain1.mark(subdomains, 1)

Which should work fine as long as no subdomains are overlapping I guess.

---

Glad to hear it worked :)

I've been poking around in dolfin's source and came across how they find a cell:

dolfin/function/function.cpp:372 (of dolfin 1.4.0)

  // Find the cell that contains x
  const double* _x = x.data();
  const Point point(mesh.geometry().dim(), _x);

  // FIXME: Testing
  int ID = 0;

  // Get index of first cell containing point
  unsigned int id
    = mesh.bounding_box_tree()->compute_first_entity_collision(point);

  // If not found, use the closest cell
  if (id == std::numeric_limits<unsigned int>::max())
  {
    if (allow_extrapolation)
    {
      id = mesh.bounding_box_tree()->compute_closest_entity(point).first;
      cout << "Extrapolating function value at x = " << point
           << " (not inside domain)." << endl;
    }
    else
    {
      dolfin_error("Function.cpp",
                   "evaluate function at point",
                   "The point is not inside the domain. Consider setting \"allow_extrapolation\" to allow extrapolation");
    }
  }

I imagine it's a faster search algorithm than simple iteration, might be worth a shot if performance is a concern.

I've noticed I have to allow extrapolation in the cases where I have ghost cells (those lovely overlapping cells) :)

---

Hello wensby and Charles,

I was wondering whether you guys have implemented these two methods, and could tell a bit about how slow/fast this performs.

I have used the method by wensby above in a vascular flow simulation, the mesh is not very heavy, and yet the method by wensby above turns out to be quite slow.

Also, Charles, how would I go about trying to implement something like what you have suggested above in python FEniCS ?

Thanks