Fix consistent restriction of one-sided interior facet integrals

cpp/dolfinx/fem/Form.h

@@ -107,6 +107,24 @@ class Form
     if (!_mesh)
       throw std::runtime_error("No mesh could be associated with the Form.");
 
+    // Extract mesh tag for cells (used for consistent restriction of interior
+    // facet integrals)
+    std::shared_ptr<const mesh::MeshTags<int>> cell_marker = nullptr;
+    for (auto& integral_type : integrals)
+    {
+      const IntegralType type = integral_type.first;
+      switch (type)
+      {
+      case IntegralType::cell:
+        if (integral_type.second.second)
+        {
+          cell_marker = std::make_shared<const mesh::MeshTags<int>>(
+              *integral_type.second.second);
+        }
+        break;
+      }
+    }
+
     // Store kernels, looping over integrals by domain type (dimension)
     for (auto& integral_type : integrals)
     {
@@ -137,7 +155,7 @@ class Form
       if (integral_type.second.second)
       {
         assert(_mesh == integral_type.second.second->mesh());
-        set_domains(type, *integral_type.second.second);
+        set_domains(type, *integral_type.second.second, cell_marker);
       }
     }
 
@@ -463,7 +481,7 @@ class Form
                                                       std::int32_t, int>>>>&
           integrals,
       const iterator& tagged_facets_begin, const iterator& tagged_facets_end,
-      const std::vector<int>& tags)
+      const std::vector<int>& tags, const xtl::span<const int>& cell_vals)
   {
     int tdim = topology.dim();
     auto f_to_c = topology.connectivity(tdim - 1, tdim);
@@ -479,8 +497,19 @@ class Form
         {
           std::array<std::pair<std::int32_t, int>, 2> pairs
               = get_cell_local_facet_pairs<2>(*f, f_to_c->links(*f), *c_to_f);
-          it->second.second.emplace_back(pairs[0].first, pairs[0].second,
-                                         pairs[1].first, pairs[1].second);
+
+          // Consistently order integral pairs such that "+" restrictions are
+          // taken on the side with the largest mesh tag value
+          if (cell_vals[pairs[0].first] < cell_vals[pairs[1].first])
+          {
+            it->second.second.emplace_back(pairs[1].first, pairs[1].second,
+                                           pairs[0].first, pairs[0].second);
+          }
+          else
+          {
+            it->second.second.emplace_back(pairs[0].first, pairs[0].second,
+                                           pairs[1].first, pairs[1].second);
+          }
         }
       }
     }
@@ -489,10 +518,13 @@ class Form
   // Sets the entity indices to assemble over for kernels with a domain
   // ID
   // @param[in] type Integral type
-  // @param[in] marker MeshTags with domain ID. Entities with marker 'i'
-  // will be assembled over using the kernel with ID 'i'. The MeshTags
-  // is not stored.
-  void set_domains(IntegralType type, const mesh::MeshTags<int>& marker)
+  // @param[in] marker MeshTags with domain ID. Entities with marker 'i' will be
+  // assembled over using the kernel with ID 'i'. The MeshTags is not stored.
+  // @param[in] cell_marker MeshTags for cells. Used to consistently order
+  // interior facet intgrals.
+  void set_domains(IntegralType type, const mesh::MeshTags<int>& marker,
+                   const std::shared_ptr<const mesh::MeshTags<int>>& cell_marker
+                   = nullptr)
   {
     std::shared_ptr<const mesh::Mesh> mesh = marker.mesh();
     const mesh::Topology& topology = mesh->topology();
@@ -527,9 +559,22 @@ class Form
                                    tagged_entities.cbegin(), entity_end, tags);
         break;
       case IntegralType::interior_facet:
