--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
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
)
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 &);
+ \}
+ }
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
--- /dev/null
+
+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::
projects / dealii.git / commitdiff