NonMatching::MeshClassifier: adequate treatment of face-conforming interface

[mschreter]

When the interface (given by a zero-level-set isosurface) is located at an element face, the current implementation of the NonMatching::MeshClassifier class detects both cells adjacent to the face-conforming interface as intersected. This can be seen here (white line refers to zero-level-set isosurface; gray cells are marked as intersected):

In our application, we (1) check if a cell is cut by the interface and (2) perform an integral over the interface.
However, when dealing with face-conforming interfaces, our current method results in incorrect outcomes due to inadvertently looping over the interface twice for adjacent cells. To address this issue, we've applied a workaround where we mark one of the adjacent cells as "intersected" and the other as "inside". I am wondering how the MeshClassifier should handle this corner case in a user-friendly way and would be happy to hear your opinion about it. Does it generally make sense to let the interface belong to only one cell in the face-conforming case?

FYI @nmuch @bergbauer

[tjhei]

FYI @vyushut
What did you do to deal with this, adjust the level set?

[simonsticko]

we've applied a workaround where we mark one of the adjacent cells as "intersected" and the other as "inside".

I think this is a good solution. Even if it is kind of weird. As it is right now, using NonMatching::FEValues gives the wrong result, but this fix would work.

In your workaround, what convention did you choose for which cell is intersected?

[mschreter]

Currently, we use the following procedure:

1. We loop over the faces and evaluate the level set values at the support points of each face. If these are all zero for a face, then we have a face conforming intersection.
2. In this case, if at least one of the other faces is marked as inside (from the MeshClassifier), we decide to mark this cell as inside.

For now, we have not modified the MeshClassifier class itself. Instead, we first set the active FE index according to the default suggestion of the MeshClassifier and then modify it according to the steps above. The corresponding code section looks as follows:

    // definition of lambda function for setting either the active or the future FE index
    // for the current cell
    const auto set_fe_index_lambda = [&](const auto &cell, const auto active_fe_index) {
      if (is_future)
        cell->set_future_fe_index(active_fe_index);
      else
        cell->set_active_fe_index(active_fe_index);
    };

    // 2) modify intersected cells that contain a face-conforming interface
    MappingQ1<dim>            mapping;
    dealii::FEFaceValues<dim> ls_eval(mapping,
                                      ls_dof_handler.get_fe(),
                                      Quadrature<dim - 1>(
                                        ls_dof_handler.get_fe().get_unit_face_support_points()),
                                      update_values);


    for (const auto &cell : dof_handler.get_triangulation().active_cell_iterators())
      {
        if (cell->is_locally_owned() && mesh_classifier.location_to_level_set(cell) ==
                                          NonMatching::LocationToLevelSet::intersected)
          {
            bool cell_has_conforming_face = false;


            TriaIterator<DoFCellAccessor<dim, dim, false>> dof_cell(
              &dof_handler.get_triangulation(), cell->level(), cell->index(), &dof_handler);
            TriaIterator<DoFCellAccessor<dim, dim, false>> ls_dof_cell(
              &dof_handler.get_triangulation(), cell->level(), cell->index(), &ls_dof_handler);


            for (unsigned int face_index = 0; face_index < GeometryInfo<dim>::faces_per_cell;
                 ++face_index)
              {
                std::vector<double> ls_values(ls_dof_handler.get_fe().n_dofs_per_face());
                ls_eval.reinit(ls_dof_cell, face_index);
                ls_eval.get_function_values(level_set, ls_values);

                if (cell_has_conforming_face =
                      std::all_of(ls_values.begin(), ls_values.end(), [&tolerance](double v) {
                        return std::abs(v) <= tolerance;
                      }))
                  break;
              }

            // shortcut if no face-conforming interface has been detected
            if (not cell_has_conforming_face)
              continue;

            // If a face with a conforming interface is detected and there is at least one inside
            // face attached to the cell, mark this cell as "inside".
            bool all_faces_are_intersected = false;
            for (unsigned int face_index = 0; face_index < GeometryInfo<dim>::faces_per_cell;
                 ++face_index)
              {
                all_faces_are_intersected =
                  std::min(all_faces_are_intersected,
                           mesh_classifier.location_to_level_set(cell, face_index) ==
                             NonMatching::LocationToLevelSet::intersected);
                if (mesh_classifier.location_to_level_set(cell, face_index) ==
                    NonMatching::LocationToLevelSet::inside)
                  {
                    set_fe_index_lambda(dof_cell, ActiveFEIndex::inside);
                    break;
                  }
              }

            AssertThrow(!all_faces_are_intersected, ExcNotImplemented());
          }
      }

However, if you think this is useful, then I can see how we could put this into the MeshClassifier class itself.