Interior facet integrals fail when exterior facets indexed using mesh domain data

[Cian Wilson]

When I mark my exterior facets directly (using subdomain.mark_facets) in the mesh domains class my interior facet integrals are ignored by assemble in 1.0-beta2 (and dev). This is caused by both the exterior_facet_domains and interior_facet_domains MeshFunctions get extracted using facet_domains at Assembler.cpp:131 and Assembler.cpp:139 respectively. Using the previous mesh data (MeshFunction) based method I would not have had an interior_facet_domain MeshFunction associated at this point. This wouldn't be a problem except that the sparse data storage of the MeshValueCollection gets converted to dense MeshFunction data storage at MeshDomains.cpp:172 by initializing the MeshFunction to the max uint. This then falls foul of the check at Assembler.cpp:405 and none of my interior facets get assembled. Initializing the MeshFunction to 0 instead (at MeshDomains.cpp:172) fixes this but is probably not always safe?

The main.cpp below (compiled with a demo CMakeLists.txt) and the attached CourantNumber.ufl illustrate this behaviour. The field u should evaluate to 10 everywhere but since the interior facets are being ignored the minimum value is 0.0. The old behaviour can be regained using 1.0-beta and commenting in and out the highlighted lines.

main.cpp:
------------------------------------------------------
#include <dolfin.h>
#include "CourantNumber.h"

using namespace dolfin;

class Side : public SubDomain
{
    public:
    Side(const double &side) : side_(side)
    {
    }

    bool inside(const Array<double>& x, bool on_boundary) const
    {
      return (std::fabs(x[0] - side_) < DOLFIN_EPS && on_boundary);
    }

  private:

    double side_;
};

int main()
{
  UnitInterval mesh(10);

  Side left(0.0);
  Side right(1.0);

  /* begin 1.0-beta2 (and dev) */
  left.mark_facets(mesh, 1);
  right.mark_facets(mesh, 2);
  /* end 1.0-beta2 */

  /* begin 1.0-beta */
// boost::shared_ptr< MeshFunction< dolfin::uint > > edgeids = mesh.data().create_mesh_function("exterior_facet_domains", 0);
// (*edgeids).set_all(0);
// left.mark(*edgeids, 1);
// right.mark(*edgeids, 2);
  /* end 1.0-beta */

  CourantNumber::FunctionSpace V(mesh);

  CourantNumber::BilinearForm a(V,V);
  CourantNumber::LinearForm L(V);

  Constant dt(1.0);
  Constant v(-1.0);
  L.dt = dt;
  L.v = v;

  Function u(V);

  solve(a == L, u);

  std::cout << "max = " << u.vector().max() << ", min = " << u.vector().min() << std::endl;

  File file("courant.pvd");
  file << u;

  plot(u);

  assert(abs(u.vector().min()-10.0)<1.e-6);

  return 0;
}
----------------------------------------------

[Garth Wells (garth-wells)]

Assemblers now only use domain data attached to forms.