+      {
+        // Unpack cell tags for consistent orientation of "+" and "-"
+        // restrictions
+        const std::int32_t num_cells = topology.index_map(tdim)->size_local()
+                                       + topology.index_map(tdim)->size_local();
+        std::vector<int> cell_vals(num_cells, 0);
+        if (cell_marker)
+        {
+          for (std::size_t i = 0; i < cell_marker->indices().size(); i++)
+            cell_vals[i] = cell_marker->values()[i];
+        }
         set_interior_facet_domains(topology, _interior_facet_integrals,
-                                   tagged_entities.cbegin(), entity_end, tags);
-        break;
+                                   tagged_entities.cbegin(), entity_end, tags,
+                                   cell_vals);
+      }
+      break;
       default:
         throw std::runtime_error(
             "Cannot set domains. Integral type not supported.");

python/test/unit/fem/test_assemble_domains.py

@@ -1,4 +1,4 @@
-# Copyright (C) 2019 Chris Richardson
+# Copyright (C) 2019-2022 Chris Richardson, Jørgen S. Dokken
 #
 # This file is part of DOLFINx (https://www.fenicsproject.org)
 #
@@ -14,8 +14,10 @@
                          dirichletbc, form)
 from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
                                set_bc)
-from dolfinx.mesh import (GhostMode, Mesh, create_unit_square, meshtags, meshtags_from_entities,
-                          locate_entities_boundary)
+from dolfinx.mesh import (GhostMode, Mesh, create_unit_cube,
+                          create_unit_square, locate_entities,
+                          locate_entities_boundary, meshtags,
+                          meshtags_from_entities)
 
 from mpi4py import MPI
 from petsc4py import PETSc
@@ -221,3 +223,51 @@ def test_additivity(mode):
     assert (J1 + J3) == pytest.approx(J13)
     assert (J2 + J3) == pytest.approx(J23)
     assert (J1 + J2 + J3) == pytest.approx(J123)
+
+
+def test_interior_facet_restriction():
+    """ Test that adding a cell marker to interior facet integral orients the '+' and '-' restrictions"""
+    mesh = create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)
+    tdim = mesh.topology.dim
+    fdim = tdim - 1
+
+    plus_restriction, minus_restriction, facet_marker = 3, 5, 1
+    assert(plus_restriction < minus_restriction)
+    minus_cells = locate_entities(mesh, tdim, lambda x: x[1] <= 0.5 + 1e-16)
+    cell_map = mesh.topology.index_map(tdim)
+    all_cells = np.arange(cell_map.size_local + cell_map.num_ghosts, dtype=np.int32)
+
+    # All cells with x[1]>0.5 is tagged with plus restriction
+    cell_marker = np.full(len(all_cells), plus_restriction, dtype=np.int32)
+    # All cells with x[1]<=0.5 tagged with minus restriction
+    cell_marker[minus_cells] = 1
+    ct = meshtags(mesh, tdim, all_cells, cell_marker)
+
+    # Tag all facets on surface between tags with facet marker
+    midline_facets = locate_entities(mesh, fdim, lambda x: np.isclose(x[1], 0.5))
+    argsort_f = np.argsort(midline_facets)
+    ft = meshtags(mesh, fdim, midline_facets[argsort_f],
+                  np.full(len(midline_facets), facet_marker, dtype=np.int32))
+
+    # Interpolate non-zero value on minus restriction
+    V = FunctionSpace(mesh, ("DG", 0))
+    u = Function(V)
+    u.x.array[:] = 0
+    u.interpolate(lambda x: np.ones(x.shape[1]), minus_cells)
+    u.x.scatter_forward()
+
+    dS = ufl.Measure("dS", domain=mesh, subdomain_data=ft)
+    dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct)
+    zero = Constant(mesh, PETSc.ScalarType(0))
+
+    # "+"" should be on side with larger value, i.e. exact value should be 0
+    plus_form = form(u("+") * dS(1) + zero * dx(88))
+    local_val = assemble_scalar(plus_form)
+    global_val = mesh.comm.allreduce(local_val, op=MPI.SUM)
+    assert np.isclose(global_val, 0)
+
+    # "-" should point out of domain with smaller value, exact value is 1
+    minus_form = form(u("-") * dS(1) + zero * dx(88))
+    local_val = assemble_scalar(minus_form)
+    global_val = mesh.comm.allreduce(local_val, op=MPI.SUM)
+    assert np.isclose(global_val, 1)