mesh_classifier.cc

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2021 - 2025 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#include <deal.II/base/quadrature.h>
 
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
 
#include <deal.II/fe/fe_bernstein.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_q_iso_q1.h>
#include <deal.II/fe/fe_values.h>
 
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/petsc_block_vector.h>
#include <deal.II/lac/petsc_vector.h>
#include <deal.II/lac/trilinos_epetra_vector.h>
#include <deal.II/lac/trilinos_parallel_block_vector.h>
#include <deal.II/lac/trilinos_tpetra_block_vector.h>
#include <deal.II/lac/trilinos_tpetra_vector.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/vector_element_access.h>
 
#include <deal.II/non_matching/mesh_classifier.h>
 
#include <algorithm>
 
DEAL_II_NAMESPACE_OPEN
 
namespace NonMatching
{
  namespace internal
  {
    namespace MeshClassifierImplementation
    {
      DeclExceptionMsg(
        ExcReclassifyNotCalled,
        "The Triangulation has not been classified. You need to call the "
        "reclassify()-function before using this function.");
 
      DeclExceptionMsg(
        ExcTriangulationMismatch,
        "The incoming cell does not belong to the triangulation passed to "
        "the constructor.");
 
      /*
       * Return LocationToLevelSet::inside/outside if all values in the vector
       * are negative/positive.
       * Return LocationToLevelSet::intersected if the values in the vector
       * include both negative and positive numbers.
       * Return LocationToLevelSet::aligned if all values in the vector are
       * exactly zero.
       */
      template <typename VectorType>
      LocationToLevelSet
      location_from_dof_signs(const VectorType &local_levelset_values)
      {
        const auto [min_element, max_element] =
          std::minmax_element(local_levelset_values.begin(),
                              local_levelset_values.end());
 
        // Note that we actually want to compare that the values are exactly
        // floating point zero here, even if the rule of thumb is that one never
        // should.
        if (*min_element == 0.0 && *max_element == 0.0)
          return LocationToLevelSet::aligned;
        if (*max_element <= 0)
          return LocationToLevelSet::inside;
        if (0 <= *min_element)
          return LocationToLevelSet::outside;
 
        return LocationToLevelSet::intersected;
      }
 
 
 
      /*
       * Type trait to identify a LinearAlgebra::distributed::Vector on the
       * host, whose ghost entries are described by a
       * Utilities::MPI::Partitioner.
       */
      template <typename VectorType>
      struct IsHostDistributedVector : std::false_type
      {};
 
      template <typename Number>
      struct IsHostDistributedVector<
        LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>>
        : std::true_type
      {};
 
 
 
      /*
       * Check that the incoming level set vector contains all locally relevant
       * DoFs of the incoming DoFHandler, either as locally owned or as ghost
       * entries. This is required since MeshClassifier::reclassify() reads the
       * level set values on the faces of all non-artificial cells, i.e., also
       * on the ghost cells. The check is only performed for
       * LinearAlgebra::distributed::Vector, for which the ghost entries might
       * be a subset of the locally relevant DoFs, e.g., if the vector was
       * initialized by MatrixFree.
       */
      template <int dim, typename VectorType>
      void
      assert_locally_relevant_dofs_are_accessible(
        [[maybe_unused]] const DoFHandler<dim> &dof_handler,
        [[maybe_unused]] const VectorType      &level_set)
      {
        if constexpr (IsHostDistributedVector<VectorType>::value)
          {
            IndexSet required_ghost_dofs =
              DoFTools::extract_locally_relevant_dofs(dof_handler);
            required_ghost_dofs.subtract_set(
              level_set.locally_owned_elements());
 
            Assert(
              required_ghost_dofs.is_subset_of(
                level_set.get_partitioner()->ghost_indices()),
              ExcMessage(
                "The level set vector does not contain all locally relevant "
                "DoFs of the DoFHandler as locally owned or ghost entries. "
                "However, the MeshClassifier needs to access the level set "
                "values on the faces of all ghost cells. This typically "
                "happens if the vector was initialized via "
                "MatrixFree::initialize_dof_vector(), whose ghost entries only "
                "cover the DoFs of the locally owned cells. Copy the level set "
                "field into a vector whose ghost entries are given by "
                "DoFTools::extract_locally_relevant_dofs()."));
          }
      }
 
 
 
      template <int dim, typename VectorType>
      class DiscreteLevelSetDescription : public LevelSetDescription<dim>
      {
      public:
        DiscreteLevelSetDescription(const DoFHandler<dim> &dof_handler,
                                    const VectorType      &level_set);
 
