]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add class determining how cells/faces relate to a level set function 12514/head
authorSimon Sticko <simon@sticko.se>
Mon, 28 Jun 2021 09:43:29 +0000 (11:43 +0200)
committerSimon Sticko <simon@sticko.se>
Wed, 18 Aug 2021 07:01:17 +0000 (09:01 +0200)
Add a class NonMatching::MeshClassifier that determines how the sign of
a level set function vary over the active cells/faces of a
triangulation. Each cell/face are associated with one of values of enum
LocationToLevelSet. If the level set is positive/negative over the
cell/face it is associated with LocationToLevelSet::inside/outside. If
the zero contour intersects the cell/face it is associated with
LocationToLevelSet::intersected.

include/deal.II/non_matching/mesh_classifier.h [new file with mode: 0644]
source/non_matching/CMakeLists.txt
source/non_matching/mesh_classifier.cc [new file with mode: 0644]
source/non_matching/mesh_classifier.inst.in [new file with mode: 0644]
tests/non_matching/mesh_classifier.cc [new file with mode: 0644]
tests/non_matching/mesh_classifier.output [new file with mode: 0644]

diff --git a/include/deal.II/non_matching/mesh_classifier.h b/include/deal.II/non_matching/mesh_classifier.h
new file mode 100644 (file)
index 0000000..e6bae28
--- /dev/null
@@ -0,0 +1,267 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_non_matching_mesh_classifier
+#define dealii_non_matching_mesh_classifier
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+
+#include <deal.II/lac/lapack_full_matrix.h>
+
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+namespace NonMatching
+{
+  namespace internal
+  {
+    namespace MeshClassifierImplementation
+    {
+      template <int dim>
+      class LevelSetDescription;
+    } // namespace MeshClassifierImplementation
+  }   // namespace internal
+
+
+  /**
+   * Type describing how a cell or a face is located relative to the zero
+   * contour of a level set function, $\psi$. The values of the type correspond
+   * to:
+   *
+   * inside      if $\psi(x) < 0$,
+   * outside     if $\psi(x) > 0$,
+   * intersected if $\psi(x)$ varies in sign,
+   *
+   * over the cell/face. The value "unassigned" is used to describe that the
+   * location of a cell/face has not yet been determined.
+   */
+  enum class LocationToLevelSet
+  {
+    inside,
+    outside,
+    intersected,
+    unassigned
+  };
+
+
+  /**
+   * Class responsible for determining how the active cells and faces of a
+   * triangulation relate to the sign of a level set function. When calling the
+   * reclassify() function each of the active cells and faces are categorized as
+   * one of the values of LocationToLevelSet: inside, outside or intersected,
+   * depending on the sign of the level set function over the cell/face. This
+   * information is typically required in immersed/cut finite element methods,
+   * both when distributing degrees of freedom over the triangulation and when
+   * the system is assembled. The given class would then be used in the
+   * following way:
+   *
+   * @code
+   * Vector<double> &level_set = ...
+   *
+   * MeshClassifier<dim> classifier(dof_handler, level_set);
+   * classifier.reclassify();
+   *
+   * LocationToLevelSet location = classifier.location_to_level_set(cell);
+   * @endcode
+   *
+   * The level set function can either be described as a discrete function by a
+   * (DoFHandler, Vector)-pair or as a general Function. In the case of a
+   * discrete function, LocationToLevelSet for a given face is determined
+   * by looking at the local degrees of freedom on the face. Since the Lagrange
+   * basis functions are not positive definite, positive/negative definite
+   * dof values do not imply that the interpolated function is
+   * positive/negative definite. Thus, to classify a face this class internally
+   * transforms the local dofs to a basis spanning the same polynomial space
+   * over the face but where definite dof values imply a definite function.
+   * Currently, only the case of FE_Q-elements is implemented, where we
+   * internally change basis to FE_Bernstein. For cells, LocationToLevelSet is
+   * determined from the faces of the cell. That is, if all faces of the cell
+   * are inside/outside the LocationToLevelSet of the cell is set to
+   * inside/outside. LocationToLevelSet of the cell is set to intersected if at
+   * least one face is intersected or if its faces have different
+   * LocationToLevelSet. Note that, this procedure will incorrectly classify the
+   * cell as inside/outside, if the mesh refinement is so low that the whole
+   * zero-contour is contained in a single cell (so that none of its faces are
+   * intersected).
+   *
+   * When the level set function is described as a Function, the level set
+   * function is locally interpolated to an FE_Q element and we proceed in the
+   * same way as for the discrete level set function.
+   */
+  template <int dim>
+  class MeshClassifier : public Subscriptor
+  {
+  public:
+    /**
+     * Constructor. Takes a level set function described as a DoFHandler and a
+     * Vector. The triangulation attached to DoFHandler is the one that will be
+     * classified.
+     */
+    template <class VECTOR>
+    MeshClassifier(const DoFHandler<dim> &level_set_dof_handler,
+                   const VECTOR &         level_set);
+
+
+    /**
+     * Constructor. Takes the triangulation that should be classified, a
+     * level set function described as a Function, and a scalar element that we
+     * interpolate the Function to in order to classify each cell/face.
+     *
+     * @note The Function and the FiniteElement must both have a single component.
+     */
+    MeshClassifier(const Triangulation<dim> &triangulation,
+                   const Function<dim> &     level_set,
+                   const FiniteElement<dim> &element);
+
+    /**
+     * Perform the classification of the non artificial cells and faces in the
+     * triangulation.
+     */
+    void
+    reclassify();
+
+    /**
+     * Return how the incoming cell is located relative to the level set
+     * function.
+     */
+    LocationToLevelSet
+    location_to_level_set(
+      const typename Triangulation<dim>::cell_iterator &cell) const;
+
+    /**
+     * Return how a face of the incoming cell is located relative to the level
+     * set function.
+     */
+    LocationToLevelSet
+    location_to_level_set(
+      const typename Triangulation<dim>::cell_iterator &cell,
+      const unsigned int                                face_index) const;
+
+  private:
+    /**
+     * For each element in the hp::FECollection returned by
+     * level_set_description, sets up the local transformation matrices.
+     */
+    void
+    initialize();
+
+    /**
+     * Computes how the face with the given index on the incoming cell is
+     * located relative to the level set function.
+     */
+    LocationToLevelSet
+    determine_face_location_to_levelset(
+      const typename Triangulation<dim>::active_cell_iterator &cell,
+      const unsigned int                                       face_index);
+
+    /**
+     * Pointer to the triangulation that should be classified.
+     */
+    const SmartPointer<const Triangulation<dim>> triangulation;
+
+    /**
+     * Pointer to an object that describes what we need to know about the
+     * level set function. The underlying object will be of different type
+     * depending on whether the level set function is discrete
+     * (DoFHandler, Vector) or described by a Function.
+     */
+    const std::unique_ptr<
+      internal::MeshClassifierImplementation::LevelSetDescription<dim>>
+      level_set_description;
+
+    /**
+     * A vector that stores how each active cell is located relative to the
+     * level set function, based on the cells active index.
+     */
+    std::vector<LocationToLevelSet> cell_locations;
+
+    /**
+     * A vector that stores how each active face is located relative to the
+     * level set function, based on the face's global index.
+     */
+    std::vector<LocationToLevelSet> face_locations;
+
+    /**
+     * For each element in the hp::FECollection returned by the
+     * LevelSetDescription, and for each local face, this vector stores a
+     * transformation matrix to a basis where positive/negative definite
+     * face dofs implies that the underlying function is positive/negative
+     * definite over the face.
+     */
+    std::vector<
+      std::array<LAPACKFullMatrix<double>, GeometryInfo<dim>::faces_per_cell>>
+      lagrange_to_bernstein_face;
+  };
+
+
+  namespace internal
+  {
+    namespace MeshClassifierImplementation
+    {
+      /**
+       * Abstract class that describes what we need to know about the level set
+       * function independently of whether it is a Function or a
+       * (DoFHandler, Vector)-pair.
+       */
+      template <int dim>
+      class LevelSetDescription
+      {
+      public:
+        /**
+         * Destructor, declared to mark it virtual.
+         */
+        virtual ~LevelSetDescription() = default;
+
+        /**
+         * Return a collection to all the elements that are used to locally
+         * describe the level set function.
+         */
+        virtual const hp::FECollection<dim> &
+        get_fe_collection() const = 0;
+
+        /**
+         * Return the index of the element in the FECollection that we associate
+         * with the level set function on the incoming cell.
+         */
+        virtual unsigned int
+        active_fe_index(const typename Triangulation<dim>::active_cell_iterator
+                          &cell) const = 0;
+
+        /**
+         * Fill the DoF values of the associated level set representation on the
+         * face of the incoming cell into the vector provided in the last
+         * argument.
+         *
+         * @note Since this function extracts the dofs on the face of the cell,
+         * it assumes that the underlying element has face support points.
+         */
+        virtual void
+        get_local_level_set_values(
+          const typename Triangulation<dim>::active_cell_iterator &cell,
+          const unsigned int                                       face_index,
+          Vector<double> &local_dofs) = 0;
+      };
+
+    } // namespace MeshClassifierImplementation
+  }   // namespace internal
+} // namespace NonMatching
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 2d1b9a8a76f79c553ce39f552d45c2d498ceb1ce..aec821f02a52a20ad423cf9b2028727e646a26bf 100644 (file)
 INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
 
 SET(_src
+  mesh_classifier.cc
   quadrature_generator.cc
   coupling.cc
   immersed_surface_quadrature.cc
   )
 
 SET(_inst
+  mesh_classifier.inst.in
   quadrature_generator.inst.in
   coupling.inst.in
   )
