From: Maximilian Bergbauer Date: Tue, 12 Apr 2022 11:33:45 +0000 (+0200) Subject: QuadratureGenerator with discrete levelset X-Git-Tag: v9.4.0-rc1~310^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=798248098132481ba40f9ec8e43b7f66c55ba16a;p=dealii.git QuadratureGenerator with discrete levelset --- diff --git a/include/deal.II/non_matching/fe_values.h b/include/deal.II/non_matching/fe_values.h index 1225cf66da..dfa7b79b08 100644 --- a/include/deal.II/non_matching/fe_values.h +++ b/include/deal.II/non_matching/fe_values.h @@ -41,16 +41,6 @@ DEAL_II_NAMESPACE_OPEN namespace NonMatching { - namespace internal - { - namespace FEValuesImplementation - { - template - 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$: @@ -375,45 +365,9 @@ namespace NonMatching /** * Object that generates the immersed quadrature rules. */ - QuadratureGenerator quadrature_generator; - - /** - * Function that describes our level set function in reference space. - */ - const std::unique_ptr< - internal::FEValuesImplementation::CellWiseFunction> - reference_space_level_set; + DiscreteQuadratureGenerator 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 - class CellWiseFunction : public Function - { - 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::active_cell_iterator &cell) = 0; - }; - } // namespace FEValuesImplementation - } // namespace internal - } // namespace NonMatching DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/non_matching/quadrature_generator.h b/include/deal.II/non_matching/quadrature_generator.h index 97653b097a..38dfae22e0 100644 --- a/include/deal.II/non_matching/quadrature_generator.h +++ b/include/deal.II/non_matching/quadrature_generator.h @@ -24,6 +24,10 @@ #include #include +#include + +#include + #include #include @@ -39,8 +43,15 @@ namespace NonMatching { template class QGenerator; - } - } // namespace internal + } // namespace QuadratureGeneratorImplementation + + + namespace DiscreteQuadratureGeneratorImplementation + { + template + class CellWiseFunction; + } // namespace DiscreteQuadratureGeneratorImplementation + } // namespace internal /** @@ -459,6 +470,118 @@ namespace NonMatching }; + /** + * 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 + class DiscreteQuadratureGenerator : public QuadratureGenerator + { + 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 + DiscreteQuadratureGenerator( + const hp::QCollection<1> &quadratures1D, + const DoFHandler & 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 + void + generate( + const TriaIterator> &cell); + + private: + /** + * Function that describes our level set function in reference space. + */ + std::unique_ptr> + 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 + class DiscreteFaceQuadratureGenerator : public FaceQuadratureGenerator + { + 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 + DiscreteFaceQuadratureGenerator( + const hp::QCollection<1> &quadratures1D, + const DoFHandler & 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 + void + generate( + const TriaIterator> &cell, + const unsigned int face_index); + + private: + /** + * Function that describes our level set function in reference space. + */ + std::unique_ptr> + reference_space_level_set; + }; + + namespace internal { namespace QuadratureGeneratorImplementation @@ -1247,6 +1370,33 @@ namespace NonMatching const std::vector> &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 + class CellWiseFunction : public Function + { + 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::active_cell_iterator &cell) = 0; + }; + + } // namespace DiscreteQuadratureGeneratorImplementation } // namespace internal } // namespace NonMatching diff --git a/source/non_matching/fe_values.cc b/source/non_matching/fe_values.cc index d4e0bcd40a..75d82d2526 100644 --- a/source/non_matching/fe_values.cc +++ b/source/non_matching/fe_values.cc @@ -43,332 +43,6 @@ 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 > - class RefSpaceFEFieldFunction : public CellWiseFunction - { - 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 &dof_handler, - const VectorType & dof_values); - - /** - * @copydoc CellWiseFunction::set_active_cell() - */ - void - set_active_cell(const typename Triangulation::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 & 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 & 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 & 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> dof_handler; - - /** - * Pointer to the vector of solution coefficients passed to the - * constructor. - */ - const SmartPointer global_dof_values; - - /** - * Pointer to the element associated with the cell in the last call to - * set_active_cell(). - */ - SmartPointer> element; - - /** - * DOF-indices of the cell in the last call to set_active_cell() - */ - std::vector local_dof_indices; - - /** - * Local solution values of the cell in the last call to - * set_active_cell() - */ - std::vector 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> poly; - - /** - * Renumbering for the tensor-product evaluator in the fast path. - */ - std::vector renumber; - - /** - * Check whether the shape functions are linear. - */ - bool polynomials_are_hat_functions; - }; - - - - template - RefSpaceFEFieldFunction::RefSpaceFEFieldFunction( - const DoFHandler &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 - void - RefSpaceFEFieldFunction::set_active_cell( - const typename Triangulation::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::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 - 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::get( - *global_dof_values, local_dof_indices[i]); - } - - - - template - bool - RefSpaceFEFieldFunction::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 - double - RefSpaceFEFieldFunction::value( - const Point & 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 - Tensor<1, dim> - RefSpaceFEFieldFunction::gradient( - const Point & 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 - SymmetricTensor<2, dim> - RefSpaceFEFieldFunction::hessian( - const Point & 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 template @@ -384,12 +58,10 @@ namespace NonMatching , 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>( - 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. @@ -418,12 +90,10 @@ namespace NonMatching , 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>( - dof_handler, - level_set)) + , quadrature_generator(q_collection_1D, + dof_handler, + level_set, + additional_data) { initialize(q_collection); } @@ -513,10 +183,7 @@ namespace NonMatching 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 unit_box = create_unit_bounding_box(); - quadrature_generator.generate(*reference_space_level_set, unit_box); + quadrature_generator.generate(cell); const Quadrature &inside_quadrature = quadrature_generator.get_inside_quadrature(); diff --git a/source/non_matching/fe_values.inst.in b/source/non_matching/fe_values.inst.in index e94c54ec6f..e9cdf571c4 100644 --- a/source/non_matching/fe_values.inst.in +++ b/source/non_matching/fe_values.inst.in @@ -21,9 +21,6 @@ for (deal_II_dimension : DIMENSIONS) for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS) { - template class internal::FEValuesImplementation:: - RefSpaceFEFieldFunction; - template FEValues::FEValues( const hp::MappingCollection &, const hp::FECollection &, diff --git a/source/non_matching/quadrature_generator.cc b/source/non_matching/quadrature_generator.cc index da551e41a1..1f5ba03a16 100644 --- a/source/non_matching/quadrature_generator.cc +++ b/source/non_matching/quadrature_generator.cc @@ -17,6 +17,20 @@ #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + #include #include @@ -1251,7 +1265,334 @@ namespace NonMatching 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 > + class RefSpaceFEFieldFunction : public CellWiseFunction + { + 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 &dof_handler, + const VectorType & dof_values); + + /** + * @copydoc CellWiseFunction::set_active_cell() + */ + void + set_active_cell(const typename Triangulation::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 & 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 & 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 & 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> dof_handler; + + /** + * Pointer to the vector of solution coefficients passed to the + * constructor. + */ + const SmartPointer global_dof_values; + + /** + * Pointer to the element associated with the cell in the last call to + * set_active_cell(). + */ + SmartPointer> element; + + /** + * DOF-indices of the cell in the last call to set_active_cell() + */ + std::vector local_dof_indices; + + /** + * Local solution values of the cell in the last call to + * set_active_cell() + */ + std::vector 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> poly; + + /** + * Renumbering for the tensor-product evaluator in the fast path. + */ + std::vector renumber; + + /** + * Check whether the shape functions are linear. + */ + bool polynomials_are_hat_functions; + }; + + + + template + RefSpaceFEFieldFunction::RefSpaceFEFieldFunction( + const DoFHandler &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 + void + RefSpaceFEFieldFunction::set_active_cell( + const typename Triangulation::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::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 + 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::get( + *global_dof_values, local_dof_indices[i]); + } + + + + template + bool + RefSpaceFEFieldFunction::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 + double + RefSpaceFEFieldFunction::value( + const Point & 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 + Tensor<1, dim> + RefSpaceFEFieldFunction::gradient( + const Point & 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 + SymmetricTensor<2, dim> + RefSpaceFEFieldFunction::hessian( + const Point & 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 @@ -1525,6 +1866,69 @@ namespace NonMatching { return surface_quadrature; } + + + + template + template + DiscreteQuadratureGenerator::DiscreteQuadratureGenerator( + const hp::QCollection<1> &quadratures1D, + const DoFHandler & dof_handler, + const VectorType & level_set, + const AdditionalData & additional_data) + : QuadratureGenerator(quadratures1D, additional_data) + , reference_space_level_set( + std::make_unique>( + dof_handler, + level_set)) + {} + + + + template + template + void + DiscreteQuadratureGenerator::generate( + const TriaIterator> &cell) + { + reference_space_level_set->set_active_cell(cell); + const BoundingBox unit_box = create_unit_bounding_box(); + QuadratureGenerator::generate(*reference_space_level_set, unit_box); + } + + + + template + template + DiscreteFaceQuadratureGenerator::DiscreteFaceQuadratureGenerator( + const hp::QCollection<1> &quadratures1D, + const DoFHandler & dof_handler, + const VectorType & level_set, + const AdditionalData & additional_data) + : FaceQuadratureGenerator(quadratures1D, additional_data) + , reference_space_level_set( + std::make_unique>( + dof_handler, + level_set)) + {} + + + + template + template + void + DiscreteFaceQuadratureGenerator::generate( + const TriaIterator> &cell, + const unsigned int face_index) + { + reference_space_level_set->set_active_cell(cell); + const BoundingBox unit_box = create_unit_bounding_box(); + FaceQuadratureGenerator::generate(*reference_space_level_set, + unit_box, + face_index); + } } // namespace NonMatching #include "quadrature_generator.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/source/non_matching/quadrature_generator.inst.in b/source/non_matching/quadrature_generator.inst.in index 3b6d03c08a..d956a98bd1 100644 --- a/source/non_matching/quadrature_generator.inst.in +++ b/source/non_matching/quadrature_generator.inst.in @@ -19,10 +19,12 @@ for (deal_II_dimension : DIMENSIONS) namespace NonMatching \{ template class QuadratureGenerator; + template class DiscreteQuadratureGenerator; #if 1 < deal_II_dimension template class FaceQuadratureGenerator; #endif + template class DiscreteFaceQuadratureGenerator; namespace internal \{ @@ -76,3 +78,33 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS) UpThroughDimensionCreator; #endif } + +for (VEC : REAL_VECTOR_TYPES; deal_II_dimension : DIMENSIONS) + { + template NonMatching::DiscreteQuadratureGenerator:: + DiscreteQuadratureGenerator(const hp::QCollection<1> &, + const DoFHandler &, + const VEC &, + const AdditionalData &); + + template NonMatching::DiscreteFaceQuadratureGenerator:: + DiscreteFaceQuadratureGenerator(const hp::QCollection<1> &, + const DoFHandler &, + const VEC &, + const AdditionalData &); + } + +// Template generate function +for (deal_II_dimension : DIMENSIONS; lda : BOOL) + { + template void + NonMatching::DiscreteQuadratureGenerator::generate( + const TriaIterator< + DoFCellAccessor> &); + + template void + NonMatching::DiscreteFaceQuadratureGenerator::generate( + const TriaIterator< + DoFCellAccessor> &, + const unsigned int); + } diff --git a/tests/non_matching/discrete_quadrature_generator.cc b/tests/non_matching/discrete_quadrature_generator.cc new file mode 100644 index 0000000000..8501821fd5 --- /dev/null +++ b/tests/non_matching/discrete_quadrature_generator.cc @@ -0,0 +1,191 @@ +// --------------------------------------------------------------------- +// +// 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 +#include +#include + +#include +#include + +#include + +#include +#include + +#include + +#include "../tests.h" + +#include "quadrature_printing.h" + +using namespace dealii; + +template +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 triangulation; + hp::FECollection fe_collection; + DoFHandler dof_handler; + + hp::MappingCollection mapping_collection; + hp::QCollection q_collection; + hp::QCollection<1> q_collection1D; + + Vector level_set; + NonMatching::MeshClassifier mesh_classifier; +}; + + + +template +Test::Test() + : dof_handler(triangulation) + , mesh_classifier(dof_handler, level_set) +{ + fe_collection.push_back(FE_Q(1)); + mapping_collection.push_back(MappingCartesian()); + const unsigned int n_quadrature_points = 1; + q_collection.push_back(QGauss(n_quadrature_points)); + q_collection1D.push_back(QGauss<1>(n_quadrature_points)); +} + + + +template +void +Test::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 +void +Test::setup_mesh() +{ + const Point lower_left; + Point 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 +void +Test::setup_discrete_level_set() +{ + Point point_on_zero_contour; + point_on_zero_contour[0] = 1.5; + const Functions::LevelSet::Plane analytical_levelset( + point_on_zero_contour, Point::unit_vector(0)); + + level_set.reinit(dof_handler.n_dofs()); + VectorTools::interpolate(dof_handler, analytical_levelset, level_set); +} + + + +template +void +Test::test_discrete_quadrature_generator() +{ + NonMatching::DiscreteQuadratureGenerator 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 +void +Test::test_discrete_face_quadrature_generator() +{ + NonMatching::DiscreteFaceQuadratureGenerator 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 +void +run_test() +{ + deallog << "dim = " << dim << std::endl; + Test test; + test.run(); + deallog << std::endl; +} + + + +int +main() +{ + initlog(); + + run_test<1>(); + run_test<2>(); + run_test<3>(); +} diff --git a/tests/non_matching/discrete_quadrature_generator.output b/tests/non_matching/discrete_quadrature_generator.output new file mode 100644 index 0000000000..4dfaa1998f --- /dev/null +++ b/tests/non_matching/discrete_quadrature_generator.output @@ -0,0 +1,40 @@ + +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::