        const hp::FECollection<dim> &
        get_fe_collection() const override;
 
        unsigned int
        active_fe_index(const typename Triangulation<dim>::active_cell_iterator
                          &cell) const override;
 
        void
        get_local_level_set_values(
          const typename Triangulation<dim>::active_cell_iterator &cell,
          const unsigned int                                       face_index,
          Vector<double> &local_levelset_values) override;
 
        void
        check_consistency() const override;
 
      private:
        const ObserverPointer<const DoFHandler<dim>> dof_handler;
 
        const ObserverPointer<const VectorType> level_set;
      };
 
 
 
      template <int dim, typename VectorType>
      DiscreteLevelSetDescription<dim, VectorType>::DiscreteLevelSetDescription(
        const DoFHandler<dim> &dof_handler,
        const VectorType      &level_set)
        : dof_handler(&dof_handler)
        , level_set(&level_set)
      {}
 
 
 
      template <int dim, typename VectorType>
      const hp::FECollection<dim> &
      DiscreteLevelSetDescription<dim, VectorType>::get_fe_collection() const
      {
        return dof_handler->get_fe_collection();
      }
 
 
 
      template <int dim, typename VectorType>
      void
      DiscreteLevelSetDescription<dim, VectorType>::get_local_level_set_values(
        const typename Triangulation<dim>::active_cell_iterator &cell,
        const unsigned int                                       face_index,
        Vector<double> &local_levelset_values)
      {
        const auto cell_with_dofs = cell->as_dof_handler_iterator(*dof_handler);
 
        const unsigned int n_dofs_per_face =
          dof_handler->get_fe().n_dofs_per_face();
        std::vector<types::global_dof_index> dof_indices(n_dofs_per_face);
        cell_with_dofs->face(face_index)->get_dof_indices(dof_indices);
 
        local_levelset_values.reinit(dof_indices.size());
 
        for (unsigned int i = 0; i < dof_indices.size(); i++)
          local_levelset_values[i] =
            ::internal::ElementAccess<VectorType>::get(*level_set,
                                                             dof_indices[i]);
      }
 
 
 
      template <int dim, typename VectorType>
      void
      DiscreteLevelSetDescription<dim, VectorType>::check_consistency() const
      {
        if constexpr (running_in_debug_mode())
          assert_locally_relevant_dofs_are_accessible(*dof_handler, *level_set);
      }
 
 
      template <int dim, typename VectorType>
      unsigned int
      DiscreteLevelSetDescription<dim, VectorType>::active_fe_index(
        const typename Triangulation<dim>::active_cell_iterator &cell) const
      {
        const auto cell_with_dofs = cell->as_dof_handler_iterator(*dof_handler);
 
        return cell_with_dofs->active_fe_index();
      }
 
 
      template <int dim>
      class AnalyticLevelSetDescription : public LevelSetDescription<dim>
      {
      public:
        AnalyticLevelSetDescription(const Function<dim>      &level_set,
                                    const FiniteElement<dim> &element);
 
        const hp::FECollection<dim> &
        get_fe_collection() const override;
 
        unsigned int
        active_fe_index(const typename Triangulation<dim>::active_cell_iterator
                          &cell) const override;
 
        void
        get_local_level_set_values(
          const typename Triangulation<dim>::active_cell_iterator &cell,
          const unsigned int                                       face_index,
          Vector<double> &local_levelset_values) override;
 
      private:
        const ObserverPointer<const Function<dim>> level_set;
 
        const hp::FECollection<dim> fe_collection;
 
        FEFaceValues<dim> fe_face_values;
      };
 
 
 
      template <int dim>
      AnalyticLevelSetDescription<dim>::AnalyticLevelSetDescription(
        const Function<dim>      &level_set,
        const FiniteElement<dim> &element)
        : level_set(&level_set)
        , fe_collection(element)
        , fe_face_values(element,
                         Quadrature<dim - 1>(
                           element.get_unit_face_support_points()),
                         update_quadrature_points)
      {}
 
 
 
      template <int dim>
      void
      AnalyticLevelSetDescription<dim>::get_local_level_set_values(
        const typename Triangulation<dim>::active_cell_iterator &cell,
        const unsigned int                                       face_index,
        Vector<double> &local_levelset_values)
      {
        AssertDimension(local_levelset_values.size(),
                        fe_face_values.n_quadrature_points);
 
        fe_face_values.reinit(cell, face_index);
        const std::vector<Point<dim>> &points =
          fe_face_values.get_quadrature_points();
 
        for (unsigned int i = 0; i < points.size(); i++)
          local_levelset_values[i] = level_set->value(points[i]);
      }
 
 
 