diff --git a/source/non_matching/mesh_classifier.cc b/source/non_matching/mesh_classifier.cc
new file mode 100644 (file)
index 0000000..4b8b994
--- /dev/null
@@ -0,0 +1,479 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/quadrature.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+
+#include <deal.II/fe/fe_bernstein.h>
+#include <deal.II/fe/fe_q.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/la_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_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
+    {
+      /**
+       * Return LocationToLevelSet::inside/outside if all values in incoming
+       * vector are negative/positive, otherwise return
+       * LocationToLevelSet::intersected.
+       */
+      template <class VECTOR>
+      LocationToLevelSet
+      location_from_dof_signs(const VECTOR &local_levelset_values)
+      {
+        const auto min_max_element =
+          std::minmax_element(local_levelset_values.begin(),
+                              local_levelset_values.end());
+
+        if (*min_max_element.second < 0)
+          return LocationToLevelSet::inside;
+        if (0 < *min_max_element.first)
+          return LocationToLevelSet::outside;
+
+        return LocationToLevelSet::intersected;
+      }
+
+
+
+      /**
+       * The concrete LevelSetDescription used when the level set function is
+       * described as a (DoFHandler, Vector)-pair.
+       */
+      template <int dim, class VECTOR>
+      class DiscreteLevelSetDescription : public LevelSetDescription<dim>
+      {
+      public:
+        /**
+         * Constructor.
+         */
+        DiscreteLevelSetDescription(const DoFHandler<dim> &dof_handler,
+                                    const VECTOR &         level_set);
+
+        /**
+         * Return the FECollection of the DoFHandler passed to the constructor.
+         */
+        const hp::FECollection<dim> &
+        get_fe_collection() const override;
+
+        /**
+         * Return the active_fe_index of the DoFCellAccessor associated with the
+         * DoFHandler and the the incoming cell in the triangulation.
+         */
+        unsigned int
+        active_fe_index(const typename Triangulation<dim>::active_cell_iterator
+                          &cell) const override;
+
+        /**
+         * Writes the local face dofs of the global level set vector to
+         * @p local_levelset_values.
+         */
+        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:
+        /**
+         * Pointer to the DoFHandler associated with the level set function.
+         */
+        const SmartPointer<const DoFHandler<dim>> dof_handler;
+
+        /**
+         * Pointer to the vector containing the level set function's global dof
+         * values.
+         */
+        const SmartPointer<const VECTOR> level_set;
+      };
+
+
+
+      template <int dim, class VECTOR>
+      DiscreteLevelSetDescription<dim, VECTOR>::DiscreteLevelSetDescription(
+        const DoFHandler<dim> &dof_handler,
+        const VECTOR &         level_set)
+        : dof_handler(&dof_handler)
+        , level_set(&level_set)
+      {}
+
+
+
+      template <int dim, class VECTOR>
+      const hp::FECollection<dim> &
+      DiscreteLevelSetDescription<dim, VECTOR>::get_fe_collection() const
+      {
+        return dof_handler->get_fe_collection();
+      }
+
+
+
+      template <int dim, class VECTOR>
+      void
+      DiscreteLevelSetDescription<dim, VECTOR>::get_local_level_set_values(
+        const typename Triangulation<dim>::active_cell_iterator &cell,
+        const unsigned int                                       face_index,
+        Vector<double> &local_levelset_values)
+      {
+        typename DoFHandler<dim>::active_cell_iterator cell_with_dofs(
+          &dof_handler->get_triangulation(),
+          cell->level(),
+          cell->index(),
+          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] =
+            dealii::internal::ElementAccess<VECTOR>::get(*level_set,
+                                                         dof_indices[i]);
+      }
+
+
+
+      template <int dim, class VECTOR>
+      unsigned int
+      DiscreteLevelSetDescription<dim, VECTOR>::active_fe_index(
+        const typename Triangulation<dim>::active_cell_iterator &cell) const
+      {
+        typename DoFHandler<dim>::active_cell_iterator cell_with_dofs(
+          &dof_handler->get_triangulation(),
+          cell->level(),
+          cell->index(),
+          dof_handler);
+
+        return cell_with_dofs->active_fe_index();
+      }
+
+
+      /**
+       * The concrete LevelSetDescription used when the level set function is
+       * described by a Function.
+       */
+      template <int dim>
+      class AnalyticLevelSetDescription : public LevelSetDescription<dim>
+      {
+      public:
+        /**
+         * Constructor. Takes the Function that describes the geometry and the
+         * element that this function should be interpolated to.
+         */
+        AnalyticLevelSetDescription(const Function<dim> &     level_set,
+                                    const FiniteElement<dim> &element);
+
+        /**
+         * Returns the finite element passed to the constructor wrapped in a
+         * collection.
+         */
+        const hp::FECollection<dim> &
+        get_fe_collection() const override;
+
+        /**
+         * Returns 0, since there is always a single element in the
+         * FECollection.
+         */
+        unsigned int
+        active_fe_index(const typename Triangulation<dim>::active_cell_iterator
+                          &cell) const override;
+
+        /**
+         * Return the level set function evaluated at the real space face
+         * support points of the finite element passed to the constructor.
+         */
+        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:
+        /**
+         * Pointer to the level set function.
+         */
+        const SmartPointer<const Function<dim>> level_set;
+
+        /**
+         * Collection containing the single element which we locally interpolate
+         * the level set function to.
+         */
+        const hp::FECollection<dim> fe_collection;
+
+        /**
+         * FEFaceValues object used to transform the support points on a face to
+         * real space.
+         */
+        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 <class VECTOR>
+  MeshClassifier<dim>::MeshClassifier(const DoFHandler<dim> &dof_handler,
+                                      const VECTOR &         level_set)
+    : triangulation(&dof_handler.get_triangulation())
+    , level_set_description(
+        std::make_unique<internal::MeshClassifierImplementation::
+                           DiscreteLevelSetDescription<dim, VECTOR>>(
+          dof_handler,
+          level_set))
+  {
+    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."));
+      }
+  }
+
+
+
+  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()
+  {
+    initialize();
+    cell_locations.resize(triangulation->n_active_cells(),
+                          LocationToLevelSet::unassigned);
+    face_locations.resize(triangulation->n_raw_faces(),
+                          LocationToLevelSet::unassigned);
+
+    // 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())
+        {
+          const LocationToLevelSet face0_location =
+            determine_face_location_to_levelset(cell, 0);
+
+          face_locations[cell->face(0)->index()] = face0_location;
+          LocationToLevelSet cell_location       = face0_location;
+
+          for (unsigned int f = 1; 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;
+
+              if (face_location != face0_location)
+                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);
+
+    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
+  {
+    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);
+
+    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<dim> *fe_q = dynamic_cast<const FE_Q<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 "mesh_classifier.inst"
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/non_matching/mesh_classifier.inst.in b/source/non_matching/mesh_classifier.inst.in
new file mode 100644 (file)
index 0000000..1505a0f
--- /dev/null
@@ -0,0 +1,32 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2010 - 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+for (deal_II_dimension : DIMENSIONS)
+  {
+    namespace NonMatching
+    \{
+      template class MeshClassifier<deal_II_dimension>;
+    \}
+  }
+
+for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS)
+  {
+    namespace NonMatching
+    \{
+      template MeshClassifier<deal_II_dimension>::MeshClassifier(
+        const DoFHandler<deal_II_dimension> &,
+        const VEC &);
+    \}
+  }
diff --git a/tests/non_matching/mesh_classifier.cc b/tests/non_matching/mesh_classifier.cc
new file mode 100644 (file)
index 0000000..09d2b7d
--- /dev/null
@@ -0,0 +1,341 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Test the MeshClassifier class by classifying a triangulation with a single
+// cell with a few different level set functions.
+
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/function_level_set.h>
+#include <deal.II/base/point.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/non_matching/mesh_classifier.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+// Return a string with the same name as the incoming value of
+// LocationToLevelSet.
+std::string
+location_to_string(const NonMatching::LocationToLevelSet location)
+{
+  std::string name;
+  switch (location)
+    {
+      case NonMatching::LocationToLevelSet::inside:
+        name = "inside";
+        break;
+      case NonMatching::LocationToLevelSet::outside:
+        name = "outside";
+        break;
+      case NonMatching::LocationToLevelSet::intersected:
+        name = "intersected";
+        break;
+      case NonMatching::LocationToLevelSet::unassigned:
+        name = "unassigned";
+        break;
+      default:
+        AssertThrow(false, ExcInternalError());
+    }
+  return name;
+}
+
+
+
+// Print LocationToLevelSet (as determined by the incoming MeshClassifier) to
+// deallog, for the incoming cell and all of its faces.
+template <int dim>
+void
+print_cell_and_face_locations(
+  const NonMatching::MeshClassifier<dim> &                 classifier,
+  const typename Triangulation<dim>::active_cell_iterator &cell)
+{
+  const NonMatching::LocationToLevelSet cell_position =
+    classifier.location_to_level_set(cell);
+  deallog << "cell " << location_to_string(cell_position) << std::endl;
+
+  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+    {
+      const NonMatching::LocationToLevelSet cell_position =
+        classifier.location_to_level_set(cell, i);
+      deallog << "face " << i << " " << location_to_string(cell_position)
+              << std::endl;
+    }
+}
+
+
+
+// Test the version MeshClassifier that takes a Vector and a DoFHandler.
+//
+// Set up single cell triangulation over [-1, 1]^dim and a DoFHandler using
+// FE_Q<1>(1). Interpolate the incoming level set function to the discrete
+// space. Classify the cell and its faces and print the result to deallog.
+template <int dim>
+void
+classify_with_discrete_level_set(const Function<dim> &level_set)
+{
+  deallog << "discrete:" << std::endl;
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube(triangulation, -1, 1);
+
+  const FE_Q<dim> element(1);
+  DoFHandler<dim> dof_handler(triangulation);
+  dof_handler.distribute_dofs(element);
+
+  Vector<double> discrete_level_set;
+  discrete_level_set.reinit(dof_handler.n_dofs());
+  VectorTools::interpolate(dof_handler, level_set, discrete_level_set);
+
+  NonMatching::MeshClassifier<dim> classifier(dof_handler, discrete_level_set);
+  classifier.reclassify();
+
+  const typename Triangulation<dim>::active_cell_iterator cell =
+    triangulation.begin_active();
+  print_cell_and_face_locations(classifier, cell);
+  deallog << std::endl;
+}
+
+
+
+// Test the version of MeshClassifier that takes a Function.
+//
+// Set up single cell triangulation over [-1, 1]^dim. Classify the cell and its
+// faces and print the result to deallog.
+template <int dim>
+void
+classify_with_analytic_level_set(const Function<dim> &level_set)
+{
+  deallog << "analytic:" << std::endl;
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube(triangulation, -1, 1);
+
+  const FE_Q<dim> element(1);
+
+  NonMatching::MeshClassifier<dim> classifier(triangulation,
+                                              level_set,
+                                              element);
+  classifier.reclassify();
+
+  const typename Triangulation<dim>::active_cell_iterator cell =
+    triangulation.begin_active();
+  print_cell_and_face_locations(classifier, cell);
+  deallog << std::endl;
+}
+
+
+
+// Test both the version of MeshClassifier that takes a Vector and a DoFHandler,
+// and the version that takes a Function.
+template <int dim>
+void
+classify_with_discrete_and_analytic_level_set(const Function<dim> &level_set)
+{
+  deallog << std::endl;
+
+  classify_with_discrete_level_set(level_set);
+
+  classify_with_analytic_level_set(level_set);
+}
+
+
+
+// Test MeshClassifier with a level set function that is constant and positive.
+template <int dim>
+void
+test_positive_function()
+{
+  deallog << "test_positive_function" << std::endl;
+
+  const Functions::ConstantFunction<dim> level_set(1);
+
+  classify_with_discrete_and_analytic_level_set(level_set);
+}
+
+
+
+// Test MeshClassifier with a level set function that is constant and negative.
+template <int dim>
+void
+test_negative_function()
+{
+  deallog << "test_negative_function" << std::endl;
+
+  const Functions::ConstantFunction<dim> level_set(-1);
+
+  classify_with_discrete_and_analytic_level_set(level_set);
+}
+
+
+
+// Test MeshClassifier with a level set function corresponding to the plane
+// (x = 0) intersecting the hypercube [-1, 1]^dim.
+template <int dim>
+void
+test_intersection_x_eq_0_plane()
+{
+  deallog << "test_intersection_x_eq_0_plane" << std::endl;
+
+  Tensor<1, dim> plane_normal;
+  plane_normal[0] = 1;
+  const Point<dim> origo;
+
+  const Functions::LevelSet::Plane<dim> level_set(origo, plane_normal);
+
+  classify_with_discrete_and_analytic_level_set(level_set);
+}
+
+
+
+// Set up local level set coefficients for an Q2 element such that all
+// coefficients are positive but the cell is still intersected.
+template <int dim>
+void
+setup_intersected_Q2_positive_coefficients(Vector<double> &level_set){
+  // This test case only makes sense in 2D and 3D, specialize for these below
+  // and do nothing by default.
+};
+
+
+
+template <>
+void
+setup_intersected_Q2_positive_coefficients<2>(Vector<double> &level_set)
+{
+  const double delta = 1e-5;
+
+  // Listing entries lexiographically.
+  level_set(0) = 100 * delta;
+  level_set(6) = delta;
+  level_set(1) = delta;
+  level_set(4) = .5;
+  level_set(8) = .5;
+  level_set(5) = .5;
+  level_set(2) = 1;
+  level_set(7) = 1;
+  level_set(3) = 1;
+}
+
+
+
+template <>
+void
+setup_intersected_Q2_positive_coefficients<3>(Vector<double> &level_set)
+{
+  const double delta = 1e-5;
+
+  level_set(0)  = 100 * delta;
+  level_set(10) = 100 * delta;
+  level_set(1)  = 100 * delta;
+
+  level_set(8)  = 100 * delta;
+  level_set(24) = delta;
+  level_set(9)  = delta;
+
+  level_set(2)  = 100 * delta;
+  level_set(11) = 100 * delta;
+  level_set(3)  = 100 * delta;
+
+  level_set(16) = .5;
+  level_set(22) = .5;
+  level_set(17) = .5;
+
+  level_set(20) = .5;
+  level_set(26) = .5;
+  level_set(21) = .5;
+
+  level_set(18) = .5;
+  level_set(23) = .5;
+  level_set(19) = .5;
+
+  level_set(4)  = 1;
+  level_set(14) = 1;
+  level_set(5)  = 1;
+
+  level_set(12) = 1;
+  level_set(25) = 1;
+  level_set(13) = 1;
+
+  level_set(6)  = 1;
+  level_set(15) = 1;
+  level_set(7)  = 1;
+}
+
+
+
+// Test the case that all the Lagrange coefficients of a Q2 element are positive
+// but the cell is still intersected. This can happen because the Lagrange shape
+// functions are negative between the support points.
+template <int dim>
+void
+test_lagrange_coefficents_positive()
+{
+  deallog << "test_lagrange_coefficents_positive" << std::endl;
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube(triangulation);
+
+  const FE_Q<dim> element(2);
+
+  DoFHandler<dim> dof_handler(triangulation);
+  dof_handler.distribute_dofs(element);
+
+  Vector<double> level_set(element.dofs_per_cell);
+  setup_intersected_Q2_positive_coefficients<dim>(level_set);
+
+  NonMatching::MeshClassifier<dim> classifier(dof_handler, level_set);
+  classifier.reclassify();
+
+  const typename Triangulation<dim>::active_cell_iterator cell =
+    triangulation.begin_active();
+  print_cell_and_face_locations(classifier, cell);
+  deallog << std::endl;
+}
+
+
+
+template <int dim>
+void
+run_test()
+{
+  deallog << "dim = " << dim << std::endl;
+
+  test_negative_function<dim>();
+  test_positive_function<dim>();
+  test_intersection_x_eq_0_plane<dim>();
+  // This test doesn't make sense in 1D.
+  if (dim != 1)
+    test_lagrange_coefficents_positive<dim>();
+}
+
+
+
+int
+main()
+{
+  initlog();
+  run_test<1>();
+  run_test<2>();
+  run_test<3>();
+}
diff --git a/tests/non_matching/mesh_classifier.output b/tests/non_matching/mesh_classifier.output
new file mode 100644 (file)
index 0000000..a29c6fa
--- /dev/null
@@ -0,0 +1,164 @@
+
+DEAL::dim = 1
+DEAL::test_negative_function
+DEAL::
+DEAL::discrete:
+DEAL::cell inside
+DEAL::face 0 inside
+DEAL::face 1 inside
+DEAL::
+DEAL::analytic:
+DEAL::cell inside
+DEAL::face 0 inside
+DEAL::face 1 inside
+DEAL::
+DEAL::test_positive_function
+DEAL::
+DEAL::discrete:
+DEAL::cell outside
+DEAL::face 0 outside
+DEAL::face 1 outside
+DEAL::
+DEAL::analytic:
+DEAL::cell outside
+DEAL::face 0 outside
+DEAL::face 1 outside
+DEAL::
+DEAL::test_intersection_x_eq_0_plane
+DEAL::
+DEAL::discrete:
+DEAL::cell intersected
+DEAL::face 0 inside
+DEAL::face 1 outside
+DEAL::
+DEAL::analytic:
+DEAL::cell intersected
+DEAL::face 0 inside
+DEAL::face 1 outside
+DEAL::
+DEAL::dim = 2
+DEAL::test_negative_function
+DEAL::
+DEAL::discrete:
+DEAL::cell inside
+DEAL::face 0 inside
+DEAL::face 1 inside
+DEAL::face 2 inside
+DEAL::face 3 inside
+DEAL::
+DEAL::analytic:
+DEAL::cell inside
+DEAL::face 0 inside
+DEAL::face 1 inside
+DEAL::face 2 inside
+DEAL::face 3 inside
+DEAL::
+DEAL::test_positive_function
+DEAL::
+DEAL::discrete:
+DEAL::cell outside
+DEAL::face 0 outside
+DEAL::face 1 outside
+DEAL::face 2 outside
+DEAL::face 3 outside
+DEAL::
+DEAL::analytic:
+DEAL::cell outside
+DEAL::face 0 outside
+DEAL::face 1 outside
+DEAL::face 2 outside
+DEAL::face 3 outside
+DEAL::
+DEAL::test_intersection_x_eq_0_plane
+DEAL::
+DEAL::discrete:
+DEAL::cell intersected
+DEAL::face 0 inside
+DEAL::face 1 outside
+DEAL::face 2 intersected
+DEAL::face 3 intersected
+DEAL::
+DEAL::analytic:
+DEAL::cell intersected
+DEAL::face 0 inside
+DEAL::face 1 outside
+DEAL::face 2 intersected
+DEAL::face 3 intersected
+DEAL::
+DEAL::test_lagrange_coefficents_positive
+DEAL::cell intersected
+DEAL::face 0 outside
+DEAL::face 1 outside
+DEAL::face 2 intersected
+DEAL::face 3 outside
+DEAL::
+DEAL::dim = 3
+DEAL::test_negative_function
+DEAL::
+DEAL::discrete:
+DEAL::cell inside
+DEAL::face 0 inside
+DEAL::face 1 inside
+DEAL::face 2 inside
+DEAL::face 3 inside
+DEAL::face 4 inside
+DEAL::face 5 inside
+DEAL::
+DEAL::analytic:
+DEAL::cell inside
+DEAL::face 0 inside
+DEAL::face 1 inside
+DEAL::face 2 inside
+DEAL::face 3 inside
+DEAL::face 4 inside
+DEAL::face 5 inside
+DEAL::
+DEAL::test_positive_function
+DEAL::
+DEAL::discrete:
+DEAL::cell outside
+DEAL::face 0 outside
+DEAL::face 1 outside
+DEAL::face 2 outside
+DEAL::face 3 outside
+DEAL::face 4 outside
+DEAL::face 5 outside
+DEAL::
+DEAL::analytic:
+DEAL::cell outside
+DEAL::face 0 outside
+DEAL::face 1 outside
+DEAL::face 2 outside
+DEAL::face 3 outside
+DEAL::face 4 outside
+DEAL::face 5 outside
+DEAL::
+DEAL::test_intersection_x_eq_0_plane
+DEAL::
+DEAL::discrete:
+DEAL::cell intersected
+DEAL::face 0 inside
+DEAL::face 1 outside
+DEAL::face 2 intersected
+DEAL::face 3 intersected
+DEAL::face 4 intersected
+DEAL::face 5 intersected
+DEAL::
+DEAL::analytic:
+DEAL::cell intersected
+DEAL::face 0 inside
+DEAL::face 1 outside
+DEAL::face 2 intersected
+DEAL::face 3 intersected
+DEAL::face 4 intersected
+DEAL::face 5 intersected
+DEAL::
+DEAL::test_lagrange_coefficents_positive
+DEAL::cell intersected
+DEAL::face 0 outside
+DEAL::face 1 intersected
+DEAL::face 2 outside
+DEAL::face 3 outside
+DEAL::face 4 intersected
+DEAL::face 5 outside
+DEAL::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.