namespace NonMatching
{
- namespace internal
- {
- namespace FEValuesImplementation
- {
- template <int dim>
- class CellWiseFunction;
- }
- } // namespace internal
-
-
/**
* Struct storing UpdateFlags for the 3 regions of a cell, $K$, that is
* defined by the sign of a level set function, $\psi$:
/**
* Object that generates the immersed quadrature rules.
*/
- QuadratureGenerator<dim> quadrature_generator;
-
- /**
- * Function that describes our level set function in reference space.
- */
- const std::unique_ptr<
- internal::FEValuesImplementation::CellWiseFunction<dim>>
- reference_space_level_set;
+ DiscreteQuadratureGenerator<dim> quadrature_generator;
};
-
- namespace internal
- {
- namespace FEValuesImplementation
- {
- /**
- * Interface for a scalar Function which has a
- * set_active_cell(..)-function. That is, a function which we in some way
- * need to associate with a given cell in order to evaluate.
- */
- template <int dim>
- class CellWiseFunction : public Function<dim>
- {
- public:
- /**
- * Destructor. Declared to make it virtual.
- */
- virtual ~CellWiseFunction() = default;
-
- /**
- * Set the cell that the function should be evaluated on.
- */
- virtual void
- set_active_cell(
- const typename Triangulation<dim>::active_cell_iterator &cell) = 0;
- };
- } // namespace FEValuesImplementation
- } // namespace internal
-
} // namespace NonMatching
DEAL_II_NAMESPACE_CLOSE
#include <deal.II/base/quadrature.h>
#include <deal.II/base/std_cxx17/optional.h>
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/grid/tria.h>
+
#include <deal.II/hp/q_collection.h>
#include <deal.II/non_matching/immersed_surface_quadrature.h>
{
template <int dim, int spacedim>
class QGenerator;
- }
- } // namespace internal
+ } // namespace QuadratureGeneratorImplementation
+
+
+ namespace DiscreteQuadratureGeneratorImplementation
+ {
+ template <int dim>
+ class CellWiseFunction;
+ } // namespace DiscreteQuadratureGeneratorImplementation
+ } // namespace internal
/**
};
+ /**
+ * This class generates the same type of immersed quadrature rules as those
+ * described in the QuadratureGenerator class. The difference is that this
+ * class handles the case when the the domain is a discrete level set
+ * function, i.e., when the level set function is described as a
+ * (DoFHandler, Vector) pair. The generate()-function of this class takes a
+ * cell in real space and constructs the immersed quadrature rules in
+ * reference space over this cell. These quadrature rules can then be obtained
+ * with one of the functions:
+ * get_inside_quadrature(),
+ * get_outside_quadrature(), and
+ * get_surface_quadrature().
+ *
+ * Internally, the quadrature generation is done by transforming the discrete
+ * level set function from real space to reference space and using the same
+ * algorithm as in the QuadratureGenerator class.
+ */
+ template <int dim>
+ class DiscreteQuadratureGenerator : public QuadratureGenerator<dim>
+ {
+ public:
+ using AdditionalData = AdditionalQGeneratorData;
+
+ /**
+ * Constructor, the discrete level set function is described by the
+ * incoming DoFHandler and Vector. Pointers to these are stored
+ * internally, so they must have a longer lifetime than the created this
+ * class. The hp::QCollection<1> and AdditionalData is passed to the
+ * QuadratureGenerator class.
+ */
+ template <class VectorType>
+ DiscreteQuadratureGenerator(
+ const hp::QCollection<1> &quadratures1D,
+ const DoFHandler<dim> & dof_handler,
+ const VectorType & level_set,
+ const AdditionalData & additional_data = AdditionalData());
+
+ /**
+ * Construct immersed quadratures rules based on the discrete level
+ * set vector over the incoming cell.
+ */
+ template <bool level_dof_access>
+ void
+ generate(
+ const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell);
+
+ private:
+ /**
+ * Function that describes our level set function in reference space.
+ */
+ std::unique_ptr<internal::DiscreteQuadratureGeneratorImplementation::
+ CellWiseFunction<dim>>
+ reference_space_level_set;
+ };
+
+ /**
+ * This class generates the same type of immersed quadrature rules as those
+ * described in the FaceQuadratureGenerator class. The difference is that this
+ * class handles the case when the the domain is a discrete level set
+ * function, i.e., when the level set function is described as a
+ * (DoFHandler, Vector) pair. The generate()-function of this class takes a
+ * cell in real space plus the respective face index and constructs the
+ * immersed quadrature rules in reference space over this face. These
+ * quadrature rules can then be obtained with one of the functions:
+ * get_inside_quadrature(),
+ * get_outside_quadrature(), and
+ * get_surface_quadrature().
+ *
+ * Internally, the quadrature generation is done by transforming the discrete
+ * level set function from real space to reference space and using the same
+ * algorithm as in the FaceQuadratureGenerator class.
+ */
+ template <int dim>
+ class DiscreteFaceQuadratureGenerator : public FaceQuadratureGenerator<dim>
+ {
+ public:
+ using AdditionalData = AdditionalQGeneratorData;
+
+ /**
+ * Constructor, the discrete level set function is described by the
+ * incoming DoFHandler and Vector. Pointers to these are stored
+ * internally, so they must have a longer lifetime than the created this
+ * class. The hp::QCollection<1> and AdditionalData is passed to the
+ * QuadratureGenerator class.
+ */
+ template <class VectorType>
+ DiscreteFaceQuadratureGenerator(
+ const hp::QCollection<1> &quadratures1D,
+ const DoFHandler<dim> & dof_handler,
+ const VectorType & level_set,
+ const AdditionalData & additional_data = AdditionalData());
+
+ /**
+ * Construct immersed quadratures rules based on the discrete level
+ * set vector over the incoming face described by cell and face index.
+ */
+ template <bool level_dof_access>
+ void
+ generate(
+ const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell,
+ const unsigned int face_index);
+
+ private:
+ /**
+ * Function that describes our level set function in reference space.
+ */
+ std::unique_ptr<internal::DiscreteQuadratureGeneratorImplementation::
+ CellWiseFunction<dim>>
+ reference_space_level_set;
+ };
+
+
namespace internal
{
namespace QuadratureGeneratorImplementation
const std::vector<FunctionBounds<dim>> &all_function_bounds);
} // namespace QuadratureGeneratorImplementation
+
+
+ namespace DiscreteQuadratureGeneratorImplementation
+ {
+ /**
+ * Interface for a scalar Function which has a
+ * set_active_cell(..)-function. That is, a function which we in some way
+ * need to associate with a given cell in order to evaluate.
+ */
+ template <int dim>
+ class CellWiseFunction : public Function<dim>
+ {
+ public:
+ /**
+ * Destructor. Declared to make it virtual.
+ */
+ virtual ~CellWiseFunction() = default;
+
+ /**
+ * Set the cell that the function should be evaluated on.
+ */
+ virtual void
+ set_active_cell(
+ const typename Triangulation<dim>::active_cell_iterator &cell) = 0;
+ };
+
+ } // namespace DiscreteQuadratureGeneratorImplementation
} // namespace internal
} // namespace NonMatching
{}
- namespace internal
- {
- namespace FEValuesImplementation
- {
- DeclExceptionMsg(
- ExcCellNotSet,
- "The set_active_cell function has to be called before calling this function.");
-
-
- /**
- * This class evaluates a function defined by a solution vector and a
- * DoFHandler transformed to reference space. To be precise, if we let
- * $\hat{x}$ be a point on the reference cell, this class implements the
- * function
- *
- * $\hat{f}(\hat{x}) = \sum_{j=0}^{n-1} f_j \hat{\phi}_j(\hat{x})$,
- *
- * where $f_j$ are the local solution values and $\hat{\phi}_j(\hat(x))$
- * are the local reference space shape functions. The gradient and Hessian
- * of this function are thus derivatives with respect to the reference
- * space coordinates, $\hat{x}_0, \hat{x}_1, \ldots$.
- *
- * Note that this class is similar to FEFieldFunction, but that
- * FEFieldFunction implements the following function on a given cell, $K$,
- *
- * $f(x) = \sum_{j=0}^{n-1} f_j \hat{\phi}_j(F_K^{-1}(x))$,
- *
- * which has the same coefficients but uses real space basis functions.
- * Here, $F_K$ is the mapping from the reference cell to the real cell.
- *
- * Before calling the value/gradient/hessian function, the set_active_cell
- * function must be called to specify which cell the function should be
- * evaluated on.
- */
- template <int dim, class VectorType = Vector<double>>
- class RefSpaceFEFieldFunction : public CellWiseFunction<dim>
- {
- public:
- /**
- * Constructor. Takes the solution vector and the associated DoFHandler.
- *
- * Pointers to the input arguments are stored internally, so they must
- * have a longer lifetime than the created RefSpaceFEFieldFunction
- * object.
- */
- RefSpaceFEFieldFunction(const DoFHandler<dim> &dof_handler,
- const VectorType & dof_values);
-
- /**
- * @copydoc CellWiseFunction::set_active_cell()
- */
- void
- set_active_cell(const typename Triangulation<dim>::active_cell_iterator
- &cell) override;
-
- /**
- * @copydoc Function::value()
- *
- * @note The set_active_cell function must be called before this function.
- * The incoming point should be on the reference cell, but this is not
- * checked.
- */
- double
- value(const Point<dim> & point,
- const unsigned int component = 0) const override;
-
- /**
- * @copydoc Function::gradient()
- *
- * @note The set_active_cell function must be called before this function.
- * The incoming point should be on the reference cell, but this is not
- * checked.
- */
- Tensor<1, dim>
- gradient(const Point<dim> & point,
- const unsigned int component = 0) const override;
-
- /**
- * @copydoc Function::hessian()
- *
- * @note The set_active_cell function must be called before this function.
- * The incoming point should be on the reference cell, but this is not
- * checked.
- */
- SymmetricTensor<2, dim>
- hessian(const Point<dim> & point,
- const unsigned int component = 0) const override;
-
- private:
- /**
- * Return whether the set_active_cell function has been called.
- */
- bool
- cell_is_set() const;
-
- /**
- * Pointer to the DoFHandler passed to the constructor.
- */
- const SmartPointer<const DoFHandler<dim>> dof_handler;
-
- /**
- * Pointer to the vector of solution coefficients passed to the
- * constructor.
- */
- const SmartPointer<const VectorType> global_dof_values;
-
- /**
- * Pointer to the element associated with the cell in the last call to
- * set_active_cell().
- */
- SmartPointer<const FiniteElement<dim>> element;
-
- /**
- * DOF-indices of the cell in the last call to set_active_cell()
- */
- std::vector<types::global_dof_index> local_dof_indices;
-
- /**
- * Local solution values of the cell in the last call to
- * set_active_cell()
- */
- std::vector<typename VectorType::value_type> local_dof_values;
-
- /**
- * Description of the 1D polynomial basis for tensor product elements
- * used for the fast path of this class using tensor product
- * evaluators.
- */
- std::vector<Polynomials::Polynomial<double>> poly;
-
- /**
- * Renumbering for the tensor-product evaluator in the fast path.
- */
- std::vector<unsigned int> renumber;
-
- /**
- * Check whether the shape functions are linear.
- */
- bool polynomials_are_hat_functions;
- };
-
-
-
- template <int dim, class VectorType>
- RefSpaceFEFieldFunction<dim, VectorType>::RefSpaceFEFieldFunction(
- const DoFHandler<dim> &dof_handler,
- const VectorType & dof_values)
- : dof_handler(&dof_handler)
- , global_dof_values(&dof_values)
- {
- Assert(dof_handler.n_dofs() == dof_values.size(),
- ExcDimensionMismatch(dof_handler.n_dofs(), dof_values.size()));
- }
-
-
-
- template <int dim, class VectorType>
- void
- RefSpaceFEFieldFunction<dim, VectorType>::set_active_cell(
- const typename Triangulation<dim>::active_cell_iterator &cell)
- {
- Assert(
- &cell->get_triangulation() == &dof_handler->get_triangulation(),
- ExcMessage(
- "The incoming cell must belong to the triangulation associated with "
- "the DoFHandler passed to the constructor."));
-
- const typename DoFHandler<dim>::active_cell_iterator dof_handler_cell(
- &dof_handler->get_triangulation(),
- cell->level(),
- cell->index(),
- dof_handler);
-
- // Save the element and the local dof values, since this is what we need
- // to evaluate the function.
-
- // Check if we can use the fast path. In case we have a different
- // element from the one used before we need to set up the data
- // structures again.
- if (element != &dof_handler_cell->get_fe())
- {
- poly.clear();
- element = &dof_handler_cell->get_fe();
-
- if (element->n_base_elements() == 1 &&
- dealii::internal::FEPointEvaluation::is_fast_path_supported(
- *element, 0))
- {
- dealii::internal::MatrixFreeFunctions::ShapeInfo<double>
- shape_info;
-
- shape_info.reinit(QMidpoint<1>(), *element, 0);
- renumber = shape_info.lexicographic_numbering;
- poly =
- dealii::internal::FEPointEvaluation::get_polynomial_space(
- element->base_element(0));
-
- polynomials_are_hat_functions =
- (poly.size() == 2 && poly[0].value(0.) == 1. &&
- poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
- poly[1].value(1.) == 1.);
- }
- }
- else
- element = &dof_handler_cell->get_fe();
-
- local_dof_indices.resize(element->dofs_per_cell);
- dof_handler_cell->get_dof_indices(local_dof_indices);
-
- local_dof_values.resize(element->dofs_per_cell);
-
- for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
- local_dof_values[i] =
- dealii::internal::ElementAccess<VectorType>::get(
- *global_dof_values, local_dof_indices[i]);
- }
-
-
-
- template <int dim, class VectorType>
- bool
- RefSpaceFEFieldFunction<dim, VectorType>::cell_is_set() const
- {
- // If set cell hasn't been called the size of local_dof_values will be
- // zero.
- return local_dof_values.size() > 0;
- }
-
-
-
- template <int dim, class VectorType>
- double
- RefSpaceFEFieldFunction<dim, VectorType>::value(
- const Point<dim> & point,
- const unsigned int component) const
- {
- AssertIndexRange(component, this->n_components);
- Assert(cell_is_set(), ExcCellNotSet());
-
- if (!poly.empty() && component == 0)
- {
- // TODO: this could be extended to a component that is not zero
- return dealii::internal::evaluate_tensor_product_value_and_gradient(
- poly,
- local_dof_values,
- point,
- polynomials_are_hat_functions,
- renumber)
- .first;
- }
- else
- {
- double value = 0;
- for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
- value += local_dof_values[i] *
- element->shape_value_component(i, point, component);
-
- return value;
- }
- }
-
-
-
- template <int dim, class VectorType>
- Tensor<1, dim>
- RefSpaceFEFieldFunction<dim, VectorType>::gradient(
- const Point<dim> & point,
- const unsigned int component) const
- {
- AssertIndexRange(component, this->n_components);
- Assert(cell_is_set(), ExcCellNotSet());
-
- if (!poly.empty() && component == 0)
- {
- // TODO: this could be extended to a component that is not zero
- return dealii::internal::evaluate_tensor_product_value_and_gradient(
- poly,
- local_dof_values,
- point,
- polynomials_are_hat_functions,
- renumber)
- .second;
- }
- else
- {
- Tensor<1, dim> gradient;
- for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
- gradient += local_dof_values[i] *
- element->shape_grad_component(i, point, component);
-
- return gradient;
- }
- }
-
-
-
- template <int dim, class VectorType>
- SymmetricTensor<2, dim>
- RefSpaceFEFieldFunction<dim, VectorType>::hessian(
- const Point<dim> & point,
- const unsigned int component) const
- {
- AssertIndexRange(component, this->n_components);
- Assert(cell_is_set(), ExcCellNotSet());
-
- if (!poly.empty() && component == 0)
- {
- // TODO: this could be extended to a component that is not zero
- return dealii::internal::evaluate_tensor_product_hessian(
- poly, local_dof_values, point, renumber);
- }
- else
- {
- Tensor<2, dim> hessian;
- for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
- hessian +=
- local_dof_values[i] *
- element->shape_grad_grad_component(i, point, component);
-
- return symmetrize(hessian);
- }
- }
- } // namespace FEValuesImplementation
- } // namespace internal
-
-
template <int dim>
template <class VectorType>
, q_collection_1D(quadrature)
, region_update_flags(region_update_flags)
, mesh_classifier(&mesh_classifier)
- , quadrature_generator(q_collection_1D, additional_data)
- , reference_space_level_set(
- std::make_unique<internal::FEValuesImplementation::
- RefSpaceFEFieldFunction<dim, VectorType>>(
- dof_handler,
- level_set))
+ , quadrature_generator(q_collection_1D,
+ dof_handler,
+ level_set,
+ additional_data)
{
// Tensor products of each quadrature in q_collection_1D. Used on the
// non-intersected cells.
, q_collection_1D(q_collection_1D)
, region_update_flags(region_update_flags)
, mesh_classifier(&mesh_classifier)
- , quadrature_generator(q_collection_1D, additional_data)
- , reference_space_level_set(
- std::make_unique<internal::FEValuesImplementation::
- RefSpaceFEFieldFunction<dim, VectorType>>(
- dof_handler,
- level_set))
+ , quadrature_generator(q_collection_1D,
+ dof_handler,
+ level_set,
+ additional_data)
{
initialize(q_collection);
}
const unsigned int q1D_index =
q_collection_1D.size() > 1 ? active_fe_index : 0;
quadrature_generator.set_1D_quadrature(q1D_index);
-
- reference_space_level_set->set_active_cell(cell);
- const BoundingBox<dim> unit_box = create_unit_bounding_box<dim>();
- quadrature_generator.generate(*reference_space_level_set, unit_box);
+ quadrature_generator.generate(cell);
const Quadrature<dim> &inside_quadrature =
quadrature_generator.get_inside_quadrature();
for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS)
{
- template class internal::FEValuesImplementation::
- RefSpaceFEFieldFunction<deal_II_dimension, VEC>;
-
template FEValues<deal_II_dimension>::FEValues(
const hp::MappingCollection<deal_II_dimension> &,
const hp::FECollection<deal_II_dimension> &,
#include <deal.II/grid/reference_cell.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/matrix_free/fe_point_evaluation.h>
+
#include <deal.II/non_matching/quadrature_generator.h>
#include <boost/math/special_functions/relative_difference.hpp>
AssertIndexRange(q_index, this->q_collection1D->size());
this->q_index = q_index;
}
+
+
+
} // namespace QuadratureGeneratorImplementation
+
+
+
+ namespace DiscreteQuadratureGeneratorImplementation
+ {
+ DeclExceptionMsg(
+ ExcCellNotSet,
+ "The set_active_cell function has to be called before calling this function.");
+
+
+ /**
+ * This class evaluates a function defined by a solution vector and a
+ * DoFHandler transformed to reference space. To be precise, if we let
+ * $\hat{x}$ be a point on the reference cell, this class implements the
+ * function
+ *
+ * $\hat{f}(\hat{x}) = \sum_{j=0}^{n-1} f_j \hat{\phi}_j(\hat{x})$,
+ *
+ * where $f_j$ are the local solution values and $\hat{\phi}_j(\hat(x))$
+ * are the local reference space shape functions. The gradient and Hessian
+ * of this function are thus derivatives with respect to the reference
+ * space coordinates, $\hat{x}_0, \hat{x}_1, \ldots$.
+ *
+ * Note that this class is similar to FEFieldFunction, but that
+ * FEFieldFunction implements the following function on a given cell, $K$,
+ *
+ * $f(x) = \sum_{j=0}^{n-1} f_j \hat{\phi}_j(F_K^{-1}(x))$,
+ *
+ * which has the same coefficients but uses real space basis functions.
+ * Here, $F_K$ is the mapping from the reference cell to the real cell.
+ *
+ * Before calling the value/gradient/hessian function, the set_active_cell
+ * function must be called to specify which cell the function should be
+ * evaluated on.
+ */
+ template <int dim, class VectorType = Vector<double>>
+ class RefSpaceFEFieldFunction : public CellWiseFunction<dim>
+ {
+ public:
+ /**
+ * Constructor. Takes the solution vector and the associated DoFHandler.
+ *
+ * Pointers to the input arguments are stored internally, so they must
+ * have a longer lifetime than the created RefSpaceFEFieldFunction
+ * object.
+ */
+ RefSpaceFEFieldFunction(const DoFHandler<dim> &dof_handler,
+ const VectorType & dof_values);
+
+ /**
+ * @copydoc CellWiseFunction::set_active_cell()
+ */
+ void
+ set_active_cell(const typename Triangulation<dim>::active_cell_iterator
+ &cell) override;
+
+ /**
+ * @copydoc Function::value()
+ *
+ * @note The set_active_cell function must be called before this function.
+ * The incoming point should be on the reference cell, but this is not
+ * checked.
+ */
+ double
+ value(const Point<dim> & point,
+ const unsigned int component = 0) const override;
+
+ /**
+ * @copydoc Function::gradient()
+ *
+ * @note The set_active_cell function must be called before this function.
+ * The incoming point should be on the reference cell, but this is not
+ * checked.
+ */
+ Tensor<1, dim>
+ gradient(const Point<dim> & point,
+ const unsigned int component = 0) const override;
+
+ /**
+ * @copydoc Function::hessian()
+ *
+ * @note The set_active_cell function must be called before this function.
+ * The incoming point should be on the reference cell, but this is not
+ * checked.
+ */
+ SymmetricTensor<2, dim>
+ hessian(const Point<dim> & point,
+ const unsigned int component = 0) const override;
+
+ private:
+ /**
+ * Return whether the set_active_cell function has been called.
+ */
+ bool
+ cell_is_set() const;
+
+ /**
+ * Pointer to the DoFHandler passed to the constructor.
+ */
+ const SmartPointer<const DoFHandler<dim>> dof_handler;
+
+ /**
+ * Pointer to the vector of solution coefficients passed to the
+ * constructor.
+ */
+ const SmartPointer<const VectorType> global_dof_values;
+
+ /**
+ * Pointer to the element associated with the cell in the last call to
+ * set_active_cell().
+ */
+ SmartPointer<const FiniteElement<dim>> element;
+
+ /**
+ * DOF-indices of the cell in the last call to set_active_cell()
+ */
+ std::vector<types::global_dof_index> local_dof_indices;
+
+ /**
+ * Local solution values of the cell in the last call to
+ * set_active_cell()
+ */
+ std::vector<typename VectorType::value_type> local_dof_values;
+
+ /**
+ * Description of the 1D polynomial basis for tensor product elements
+ * used for the fast path of this class using tensor product
+ * evaluators.
+ */
+ std::vector<Polynomials::Polynomial<double>> poly;
+
+ /**
+ * Renumbering for the tensor-product evaluator in the fast path.
+ */
+ std::vector<unsigned int> renumber;
+
+ /**
+ * Check whether the shape functions are linear.
+ */
+ bool polynomials_are_hat_functions;
+ };
+
+
+
+ template <int dim, class VectorType>
+ RefSpaceFEFieldFunction<dim, VectorType>::RefSpaceFEFieldFunction(
+ const DoFHandler<dim> &dof_handler,
+ const VectorType & dof_values)
+ : dof_handler(&dof_handler)
+ , global_dof_values(&dof_values)
+ {
+ Assert(dof_handler.n_dofs() == dof_values.size(),
+ ExcDimensionMismatch(dof_handler.n_dofs(), dof_values.size()));
+ }
+
+
+
+ template <int dim, class VectorType>
+ void
+ RefSpaceFEFieldFunction<dim, VectorType>::set_active_cell(
+ const typename Triangulation<dim>::active_cell_iterator &cell)
+ {
+ Assert(
+ &cell->get_triangulation() == &dof_handler->get_triangulation(),
+ ExcMessage(
+ "The incoming cell must belong to the triangulation associated with "
+ "the DoFHandler passed to the constructor."));
+
+ const typename DoFHandler<dim>::active_cell_iterator dof_handler_cell(
+ &dof_handler->get_triangulation(),
+ cell->level(),
+ cell->index(),
+ dof_handler);
+
+ // Save the element and the local dof values, since this is what we need
+ // to evaluate the function.
+
+ // Check if we can use the fast path. In case we have a different
+ // element from the one used before we need to set up the data
+ // structures again.
+ if (element != &dof_handler_cell->get_fe())
+ {
+ poly.clear();
+ element = &dof_handler_cell->get_fe();
+
+ if (element->n_base_elements() == 1 &&
+ dealii::internal::FEPointEvaluation::is_fast_path_supported(
+ *element, 0))
+ {
+ dealii::internal::MatrixFreeFunctions::ShapeInfo<double>
+ shape_info;
+
+ shape_info.reinit(QMidpoint<1>(), *element, 0);
+ renumber = shape_info.lexicographic_numbering;
+ poly =
+ dealii::internal::FEPointEvaluation::get_polynomial_space(
+ element->base_element(0));
+
+ polynomials_are_hat_functions =
+ (poly.size() == 2 && poly[0].value(0.) == 1. &&
+ poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
+ poly[1].value(1.) == 1.);
+ }
+ }
+ else
+ element = &dof_handler_cell->get_fe();
+
+ local_dof_indices.resize(element->dofs_per_cell);
+ dof_handler_cell->get_dof_indices(local_dof_indices);
+
+ local_dof_values.resize(element->dofs_per_cell);
+
+ for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
+ local_dof_values[i] =
+ dealii::internal::ElementAccess<VectorType>::get(
+ *global_dof_values, local_dof_indices[i]);
+ }
+
+
+
+ template <int dim, class VectorType>
+ bool
+ RefSpaceFEFieldFunction<dim, VectorType>::cell_is_set() const
+ {
+ // If set cell hasn't been called the size of local_dof_values will be
+ // zero.
+ return local_dof_values.size() > 0;
+ }
+
+
+
+ template <int dim, class VectorType>
+ double
+ RefSpaceFEFieldFunction<dim, VectorType>::value(
+ const Point<dim> & point,
+ const unsigned int component) const
+ {
+ AssertIndexRange(component, this->n_components);
+ Assert(cell_is_set(), ExcCellNotSet());
+
+ if (!poly.empty() && component == 0)
+ {
+ // TODO: this could be extended to a component that is not zero
+ return dealii::internal::evaluate_tensor_product_value_and_gradient(
+ poly,
+ local_dof_values,
+ point,
+ polynomials_are_hat_functions,
+ renumber)
+ .first;
+ }
+ else
+ {
+ double value = 0;
+ for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
+ value += local_dof_values[i] *
+ element->shape_value_component(i, point, component);
+
+ return value;
+ }
+ }
+
+
+
+ template <int dim, class VectorType>
+ Tensor<1, dim>
+ RefSpaceFEFieldFunction<dim, VectorType>::gradient(
+ const Point<dim> & point,
+ const unsigned int component) const
+ {
+ AssertIndexRange(component, this->n_components);
+ Assert(cell_is_set(), ExcCellNotSet());
+
+ if (!poly.empty() && component == 0)
+ {
+ // TODO: this could be extended to a component that is not zero
+ return dealii::internal::evaluate_tensor_product_value_and_gradient(
+ poly,
+ local_dof_values,
+ point,
+ polynomials_are_hat_functions,
+ renumber)
+ .second;
+ }
+ else
+ {
+ Tensor<1, dim> gradient;
+ for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
+ gradient += local_dof_values[i] *
+ element->shape_grad_component(i, point, component);
+
+ return gradient;
+ }
+ }
+
+
+
+ template <int dim, class VectorType>
+ SymmetricTensor<2, dim>
+ RefSpaceFEFieldFunction<dim, VectorType>::hessian(
+ const Point<dim> & point,
+ const unsigned int component) const
+ {
+ AssertIndexRange(component, this->n_components);
+ Assert(cell_is_set(), ExcCellNotSet());
+
+ if (!poly.empty() && component == 0)
+ {
+ // TODO: this could be extended to a component that is not zero
+ return dealii::internal::evaluate_tensor_product_hessian(
+ poly, local_dof_values, point, renumber);
+ }
+ else
+ {
+ Tensor<2, dim> hessian;
+ for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
+ hessian +=
+ local_dof_values[i] *
+ element->shape_grad_grad_component(i, point, component);
+
+ return symmetrize(hessian);
+ }
+ }
+ } // namespace DiscreteQuadratureGeneratorImplementation
} // namespace internal
{
return surface_quadrature;
}
+
+
+
+ template <int dim>
+ template <class VectorType>
+ DiscreteQuadratureGenerator<dim>::DiscreteQuadratureGenerator(
+ const hp::QCollection<1> &quadratures1D,
+ const DoFHandler<dim> & dof_handler,
+ const VectorType & level_set,
+ const AdditionalData & additional_data)
+ : QuadratureGenerator<dim>(quadratures1D, additional_data)
+ , reference_space_level_set(
+ std::make_unique<internal::DiscreteQuadratureGeneratorImplementation::
+ RefSpaceFEFieldFunction<dim, VectorType>>(
+ dof_handler,
+ level_set))
+ {}
+
+
+
+ template <int dim>
+ template <bool level_dof_access>
+ void
+ DiscreteQuadratureGenerator<dim>::generate(
+ const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell)
+ {
+ reference_space_level_set->set_active_cell(cell);
+ const BoundingBox<dim> unit_box = create_unit_bounding_box<dim>();
+ QuadratureGenerator<dim>::generate(*reference_space_level_set, unit_box);
+ }
+
+
+
+ template <int dim>
+ template <class VectorType>
+ DiscreteFaceQuadratureGenerator<dim>::DiscreteFaceQuadratureGenerator(
+ const hp::QCollection<1> &quadratures1D,
+ const DoFHandler<dim> & dof_handler,
+ const VectorType & level_set,
+ const AdditionalData & additional_data)
+ : FaceQuadratureGenerator<dim>(quadratures1D, additional_data)
+ , reference_space_level_set(
+ std::make_unique<internal::DiscreteQuadratureGeneratorImplementation::
+ RefSpaceFEFieldFunction<dim, VectorType>>(
+ dof_handler,
+ level_set))
+ {}
+
+
+
+ template <int dim>
+ template <bool level_dof_access>
+ void
+ DiscreteFaceQuadratureGenerator<dim>::generate(
+ const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell,
+ const unsigned int face_index)
+ {
+ reference_space_level_set->set_active_cell(cell);
+ const BoundingBox<dim> unit_box = create_unit_bounding_box<dim>();
+ FaceQuadratureGenerator<dim>::generate(*reference_space_level_set,
+ unit_box,
+ face_index);
+ }
} // namespace NonMatching
#include "quadrature_generator.inst"
DEAL_II_NAMESPACE_CLOSE
namespace NonMatching
\{
template class QuadratureGenerator<deal_II_dimension>;
+ template class DiscreteQuadratureGenerator<deal_II_dimension>;
#if 1 < deal_II_dimension
template class FaceQuadratureGenerator<deal_II_dimension>;
#endif
+ template class DiscreteFaceQuadratureGenerator<deal_II_dimension>;
namespace internal
\{
UpThroughDimensionCreator<deal_II_dimension, deal_II_space_dimension>;
#endif
}
+
+for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS)
+ {
+ template NonMatching::DiscreteQuadratureGenerator<deal_II_dimension>::
+ DiscreteQuadratureGenerator(const hp::QCollection<1> &,
+ const DoFHandler<deal_II_dimension> &,
+ const VEC &,
+ const AdditionalData &);
+
+ template NonMatching::DiscreteFaceQuadratureGenerator<deal_II_dimension>::
+ DiscreteFaceQuadratureGenerator(const hp::QCollection<1> &,
+ const DoFHandler<deal_II_dimension> &,
+ const VEC &,
+ const AdditionalData &);
+ }
+
+// Template generate function
+for (deal_II_dimension : DIMENSIONS; lda : BOOL)
+ {
+ template void
+ NonMatching::DiscreteQuadratureGenerator<deal_II_dimension>::generate(
+ const TriaIterator<
+ DoFCellAccessor<deal_II_dimension, deal_II_dimension, lda>> &);
+
+ template void
+ NonMatching::DiscreteFaceQuadratureGenerator<deal_II_dimension>::generate(
+ const TriaIterator<
+ DoFCellAccessor<deal_II_dimension, deal_II_dimension, lda>> &,
+ const unsigned int);
+ }
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/function_level_set.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_cartesian.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/non_matching/mesh_classifier.h>
+#include <deal.II/non_matching/quadrature_generator.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+#include "quadrature_printing.h"
+
+using namespace dealii;
+
+template <int dim>
+class Test
+{
+public:
+ Test();
+
+ void
+ run();
+
+private:
+ void
+ setup_mesh();
+
+ // Setup a discrete level set function corresponding to
+ // $\psi(x) = (x_0 - 1.5) = 0$
+ void
+ setup_discrete_level_set();
+
+ void
+ test_discrete_quadrature_generator();
+
+ void
+ test_discrete_face_quadrature_generator();
+
+ Triangulation<dim> triangulation;
+ hp::FECollection<dim> fe_collection;
+ DoFHandler<dim> dof_handler;
+
+ hp::MappingCollection<dim> mapping_collection;
+ hp::QCollection<dim> q_collection;
+ hp::QCollection<1> q_collection1D;
+
+ Vector<double> level_set;
+ NonMatching::MeshClassifier<dim> mesh_classifier;
+};
+
+
+
+template <int dim>
+Test<dim>::Test()
+ : dof_handler(triangulation)
+ , mesh_classifier(dof_handler, level_set)
+{
+ fe_collection.push_back(FE_Q<dim>(1));
+ mapping_collection.push_back(MappingCartesian<dim>());
+ const unsigned int n_quadrature_points = 1;
+ q_collection.push_back(QGauss<dim>(n_quadrature_points));
+ q_collection1D.push_back(QGauss<1>(n_quadrature_points));
+}
+
+
+
+template <int dim>
+void
+Test<dim>::run()
+{
+ setup_mesh();
+ dof_handler.distribute_dofs(fe_collection);
+ setup_discrete_level_set();
+ mesh_classifier.reclassify();
+ test_discrete_quadrature_generator();
+ test_discrete_face_quadrature_generator();
+}
+
+
+
+template <int dim>
+void
+Test<dim>::setup_mesh()
+{
+ const Point<dim> lower_left;
+ Point<dim> upper_right;
+ upper_right[0] = 2.;
+
+ for (unsigned int d = 1; d < dim; ++d)
+ {
+ upper_right[d] = 1.;
+ }
+
+ GridGenerator::hyper_rectangle(triangulation, lower_left, upper_right);
+}
+
+
+
+template <int dim>
+void
+Test<dim>::setup_discrete_level_set()
+{
+ Point<dim> point_on_zero_contour;
+ point_on_zero_contour[0] = 1.5;
+ const Functions::LevelSet::Plane<dim> analytical_levelset(
+ point_on_zero_contour, Point<dim>::unit_vector(0));
+
+ level_set.reinit(dof_handler.n_dofs());
+ VectorTools::interpolate(dof_handler, analytical_levelset, level_set);
+}
+
+
+
+template <int dim>
+void
+Test<dim>::test_discrete_quadrature_generator()
+{
+ NonMatching::DiscreteQuadratureGenerator<dim> quadrature_generator(
+ q_collection1D, dof_handler, level_set);
+ quadrature_generator.generate(dof_handler.begin_active());
+
+ print_quadrature(quadrature_generator.get_inside_quadrature());
+ print_quadrature(quadrature_generator.get_outside_quadrature());
+ print_surface_quadrature(quadrature_generator.get_surface_quadrature());
+}
+
+
+
+template <int dim>
+void
+Test<dim>::test_discrete_face_quadrature_generator()
+{
+ NonMatching::DiscreteFaceQuadratureGenerator<dim> face_quadrature_generator(
+ q_collection1D, dof_handler, level_set);
+
+ const auto &cell = dof_handler.begin_active();
+ for (const auto f : cell->face_indices())
+ {
+ face_quadrature_generator.generate(cell, f);
+
+ print_quadrature(face_quadrature_generator.get_inside_quadrature());
+ print_quadrature(face_quadrature_generator.get_outside_quadrature());
+ print_surface_quadrature(
+ face_quadrature_generator.get_surface_quadrature());
+ }
+}
+
+
+
+template <int dim>
+void
+run_test()
+{
+ deallog << "dim = " << dim << std::endl;
+ Test<dim> test;
+ test.run();
+ deallog << std::endl;
+}
+
+
+
+int
+main()
+{
+ initlog();
+
+ run_test<1>();
+ run_test<2>();
+ run_test<3>();
+}
--- /dev/null
+
+DEAL::dim = 1
+DEAL::0.375000, 0.750000
+DEAL::0.875000, 0.250000
+DEAL::0.750000, 1.00000, 1.00000
+DEAL::1.00000
+DEAL::1.00000
+DEAL::
+DEAL::dim = 2
+DEAL::0.375000, 0.500000, 0.750000
+DEAL::0.875000, 0.500000, 0.250000
+DEAL::0.750000, 0.500000, 1.00000, 1.00000, 0.00000
+DEAL::0.500000, 1.00000
+DEAL::0.500000, 1.00000
+DEAL::0.375000, 0.750000
+DEAL::0.875000, 0.250000
+DEAL::0.750000, 1.00000, 1.00000, 0.00000
+DEAL::0.375000, 0.750000
+DEAL::0.875000, 0.250000
+DEAL::0.750000, 1.00000, 1.00000, 0.00000
+DEAL::
+DEAL::dim = 3
+DEAL::0.375000, 0.500000, 0.500000, 0.750000
+DEAL::0.875000, 0.500000, 0.500000, 0.250000
+DEAL::0.750000, 0.500000, 0.500000, 1.00000, 1.00000, 0.00000, 0.00000
+DEAL::0.500000, 0.500000, 1.00000
+DEAL::0.500000, 0.500000, 1.00000
+DEAL::0.500000, 0.375000, 0.750000
+DEAL::0.500000, 0.875000, 0.250000
+DEAL::0.500000, 0.750000, 1.00000, 1.00000, 0.00000, 0.00000
+DEAL::0.500000, 0.375000, 0.750000
+DEAL::0.500000, 0.875000, 0.250000
+DEAL::0.500000, 0.750000, 1.00000, 1.00000, 0.00000, 0.00000
+DEAL::0.375000, 0.500000, 0.750000
+DEAL::0.875000, 0.500000, 0.250000
+DEAL::0.750000, 0.500000, 1.00000, 1.00000, 0.00000, 0.00000
+DEAL::0.375000, 0.500000, 0.750000
+DEAL::0.875000, 0.500000, 0.250000
+DEAL::0.750000, 0.500000, 1.00000, 1.00000, 0.00000, 0.00000
+DEAL::
projects / dealii.git / commitdiff