      template <int dim>
      const hp::FECollection<dim> &
      AnalyticLevelSetDescription<dim>::get_fe_collection() const
      {
        return fe_collection;
      }
 
 
 
      template <int dim>
      unsigned int
      AnalyticLevelSetDescription<dim>::active_fe_index(
        const typename Triangulation<dim>::active_cell_iterator &) const
      {
        return 0;
      }
    } // namespace MeshClassifierImplementation
  }   // namespace internal
 
 
 
  template <int dim>
  template <typename VectorType>
  MeshClassifier<dim>::MeshClassifier(const DoFHandler<dim> &dof_handler,
                                      const VectorType      &level_set)
    : triangulation(&dof_handler.get_triangulation())
    , level_set_description(
        std::make_unique<internal::MeshClassifierImplementation::
                           DiscreteLevelSetDescription<dim, VectorType>>(
          dof_handler,
          level_set))
  {
#ifdef DEAL_II_WITH_LAPACK
    const hp::FECollection<dim> &fe_collection =
      dof_handler.get_fe_collection();
    for (unsigned int i = 0; i < fe_collection.size(); i++)
      {
        // The level set function must be scalar.
        AssertDimension(fe_collection[i].n_components(), 1);
 
        Assert(fe_collection[i].has_face_support_points(),
               ExcMessage(
                 "The elements in the FECollection of the incoming DoFHandler "
                 "must have face support points."));
      }
#else
    AssertThrow(false, ExcNeedsLAPACK());
#endif
  }
 
 
 
  template <int dim>
  MeshClassifier<dim>::MeshClassifier(const Triangulation<dim> &triangulation,
                                      const Function<dim>      &level_set,
                                      const FiniteElement<dim> &element)
    : triangulation(&triangulation)
    , level_set_description(
        std::make_unique<internal::MeshClassifierImplementation::
                           AnalyticLevelSetDescription<dim>>(level_set,
                                                             element))
  {
    // The level set function must be scalar.
    AssertDimension(level_set.n_components, 1);
    AssertDimension(element.n_components(), 1);
  }
 
 
 
  template <int dim>
  void
  MeshClassifier<dim>::reclassify()
  {
    level_set_description->check_consistency();
    initialize();
    cell_locations.assign(triangulation->n_active_cells(),
                          LocationToLevelSet::unassigned);
    face_locations.assign(triangulation->n_raw_faces(),
                          LocationToLevelSet::unassigned);
 
    // Returns weather the incoming LocationToLevelSet is in the incoming set.
    // This lambda can be factored away once C++20 is enabled, since set then
    // has a contains function.
    const auto contains =
      [](const std::set<LocationToLevelSet> &local_face_locations,
         const LocationToLevelSet           &location) {
        return local_face_locations.count(location) > 0;
      };
 
    // Loop over all cells and determine the location of all non artificial
    // cells and faces.
    for (const auto &cell : triangulation->active_cell_iterators())
      if (!cell->is_artificial())
        {
          std::set<LocationToLevelSet> local_face_locations;
 
          for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
            {
              const LocationToLevelSet face_location =
                determine_face_location_to_levelset(cell, f);
 
              face_locations[cell->face(f)->index()] = face_location;
              local_face_locations.insert(face_location);
            }
 
          LocationToLevelSet cell_location = LocationToLevelSet::unassigned;
 
          const bool all_faces_have_same_location =
            local_face_locations.size() == 1;
 
          if (all_faces_have_same_location)
            {
              cell_location = *local_face_locations.cbegin();
            }
          else if (contains(local_face_locations,
                            LocationToLevelSet::intersected))
            {
              cell_location = LocationToLevelSet::intersected;
            }
          else if (contains(local_face_locations, LocationToLevelSet::aligned))
            {
              Assert(local_face_locations.size() == 2, ExcNotImplemented());
 
              if (contains(local_face_locations, LocationToLevelSet::outside))
                cell_location = LocationToLevelSet::outside;
              else // contains LocationToLevelSet::inside
                cell_location = LocationToLevelSet::intersected;
            }
          else
            {
              // In 2D, a cell can be intersected without any faces being
              // intersected, if the zero contour cuts the cell diagonally
              // through the vertices. In this case, two faces are inside and
              // two are outside.
              Assert(local_face_locations.size() == 2, ExcNotImplemented());
              cell_location = LocationToLevelSet::intersected;
            }
 
 
          cell_locations[cell->active_cell_index()] = cell_location;
        }
  }
 
 
 
