From 3fa370ba6b8f11bcd473c904f99fdf223011df23 Mon Sep 17 00:00:00 2001 From: Simon Sticko Date: Mon, 28 Jun 2021 11:43:29 +0200 Subject: [PATCH] Add class determining how cells/faces relate to a level set function 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. --- .../deal.II/non_matching/mesh_classifier.h | 267 ++++++++++ source/non_matching/CMakeLists.txt | 2 + source/non_matching/mesh_classifier.cc | 479 ++++++++++++++++++ source/non_matching/mesh_classifier.inst.in | 32 ++ tests/non_matching/mesh_classifier.cc | 341 +++++++++++++ tests/non_matching/mesh_classifier.output | 164 ++++++ 6 files changed, 1285 insertions(+) create mode 100644 include/deal.II/non_matching/mesh_classifier.h create mode 100644 source/non_matching/mesh_classifier.cc create mode 100644 source/non_matching/mesh_classifier.inst.in create mode 100644 tests/non_matching/mesh_classifier.cc create mode 100644 tests/non_matching/mesh_classifier.output diff --git a/include/deal.II/non_matching/mesh_classifier.h b/include/deal.II/non_matching/mesh_classifier.h new file mode 100644 index 0000000000..e6bae28b96 --- /dev/null +++ b/include/deal.II/non_matching/mesh_classifier.h @@ -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 + +#include + +#include +#include + +#include + +#include + +DEAL_II_NAMESPACE_OPEN +namespace NonMatching +{ + namespace internal + { + namespace MeshClassifierImplementation + { + template + 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 &level_set = ... + * + * MeshClassifier 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 + 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 + MeshClassifier(const DoFHandler &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 &triangulation, + const Function & level_set, + const FiniteElement &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::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::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::active_cell_iterator &cell, + const unsigned int face_index); + + /** + * Pointer to the triangulation that should be classified. + */ + const SmartPointer> 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> + 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 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 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, GeometryInfo::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 + 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 & + 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::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::active_cell_iterator &cell, + const unsigned int face_index, + Vector &local_dofs) = 0; + }; + + } // namespace MeshClassifierImplementation + } // namespace internal +} // namespace NonMatching +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/non_matching/CMakeLists.txt b/source/non_matching/CMakeLists.txt index 2d1b9a8a76..aec821f02a 100644 --- a/source/non_matching/CMakeLists.txt +++ b/source/non_matching/CMakeLists.txt @@ -16,12 +16,14 @@ 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 index 0000000000..4b8b9941d3 --- /dev/null +++ b/source/non_matching/mesh_classifier.cc @@ -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 + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +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 + 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 + class DiscreteLevelSetDescription : public LevelSetDescription + { + public: + /** + * Constructor. + */ + DiscreteLevelSetDescription(const DoFHandler &dof_handler, + const VECTOR & level_set); + + /** + * Return the FECollection of the DoFHandler passed to the constructor. + */ + const hp::FECollection & + 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::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::active_cell_iterator &cell, + const unsigned int face_index, + Vector &local_levelset_values) override; + + private: + /** + * Pointer to the DoFHandler associated with the level set function. + */ + const SmartPointer> dof_handler; + + /** + * Pointer to the vector containing the level set function's global dof + * values. + */ + const SmartPointer level_set; + }; + + + + template + DiscreteLevelSetDescription::DiscreteLevelSetDescription( + const DoFHandler &dof_handler, + const VECTOR & level_set) + : dof_handler(&dof_handler) + , level_set(&level_set) + {} + + + + template + const hp::FECollection & + DiscreteLevelSetDescription::get_fe_collection() const + { + return dof_handler->get_fe_collection(); + } + + + + template + void + DiscreteLevelSetDescription::get_local_level_set_values( + const typename Triangulation::active_cell_iterator &cell, + const unsigned int face_index, + Vector &local_levelset_values) + { + typename DoFHandler::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 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::get(*level_set, + dof_indices[i]); + } + + + + template + unsigned int + DiscreteLevelSetDescription::active_fe_index( + const typename Triangulation::active_cell_iterator &cell) const + { + typename DoFHandler::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 + class AnalyticLevelSetDescription : public LevelSetDescription + { + public: + /** + * Constructor. Takes the Function that describes the geometry and the + * element that this function should be interpolated to. + */ + AnalyticLevelSetDescription(const Function & level_set, + const FiniteElement &element); + + /** + * Returns the finite element passed to the constructor wrapped in a + * collection. + */ + const hp::FECollection & + 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::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::active_cell_iterator &cell, + const unsigned int face_index, + Vector &local_levelset_values) override; + + private: + /** + * Pointer to the level set function. + */ + const SmartPointer> level_set; + + /** + * Collection containing the single element which we locally interpolate + * the level set function to. + */ + const hp::FECollection fe_collection; + + /** + * FEFaceValues object used to transform the support points on a face to + * real space. + */ + FEFaceValues fe_face_values; + }; + + + + template + AnalyticLevelSetDescription::AnalyticLevelSetDescription( + const Function & level_set, + const FiniteElement &element) + : level_set(&level_set) + , fe_collection(element) + , fe_face_values(element, + Quadrature( + element.get_unit_face_support_points()), + update_quadrature_points) + {} + + + + template + void + AnalyticLevelSetDescription::get_local_level_set_values( + const typename Triangulation::active_cell_iterator &cell, + const unsigned int face_index, + Vector &local_levelset_values) + { + AssertDimension(local_levelset_values.size(), + fe_face_values.n_quadrature_points); + + fe_face_values.reinit(cell, face_index); + const std::vector> &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 + const hp::FECollection & + AnalyticLevelSetDescription::get_fe_collection() const + { + return fe_collection; + } + + + + template + unsigned int + AnalyticLevelSetDescription::active_fe_index( + const typename Triangulation::active_cell_iterator &) const + { + return 0; + } + } // namespace MeshClassifierImplementation + } // namespace internal + + + + template + template + MeshClassifier::MeshClassifier(const DoFHandler &dof_handler, + const VECTOR & level_set) + : triangulation(&dof_handler.get_triangulation()) + , level_set_description( + std::make_unique>( + dof_handler, + level_set)) + { + const hp::FECollection &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 + MeshClassifier::MeshClassifier(const Triangulation &triangulation, + const Function & level_set, + const FiniteElement &element) + : triangulation(&triangulation) + , level_set_description( + std::make_unique>(level_set, + element)) + { + // The level set function must be scalar. + AssertDimension(level_set.n_components, 1); + AssertDimension(element.n_components(), 1); + } + + + + template + void + MeshClassifier::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::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 + LocationToLevelSet + MeshClassifier::determine_face_location_to_levelset( + const typename Triangulation::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 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 + LocationToLevelSet + MeshClassifier::location_to_level_set( + const typename Triangulation::cell_iterator &cell) const + { + return cell_locations.at(cell->active_cell_index()); + } + + + + template + LocationToLevelSet + MeshClassifier::location_to_level_set( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_index) const + { + AssertIndexRange(face_index, GeometryInfo::faces_per_cell); + + return face_locations.at(cell->face(face_index)->index()); + } + + + + template + void + MeshClassifier::initialize() + { + const hp::FECollection &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 &element = fe_collection[i]; + const FE_Q *fe_q = dynamic_cast *>(&element); + Assert(fe_q != nullptr, ExcNotImplemented()); + + const FE_Bernstein fe_bernstein(fe_q->get_degree()); + + const unsigned int dofs_per_face = fe_q->dofs_per_face; + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; f++) + { + FullMatrix 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 index 0000000000..1505a0f8b5 --- /dev/null +++ b/source/non_matching/mesh_classifier.inst.in @@ -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; + \} + } + +for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS) + { + namespace NonMatching + \{ + template MeshClassifier::MeshClassifier( + const DoFHandler &, + const VEC &); + \} + } diff --git a/tests/non_matching/mesh_classifier.cc b/tests/non_matching/mesh_classifier.cc new file mode 100644 index 0000000000..09d2b7d475 --- /dev/null +++ b/tests/non_matching/mesh_classifier.cc @@ -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 +#include +#include +#include + +#include + +#include + +#include +#include + +#include + +#include + +#include + +#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 +void +print_cell_and_face_locations( + const NonMatching::MeshClassifier & classifier, + const typename Triangulation::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::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 +void +classify_with_discrete_level_set(const Function &level_set) +{ + deallog << "discrete:" << std::endl; + Triangulation triangulation; + GridGenerator::hyper_cube(triangulation, -1, 1); + + const FE_Q element(1); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(element); + + Vector discrete_level_set; + discrete_level_set.reinit(dof_handler.n_dofs()); + VectorTools::interpolate(dof_handler, level_set, discrete_level_set); + + NonMatching::MeshClassifier classifier(dof_handler, discrete_level_set); + classifier.reclassify(); + + const typename Triangulation::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 +void +classify_with_analytic_level_set(const Function &level_set) +{ + deallog << "analytic:" << std::endl; + Triangulation triangulation; + GridGenerator::hyper_cube(triangulation, -1, 1); + + const FE_Q element(1); + + NonMatching::MeshClassifier classifier(triangulation, + level_set, + element); + classifier.reclassify(); + + const typename Triangulation::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 +void +classify_with_discrete_and_analytic_level_set(const Function &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 +void +test_positive_function() +{ + deallog << "test_positive_function" << std::endl; + + const Functions::ConstantFunction 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 +void +test_negative_function() +{ + deallog << "test_negative_function" << std::endl; + + const Functions::ConstantFunction 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 +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 origo; + + const Functions::LevelSet::Plane 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 +void +setup_intersected_Q2_positive_coefficients(Vector &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 &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 &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 +void +test_lagrange_coefficents_positive() +{ + deallog << "test_lagrange_coefficents_positive" << std::endl; + Triangulation triangulation; + GridGenerator::hyper_cube(triangulation); + + const FE_Q element(2); + + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(element); + + Vector level_set(element.dofs_per_cell); + setup_intersected_Q2_positive_coefficients(level_set); + + NonMatching::MeshClassifier classifier(dof_handler, level_set); + classifier.reclassify(); + + const typename Triangulation::active_cell_iterator cell = + triangulation.begin_active(); + print_cell_and_face_locations(classifier, cell); + deallog << std::endl; +} + + + +template +void +run_test() +{ + deallog << "dim = " << dim << std::endl; + + test_negative_function(); + test_positive_function(); + test_intersection_x_eq_0_plane(); + // This test doesn't make sense in 1D. + if (dim != 1) + test_lagrange_coefficents_positive(); +} + + + +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 index 0000000000..a29c6fa6be --- /dev/null +++ b/tests/non_matching/mesh_classifier.output @@ -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:: -- 2.39.5