  template <int dim>
  LocationToLevelSet
  MeshClassifier<dim>::determine_face_location_to_levelset(
    const typename Triangulation<dim>::active_cell_iterator &cell,
    const unsigned int                                       face_index)
  {
    // The location of the face might already be computed on the neighboring
    // cell. If this is the case we just return the value.
    const LocationToLevelSet location =
      face_locations.at(cell->face(face_index)->index());
    if (location != LocationToLevelSet::unassigned)
      return location;
 
    // Determine the location by changing basis to FE_Bernstein and checking
    // the signs of the dofs.
    const unsigned int fe_index = level_set_description->active_fe_index(cell);
    const unsigned int n_local_dofs =
      lagrange_to_bernstein_face[fe_index][face_index].m();
 
    Vector<double> local_levelset_values(n_local_dofs);
    level_set_description->get_local_level_set_values(cell,
                                                      face_index,
                                                      local_levelset_values);
 
    const FiniteElement<dim> &fe =
      level_set_description->get_fe_collection()[fe_index];
 
    const FE_Q_iso_Q1<dim> *fe_q_iso_q1 =
      dynamic_cast<const FE_Q_iso_Q1<dim> *>(&fe);
 
    const FE_Poly<dim> *fe_poly = dynamic_cast<const FE_Poly<dim> *>(&fe);
 
    const bool is_linear = fe_q_iso_q1 != nullptr ||
                           (fe_poly != nullptr && fe_poly->get_degree() == 1);
 
    // shortcut for linear elements
    if (is_linear)
      {
        return internal::MeshClassifierImplementation::location_from_dof_signs(
          local_levelset_values);
      }
 
    lagrange_to_bernstein_face[fe_index][face_index].solve(
      local_levelset_values);
 
    return internal::MeshClassifierImplementation::location_from_dof_signs(
      local_levelset_values);
  }
 
 
 
  template <int dim>
  LocationToLevelSet
  MeshClassifier<dim>::location_to_level_set(
    const typename Triangulation<dim>::cell_iterator &cell) const
  {
    Assert(cell_locations.size() == triangulation->n_active_cells(),
           internal::MeshClassifierImplementation::ExcReclassifyNotCalled());
    Assert(&cell->get_triangulation() == triangulation,
           internal::MeshClassifierImplementation::ExcTriangulationMismatch());
 
    return cell_locations.at(cell->active_cell_index());
  }
 
 
 
  template <int dim>
  LocationToLevelSet
  MeshClassifier<dim>::location_to_level_set(
    const typename Triangulation<dim>::cell_iterator &cell,
    const unsigned int                                face_index) const
  {
    AssertIndexRange(face_index, GeometryInfo<dim>::faces_per_cell);
    Assert(face_locations.size() == triangulation->n_raw_faces(),
           internal::MeshClassifierImplementation::ExcReclassifyNotCalled());
    Assert(&cell->get_triangulation() == triangulation,
           internal::MeshClassifierImplementation::ExcTriangulationMismatch());
 
    return face_locations.at(cell->face(face_index)->index());
  }
 
 
 
  template <int dim>
  void
  MeshClassifier<dim>::initialize()
  {
    const hp::FECollection<dim> &fe_collection =
      level_set_description->get_fe_collection();
 
    // The level set function must be scalar.
    AssertDimension(fe_collection.n_components(), 1);
 
    lagrange_to_bernstein_face.resize(fe_collection.size());
 
    for (unsigned int i = 0; i < fe_collection.size(); i++)
      {
        const FiniteElement<dim> &element = fe_collection[i];
        const FE_Q_Base<dim>     *fe_q =
          dynamic_cast<const FE_Q_Base<dim> *>(&element);
        Assert(fe_q != nullptr, ExcNotImplemented());
 
        const FE_Bernstein<dim> fe_bernstein(fe_q->get_degree());
 
        const unsigned int dofs_per_face = fe_q->dofs_per_face;
        for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; f++)
          {
            FullMatrix<double> face_interpolation_matrix(dofs_per_face,
                                                         dofs_per_face);
 
            fe_bernstein.get_face_interpolation_matrix(
              *fe_q, face_interpolation_matrix, f);
            lagrange_to_bernstein_face[i][f].reinit(dofs_per_face);
            lagrange_to_bernstein_face[i][f] = face_interpolation_matrix;
            lagrange_to_bernstein_face[i][f].compute_lu_factorization();
          }
      }
  }
 
} // namespace NonMatching
 
#include "non_matching/mesh_classifier.inst"
 
DEAL_II_NAMESPACE_CLOSE