From 6071656d41e9ffbfd9271f5a4500e6070b2e5264 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 13 Feb 2021 21:50:35 +0100 Subject: [PATCH] Remove simplex folders --- .../polynomials_barycentric.h} | 0 include/deal.II/base/polynomials_pyramid.h | 151 ++++ .../polynomials_wedge.h} | 151 +--- include/deal.II/base/quadrature_lib.h | 92 ++ include/deal.II/fe/fe_pyramid_p.h | 141 +++ include/deal.II/fe/fe_simplex_p.h | 212 +++++ include/deal.II/fe/fe_simplex_p_bubbles.h | 119 +++ include/deal.II/fe/fe_wedge_p.h | 141 +++ include/deal.II/grid/grid_generator.h | 46 +- .../matrix_free/matrix_free.templates.h | 3 +- .../matrix_free/shape_info.templates.h | 6 +- include/deal.II/matrix_free/util.h | 3 +- .../numerics/data_out_dof_data.templates.h | 4 +- include/deal.II/simplex/fe_lib.h | 517 ----------- include/deal.II/simplex/grid_generator.h | 80 -- include/deal.II/simplex/quadrature_lib.h | 118 --- source/CMakeLists.txt | 1 - source/base/CMakeLists.txt | 3 + .../polynomials_barycentric.cc} | 2 +- .../polynomials_pyramid.cc} | 229 +---- source/base/polynomials_wedge.cc | 255 ++++++ source/base/qprojector.cc | 3 +- source/base/quadrature_lib.cc | 457 ++++++++++ source/fe/CMakeLists.txt | 10 +- source/fe/fe_dgq.cc | 6 +- source/fe/fe_pyramid_p.cc | 315 +++++++ .../fe_pyramid_p.inst.in} | 11 - source/fe/fe_q.cc | 6 +- source/fe/fe_q_base.cc | 6 +- .../{simplex/fe_lib.cc => fe/fe_simplex_p.cc} | 838 +----------------- .../fe_simplex_p.inst.in} | 20 +- source/fe/fe_simplex_p_bubbles.cc | 390 ++++++++ source/fe/fe_simplex_p_bubbles.inst.in | 24 + source/fe/fe_wedge_p.cc | 325 +++++++ source/fe/fe_wedge_p.inst.in | 28 + source/grid/grid_generator.cc | 190 ++++ source/grid/grid_generator.inst.in | 17 + source/grid/reference_cell.cc | 9 +- source/simplex/CMakeLists.txt | 48 - source/simplex/grid_generator.cc | 214 ----- source/simplex/quadrature_lib.cc | 485 ---------- tests/matrix_free/compute_diagonal_05.cc | 7 +- .../multigrid_util.h | 5 +- tests/simplex/barycentric_01.cc | 8 +- tests/simplex/cell_measure_01.cc | 6 +- tests/simplex/compute_point_locations_01.cc | 11 +- .../simplex/compute_projection_matrices_01.cc | 12 +- tests/simplex/conv_hex_to_simplex.cc | 2 - tests/simplex/data_out_write_hdf5_01.cc | 8 +- tests/simplex/data_out_write_hdf5_02.cc | 8 +- tests/simplex/data_out_write_vtk_01.cc | 11 +- tests/simplex/data_out_write_vtk_02.cc | 11 +- tests/simplex/data_out_write_vtu_01.cc | 8 +- tests/simplex/extract_boundary_dofs.cc | 8 +- tests/simplex/fe_lib_01.cc | 5 +- tests/simplex/fe_lib_02.cc | 8 +- tests/simplex/fe_p_bubbles_01.cc | 8 +- tests/simplex/fe_p_bubbles_02.cc | 8 +- tests/simplex/get_projection_matrix_01.cc | 12 +- tests/simplex/hanging_nodes_01.cc | 8 +- tests/simplex/hanging_nodes_02.cc | 8 +- tests/simplex/hanging_nodes_03.cc | 8 +- tests/simplex/hanging_nodes_hybrid_01.cc | 6 +- tests/simplex/hanging_nodes_hybrid_02.cc | 6 +- tests/simplex/hanging_nodes_hybrid_03.cc | 6 +- tests/simplex/mapping_fe_01.cc | 9 +- tests/simplex/mapping_fe_fields_01.cc | 11 +- tests/simplex/mapping_fe_fields_02.cc | 11 +- tests/simplex/mapping_transformations_01.cc | 6 +- tests/simplex/matrix_free_01.cc | 10 +- tests/simplex/matrix_free_02.cc | 10 +- tests/simplex/matrix_free_03.cc | 10 +- tests/simplex/matrix_free_04.cc | 10 +- tests/simplex/matrix_free_shape_info_01.cc | 11 +- tests/simplex/poisson_01.cc | 9 +- tests/simplex/poisson_02.cc | 8 +- tests/simplex/polynomials_01.cc | 3 +- tests/simplex/q_projection_01.cc | 3 +- tests/simplex/q_witherden_vincent_01.cc | 3 +- tests/simplex/quadrature_lib_01.cc | 2 +- tests/simplex/reference_cell_kind_01.cc | 4 +- tests/simplex/refinement_01.cc | 3 +- tests/simplex/refinement_02.cc | 3 +- tests/simplex/refinement_03.cc | 3 +- tests/simplex/step-01.cc | 2 - tests/simplex/step-02.cc | 7 +- tests/simplex/step-03.cc | 8 +- tests/simplex/step-04.cc | 8 +- tests/simplex/step-06.cc | 7 +- tests/simplex/step-07.cc | 7 +- tests/simplex/step-08.cc | 8 +- tests/simplex/step-12.cc | 8 +- tests/simplex/step-12a.cc | 12 +- tests/simplex/step-18.cc | 10 +- tests/simplex/step-20.cc | 8 +- tests/simplex/step-23.cc | 10 +- tests/simplex/step-38.cc | 11 +- tests/simplex/step-55.cc | 9 +- tests/simplex/step-68.cc | 9 +- tests/simplex/step-74.cc | 7 +- tests/simplex/unit_tangential_vectors_01.cc | 4 +- tests/simplex/variable_face_quadratures_03.cc | 11 +- tests/simplex/wedge_01.cc | 7 +- 103 files changed, 3241 insertions(+), 2905 deletions(-) rename include/deal.II/{simplex/barycentric_polynomials.h => base/polynomials_barycentric.h} (100%) create mode 100644 include/deal.II/base/polynomials_pyramid.h rename include/deal.II/{simplex/polynomials.h => base/polynomials_wedge.h} (54%) create mode 100644 include/deal.II/fe/fe_pyramid_p.h create mode 100644 include/deal.II/fe/fe_simplex_p.h create mode 100644 include/deal.II/fe/fe_simplex_p_bubbles.h create mode 100644 include/deal.II/fe/fe_wedge_p.h delete mode 100644 include/deal.II/simplex/fe_lib.h delete mode 100644 include/deal.II/simplex/grid_generator.h delete mode 100644 include/deal.II/simplex/quadrature_lib.h rename source/{simplex/barycentric_polynomials.cc => base/polynomials_barycentric.cc} (99%) rename source/{simplex/polynomials.cc => base/polynomials_pyramid.cc} (53%) create mode 100644 source/base/polynomials_wedge.cc create mode 100644 source/fe/fe_pyramid_p.cc rename source/{simplex/fe_lib.inst.in => fe/fe_pyramid_p.inst.in} (61%) rename source/{simplex/fe_lib.cc => fe/fe_simplex_p.cc} (52%) rename source/{simplex/grid_generator.inst.in => fe/fe_simplex_p.inst.in} (51%) create mode 100644 source/fe/fe_simplex_p_bubbles.cc create mode 100644 source/fe/fe_simplex_p_bubbles.inst.in create mode 100644 source/fe/fe_wedge_p.cc create mode 100644 source/fe/fe_wedge_p.inst.in delete mode 100644 source/simplex/CMakeLists.txt delete mode 100644 source/simplex/grid_generator.cc delete mode 100644 source/simplex/quadrature_lib.cc diff --git a/include/deal.II/simplex/barycentric_polynomials.h b/include/deal.II/base/polynomials_barycentric.h similarity index 100% rename from include/deal.II/simplex/barycentric_polynomials.h rename to include/deal.II/base/polynomials_barycentric.h diff --git a/include/deal.II/base/polynomials_pyramid.h b/include/deal.II/base/polynomials_pyramid.h new file mode 100644 index 0000000000..31421e95f7 --- /dev/null +++ b/include/deal.II/base/polynomials_pyramid.h @@ -0,0 +1,151 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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_polynomials_pyramid_h +#define dealii_polynomials_pyramid_h + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +/** + * A namespace for functions and classes that provide support for simplex + * reference cell entities, i.e., triangles and tetrahedrons. + * + * @ingroup simplex + */ +namespace Simplex +{ + /** + * Polynomials defined on pyramid entities. This class is basis of + * Simplex::FE_PyramidP. + */ + template + class ScalarPyramidPolynomial : public ScalarPolynomialsBase + { + public: + /** + * Make the dimension available to the outside. + */ + static const unsigned int dimension = dim; + + /* + * Constructor taking the polynomial @p degree as input. + * + * @note Currently, only linear polynomials (degree=1) are implemented. + */ + ScalarPyramidPolynomial(const unsigned int degree); + + /** + * @copydoc ScalarPolynomialsBase::evaluate() + * + * @note Currently, only the vectors @p values and @p grads are filled. + */ + void + evaluate(const Point & unit_point, + std::vector & values, + std::vector> &grads, + std::vector> &grad_grads, + std::vector> &third_derivatives, + std::vector> &fourth_derivatives) const override; + + double + compute_value(const unsigned int i, const Point &p) const override; + + /** + * @copydoc ScalarPolynomialsBase::compute_derivative() + * + * @note Currently, only implemented for first derivative. + */ + template + Tensor + compute_derivative(const unsigned int i, const Point &p) const; + + Tensor<1, dim> + compute_1st_derivative(const unsigned int i, + const Point & p) const override; + + Tensor<2, dim> + compute_2nd_derivative(const unsigned int i, + const Point & p) const override; + + /** + * @copydoc ScalarPolynomialsBase::compute_3rd_derivative() + * + * @note Not implemented yet. + */ + Tensor<3, dim> + compute_3rd_derivative(const unsigned int i, + const Point & p) const override; + + /** + * @copydoc ScalarPolynomialsBase::compute_4th_derivative() + * + * @note Not implemented yet. + */ + Tensor<4, dim> + compute_4th_derivative(const unsigned int i, + const Point & p) const override; + + /** + * @copydoc ScalarPolynomialsBase::compute_grad() + * + * @note Not implemented yet. + */ + Tensor<1, dim> + compute_grad(const unsigned int i, const Point &p) const override; + + /** + * @copydoc ScalarPolynomialsBase::compute_grad_grad() + * + * @note Not implemented yet. + */ + Tensor<2, dim> + compute_grad_grad(const unsigned int i, const Point &p) const override; + + std::string + name() const override; + + virtual std::unique_ptr> + clone() const override; + }; + + + + template + template + Tensor + ScalarPyramidPolynomial::compute_derivative(const unsigned int i, + const Point & p) const + { + Tensor der; + + Assert(order == 1, ExcNotImplemented()); + const auto grad = compute_grad(i, p); + + for (unsigned int i = 0; i < dim; i++) + der[i] = grad[i]; + + return der; + } + +} // namespace Simplex + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/simplex/polynomials.h b/include/deal.II/base/polynomials_wedge.h similarity index 54% rename from include/deal.II/simplex/polynomials.h rename to include/deal.II/base/polynomials_wedge.h index 558e1f97a4..880660f27d 100644 --- a/include/deal.II/simplex/polynomials.h +++ b/include/deal.II/base/polynomials_wedge.h @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2020 by the deal.II authors +// Copyright (C) 2021 by the deal.II authors // // This file is part of the deal.II library. // @@ -14,15 +14,14 @@ // --------------------------------------------------------------------- -#ifndef dealii_simplex_polynomials_h -#define dealii_simplex_polynomials_h +#ifndef dealii_base_polynomials_wedge_h +#define dealii_base_polynomials_wedge_h #include +#include #include -#include - DEAL_II_NAMESPACE_OPEN /** @@ -72,9 +71,6 @@ namespace Simplex std::vector> &third_derivatives, std::vector> &fourth_derivatives) const override; - /** - * @copydoc ScalarPolynomialsBase::compute_value() - */ double compute_value(const unsigned int i, const Point &p) const override; @@ -87,9 +83,6 @@ namespace Simplex Tensor compute_derivative(const unsigned int i, const Point &p) const; - /** - * @copydoc ScalarPolynomialsBase::compute_1st_derivative() - */ Tensor<1, dim> compute_1st_derivative(const unsigned int i, const Point & p) const override; @@ -137,15 +130,9 @@ namespace Simplex Tensor<2, dim> compute_grad_grad(const unsigned int i, const Point &p) const override; - /** - * @copydoc ScalarPolynomialsBase::name() - */ std::string name() const override; - /** - * @copydoc ScalarPolynomialsBase::clone() - */ virtual std::unique_ptr> clone() const override; @@ -163,117 +150,6 @@ namespace Simplex - /** - * Polynomials defined on pyramid entities. This class is basis of - * Simplex::FE_PyramidP. - */ - template - class ScalarPyramidPolynomial : public ScalarPolynomialsBase - { - public: - /** - * Make the dimension available to the outside. - */ - static const unsigned int dimension = dim; - - /* - * Constructor taking the polynomial @p degree as input. - * - * @note Currently, only linear polynomials (degree=1) are implemented. - */ - ScalarPyramidPolynomial(const unsigned int degree); - - /** - * @copydoc ScalarPolynomialsBase::evaluate() - * - * @note Currently, only the vectors @p values and @p grads are filled. - */ - void - evaluate(const Point & unit_point, - std::vector & values, - std::vector> &grads, - std::vector> &grad_grads, - std::vector> &third_derivatives, - std::vector> &fourth_derivatives) const override; - - /** - * @copydoc ScalarPolynomialsBase::compute_value() - */ - double - compute_value(const unsigned int i, const Point &p) const override; - - /** - * @copydoc ScalarPolynomialsBase::compute_derivative() - * - * @note Currently, only implemented for first derivative. - */ - template - Tensor - compute_derivative(const unsigned int i, const Point &p) const; - - /** - * @copydoc ScalarPolynomialsBase::compute_1st_derivative() - */ - Tensor<1, dim> - compute_1st_derivative(const unsigned int i, - const Point & p) const override; - - /** - * @copydoc ScalarPolynomialsBase::compute_2nd_derivative() - */ - Tensor<2, dim> - compute_2nd_derivative(const unsigned int i, - const Point & p) const override; - - /** - * @copydoc ScalarPolynomialsBase::compute_3rd_derivative() - * - * @note Not implemented yet. - */ - Tensor<3, dim> - compute_3rd_derivative(const unsigned int i, - const Point & p) const override; - - /** - * @copydoc ScalarPolynomialsBase::compute_4th_derivative() - * - * @note Not implemented yet. - */ - Tensor<4, dim> - compute_4th_derivative(const unsigned int i, - const Point & p) const override; - - /** - * @copydoc ScalarPolynomialsBase::compute_grad() - * - * @note Not implemented yet. - */ - Tensor<1, dim> - compute_grad(const unsigned int i, const Point &p) const override; - - /** - * @copydoc ScalarPolynomialsBase::compute_grad_grad() - * - * @note Not implemented yet. - */ - Tensor<2, dim> - compute_grad_grad(const unsigned int i, const Point &p) const override; - - /** - * @copydoc ScalarPolynomialsBase::name() - */ - std::string - name() const override; - - /** - * @copydoc ScalarPolynomialsBase::clone() - */ - virtual std::unique_ptr> - clone() const override; - }; - - - template template Tensor @@ -291,25 +167,6 @@ namespace Simplex return der; } - - - template - template - Tensor - ScalarPyramidPolynomial::compute_derivative(const unsigned int i, - const Point & p) const - { - Tensor der; - - AssertDimension(order, 1); - const auto grad = compute_grad(i, p); - - for (unsigned int i = 0; i < dim; i++) - der[i] = grad[i]; - - return der; - } - } // namespace Simplex DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h index 3a832504d6..890440acd4 100644 --- a/include/deal.II/base/quadrature_lib.h +++ b/include/deal.II/base/quadrature_lib.h @@ -792,6 +792,98 @@ public: QSplit(const QSimplex &base, const Point &split_point); }; +namespace Simplex +{ + /** + * Integration rule for simplex entities. + * + * Users specify a number `n_points_1D` as an indication of what polynomial + * degree to be integrated exactly, similarly to the number of points in a + * QGauss quadrature object, even though the present quadrature formula is not + * a tensor product. The given value is translated for n_points_1D=1,2,3,4 to + * following number of quadrature points for 2D and 3D: + * - 2D: 1, 3, 7, 15 + * - 3D: 1, 4, 10, 35 + * + * For 1D, the quadrature rule degenerates to a + * `dealii::QGauss<1>(n_points_1D)`. + * + * @ingroup simplex + */ + template + class QGauss : public QSimplex + { + public: + /** + * Constructor taking the number of quadrature points in 1D direction + * @p n_points_1D. + */ + explicit QGauss(const unsigned int n_points_1D); + }; + + /** + * Witherden-Vincent rules for simplex entities. + * + * Like QGauss, users should specify a number `n_points_1D` as an indication + * of what polynomial degree to be integrated exactly (e.g., for $n$ points, + * the rule can integrate polynomials of degree $2 n - 1$ exactly). The given + * value for n_points_1D = 1, 2, 3, 4, 5 results in the following number of + * quadrature points in 2D and 3D: + * - 2D: 1, 6, 7, 15, 19 + * - 3D: 1, 8, 14, 35, 59 + * + * For 1D, the quadrature rule degenerates to a + * `dealii::QGauss<1>(n_points_1D)`. + * + * These rules match the ones listed for Witherden-Vincent in the quadpy + * @cite quadpy library and were first described in + * @cite witherden2015identification. + * + * @ingroup simplex + */ + template + class QWitherdenVincent : public QSimplex + { + public: + /** + * Constructor taking the number of quadrature points in 1D direction + * @p n_points_1D. + */ + explicit QWitherdenVincent(const unsigned int n_points_1D); + }; + + /** + * Integration rule for wedge entities. + */ + template + class QGaussWedge : public Quadrature + { + public: + /** + * Users specify a number `n_points_1D` as an indication of what polynomial + * degree to be integrated exactly. For details, see the comments of + * Simplex::QGauss. + */ + explicit QGaussWedge(const unsigned int n_points_1D); + }; + + /** + * Integration rule for pyramid entities. + */ + template + class QGaussPyramid : public Quadrature + { + public: + /** + * Users specify a number `n_points_1D` as an indication of what polynomial + * degree to be integrated exactly. For details, see the comments of + * Simplex::QGauss. + */ + explicit QGaussPyramid(const unsigned int n_points_1D); + }; + +} // namespace Simplex + /*@}*/ /* -------------- declaration of explicit specializations ------------- */ diff --git a/include/deal.II/fe/fe_pyramid_p.h b/include/deal.II/fe/fe_pyramid_p.h new file mode 100644 index 0000000000..0079600da6 --- /dev/null +++ b/include/deal.II/fe/fe_pyramid_p.h @@ -0,0 +1,141 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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_fe_fe_p_pyramids_h +#define dealii_fe_fe_p_pyramids_h + +#include + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Simplex +{ + /** + * Base class of FE_PyramidP and FE_PyramidDGP. + * + * @note Only implemented for 3D. + * + * @ingroup simplex + */ + template + class FE_Pyramid : public dealii::FE_Poly + { + public: + /** + * Constructor. + */ + FE_Pyramid(const unsigned int degree, + const internal::GenericDoFsPerObject & dpos, + const typename FiniteElementData::Conformity conformity); + }; + + /** + * Implementation of a scalar Lagrange finite element on a pyramid that yields + * the finite element space of continuous, piecewise polynomials of + * degree $k$. + * + * @ingroup simplex + */ + template + class FE_PyramidP : public FE_Pyramid + { + public: + /** + * Constructor. + */ + FE_PyramidP(const unsigned int degree); + + /** + * @copydoc dealii::FiniteElement::clone() + */ + std::unique_ptr> + clone() const override; + + /** + * Return a string that uniquely identifies a finite element. This class + * returns Simplex::FE_PyramidP(degree), with @p dim and @p degree + * replaced by appropriate values. + */ + std::string + get_name() const override; + + /** + * @copydoc dealii::FiniteElement::compare_for_domination() + */ + FiniteElementDomination::Domination + compare_for_domination(const FiniteElement &fe_other, + const unsigned int codim) const override; + + /** + * @copydoc dealii::FiniteElement::hp_vertex_dof_identities() + */ + std::vector> + hp_vertex_dof_identities( + const FiniteElement &fe_other) const override; + + /** + * @copydoc dealii::FiniteElement::hp_line_dof_identities() + */ + std::vector> + hp_line_dof_identities( + const FiniteElement &fe_other) const override; + + /** + * @copydoc dealii::FiniteElement::hp_quad_dof_identities() + */ + std::vector> + hp_quad_dof_identities(const FiniteElement &fe_other, + const unsigned int face_no = 0) const override; + }; + + /** + * Implementation of a scalar Lagrange finite element on a pyramid that yields + * the finite element space of discontinuous, piecewise polynomials of + * degree $k$. + * + * @ingroup simplex + */ + template + class FE_PyramidDGP : public FE_Pyramid + { + public: + /** + * Constructor. + */ + FE_PyramidDGP(const unsigned int degree); + + /** + * @copydoc dealii::FiniteElement::clone() + */ + std::unique_ptr> + clone() const override; + + /** + * Return a string that uniquely identifies a finite element. This class + * returns Simplex::FE_PyramidDGP(degree), with @p dim and @p degree + * replaced by appropriate values. + */ + std::string + get_name() const override; + }; +} // namespace Simplex + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/fe/fe_simplex_p.h b/include/deal.II/fe/fe_simplex_p.h new file mode 100644 index 0000000000..8105d22376 --- /dev/null +++ b/include/deal.II/fe/fe_simplex_p.h @@ -0,0 +1,212 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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_fe_fe_p_h +#define dealii_fe_fe_p_h + +#include + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Simplex +{ + /** + * Base class of FE_P and FE_DGP. + * + * @note Only implemented for 2D and 3D. + * + * @ingroup simplex + */ + template + class FE_Poly : public dealii::FE_Poly + { + public: + /** + * Constructor. + */ + FE_Poly(const unsigned int degree, + const std::vector & dpo_vector, + const typename FiniteElementData::Conformity conformity); + + /** + * Return a list of constant modes of the element. For this element, the + * list consists of true arguments for all components. + */ + virtual std::pair, std::vector> + get_constant_modes() const override; + + /** + * @copydoc dealii::FiniteElement::get_prolongation_matrix() + * + * @note Only implemented for RefinementCase::isotropic_refinement. + */ + virtual const FullMatrix & + get_prolongation_matrix( + const unsigned int child, + const RefinementCase &refinement_case = + RefinementCase::isotropic_refinement) const override; + + /** + * @copydoc dealii::FiniteElement::get_face_interpolation_matrix() + */ + void + get_face_interpolation_matrix(const FiniteElement &source_fe, + FullMatrix &interpolation_matrix, + const unsigned int face_no) const override; + + /** + * @copydoc dealii::FiniteElement::get_subface_interpolation_matrix() + */ + void + get_subface_interpolation_matrix( + const FiniteElement &x_source_fe, + const unsigned int subface, + FullMatrix & interpolation_matrix, + const unsigned int face_no) const override; + + /** + * @copydoc dealii::FiniteElement::hp_constraints_are_implemented() + */ + bool + hp_constraints_are_implemented() const override; + + /** + * @copydoc dealii::FiniteElement::convert_generalized_support_point_values_to_dof_values() + */ + virtual void + convert_generalized_support_point_values_to_dof_values( + const std::vector> &support_point_values, + std::vector & nodal_values) const override; + + mutable Threads::Mutex mutex; + }; + + + + /** + * Implementation of a scalar Lagrange finite element $P_k$ that yields + * the finite element space of continuous, piecewise polynomials of + * degree $k$. + * + * @ingroup simplex + */ + template + class FE_P : public FE_Poly + { + public: + /** + * Constructor. + */ + FE_P(const unsigned int degree); + + /** + * @copydoc dealii::FiniteElement::clone() + */ + std::unique_ptr> + clone() const override; + + /** + * Return a string that uniquely identifies a finite element. This class + * returns Simplex::FE_P(degree), with @p dim and @p degree + * replaced by appropriate values. + */ + std::string + get_name() const override; + + /** + * @copydoc dealii::FiniteElement::compare_for_domination() + */ + FiniteElementDomination::Domination + compare_for_domination(const FiniteElement &fe_other, + const unsigned int codim) const override; + + /** + * @copydoc dealii::FiniteElement::hp_vertex_dof_identities() + */ + std::vector> + hp_vertex_dof_identities( + const FiniteElement &fe_other) const override; + + /** + * @copydoc dealii::FiniteElement::hp_line_dof_identities() + */ + std::vector> + hp_line_dof_identities( + const FiniteElement &fe_other) const override; + }; + + + + /** + * Implementation of a scalar discontinuous Lagrange finite element + * $P_k$, sometimes denoted as $P_{-k}$, that yields the finite + * element space of discontinuous, piecewise polynomials of degree + * $k$. + * + * @ingroup simplex + */ + template + class FE_DGP : public FE_Poly + { + public: + /** + * Constructor. + */ + FE_DGP(const unsigned int degree); + + /** + * @copydoc dealii::FiniteElement::clone() + */ + std::unique_ptr> + clone() const override; + + /** + * Return a string that uniquely identifies a finite element. This class + * returns Simplex::FE_DGP(degree), with @p dim and @p degree + * replaced by appropriate values. + */ + std::string + get_name() const override; + + /** + * @copydoc dealii::FiniteElement::compare_for_domination() + */ + FiniteElementDomination::Domination + compare_for_domination(const FiniteElement &fe_other, + const unsigned int codim) const override; + + /** + * @copydoc dealii::FiniteElement::hp_vertex_dof_identities() + */ + std::vector> + hp_vertex_dof_identities( + const FiniteElement &fe_other) const override; + + /** + * @copydoc dealii::FiniteElement::hp_line_dof_identities() + */ + std::vector> + hp_line_dof_identities( + const FiniteElement &fe_other) const override; + }; +} // namespace Simplex + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/fe/fe_simplex_p_bubbles.h b/include/deal.II/fe/fe_simplex_p_bubbles.h new file mode 100644 index 0000000000..d6cdb8abe8 --- /dev/null +++ b/include/deal.II/fe/fe_simplex_p_bubbles.h @@ -0,0 +1,119 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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_fe_fe_p_bubbles_h +#define dealii_fe_fe_p_bubbles_h + +#include + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Simplex +{ + /** + * @brief Enriched version of FE_P that can be used with nodal quadrature. + * + * Many explicit time integration schemes require solving a mass matrix at + * each time step. There are various ways around this requirement - for + * example, step-48 replaces the mass matrix with a diagonal approximation, + * which makes the solution step trivial. In step-48, and also commonly for + * tensor-product elements, this is done by computing the mass matrix with a + * lower-order quadrature point based on the nodes of the finite element + * (i.e., the nodal quadrature rule one obtains by using the shape functions + * as an interpolatory basis). + * + * A major drawback of standard simplex-based finite elements is that they + * cannot be used with nodal quadrature since some of the quadrature weights + * end up being either zero or negative, resulting in either an unsolvable or + * unstable approximation to the mass matrix. For example: the shape functions + * of FE_P<2>(2) with support points at vertices have mean values of zero so + * that element cannot be used with mass lumping. + + * This element avoids this issue by replacing the shape functions of FE_P + * with an augmented space amendable to the construction of nodal quadrature + * rules. For example, on the triangle a single basis function is added + * corresponding to interpolation at the centroid (and all other basis + * functions are updated to preserve the partition of unity property). This + * results in shape functions with positive means (i.e., a valid nodal + * quadrature formula). Similarly, in 3D, the polynomial space of FE_P<3>(2) + * is enriched with five additional degrees of freedom (where four have + * support points at face centroids and one has a support point at the + * centroid) to enable construction of valid nodal quadrature rule. + * + * Since this FE space includes bubbles (i.e., extra functions which are + * nonzero only on element interiors), the polynomial degrees of the component + * basis functions are higher than the actual approximation degree of the + * element. For example, with a constructor argument degree = 2 + * in 3D, the polynomials are in fact cubic (degree 3) but the order of the + * approximation is the same as if we were using quadratic (degree 2) finite + * elements. + * + * The 2D quadratic element was first described in @cite fried1975finite. The + * 3D quadratic element implemented here was first described in + * @cite Geevers_2018. Higher degree elements amendable to lumping exist but + * are not yet implemented in this class. + */ + template + class FE_P_Bubbles : public dealii::FE_Poly + { + public: + /** + * Constructor, taking the approximation degree as an argument. The + * polynomial space is typically one degree higher than the approximation + * space for this element: see the general documentation of this class for + * more information. + * + * @note For degree == 1 this element is equivalent to FE_P(1). + */ + FE_P_Bubbles(const unsigned int degree); + + /** + * @copydoc dealii::FiniteElement::clone() + */ + virtual std::unique_ptr> + clone() const override; + + /** + * Return a string that uniquely identifies a finite element. This class + * returns Simplex::FE_F_Bubbles(degree), with + * @p dim, @p spacedim, and @p degree replaced by appropriate values. As + * usual, @p spacedim is omitted in the codimension zero case. + */ + virtual std::string + get_name() const override; + + /** + * @copydoc dealii::FiniteElement::convert_generalized_support_point_values_to_dof_values() + */ + virtual void + convert_generalized_support_point_values_to_dof_values( + const std::vector> &support_point_values, + std::vector & nodal_values) const override; + + protected: + /** + * Degree of the approximation (i.e., the constructor argument). + */ + unsigned int approximation_degree; + }; +} // namespace Simplex + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/fe/fe_wedge_p.h b/include/deal.II/fe/fe_wedge_p.h new file mode 100644 index 0000000000..3a01b2e7c3 --- /dev/null +++ b/include/deal.II/fe/fe_wedge_p.h @@ -0,0 +1,141 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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_fe_fe_p_wedge_h +#define dealii_fe_fe_p_wedge_h + +#include + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Simplex +{ + /** + * Base class of FE_WedgeP and FE_WedgeDGP. + * + * @note Only implemented for 3D. + * + * @ingroup simplex + */ + template + class FE_Wedge : public dealii::FE_Poly + { + public: + /** + * Constructor. + */ + FE_Wedge(const unsigned int degree, + const internal::GenericDoFsPerObject & dpos, + const typename FiniteElementData::Conformity conformity); + }; + + /** + * Implementation of a scalar Lagrange finite element on a wedge that yields + * the finite element space of continuous, piecewise polynomials of + * degree $k$. + * + * @ingroup simplex + */ + template + class FE_WedgeP : public FE_Wedge + { + public: + /** + * Constructor. + */ + FE_WedgeP(const unsigned int degree); + + /** + * @copydoc dealii::FiniteElement::clone() + */ + std::unique_ptr> + clone() const override; + + /** + * Return a string that uniquely identifies a finite element. This class + * returns Simplex::FE_WedgeP(degree), with @p dim and @p degree + * replaced by appropriate values. + */ + std::string + get_name() const override; + + /** + * @copydoc dealii::FiniteElement::compare_for_domination() + */ + FiniteElementDomination::Domination + compare_for_domination(const FiniteElement &fe_other, + const unsigned int codim) const override; + + /** + * @copydoc dealii::FiniteElement::hp_vertex_dof_identities() + */ + std::vector> + hp_vertex_dof_identities( + const FiniteElement &fe_other) const override; + + /** + * @copydoc dealii::FiniteElement::hp_line_dof_identities() + */ + std::vector> + hp_line_dof_identities( + const FiniteElement &fe_other) const override; + + /** + * @copydoc dealii::FiniteElement::hp_quad_dof_identities() + */ + std::vector> + hp_quad_dof_identities(const FiniteElement &fe_other, + const unsigned int face_no = 0) const override; + }; + + /** + * Implementation of a scalar Lagrange finite element on a wedge that yields + * the finite element space of discontinuous, piecewise polynomials of + * degree $k$. + * + * @ingroup simplex + */ + template + class FE_WedgeDGP : public FE_Wedge + { + public: + /** + * Constructor. + */ + FE_WedgeDGP(const unsigned int degree); + + /** + * @copydoc dealii::FiniteElement::clone() + */ + std::unique_ptr> + clone() const override; + + /** + * Return a string that uniquely identifies a finite element. This class + * returns Simplex::FE_WedgeDGP(degree), with @p dim and @p degree + * replaced by appropriate values. + */ + std::string + get_name() const override; + }; +} // namespace Simplex + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index f301aac359..aaf94235ad 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -2146,8 +2146,6 @@ namespace GridGenerator add_parameters(ParameterHandler &prm); }; - - /** * Initialize the given triangulation with a flow field around an airfoil, * i.e., a mesh of C-Type approximating Joukowski or NACA (4 digit) @@ -2207,6 +2205,50 @@ namespace GridGenerator } // namespace Airfoil + /** + * Create a coordinate-parallel brick from the two diagonally opposite + * corner points @p p1 and @p p2. The number of vertices in coordinate + * direction @p i is given by repetitions[i]+1. + * + * @note This function connects internally 4/8 vertices to + * quadrilateral/hexahedral cells and subdivides these into 2/5 + * triangular/tetrahedral cells. + * + * @note Currently, this function only works for `dim==spacedim`. + * + * @ingroup simplex + */ + template + void + subdivided_hyper_rectangle_with_simplices( + Triangulation & tria, + const std::vector &repetitions, + const Point & p1, + const Point & p2, + const bool colorize = false); + + /** + * Initialize the given triangulation with a hypercube (square in 2D and + * cube in 3D) consisting of @p repetitions cells in each direction. + * The hypercube volume is the tensor product interval + * $[left,right]^{\text{dim}}$ in the present number of dimensions, where + * the limits are given as arguments. They default to zero and unity, then + * producing the unit hypercube. + * + * @note This function connects internally 4/8 vertices to + * quadrilateral/hexahedral cells and subdivides these into 2/5 + * triangular/tetrahedral cells. + * + * @ingroup simplex + */ + template + void + subdivided_hyper_cube_with_simplices(Triangulation &tria, + const unsigned int repetitions, + const double p1 = 0.0, + const double p2 = 1.0, + const bool colorize = false); + ///@} /** diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index df3f382117..1e0d715dd4 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -23,6 +23,7 @@ #include #include #include +#include #include #include @@ -41,8 +42,6 @@ #include #include -#include - #ifdef DEAL_II_WITH_TBB # include diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 4298f7a9ba..6cf1cb2125 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -29,7 +29,11 @@ #include #include #include +#include #include +#include +#include +#include #include @@ -38,8 +42,6 @@ #include #include -#include - DEAL_II_NAMESPACE_OPEN diff --git a/include/deal.II/matrix_free/util.h b/include/deal.II/matrix_free/util.h index 495c6a6e6d..dec9542a16 100644 --- a/include/deal.II/matrix_free/util.h +++ b/include/deal.II/matrix_free/util.h @@ -21,13 +21,12 @@ #include #include +#include #include #include -#include - DEAL_II_NAMESPACE_OPEN diff --git a/include/deal.II/numerics/data_out_dof_data.templates.h b/include/deal.II/numerics/data_out_dof_data.templates.h index bf8962e282..4222b7f112 100644 --- a/include/deal.II/numerics/data_out_dof_data.templates.h +++ b/include/deal.II/numerics/data_out_dof_data.templates.h @@ -30,7 +30,9 @@ #include #include +#include #include +#include #include #include @@ -46,8 +48,6 @@ #include #include -#include - #include #include #include diff --git a/include/deal.II/simplex/fe_lib.h b/include/deal.II/simplex/fe_lib.h deleted file mode 100644 index 7cc78009fe..0000000000 --- a/include/deal.II/simplex/fe_lib.h +++ /dev/null @@ -1,517 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2020 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_simplex_fe_lib_h -#define dealii_simplex_fe_lib_h - -#include - -#include - -#include - -DEAL_II_NAMESPACE_OPEN - -namespace Simplex -{ - /** - * Base class of FE_P and FE_DGP. - * - * @note Only implemented for 2D and 3D. - * - * @ingroup simplex - */ - template - class FE_Poly : public dealii::FE_Poly - { - public: - /** - * Constructor. - */ - FE_Poly(const unsigned int degree, - const std::vector & dpo_vector, - const typename FiniteElementData::Conformity conformity); - - /** - * Return a list of constant modes of the element. For this element, the - * list consists of true arguments for all components. - */ - virtual std::pair, std::vector> - get_constant_modes() const override; - - /** - * @copydoc dealii::FiniteElement::get_prolongation_matrix() - * - * @note Only implemented for RefinementCase::isotropic_refinement. - */ - virtual const FullMatrix & - get_prolongation_matrix( - const unsigned int child, - const RefinementCase &refinement_case = - RefinementCase::isotropic_refinement) const override; - - /** - * @copydoc dealii::FiniteElement::get_face_interpolation_matrix() - */ - void - get_face_interpolation_matrix(const FiniteElement &source_fe, - FullMatrix &interpolation_matrix, - const unsigned int face_no) const override; - - /** - * @copydoc dealii::FiniteElement::get_subface_interpolation_matrix() - */ - void - get_subface_interpolation_matrix( - const FiniteElement &x_source_fe, - const unsigned int subface, - FullMatrix & interpolation_matrix, - const unsigned int face_no) const override; - - /** - * @copydoc dealii::FiniteElement::hp_constraints_are_implemented() - */ - bool - hp_constraints_are_implemented() const override; - - /** - * @copydoc dealii::FiniteElement::convert_generalized_support_point_values_to_dof_values() - */ - virtual void - convert_generalized_support_point_values_to_dof_values( - const std::vector> &support_point_values, - std::vector & nodal_values) const override; - - mutable Threads::Mutex mutex; - }; - - - - /** - * Implementation of a scalar Lagrange finite element $P_k$ that yields - * the finite element space of continuous, piecewise polynomials of - * degree $k$. - * - * @ingroup simplex - */ - template - class FE_P : public FE_Poly - { - public: - /** - * Constructor. - */ - FE_P(const unsigned int degree); - - /** - * @copydoc dealii::FiniteElement::clone() - */ - std::unique_ptr> - clone() const override; - - /** - * Return a string that uniquely identifies a finite element. This class - * returns Simplex::FE_P(degree), with @p dim and @p degree - * replaced by appropriate values. - */ - std::string - get_name() const override; - - /** - * @copydoc dealii::FiniteElement::compare_for_domination() - */ - FiniteElementDomination::Domination - compare_for_domination(const FiniteElement &fe_other, - const unsigned int codim) const override; - - /** - * @copydoc dealii::FiniteElement::hp_vertex_dof_identities() - */ - std::vector> - hp_vertex_dof_identities( - const FiniteElement &fe_other) const override; - - /** - * @copydoc dealii::FiniteElement::hp_line_dof_identities() - */ - std::vector> - hp_line_dof_identities( - const FiniteElement &fe_other) const override; - }; - - - - /** - * Implementation of a scalar discontinuous Lagrange finite element - * $P_k$, sometimes denoted as $P_{-k}$, that yields the finite - * element space of discontinuous, piecewise polynomials of degree - * $k$. - * - * @ingroup simplex - */ - template - class FE_DGP : public FE_Poly - { - public: - /** - * Constructor. - */ - FE_DGP(const unsigned int degree); - - /** - * @copydoc dealii::FiniteElement::clone() - */ - std::unique_ptr> - clone() const override; - - /** - * Return a string that uniquely identifies a finite element. This class - * returns Simplex::FE_DGP(degree), with @p dim and @p degree - * replaced by appropriate values. - */ - std::string - get_name() const override; - - /** - * @copydoc dealii::FiniteElement::compare_for_domination() - */ - FiniteElementDomination::Domination - compare_for_domination(const FiniteElement &fe_other, - const unsigned int codim) const override; - - /** - * @copydoc dealii::FiniteElement::hp_vertex_dof_identities() - */ - std::vector> - hp_vertex_dof_identities( - const FiniteElement &fe_other) const override; - - /** - * @copydoc dealii::FiniteElement::hp_line_dof_identities() - */ - std::vector> - hp_line_dof_identities( - const FiniteElement &fe_other) const override; - }; - - /** - * Base class of FE_WedgeP and FE_WedgeDGP. - * - * @note Only implemented for 3D. - * - * @ingroup simplex - */ - template - class FE_Wedge : public dealii::FE_Poly - { - public: - /** - * Constructor. - */ - FE_Wedge(const unsigned int degree, - const internal::GenericDoFsPerObject & dpos, - const typename FiniteElementData::Conformity conformity); - }; - - /** - * Implementation of a scalar Lagrange finite element on a wedge that yields - * the finite element space of continuous, piecewise polynomials of - * degree $k$. - * - * @ingroup simplex - */ - template - class FE_WedgeP : public FE_Wedge - { - public: - /** - * Constructor. - */ - FE_WedgeP(const unsigned int degree); - - /** - * @copydoc dealii::FiniteElement::clone() - */ - std::unique_ptr> - clone() const override; - - /** - * Return a string that uniquely identifies a finite element. This class - * returns Simplex::FE_WedgeP(degree), with @p dim and @p degree - * replaced by appropriate values. - */ - std::string - get_name() const override; - - /** - * @copydoc dealii::FiniteElement::compare_for_domination() - */ - FiniteElementDomination::Domination - compare_for_domination(const FiniteElement &fe_other, - const unsigned int codim) const override; - - /** - * @copydoc dealii::FiniteElement::hp_vertex_dof_identities() - */ - std::vector> - hp_vertex_dof_identities( - const FiniteElement &fe_other) const override; - - /** - * @copydoc dealii::FiniteElement::hp_line_dof_identities() - */ - std::vector> - hp_line_dof_identities( - const FiniteElement &fe_other) const override; - - /** - * @copydoc dealii::FiniteElement::hp_quad_dof_identities() - */ - std::vector> - hp_quad_dof_identities(const FiniteElement &fe_other, - const unsigned int face_no = 0) const override; - }; - - /** - * Implementation of a scalar Lagrange finite element on a wedge that yields - * the finite element space of discontinuous, piecewise polynomials of - * degree $k$. - * - * @ingroup simplex - */ - template - class FE_WedgeDGP : public FE_Wedge - { - public: - /** - * Constructor. - */ - FE_WedgeDGP(const unsigned int degree); - - /** - * @copydoc dealii::FiniteElement::clone() - */ - std::unique_ptr> - clone() const override; - - /** - * Return a string that uniquely identifies a finite element. This class - * returns Simplex::FE_WedgeDGP(degree), with @p dim and @p degree - * replaced by appropriate values. - */ - std::string - get_name() const override; - }; - - /** - * Base class of FE_PyramidP and FE_PyramidDGP. - * - * @note Only implemented for 3D. - * - * @ingroup simplex - */ - template - class FE_Pyramid : public dealii::FE_Poly - { - public: - /** - * Constructor. - */ - FE_Pyramid(const unsigned int degree, - const internal::GenericDoFsPerObject & dpos, - const typename FiniteElementData::Conformity conformity); - }; - - /** - * Implementation of a scalar Lagrange finite element on a pyramid that yields - * the finite element space of continuous, piecewise polynomials of - * degree $k$. - * - * @ingroup simplex - */ - template - class FE_PyramidP : public FE_Pyramid - { - public: - /** - * Constructor. - */ - FE_PyramidP(const unsigned int degree); - - /** - * @copydoc dealii::FiniteElement::clone() - */ - std::unique_ptr> - clone() const override; - - /** - * Return a string that uniquely identifies a finite element. This class - * returns Simplex::FE_PyramidP(degree), with @p dim and @p degree - * replaced by appropriate values. - */ - std::string - get_name() const override; - - /** - * @copydoc dealii::FiniteElement::compare_for_domination() - */ - FiniteElementDomination::Domination - compare_for_domination(const FiniteElement &fe_other, - const unsigned int codim) const override; - - /** - * @copydoc dealii::FiniteElement::hp_vertex_dof_identities() - */ - std::vector> - hp_vertex_dof_identities( - const FiniteElement &fe_other) const override; - - /** - * @copydoc dealii::FiniteElement::hp_line_dof_identities() - */ - std::vector> - hp_line_dof_identities( - const FiniteElement &fe_other) const override; - - /** - * @copydoc dealii::FiniteElement::hp_quad_dof_identities() - */ - std::vector> - hp_quad_dof_identities(const FiniteElement &fe_other, - const unsigned int face_no = 0) const override; - }; - - /** - * Implementation of a scalar Lagrange finite element on a pyramid that yields - * the finite element space of discontinuous, piecewise polynomials of - * degree $k$. - * - * @ingroup simplex - */ - template - class FE_PyramidDGP : public FE_Pyramid - { - public: - /** - * Constructor. - */ - FE_PyramidDGP(const unsigned int degree); - - /** - * @copydoc dealii::FiniteElement::clone() - */ - std::unique_ptr> - clone() const override; - - /** - * Return a string that uniquely identifies a finite element. This class - * returns Simplex::FE_PyramidDGP(degree), with @p dim and @p degree - * replaced by appropriate values. - */ - std::string - get_name() const override; - }; - - /** - * @brief Enriched version of FE_P that can be used with nodal quadrature. - * - * Many explicit time integration schemes require solving a mass matrix at - * each time step. There are various ways around this requirement - for - * example, step-48 replaces the mass matrix with a diagonal approximation, - * which makes the solution step trivial. In step-48, and also commonly for - * tensor-product elements, this is done by computing the mass matrix with a - * lower-order quadrature point based on the nodes of the finite element - * (i.e., the nodal quadrature rule one obtains by using the shape functions - * as an interpolatory basis). - * - * A major drawback of standard simplex-based finite elements is that they - * cannot be used with nodal quadrature since some of the quadrature weights - * end up being either zero or negative, resulting in either an unsolvable or - * unstable approximation to the mass matrix. For example: the shape functions - * of FE_P<2>(2) with support points at vertices have mean values of zero so - * that element cannot be used with mass lumping. - - * This element avoids this issue by replacing the shape functions of FE_P - * with an augmented space amendable to the construction of nodal quadrature - * rules. For example, on the triangle a single basis function is added - * corresponding to interpolation at the centroid (and all other basis - * functions are updated to preserve the partition of unity property). This - * results in shape functions with positive means (i.e., a valid nodal - * quadrature formula). Similarly, in 3D, the polynomial space of FE_P<3>(2) - * is enriched with five additional degrees of freedom (where four have - * support points at face centroids and one has a support point at the - * centroid) to enable construction of valid nodal quadrature rule. - * - * Since this FE space includes bubbles (i.e., extra functions which are - * nonzero only on element interiors), the polynomial degrees of the component - * basis functions are higher than the actual approximation degree of the - * element. For example, with a constructor argument degree = 2 - * in 3D, the polynomials are in fact cubic (degree 3) but the order of the - * approximation is the same as if we were using quadratic (degree 2) finite - * elements. - * - * The 2D quadratic element was first described in @cite fried1975finite. The - * 3D quadratic element implemented here was first described in - * @cite Geevers_2018. Higher degree elements amendable to lumping exist but - * are not yet implemented in this class. - */ - template - class FE_P_Bubbles : public dealii::FE_Poly - { - public: - /** - * Constructor, taking the approximation degree as an argument. The - * polynomial space is typically one degree higher than the approximation - * space for this element: see the general documentation of this class for - * more information. - * - * @note For degree == 1 this element is equivalent to FE_P(1). - */ - FE_P_Bubbles(const unsigned int degree); - - /** - * @copydoc dealii::FiniteElement::clone() - */ - virtual std::unique_ptr> - clone() const override; - - /** - * Return a string that uniquely identifies a finite element. This class - * returns Simplex::FE_F_Bubbles(degree), with - * @p dim, @p spacedim, and @p degree replaced by appropriate values. As - * usual, @p spacedim is omitted in the codimension zero case. - */ - virtual std::string - get_name() const override; - - /** - * @copydoc dealii::FiniteElement::convert_generalized_support_point_values_to_dof_values() - */ - virtual void - convert_generalized_support_point_values_to_dof_values( - const std::vector> &support_point_values, - std::vector & nodal_values) const override; - - protected: - /** - * Degree of the approximation (i.e., the constructor argument). - */ - unsigned int approximation_degree; - }; -} // namespace Simplex - -DEAL_II_NAMESPACE_CLOSE - -#endif diff --git a/include/deal.II/simplex/grid_generator.h b/include/deal.II/simplex/grid_generator.h deleted file mode 100644 index 3eaba788c9..0000000000 --- a/include/deal.II/simplex/grid_generator.h +++ /dev/null @@ -1,80 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2020 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_simplex_grid_generator_h -#define dealii_simplex_grid_generator_h - - -#include - -#include - -#include - -#include - -DEAL_II_NAMESPACE_OPEN - -namespace GridGenerator -{ - /** - * Create a coordinate-parallel brick from the two diagonally opposite - * corner points @p p1 and @p p2. The number of vertices in coordinate - * direction @p i is given by repetitions[i]+1. - * - * @note This function connects internally 4/8 vertices to - * quadrilateral/hexahedral cells and subdivides these into 2/5 - * triangular/tetrahedral cells. - * - * @note Currently, this function only works for `dim==spacedim`. - * - * @ingroup simplex - */ - template - void - subdivided_hyper_rectangle_with_simplices( - Triangulation & tria, - const std::vector &repetitions, - const Point & p1, - const Point & p2, - const bool colorize = false); - - /** - * Initialize the given triangulation with a hypercube (square in 2D and - * cube in 3D) consisting of @p repetitions cells in each direction. - * The hypercube volume is the tensor product interval - * $[left,right]^{\text{dim}}$ in the present number of dimensions, where - * the limits are given as arguments. They default to zero and unity, then - * producing the unit hypercube. - * - * @note This function connects internally 4/8 vertices to - * quadrilateral/hexahedral cells and subdivides these into 2/5 - * triangular/tetrahedral cells. - * - * @ingroup simplex - */ - template - void - subdivided_hyper_cube_with_simplices(Triangulation &tria, - const unsigned int repetitions, - const double p1 = 0.0, - const double p2 = 1.0, - const bool colorize = false); - -} // namespace GridGenerator - -DEAL_II_NAMESPACE_CLOSE - -#endif diff --git a/include/deal.II/simplex/quadrature_lib.h b/include/deal.II/simplex/quadrature_lib.h deleted file mode 100644 index ba55488871..0000000000 --- a/include/deal.II/simplex/quadrature_lib.h +++ /dev/null @@ -1,118 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2020 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_tet_quadrature_lib_h -#define dealii_tet_quadrature_lib_h - - -#include - -#include - -DEAL_II_NAMESPACE_OPEN - -namespace Simplex -{ - /** - * Integration rule for simplex entities. - * - * Users specify a number `n_points_1D` as an indication of what polynomial - * degree to be integrated exactly, similarly to the number of points in a - * QGauss quadrature object, even though the present quadrature formula is not - * a tensor product. The given value is translated for n_points_1D=1,2,3,4 to - * following number of quadrature points for 2D and 3D: - * - 2D: 1, 3, 7, 15 - * - 3D: 1, 4, 10, 35 - * - * For 1D, the quadrature rule degenerates to a `QGauss<1>(n_points_1D)`. - * - * @ingroup simplex - */ - template - class QGauss : public QSimplex - { - public: - /** - * Constructor taking the number of quadrature points in 1D direction - * @p n_points_1D. - */ - explicit QGauss(const unsigned int n_points_1D); - }; - - /** - * Witherden-Vincent rules for simplex entities. - * - * Like QGauss, users should specify a number `n_points_1D` as an indication - * of what polynomial degree to be integrated exactly (e.g., for $n$ points, - * the rule can integrate polynomials of degree $2 n - 1$ exactly). The given - * value for n_points_1D = 1, 2, 3, 4, 5 results in the following number of - * quadrature points in 2D and 3D: - * - 2D: 1, 6, 7, 15, 19 - * - 3D: 1, 8, 14, 35, 59 - * - * For 1D, the quadrature rule degenerates to a `QGauss<1>(n_points_1D)`. - * - * These rules match the ones listed for Witherden-Vincent in the quadpy - * @cite quadpy library and were first described in - * @cite witherden2015identification. - * - * @ingroup simplex - */ - template - class QWitherdenVincent : public QSimplex - { - public: - /** - * Constructor taking the number of quadrature points in 1D direction - * @p n_points_1D. - */ - explicit QWitherdenVincent(const unsigned int n_points_1D); - }; - - /** - * Integration rule for wedge entities. - */ - template - class QGaussWedge : public Quadrature - { - public: - /** - * Users specify a number `n_points_1D` as an indication of what polynomial - * degree to be integrated exactly. For details, see the comments of - * Simplex::QGauss. - */ - explicit QGaussWedge(const unsigned int n_points_1D); - }; - - /** - * Integration rule for pyramid entities. - */ - template - class QGaussPyramid : public Quadrature - { - public: - /** - * Users specify a number `n_points_1D` as an indication of what polynomial - * degree to be integrated exactly. For details, see the comments of - * Simplex::QGauss. - */ - explicit QGaussPyramid(const unsigned int n_points_1D); - }; - -} // namespace Simplex - -DEAL_II_NAMESPACE_CLOSE - -#endif diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index a6d85d8b68..7173d25b85 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -71,7 +71,6 @@ ADD_SUBDIRECTORY(differentiation) ADD_SUBDIRECTORY(physics) ADD_SUBDIRECTORY(optimization/rol) ADD_SUBDIRECTORY(non_matching) -ADD_SUBDIRECTORY(simplex) ADD_SUBDIRECTORY(sundials) FOREACH(build ${DEAL_II_BUILD_TYPES}) diff --git a/source/base/CMakeLists.txt b/source/base/CMakeLists.txt index 2a924e7fb6..c48dce79cc 100644 --- a/source/base/CMakeLists.txt +++ b/source/base/CMakeLists.txt @@ -64,6 +64,7 @@ SET(_unity_include_src polynomials_abf.cc polynomials_adini.cc tensor_polynomials_base.cc + polynomials_barycentric.cc polynomials_bernstein.cc polynomials_bdm.cc polynomials_bernardi_raugel.cc @@ -72,9 +73,11 @@ SET(_unity_include_src polynomial_space.cc polynomials_p.cc polynomials_piecewise.cc + polynomials_pyramid.cc polynomials_rannacher_turek.cc polynomials_raviart_thomas.cc polynomials_rt_bubbles.cc + polynomials_wedge.cc process_grid.cc qprojector.cc quadrature.cc diff --git a/source/simplex/barycentric_polynomials.cc b/source/base/polynomials_barycentric.cc similarity index 99% rename from source/simplex/barycentric_polynomials.cc rename to source/base/polynomials_barycentric.cc index 8e4e5416a5..d916f9e944 100644 --- a/source/simplex/barycentric_polynomials.cc +++ b/source/base/polynomials_barycentric.cc @@ -14,7 +14,7 @@ // --------------------------------------------------------------------- -#include +#include DEAL_II_NAMESPACE_OPEN diff --git a/source/simplex/polynomials.cc b/source/base/polynomials_pyramid.cc similarity index 53% rename from source/simplex/polynomials.cc rename to source/base/polynomials_pyramid.cc index 4808d45584..9321dfb3a9 100644 --- a/source/simplex/polynomials.cc +++ b/source/base/polynomials_pyramid.cc @@ -14,10 +14,8 @@ // --------------------------------------------------------------------- -#include - -#include -#include +#include +#include DEAL_II_NAMESPACE_OPEN @@ -39,230 +37,10 @@ namespace Simplex return 0; } - - unsigned int - compute_n_polynomials_wedge(const unsigned int dim, - const unsigned int degree) - { - if (dim == 3) - { - if (degree == 1) - return 6; - if (degree == 2) - return 18; - } - - Assert(false, ExcNotImplemented()); - - return 0; - } } // namespace - template - ScalarWedgePolynomial::ScalarWedgePolynomial(const unsigned int degree) - : ScalarPolynomialsBase(degree, - compute_n_polynomials_wedge(dim, degree)) - , poly_tri(BarycentricPolynomials<2>::get_fe_p_basis(degree)) - , poly_line(BarycentricPolynomials<1>::get_fe_p_basis(degree)) - {} - - - namespace - { - /** - * Decompose the shape-function index of a linear wedge into an index - * to access the right shape function within the triangle and and within - * the line. - */ - static const constexpr ndarray wedge_table_1{ - {{{0, 0}}, {{1, 0}}, {{2, 0}}, {{0, 1}}, {{1, 1}}, {{2, 1}}}}; - - /** - * Decompose the shape-function index of a quadratic wedge into an index - * to access the right shape function within the triangle and and within - * the line. - */ - static const constexpr ndarray wedge_table_2{ - {{{0, 0}}, - {{1, 0}}, - {{2, 0}}, - {{0, 1}}, - {{1, 1}}, - {{2, 1}}, - {{3, 0}}, - {{4, 0}}, - {{5, 0}}, - {{3, 1}}, - {{4, 1}}, - {{5, 1}}, - {{0, 2}}, - {{1, 2}}, - {{2, 2}}, - {{3, 2}}, - {{4, 2}}, - {{5, 2}}}}; - } // namespace - - - - template - double - ScalarWedgePolynomial::compute_value(const unsigned int i, - const Point & p) const - { - const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i]; - - const Point<2> p_tri(p[0], p[1]); - const auto v_tri = poly_tri.compute_value(pair[0], p_tri); - - const Point<1> p_line(p[2]); - const auto v_line = poly_line.compute_value(pair[1], p_line); - - return v_tri * v_line; - } - - - - template - Tensor<1, dim> - ScalarWedgePolynomial::compute_grad(const unsigned int i, - const Point & p) const - { - const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i]; - - const Point<2> p_tri(p[0], p[1]); - const auto v_tri = poly_tri.compute_value(pair[0], p_tri); - const auto g_tri = poly_tri.compute_grad(pair[0], p_tri); - - const Point<1> p_line(p[2]); - const auto v_line = poly_line.compute_value(pair[1], p_line); - const auto g_line = poly_line.compute_grad(pair[1], p_line); - - Tensor<1, dim> grad; - grad[0] = g_tri[0] * v_line; - grad[1] = g_tri[1] * v_line; - grad[2] = v_tri * g_line[0]; - - return grad; - } - - - - template - Tensor<2, dim> - ScalarWedgePolynomial::compute_grad_grad(const unsigned int i, - const Point & p) const - { - (void)i; - (void)p; - - Assert(false, ExcNotImplemented()); - return Tensor<2, dim>(); - } - - - - template - void - ScalarWedgePolynomial::evaluate( - const Point & unit_point, - std::vector & values, - std::vector> &grads, - std::vector> &grad_grads, - std::vector> &third_derivatives, - std::vector> &fourth_derivatives) const - { - (void)grads; - (void)grad_grads; - (void)third_derivatives; - (void)fourth_derivatives; - - if (values.size() == this->n()) - for (unsigned int i = 0; i < this->n(); i++) - values[i] = compute_value(i, unit_point); - - if (grads.size() == this->n()) - for (unsigned int i = 0; i < this->n(); i++) - grads[i] = compute_grad(i, unit_point); - } - - - - template - Tensor<1, dim> - ScalarWedgePolynomial::compute_1st_derivative(const unsigned int i, - const Point & p) const - { - return compute_grad(i, p); - } - - - - template - Tensor<2, dim> - ScalarWedgePolynomial::compute_2nd_derivative(const unsigned int i, - const Point & p) const - { - (void)i; - (void)p; - - Assert(false, ExcNotImplemented()); - - return {}; - } - - - - template - Tensor<3, dim> - ScalarWedgePolynomial::compute_3rd_derivative(const unsigned int i, - const Point & p) const - { - (void)i; - (void)p; - - Assert(false, ExcNotImplemented()); - - return {}; - } - - - - template - Tensor<4, dim> - ScalarWedgePolynomial::compute_4th_derivative(const unsigned int i, - const Point & p) const - { - (void)i; - (void)p; - - Assert(false, ExcNotImplemented()); - - return {}; - } - - - - template - std::string - ScalarWedgePolynomial::name() const - { - return "ScalarWedgePolynomial"; - } - - - - template - std::unique_ptr> - ScalarWedgePolynomial::clone() const - { - return std::make_unique>(*this); - } - - - template ScalarPyramidPolynomial::ScalarPyramidPolynomial( const unsigned int degree) @@ -503,9 +281,6 @@ namespace Simplex - template class ScalarWedgePolynomial<1>; - template class ScalarWedgePolynomial<2>; - template class ScalarWedgePolynomial<3>; template class ScalarPyramidPolynomial<1>; template class ScalarPyramidPolynomial<2>; template class ScalarPyramidPolynomial<3>; diff --git a/source/base/polynomials_wedge.cc b/source/base/polynomials_wedge.cc new file mode 100644 index 0000000000..8ac20bb9de --- /dev/null +++ b/source/base/polynomials_wedge.cc @@ -0,0 +1,255 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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 + +DEAL_II_NAMESPACE_OPEN + +namespace Simplex +{ + namespace + { + unsigned int + compute_n_polynomials_wedge(const unsigned int dim, + const unsigned int degree) + { + if (dim == 3) + { + if (degree == 1) + return 6; + if (degree == 2) + return 18; + } + + Assert(false, ExcNotImplemented()); + + return 0; + } + } // namespace + + + + template + ScalarWedgePolynomial::ScalarWedgePolynomial(const unsigned int degree) + : ScalarPolynomialsBase(degree, + compute_n_polynomials_wedge(dim, degree)) + , poly_tri(BarycentricPolynomials<2>::get_fe_p_basis(degree)) + , poly_line(BarycentricPolynomials<1>::get_fe_p_basis(degree)) + {} + + + namespace + { + /** + * Decompose the shape-function index of a linear wedge into an index + * to access the right shape function within the triangle and and within + * the line. + */ + static const constexpr ndarray wedge_table_1{ + {{{0, 0}}, {{1, 0}}, {{2, 0}}, {{0, 1}}, {{1, 1}}, {{2, 1}}}}; + + /** + * Decompose the shape-function index of a quadratic wedge into an index + * to access the right shape function within the triangle and and within + * the line. + */ + static const constexpr ndarray wedge_table_2{ + {{{0, 0}}, + {{1, 0}}, + {{2, 0}}, + {{0, 1}}, + {{1, 1}}, + {{2, 1}}, + {{3, 0}}, + {{4, 0}}, + {{5, 0}}, + {{3, 1}}, + {{4, 1}}, + {{5, 1}}, + {{0, 2}}, + {{1, 2}}, + {{2, 2}}, + {{3, 2}}, + {{4, 2}}, + {{5, 2}}}}; + } // namespace + + + + template + double + ScalarWedgePolynomial::compute_value(const unsigned int i, + const Point & p) const + { + const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i]; + + const Point<2> p_tri(p[0], p[1]); + const auto v_tri = poly_tri.compute_value(pair[0], p_tri); + + const Point<1> p_line(p[2]); + const auto v_line = poly_line.compute_value(pair[1], p_line); + + return v_tri * v_line; + } + + + + template + Tensor<1, dim> + ScalarWedgePolynomial::compute_grad(const unsigned int i, + const Point & p) const + { + const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i]; + + const Point<2> p_tri(p[0], p[1]); + const auto v_tri = poly_tri.compute_value(pair[0], p_tri); + const auto g_tri = poly_tri.compute_grad(pair[0], p_tri); + + const Point<1> p_line(p[2]); + const auto v_line = poly_line.compute_value(pair[1], p_line); + const auto g_line = poly_line.compute_grad(pair[1], p_line); + + Tensor<1, dim> grad; + grad[0] = g_tri[0] * v_line; + grad[1] = g_tri[1] * v_line; + grad[2] = v_tri * g_line[0]; + + return grad; + } + + + + template + Tensor<2, dim> + ScalarWedgePolynomial::compute_grad_grad(const unsigned int i, + const Point & p) const + { + (void)i; + (void)p; + + Assert(false, ExcNotImplemented()); + return Tensor<2, dim>(); + } + + + + template + void + ScalarWedgePolynomial::evaluate( + const Point & unit_point, + std::vector & values, + std::vector> &grads, + std::vector> &grad_grads, + std::vector> &third_derivatives, + std::vector> &fourth_derivatives) const + { + (void)grads; + (void)grad_grads; + (void)third_derivatives; + (void)fourth_derivatives; + + if (values.size() == this->n()) + for (unsigned int i = 0; i < this->n(); i++) + values[i] = compute_value(i, unit_point); + + if (grads.size() == this->n()) + for (unsigned int i = 0; i < this->n(); i++) + grads[i] = compute_grad(i, unit_point); + } + + + + template + Tensor<1, dim> + ScalarWedgePolynomial::compute_1st_derivative(const unsigned int i, + const Point & p) const + { + return compute_grad(i, p); + } + + + + template + Tensor<2, dim> + ScalarWedgePolynomial::compute_2nd_derivative(const unsigned int i, + const Point & p) const + { + (void)i; + (void)p; + + Assert(false, ExcNotImplemented()); + + return {}; + } + + + + template + Tensor<3, dim> + ScalarWedgePolynomial::compute_3rd_derivative(const unsigned int i, + const Point & p) const + { + (void)i; + (void)p; + + Assert(false, ExcNotImplemented()); + + return {}; + } + + + + template + Tensor<4, dim> + ScalarWedgePolynomial::compute_4th_derivative(const unsigned int i, + const Point & p) const + { + (void)i; + (void)p; + + Assert(false, ExcNotImplemented()); + + return {}; + } + + + + template + std::string + ScalarWedgePolynomial::name() const + { + return "ScalarWedgePolynomial"; + } + + + + template + std::unique_ptr> + ScalarWedgePolynomial::clone() const + { + return std::make_unique>(*this); + } + + + + template class ScalarWedgePolynomial<1>; + template class ScalarWedgePolynomial<2>; + template class ScalarWedgePolynomial<3>; +} // namespace Simplex + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/base/qprojector.cc b/source/base/qprojector.cc index 0a310707ce..f64d67cc5a 100644 --- a/source/base/qprojector.cc +++ b/source/base/qprojector.cc @@ -15,11 +15,10 @@ #include #include +#include #include #include -#include - DEAL_II_NAMESPACE_OPEN diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc index 9be07c25bd..f84f13fb6c 100644 --- a/source/base/quadrature_lib.cc +++ b/source/base/quadrature_lib.cc @@ -1350,6 +1350,449 @@ QSplit::QSplit(const QSimplex &base, const Point &split_point) +namespace Simplex +{ + template + QGauss::QGauss(const unsigned int n_points_1D) + : QSimplex(Quadrature()) + { + // fill quadrature points and quadrature weights + if (dim == 1) + { + const dealii::QGauss quad(n_points_1D); + + this->quadrature_points = quad.get_points(); + this->weights = quad.get_weights(); + } + else if (dim == 2) + { + if (n_points_1D == 1) + { + const double p = 1.0 / 3.0; + this->quadrature_points.emplace_back(p, p); + this->weights.emplace_back(0.5); + } + else if (n_points_1D == 2) + { + const double Q23 = 2.0 / 3.0; + const double Q16 = 1.0 / 6.0; + + this->quadrature_points.emplace_back(Q23, Q16); + this->quadrature_points.emplace_back(Q16, Q23); + this->quadrature_points.emplace_back(Q16, Q16); + this->weights.emplace_back(Q16); + this->weights.emplace_back(Q16); + this->weights.emplace_back(Q16); + } + else if (n_points_1D == 3) + { + const double q12 = 0.5; + + // clang-format off + this->quadrature_points.emplace_back(0.3333333333330, 0.3333333333330); + this->quadrature_points.emplace_back(0.7974269853530, 0.1012865073230); + this->quadrature_points.emplace_back(0.1012865073230, 0.7974269853530); + this->quadrature_points.emplace_back(0.1012865073230, 0.1012865073230); + this->quadrature_points.emplace_back(0.0597158717898, 0.4701420641050); + this->quadrature_points.emplace_back(0.4701420641050, 0.0597158717898); + this->quadrature_points.emplace_back(0.4701420641050, 0.4701420641050); + // clang-format on + + this->weights.emplace_back(q12 * 0.225); + this->weights.emplace_back(q12 * 0.125939180545); + this->weights.emplace_back(q12 * 0.125939180545); + this->weights.emplace_back(q12 * 0.125939180545); + this->weights.emplace_back(q12 * 0.132394152789); + this->weights.emplace_back(q12 * 0.132394152789); + this->weights.emplace_back(q12 * 0.132394152789); + } + else if (n_points_1D == 4) + { + Quadrature::operator=(QWitherdenVincent(n_points_1D)); + } + } + else if (dim == 3) + { + if (n_points_1D == 1) + { + const double Q14 = 1.0 / 4.0; + const double Q16 = 1.0 / 6.0; + + this->quadrature_points.emplace_back(Q14, Q14, Q14); + this->weights.emplace_back(Q16); + } + else if (n_points_1D == 2) + { + const double Q124 = 1.0 / 6.0 / 4.0; + + const double palpha = (5.0 + 3.0 * sqrt(5.0)) / 20.0; + const double pbeta = (5.0 - sqrt(5.0)) / 20.0; + this->quadrature_points.emplace_back(pbeta, pbeta, pbeta); + this->quadrature_points.emplace_back(palpha, pbeta, pbeta); + this->quadrature_points.emplace_back(pbeta, palpha, pbeta); + this->quadrature_points.emplace_back(pbeta, pbeta, palpha); + this->weights.emplace_back(Q124); + this->weights.emplace_back(Q124); + this->weights.emplace_back(Q124); + this->weights.emplace_back(Q124); + } + else if (n_points_1D == 3) + { + const double Q16 = 1.0 / 6.0; + + // clang-format off + this->quadrature_points.emplace_back(0.5684305841968444, 0.1438564719343852, 0.1438564719343852); + this->quadrature_points.emplace_back(0.1438564719343852, 0.1438564719343852, 0.1438564719343852); + this->quadrature_points.emplace_back(0.1438564719343852, 0.1438564719343852, 0.5684305841968444); + this->quadrature_points.emplace_back(0.1438564719343852, 0.5684305841968444, 0.1438564719343852); + this->quadrature_points.emplace_back(0.0000000000000000, 0.5000000000000000, 0.5000000000000000); + this->quadrature_points.emplace_back(0.5000000000000000, 0.0000000000000000, 0.5000000000000000); + this->quadrature_points.emplace_back(0.5000000000000000, 0.5000000000000000, 0.0000000000000000); + this->quadrature_points.emplace_back(0.5000000000000000, 0.0000000000000000, 0.0000000000000000); + this->quadrature_points.emplace_back(0.0000000000000000, 0.5000000000000000, 0.0000000000000000); + this->quadrature_points.emplace_back(0.0000000000000000, 0.0000000000000000, 0.5000000000000000); + // clang-format on + + this->weights.emplace_back(0.2177650698804054 * Q16); + this->weights.emplace_back(0.2177650698804054 * Q16); + this->weights.emplace_back(0.2177650698804054 * Q16); + this->weights.emplace_back(0.2177650698804054 * Q16); + this->weights.emplace_back(0.0214899534130631 * Q16); + this->weights.emplace_back(0.0214899534130631 * Q16); + this->weights.emplace_back(0.0214899534130631 * Q16); + this->weights.emplace_back(0.0214899534130631 * Q16); + this->weights.emplace_back(0.0214899534130631 * Q16); + this->weights.emplace_back(0.0214899534130631 * Q16); + } + else if (n_points_1D == 4) + { + Quadrature::operator=(QWitherdenVincent(n_points_1D)); + } + } + + AssertDimension(this->quadrature_points.size(), this->weights.size()); + Assert(this->quadrature_points.size() > 0, + ExcNotImplemented( + "Simplex::QGauss is currently only implemented for " + "n_points_1D = 1, 2, 3, and 4 while you are asking for " + "n_points_1D = " + + Utilities::to_string(n_points_1D))); + } + + namespace + { + template + std::vector> + all_permutations(const std::array &b_point) + { + std::vector> output; + + // We want all possible permutations of the barycentric coordinates. + // The easiest way to get all of them is to sort the input first and + // then use next_permutation to cycle through them all. + std::array temp = b_point; + std::sort(temp.begin(), temp.end()); + do + { + output.push_back(temp); + } + while (std::next_permutation(temp.begin(), temp.end())); + + return output; + } + } // namespace + + + + template + QWitherdenVincent::QWitherdenVincent(const unsigned int n_points_1D) + : QSimplex(Quadrature()) + { + Assert(1 <= dim && dim <= 3, ExcNotImplemented()); + // Just use Gauss in 1D: this is a high-order open rule so this is a + // reasonable equivalent for generic programming. + if (dim == 1) + { + Quadrature::operator=(dealii::QGauss(n_points_1D)); + return; + } + + std::array centroid; + std::fill(centroid.begin(), centroid.end(), 1.0 / (dim + 1.0)); + std::vector>> b_point_permutations; + std::vector b_weights; + + // We can simplify the implementation of these quadrature rules + // by quite a bit by exploiting symmetry - we do essentially the + // same thing for each barycentric coordinate, so we can express + // our quadrature rule as permutations of barycentric points + // instead of writing things out explicitly. + + // Apply a Barycentric permutation where one point is different. + auto process_point_1 = [&](const double a, const double w) { + const double b = 1.0 - dim * a; + std::array b_point; + std::fill(b_point.begin(), b_point.begin() + dim, a); + b_point[dim] = b; + + b_weights.push_back(w); + b_point_permutations.push_back(all_permutations(b_point)); + }; + + // Apply a Barycentric permutation where two points (in 3D) are different. + auto process_point_2 = [&](const double a, const double w) { + Assert(dim == 3, ExcInternalError()); + const double b = (1.0 - 2.0 * a) / 2.0; + std::array b_point; + std::fill(b_point.begin(), b_point.begin() + dim - 1, a); + b_point[dim - 1] = b; + b_point[dim] = b; + + b_weights.push_back(w); + b_point_permutations.push_back(all_permutations(b_point)); + }; + + // Apply a Barycentric permutation where three (or four) points + // are different (since there are two inputs). + auto process_point_3 = [&](const double a, const double b, const double w) { + const double c = 1.0 - (dim - 1.0) * a - b; + std::array b_point; + std::fill(b_point.begin(), b_point.begin() + dim - 1, a); + b_point[dim - 1] = b; + b_point[dim] = c; + + b_weights.push_back(w); + b_point_permutations.push_back(all_permutations(b_point)); + }; + + if (n_points_1D == 1) + { + b_point_permutations.push_back({centroid}); + b_weights.push_back(1.0); + } + else if (n_points_1D == 2) + { + // This is WV-4 in 2D and WV-3 in 3D + if (dim == 2) + { + process_point_1(9.1576213509770743e-02, 1.0995174365532187e-01); + process_point_1(4.4594849091596489e-01, 2.2338158967801147e-01); + } + else if (dim == 3) + { + process_point_1(3.281633025163817e-01, 1.362178425370874e-01); + process_point_1(1.080472498984286e-01, 1.137821574629126e-01); + } + } + else if (n_points_1D == 3) + { + // This is the WV-5 rule in both 2D and 3D + if (dim == 2) + { + b_weights.push_back(0.225); + b_point_permutations.push_back({centroid}); + + process_point_1(1.0128650732345634e-01, 1.2593918054482714e-01); + process_point_1(4.7014206410511511e-01, 1.3239415278850619e-01); + } + else if (dim == 3) + { + process_point_1(3.108859192633006e-01, 1.126879257180159e-01); + process_point_1(9.273525031089125e-02, 7.349304311636196e-02); + + process_point_2(4.550370412564964e-02, 4.254602077708147e-02); + } + } + else if (n_points_1D == 4) + { + // This is the WV-7 rule in both 2D and 3D + if (dim == 2) + { + process_point_1(3.3730648554587850e-02, 1.6545050110792131e-02); + process_point_1(4.7430969250471822e-01, 7.7086646185986069e-02); + process_point_1(2.4157738259540357e-01, 1.2794417123015558e-01); + process_point_3(4.7036644652595216e-02, + 1.9868331479735168e-01, + 5.5878732903199779e-02); + } + else if (dim == 3) + { + b_point_permutations.push_back({centroid}); + b_weights.push_back(9.548528946413085e-02); + + process_point_1(3.157011497782028e-01, 4.232958120996703e-02); + process_point_2(5.048982259839635e-02, 3.189692783285758e-02); + + process_point_3(1.888338310260010e-01, + 5.751716375870000e-01, + 3.720713072833462e-02); + process_point_3(2.126547254148314e-02, + 8.108302410985486e-01, + 8.110770829903342e-03); + } + } + else if (n_points_1D == 5) + { + // This is the WV-9 rule in both 2D and 3D + if (dim == 2) + { + b_point_permutations.push_back({centroid}); + b_weights.push_back(9.7135796282798836e-02); + + process_point_1(4.4729513394452691e-02, 2.5577675658698031e-02); + process_point_1(4.8968251919873762e-01, 3.1334700227139071e-02); + process_point_1(4.3708959149293664e-01, 7.7827541004774278e-02); + process_point_1(1.8820353561903275e-01, 7.9647738927210249e-02); + + process_point_3(3.6838412054736258e-02, + 2.2196298916076568e-01, + 4.3283539377289376e-02); + } + else if (dim == 3) + { + b_point_permutations.push_back({centroid}); + b_weights.push_back(5.801054891248025e-02); + + process_point_1(6.198169755222693e-10, 6.431928175925639e-05); + process_point_1(1.607745353952616e-01, 2.317333846242546e-02); + process_point_1(3.222765218214210e-01, 2.956291233542929e-02); + process_point_1(4.510891834541358e-02, 8.063979979616182e-03); + + process_point_2(1.122965460043761e-01, 3.813408010370246e-02); + + process_point_3(4.588714487524592e-01, + 2.554579233041310e-03, + 8.384422198298552e-03); + process_point_3(3.377587068533860e-02, + 7.183503264420745e-01, + 1.023455935274533e-02); + process_point_3(1.836413698099279e-01, + 3.441591057817528e-02, + 2.052491596798814e-02); + } + } + else if (n_points_1D == 6) + { + // There is no WV-11 rule in 3D yet + if (dim == 2) + { + b_point_permutations.push_back({centroid}); + b_weights.push_back(8.5761179732224219e-02); + + process_point_1(2.8485417614371900e-02, 1.0431870512894697e-02); + process_point_1(4.9589190096589092e-01, 1.6606273054585369e-02); + process_point_1(1.0263548271224643e-01, 3.8630759237019321e-02); + process_point_1(4.3846592676435220e-01, 6.7316154079468296e-02); + process_point_1(2.1021995670317828e-01, 7.0515684111716576e-02); + + process_point_3(7.3254276860644785e-03, + 1.4932478865208237e-01, + 1.0290289572953278e-02); + process_point_3(4.6010500165429957e-02, + 2.8958112563770588e-01, + 4.0332476640500554e-02); + } + else if (dim == 3) + { + Assert(false, ExcNotImplemented()); + } + } + else + { + Assert(false, ExcNotImplemented()); + } + + Assert(b_point_permutations.size() == b_weights.size(), ExcInternalError()); + for (unsigned int permutation_n = 0; permutation_n < b_weights.size(); + ++permutation_n) + { + for (const std::array &b_point : + b_point_permutations[permutation_n]) + { + const double volume = (dim == 2 ? 1.0 / 2.0 : 1.0 / 6.0); + this->weights.emplace_back(volume * b_weights[permutation_n]); + Point c_point; + std::copy(b_point.begin(), + b_point.begin() + dim, + c_point.begin_raw()); + this->quadrature_points.emplace_back(c_point); + } + } + } + + + + template + QGaussWedge::QGaussWedge(const unsigned int n_points) + : Quadrature() + { + AssertDimension(dim, 3); + + Simplex::QGauss<2> quad_tri(n_points); + QGauss<1> quad_line(n_points); + + for (unsigned int i = 0; i < quad_line.size(); ++i) + for (unsigned int j = 0; j < quad_tri.size(); ++j) + { + this->quadrature_points.emplace_back(quad_tri.point(j)[0], + quad_tri.point(j)[1], + quad_line.point(i)[0]); + this->weights.emplace_back(quad_tri.weight(j) * quad_line.weight(i)); + } + + AssertDimension(this->quadrature_points.size(), this->weights.size()); + Assert(this->quadrature_points.size() > 0, + ExcMessage("No valid quadrature points!")); + } + + + + template + QGaussPyramid::QGaussPyramid(const unsigned int n_points_1D) + : Quadrature() + { + AssertDimension(dim, 3); + + if (n_points_1D == 1) + { + const double Q14 = 1.0 / 4.0; + const double Q43 = 4.0 / 3.0; + + this->quadrature_points.emplace_back(0, 0, Q14); + this->weights.emplace_back(Q43); + } + else if (n_points_1D == 2) + { + // clang-format off + this->quadrature_points.emplace_back(-0.26318405556971, -0.26318405556971, 0.54415184401122); + this->quadrature_points.emplace_back(-0.50661630334979, -0.50661630334979, 0.12251482265544); + this->quadrature_points.emplace_back(-0.26318405556971, +0.26318405556971, 0.54415184401122); + this->quadrature_points.emplace_back(-0.50661630334979, +0.50661630334979, 0.12251482265544); + this->quadrature_points.emplace_back(+0.26318405556971, -0.26318405556971, 0.54415184401122); + this->quadrature_points.emplace_back(+0.50661630334979, -0.50661630334979, 0.12251482265544); + this->quadrature_points.emplace_back(+0.26318405556971, +0.26318405556971, 0.54415184401122); + this->quadrature_points.emplace_back(+0.50661630334979, +0.50661630334979, 0.12251482265544); + // clang-format on + + this->weights.emplace_back(0.10078588207983); + this->weights.emplace_back(0.23254745125351); + this->weights.emplace_back(0.10078588207983); + this->weights.emplace_back(0.23254745125351); + this->weights.emplace_back(0.10078588207983); + this->weights.emplace_back(0.23254745125351); + this->weights.emplace_back(0.10078588207983); + this->weights.emplace_back(0.23254745125351); + } + + AssertDimension(this->quadrature_points.size(), this->weights.size()); + Assert(this->quadrature_points.size() > 0, + ExcMessage("No valid quadrature points!")); + } + +} // namespace Simplex + + + // explicit specialization // note that 1d formulae are specialized by implementation above template class QGauss<2>; @@ -1396,4 +1839,18 @@ template class QSplit<1>; template class QSplit<2>; template class QSplit<3>; +template class Simplex::QGauss<1>; +template class Simplex::QGauss<2>; +template class Simplex::QGauss<3>; +template class Simplex::QGaussWedge<1>; +template class Simplex::QGaussWedge<2>; +template class Simplex::QGaussWedge<3>; +template class Simplex::QGaussPyramid<1>; +template class Simplex::QGaussPyramid<2>; +template class Simplex::QGaussPyramid<3>; + +template class Simplex::QWitherdenVincent<1>; +template class Simplex::QWitherdenVincent<2>; +template class Simplex::QWitherdenVincent<3>; + DEAL_II_NAMESPACE_CLOSE diff --git a/source/fe/CMakeLists.txt b/source/fe/CMakeLists.txt index 666999b0a5..9d9fc09703 100644 --- a/source/fe/CMakeLists.txt +++ b/source/fe/CMakeLists.txt @@ -29,12 +29,14 @@ SET(_unity_include_src fe_dgp_nonparametric.cc fe_dgq.cc fe_dg_vector.cc + fe_enriched.cc fe_face.cc fe_nedelec.cc fe_nedelec_sz.cc fe_nothing.cc fe_poly.cc fe_poly_tensor.cc + fe_pyramid_p.cc fe_p1nc.cc fe_q_base.cc fe_q.cc @@ -50,8 +52,10 @@ SET(_unity_include_src fe_series_fourier.cc fe_series_legendre.cc fe_system.cc - fe_enriched.cc + fe_simplex_p.cc + fe_simplex_p_bubbles.cc fe_trace.cc + fe_wedge_p.cc mapping_c1.cc mapping_cartesian.cc mapping.cc @@ -105,6 +109,7 @@ SET(_inst fe_nothing.inst.in fe_poly.inst.in fe_poly_tensor.inst.in + fe_pyramid_p.inst.in fe_q_base.inst.in fe_q_bubbles.inst.in fe_q_dg0.inst.in @@ -117,6 +122,8 @@ SET(_inst fe_rt_bubbles.inst.in fe_series_fourier.inst.in fe_series_legendre.inst.in + fe_simplex_p.inst.in + fe_simplex_p_bubbles.inst.in fe_system.inst.in fe_enriched.inst.in fe_tools.inst.in @@ -128,6 +135,7 @@ SET(_inst fe_values.impl.1.inst.in fe_values.impl.2.inst.in fe_values.inst.in + fe_wedge_p.inst.in mapping_c1.inst.in mapping_cartesian.inst.in mapping.inst.in diff --git a/source/fe/fe_dgq.cc b/source/fe/fe_dgq.cc index f75719806e..73acc4e8d5 100644 --- a/source/fe/fe_dgq.cc +++ b/source/fe/fe_dgq.cc @@ -21,17 +21,19 @@ #include #include #include +#include #include #include #include #include #include +#include +#include #include +#include #include -#include - #include #include #include diff --git a/source/fe/fe_pyramid_p.cc b/source/fe/fe_pyramid_p.cc new file mode 100644 index 0000000000..fe286a52ed --- /dev/null +++ b/source/fe/fe_pyramid_p.cc @@ -0,0 +1,315 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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 + +DEAL_II_NAMESPACE_OPEN + +namespace Simplex +{ + namespace + { + /** + * Helper function to set up the dpo vector of FE_PyramidP for a given @p degree. + */ + internal::GenericDoFsPerObject + get_dpo_vector_fe_pyramid_p(const unsigned int degree) + { + internal::GenericDoFsPerObject dpo; + + if (degree == 1) + { + dpo.dofs_per_object_exclusive = {{1}, {0}, {0, 0, 0, 0, 0}, {0}}; + dpo.dofs_per_object_inclusive = {{1}, {2}, {4, 3, 3, 3, 3}, {5}}; + dpo.object_index = {{}, {5}, {5}, {5}}; + dpo.first_object_index_on_face = {{}, + {4, 3, 3, 3, 3}, + {4, 3, 3, 3, 3}}; + } + else + { + Assert(false, ExcNotImplemented()); + } + + return dpo; + } + + /** + * Helper function to set up the dpo vector of FE_PyramidDGP for a given @p degree. + */ + internal::GenericDoFsPerObject + get_dpo_vector_fe_pyramid_dgp(const unsigned int degree) + { + unsigned int n_dofs = 0; + + if (degree == 1) + n_dofs = 5; + else + Assert(false, ExcNotImplemented()); + + return internal::expand(3, {{0, 0, 0, n_dofs}}, ReferenceCells::Pyramid); + } + } // namespace + + + template + FE_Pyramid::FE_Pyramid( + const unsigned int degree, + const internal::GenericDoFsPerObject & dpos, + const typename FiniteElementData::Conformity conformity) + : dealii::FE_Poly( + Simplex::ScalarPyramidPolynomial(degree), + FiniteElementData(dpos, + ReferenceCells::Pyramid, + 1, + degree, + conformity), + std::vector( + FiniteElementData(dpos, ReferenceCells::Pyramid, 1, degree) + .dofs_per_cell, + true), + std::vector( + FiniteElementData(dpos, ReferenceCells::Pyramid, 1, degree) + .dofs_per_cell, + std::vector(1, true))) + { + AssertDimension(dim, 3); + + + if (degree == 1) + { + this->unit_support_points.emplace_back(-1.0, -1.0, 0.0); + this->unit_support_points.emplace_back(+1.0, -1.0, 0.0); + this->unit_support_points.emplace_back(-1.0, +1.0, 0.0); + this->unit_support_points.emplace_back(+1.0, +1.0, 0.0); + this->unit_support_points.emplace_back(+0.0, +0.0, 1.0); + } + } + + + + template + FE_PyramidP::FE_PyramidP(const unsigned int degree) + : FE_Pyramid(degree, + get_dpo_vector_fe_pyramid_p(degree), + FiniteElementData::H1) + {} + + + + template + std::unique_ptr> + FE_PyramidP::clone() const + { + return std::make_unique>(*this); + } + + + + template + std::string + FE_PyramidP::get_name() const + { + std::ostringstream namebuf; + namebuf << "FE_PyramidP<" << dim << ">(" << this->degree << ")"; + + return namebuf.str(); + } + + + + template + FiniteElementDomination::Domination + FE_PyramidP::compare_for_domination( + const FiniteElement &fe_other, + const unsigned int codim) const + { + Assert(codim <= dim, ExcImpossibleInDim(dim)); + + // vertex/line/face domination + // (if fe_other is derived from FE_DGP) + // ------------------------------------ + if (codim > 0) + if (dynamic_cast *>(&fe_other) != nullptr) + // there are no requirements between continuous and discontinuous + // elements + return FiniteElementDomination::no_requirements; + + // vertex/line/face domination + // (if fe_other is not derived from FE_DGP) + // & cell domination + // ---------------------------------------- + if (const FE_PyramidP *fe_pp_other = + dynamic_cast *>(&fe_other)) + { + if (this->degree < fe_pp_other->degree) + return FiniteElementDomination::this_element_dominates; + else if (this->degree == fe_pp_other->degree) + return FiniteElementDomination::either_element_can_dominate; + else + return FiniteElementDomination::other_element_dominates; + } + else if (const FE_P *fe_p_other = + dynamic_cast *>(&fe_other)) + { + if (this->degree < fe_p_other->degree) + return FiniteElementDomination::this_element_dominates; + else if (this->degree == fe_p_other->degree) + return FiniteElementDomination::either_element_can_dominate; + else + return FiniteElementDomination::other_element_dominates; + } + else if (const FE_Q *fe_q_other = + dynamic_cast *>(&fe_other)) + { + if (this->degree < fe_q_other->degree) + return FiniteElementDomination::this_element_dominates; + else if (this->degree == fe_q_other->degree) + return FiniteElementDomination::either_element_can_dominate; + else + return FiniteElementDomination::other_element_dominates; + } + else if (const FE_Nothing *fe_nothing = + dynamic_cast *>(&fe_other)) + { + if (fe_nothing->is_dominating()) + return FiniteElementDomination::other_element_dominates; + else + // the FE_Nothing has no degrees of freedom and it is typically used + // in a context where we don't require any continuity along the + // interface + return FiniteElementDomination::no_requirements; + } + + Assert(false, ExcNotImplemented()); + return FiniteElementDomination::neither_element_dominates; + } + + + + template + std::vector> + FE_PyramidP::hp_vertex_dof_identities( + const FiniteElement &fe_other) const + { + (void)fe_other; + + Assert((dynamic_cast *>(&fe_other)) || + (dynamic_cast *>(&fe_other)), + ExcNotImplemented()); + + return {{0, 0}}; + } + + + + template + std::vector> + FE_PyramidP::hp_line_dof_identities( + const FiniteElement &fe_other) const + { + (void)fe_other; + + Assert((dynamic_cast *>(&fe_other)) || + (dynamic_cast *>(&fe_other)), + ExcNotImplemented()); + + std::vector> result; + + for (unsigned int i = 0; i < this->degree - 1; ++i) + result.emplace_back(i, i); + + return result; + } + + + + template + std::vector> + FE_PyramidP::hp_quad_dof_identities( + const FiniteElement &fe_other, + const unsigned int face_no) const + { + (void)fe_other; + + + AssertIndexRange(face_no, 5); + + if (face_no == 0) + { + Assert((dynamic_cast *>(&fe_other)), + ExcNotImplemented()); + } + else + { + Assert((dynamic_cast *>(&fe_other)), + ExcNotImplemented()); + } + + std::vector> result; + + for (unsigned int i = 0; i < this->n_dofs_per_quad(face_no); ++i) + result.emplace_back(i, i); + + return result; + } + + + + template + FE_PyramidDGP::FE_PyramidDGP(const unsigned int degree) + : FE_Pyramid(degree, + get_dpo_vector_fe_pyramid_dgp(degree), + FiniteElementData::L2) + {} + + + + template + std::unique_ptr> + FE_PyramidDGP::clone() const + { + return std::make_unique>(*this); + } + + + + template + std::string + FE_PyramidDGP::get_name() const + { + std::ostringstream namebuf; + namebuf << "FE_PyramidDGP<" << dim << ">(" << this->degree << ")"; + + return namebuf.str(); + } +} // namespace Simplex + +// explicit instantiations +#include "fe_pyramid_p.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/simplex/fe_lib.inst.in b/source/fe/fe_pyramid_p.inst.in similarity index 61% rename from source/simplex/fe_lib.inst.in rename to source/fe/fe_pyramid_p.inst.in index f40343af26..cce82dcff2 100644 --- a/source/simplex/fe_lib.inst.in +++ b/source/fe/fe_pyramid_p.inst.in @@ -18,22 +18,11 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS) { #if deal_II_dimension <= deal_II_space_dimension - template class Simplex::FE_Poly; - template class Simplex::FE_P; - template class Simplex::FE_DGP; - template class Simplex::FE_Wedge; - template class Simplex::FE_WedgeP; - template class Simplex::FE_WedgeDGP; template class Simplex::FE_Pyramid; template class Simplex::FE_PyramidP; template class Simplex::FE_PyramidDGP; - template class Simplex::FE_P_Bubbles; #endif } diff --git a/source/fe/fe_q.cc b/source/fe/fe_q.cc index 9be753fa61..39db78bb5c 100644 --- a/source/fe/fe_q.cc +++ b/source/fe/fe_q.cc @@ -18,12 +18,14 @@ #include #include +#include #include +#include +#include +#include #include -#include - #include #include #include diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc index 64eff9f429..c663f75efe 100644 --- a/source/fe/fe_q_base.cc +++ b/source/fe/fe_q_base.cc @@ -26,10 +26,12 @@ #include #include #include +#include #include +#include +#include #include - -#include +#include #include #include diff --git a/source/simplex/fe_lib.cc b/source/fe/fe_simplex_p.cc similarity index 52% rename from source/simplex/fe_lib.cc rename to source/fe/fe_simplex_p.cc index fd51f04024..9adbe33f9e 100644 --- a/source/simplex/fe_lib.cc +++ b/source/fe/fe_simplex_p.cc @@ -15,16 +15,15 @@ #include +#include #include #include #include #include +#include #include -#include -#include - DEAL_II_NAMESPACE_OPEN namespace Simplex @@ -260,99 +259,6 @@ namespace Simplex return dpo; } - - /** - * Helper function to set up the dpo vector of FE_WedgeP for a given @p degree. - */ - internal::GenericDoFsPerObject - get_dpo_vector_fe_wedge_p(const unsigned int degree) - { - internal::GenericDoFsPerObject dpo; - - if (degree == 1) - { - dpo.dofs_per_object_exclusive = {{1}, {0}, {0, 0, 0, 0, 0}, {0}}; - dpo.dofs_per_object_inclusive = {{1}, {2}, {3, 3, 4, 4, 4}, {6}}; - dpo.object_index = {{}, {6}, {6}, {6}}; - dpo.first_object_index_on_face = {{}, - {3, 3, 4, 4, 4}, - {3, 3, 4, 4, 4}}; - } - else if (degree == 2) - { - dpo.dofs_per_object_exclusive = {{1}, {1}, {0, 0, 1, 1, 1}, {0}}; - dpo.dofs_per_object_inclusive = {{1}, {3}, {6, 6, 9, 9, 9}, {18}}; - dpo.object_index = {{}, {6}, {15, 15, 15, 16, 17}, {18}}; - dpo.first_object_index_on_face = {{}, - {3, 3, 4, 4, 4}, - {6, 6, 8, 8, 8}}; - } - else - { - Assert(false, ExcNotImplemented()); - } - - return dpo; - } - - /** - * Helper function to set up the dpo vector of FE_WedgeDGP for a given @p degree. - */ - internal::GenericDoFsPerObject - get_dpo_vector_fe_wedge_dgp(const unsigned int degree) - { - unsigned int n_dofs = 0; - - if (degree == 1) - n_dofs = 6; - else if (degree == 2) - n_dofs = 18; - else - Assert(false, ExcNotImplemented()); - - return internal::expand(3, {{0, 0, 0, n_dofs}}, ReferenceCells::Wedge); - } - - /** - * Helper function to set up the dpo vector of FE_PyramidP for a given @p degree. - */ - internal::GenericDoFsPerObject - get_dpo_vector_fe_pyramid_p(const unsigned int degree) - { - internal::GenericDoFsPerObject dpo; - - if (degree == 1) - { - dpo.dofs_per_object_exclusive = {{1}, {0}, {0, 0, 0, 0, 0}, {0}}; - dpo.dofs_per_object_inclusive = {{1}, {2}, {4, 3, 3, 3, 3}, {5}}; - dpo.object_index = {{}, {5}, {5}, {5}}; - dpo.first_object_index_on_face = {{}, - {4, 3, 3, 3, 3}, - {4, 3, 3, 3, 3}}; - } - else - { - Assert(false, ExcNotImplemented()); - } - - return dpo; - } - - /** - * Helper function to set up the dpo vector of FE_PyramidDGP for a given @p degree. - */ - internal::GenericDoFsPerObject - get_dpo_vector_fe_pyramid_dgp(const unsigned int degree) - { - unsigned int n_dofs = 0; - - if (degree == 1) - n_dofs = 5; - else - Assert(false, ExcNotImplemented()); - - return internal::expand(3, {{0, 0, 0, n_dofs}}, ReferenceCells::Pyramid); - } } // namespace @@ -949,747 +855,9 @@ namespace Simplex return {}; } - - - - template - FE_Wedge::FE_Wedge( - const unsigned int degree, - const internal::GenericDoFsPerObject & dpos, - const typename FiniteElementData::Conformity conformity) - : dealii::FE_Poly( - Simplex::ScalarWedgePolynomial(degree), - FiniteElementData(dpos, - ReferenceCells::Wedge, - 1, - degree, - conformity), - std::vector( - FiniteElementData(dpos, ReferenceCells::Wedge, 1, degree) - .dofs_per_cell, - true), - std::vector( - FiniteElementData(dpos, ReferenceCells::Wedge, 1, degree) - .dofs_per_cell, - std::vector(1, true))) - { - AssertDimension(dim, 3); - - if (degree == 1) - { - this->unit_support_points.emplace_back(0.0, 0.0, 0.0); - this->unit_support_points.emplace_back(1.0, 0.0, 0.0); - this->unit_support_points.emplace_back(0.0, 1.0, 0.0); - this->unit_support_points.emplace_back(0.0, 0.0, 1.0); - this->unit_support_points.emplace_back(1.0, 0.0, 1.0); - this->unit_support_points.emplace_back(0.0, 1.0, 1.0); - } - } - - - - template - FE_WedgeP::FE_WedgeP(const unsigned int degree) - : FE_Wedge(degree, - get_dpo_vector_fe_wedge_p(degree), - FiniteElementData::H1) - {} - - - - template - std::unique_ptr> - FE_WedgeP::clone() const - { - return std::make_unique>(*this); - } - - - - template - std::string - FE_WedgeP::get_name() const - { - std::ostringstream namebuf; - namebuf << "FE_WedgeP<" << dim << ">(" << this->degree << ")"; - - return namebuf.str(); - } - - - - template - FiniteElementDomination::Domination - FE_WedgeP::compare_for_domination( - const FiniteElement &fe_other, - const unsigned int codim) const - { - Assert(codim <= dim, ExcImpossibleInDim(dim)); - - // vertex/line/face domination - // (if fe_other is derived from FE_DGP) - // ------------------------------------ - if (codim > 0) - if (dynamic_cast *>(&fe_other) != nullptr) - // there are no requirements between continuous and discontinuous - // elements - return FiniteElementDomination::no_requirements; - - - // vertex/line/face domination - // (if fe_other is not derived from FE_DGP) - // & cell domination - // ---------------------------------------- - if (const FE_WedgeP *fe_wp_other = - dynamic_cast *>(&fe_other)) - { - if (this->degree < fe_wp_other->degree) - return FiniteElementDomination::this_element_dominates; - else if (this->degree == fe_wp_other->degree) - return FiniteElementDomination::either_element_can_dominate; - else - return FiniteElementDomination::other_element_dominates; - } - else if (const FE_P *fe_p_other = - dynamic_cast *>(&fe_other)) - { - if (this->degree < fe_p_other->degree) - return FiniteElementDomination::this_element_dominates; - else if (this->degree == fe_p_other->degree) - return FiniteElementDomination::either_element_can_dominate; - else - return FiniteElementDomination::other_element_dominates; - } - else if (const FE_Q *fe_q_other = - dynamic_cast *>(&fe_other)) - { - if (this->degree < fe_q_other->degree) - return FiniteElementDomination::this_element_dominates; - else if (this->degree == fe_q_other->degree) - return FiniteElementDomination::either_element_can_dominate; - else - return FiniteElementDomination::other_element_dominates; - } - else if (const FE_Nothing *fe_nothing = - dynamic_cast *>(&fe_other)) - { - if (fe_nothing->is_dominating()) - return FiniteElementDomination::other_element_dominates; - else - // the FE_Nothing has no degrees of freedom and it is typically used - // in a context where we don't require any continuity along the - // interface - return FiniteElementDomination::no_requirements; - } - - Assert(false, ExcNotImplemented()); - return FiniteElementDomination::neither_element_dominates; - } - - - - template - std::vector> - FE_WedgeP::hp_vertex_dof_identities( - const FiniteElement &fe_other) const - { - (void)fe_other; - - Assert((dynamic_cast *>(&fe_other)) || - (dynamic_cast *>(&fe_other)), - ExcNotImplemented()); - - return {{0, 0}}; - } - - - - template - std::vector> - FE_WedgeP::hp_line_dof_identities( - const FiniteElement &fe_other) const - { - (void)fe_other; - - Assert((dynamic_cast *>(&fe_other)) || - (dynamic_cast *>(&fe_other)), - ExcNotImplemented()); - - std::vector> result; - - for (unsigned int i = 0; i < this->degree - 1; ++i) - result.emplace_back(i, i); - - return result; - } - - - - template - std::vector> - FE_WedgeP::hp_quad_dof_identities( - const FiniteElement &fe_other, - const unsigned int face_no) const - { - (void)fe_other; - - AssertIndexRange(face_no, 5); - - if (face_no < 2) - { - Assert((dynamic_cast *>(&fe_other)), - ExcNotImplemented()); - } - else - { - Assert((dynamic_cast *>(&fe_other)), - ExcNotImplemented()); - } - - std::vector> result; - - for (unsigned int i = 0; i < this->n_dofs_per_quad(face_no); ++i) - result.emplace_back(i, i); - - return result; - } - - - - template - FE_WedgeDGP::FE_WedgeDGP(const unsigned int degree) - : FE_Wedge(degree, - get_dpo_vector_fe_wedge_dgp(degree), - FiniteElementData::L2) - {} - - - - template - std::unique_ptr> - FE_WedgeDGP::clone() const - { - return std::make_unique>(*this); - } - - - - template - std::string - FE_WedgeDGP::get_name() const - { - std::ostringstream namebuf; - namebuf << "FE_WedgeDGP<" << dim << ">(" << this->degree << ")"; - - return namebuf.str(); - } - - - - template - FE_Pyramid::FE_Pyramid( - const unsigned int degree, - const internal::GenericDoFsPerObject & dpos, - const typename FiniteElementData::Conformity conformity) - : dealii::FE_Poly( - Simplex::ScalarPyramidPolynomial(degree), - FiniteElementData(dpos, - ReferenceCells::Pyramid, - 1, - degree, - conformity), - std::vector( - FiniteElementData(dpos, ReferenceCells::Pyramid, 1, degree) - .dofs_per_cell, - true), - std::vector( - FiniteElementData(dpos, ReferenceCells::Pyramid, 1, degree) - .dofs_per_cell, - std::vector(1, true))) - { - AssertDimension(dim, 3); - - - if (degree == 1) - { - this->unit_support_points.emplace_back(-1.0, -1.0, 0.0); - this->unit_support_points.emplace_back(+1.0, -1.0, 0.0); - this->unit_support_points.emplace_back(-1.0, +1.0, 0.0); - this->unit_support_points.emplace_back(+1.0, +1.0, 0.0); - this->unit_support_points.emplace_back(+0.0, +0.0, 1.0); - } - } - - - - template - FE_PyramidP::FE_PyramidP(const unsigned int degree) - : FE_Pyramid(degree, - get_dpo_vector_fe_pyramid_p(degree), - FiniteElementData::H1) - {} - - - - template - std::unique_ptr> - FE_PyramidP::clone() const - { - return std::make_unique>(*this); - } - - - - template - std::string - FE_PyramidP::get_name() const - { - std::ostringstream namebuf; - namebuf << "FE_PyramidP<" << dim << ">(" << this->degree << ")"; - - return namebuf.str(); - } - - - - template - FiniteElementDomination::Domination - FE_PyramidP::compare_for_domination( - const FiniteElement &fe_other, - const unsigned int codim) const - { - Assert(codim <= dim, ExcImpossibleInDim(dim)); - - // vertex/line/face domination - // (if fe_other is derived from FE_DGP) - // ------------------------------------ - if (codim > 0) - if (dynamic_cast *>(&fe_other) != nullptr) - // there are no requirements between continuous and discontinuous - // elements - return FiniteElementDomination::no_requirements; - - // vertex/line/face domination - // (if fe_other is not derived from FE_DGP) - // & cell domination - // ---------------------------------------- - if (const FE_PyramidP *fe_pp_other = - dynamic_cast *>(&fe_other)) - { - if (this->degree < fe_pp_other->degree) - return FiniteElementDomination::this_element_dominates; - else if (this->degree == fe_pp_other->degree) - return FiniteElementDomination::either_element_can_dominate; - else - return FiniteElementDomination::other_element_dominates; - } - else if (const FE_P *fe_p_other = - dynamic_cast *>(&fe_other)) - { - if (this->degree < fe_p_other->degree) - return FiniteElementDomination::this_element_dominates; - else if (this->degree == fe_p_other->degree) - return FiniteElementDomination::either_element_can_dominate; - else - return FiniteElementDomination::other_element_dominates; - } - else if (const FE_Q *fe_q_other = - dynamic_cast *>(&fe_other)) - { - if (this->degree < fe_q_other->degree) - return FiniteElementDomination::this_element_dominates; - else if (this->degree == fe_q_other->degree) - return FiniteElementDomination::either_element_can_dominate; - else - return FiniteElementDomination::other_element_dominates; - } - else if (const FE_Nothing *fe_nothing = - dynamic_cast *>(&fe_other)) - { - if (fe_nothing->is_dominating()) - return FiniteElementDomination::other_element_dominates; - else - // the FE_Nothing has no degrees of freedom and it is typically used - // in a context where we don't require any continuity along the - // interface - return FiniteElementDomination::no_requirements; - } - - Assert(false, ExcNotImplemented()); - return FiniteElementDomination::neither_element_dominates; - } - - - - template - std::vector> - FE_PyramidP::hp_vertex_dof_identities( - const FiniteElement &fe_other) const - { - (void)fe_other; - - Assert((dynamic_cast *>(&fe_other)) || - (dynamic_cast *>(&fe_other)), - ExcNotImplemented()); - - return {{0, 0}}; - } - - - - template - std::vector> - FE_PyramidP::hp_line_dof_identities( - const FiniteElement &fe_other) const - { - (void)fe_other; - - Assert((dynamic_cast *>(&fe_other)) || - (dynamic_cast *>(&fe_other)), - ExcNotImplemented()); - - std::vector> result; - - for (unsigned int i = 0; i < this->degree - 1; ++i) - result.emplace_back(i, i); - - return result; - } - - - - template - std::vector> - FE_PyramidP::hp_quad_dof_identities( - const FiniteElement &fe_other, - const unsigned int face_no) const - { - (void)fe_other; - - - AssertIndexRange(face_no, 5); - - if (face_no == 0) - { - Assert((dynamic_cast *>(&fe_other)), - ExcNotImplemented()); - } - else - { - Assert((dynamic_cast *>(&fe_other)), - ExcNotImplemented()); - } - - std::vector> result; - - for (unsigned int i = 0; i < this->n_dofs_per_quad(face_no); ++i) - result.emplace_back(i, i); - - return result; - } - - - - template - FE_PyramidDGP::FE_PyramidDGP(const unsigned int degree) - : FE_Pyramid(degree, - get_dpo_vector_fe_pyramid_dgp(degree), - FiniteElementData::L2) - {} - - - - template - std::unique_ptr> - FE_PyramidDGP::clone() const - { - return std::make_unique>(*this); - } - - - - template - std::string - FE_PyramidDGP::get_name() const - { - std::ostringstream namebuf; - namebuf << "FE_PyramidDGP<" << dim << ">(" << this->degree << ")"; - - return namebuf.str(); - } - - - namespace FE_P_BubblesImplementation - { - template - std::vector - get_dpo_vector(const unsigned int degree) - { - std::vector dpo(dim + 1); - if (degree == 0) - { - dpo[dim] = 1; // single interior dof - } - else - { - Assert(degree == 1 || degree == 2, ExcNotImplemented()); - dpo[0] = 1; // vertex dofs - - if (degree == 2) - { - dpo[1] = 1; // line dofs - - if (dim > 1) - dpo[dim] = 1; // the internal bubble function - if (dim == 3) - dpo[dim - 1] = 1; // face bubble functions - } - } - - return dpo; - } - - - - template - std::vector> - unit_support_points(const unsigned int degree) - { - Assert(degree < 3, ExcNotImplemented()); - std::vector> points = unit_support_points_fe_poly(degree); - - Point centroid; - std::fill(centroid.begin_raw(), - centroid.end_raw(), - 1.0 / double(dim + 1)); - - switch (dim) - { - case 1: - // nothing more to do - return points; - case 2: - { - if (degree == 2) - points.push_back(centroid); - return points; - } - case 3: - { - if (degree == 2) - { - const double q13 = 1.0 / 3.0; - points.emplace_back(q13, q13, 0.0); - points.emplace_back(q13, 0.0, q13); - points.emplace_back(0.0, q13, q13); - points.emplace_back(q13, q13, q13); - points.push_back(centroid); - } - return points; - } - default: - Assert(false, ExcNotImplemented()); - } - return points; - } - - - - template - BarycentricPolynomials - get_basis(const unsigned int degree) - { - Point centroid; - std::fill(centroid.begin_raw(), - centroid.end_raw(), - 1.0 / double(dim + 1)); - - auto M = [](const unsigned int d) { - return BarycentricPolynomial::monomial(d); - }; - - switch (degree) - { - // we don't need to add bubbles to P0 or P1 - case 0: - case 1: - return BarycentricPolynomials::get_fe_p_basis(degree); - case 2: - { - const auto fe_p = - BarycentricPolynomials::get_fe_p_basis(degree); - // no further work is needed in 1D - if (dim == 1) - return fe_p; - - // in 2D and 3D we add a centroid bubble function - auto c_bubble = BarycentricPolynomial() + 1; - for (unsigned int d = 0; d < dim + 1; ++d) - c_bubble = c_bubble * M(d); - c_bubble = c_bubble / c_bubble.value(centroid); - - std::vector> bubble_functions; - if (dim == 2) - { - bubble_functions.push_back(c_bubble); - } - else if (dim == 3) - { - // need 'face bubble' functions in addition to the centroid. - // Furthermore we need to subtract them off from the other - // functions so that we end up with an interpolatory basis - auto b0 = 27 * M(0) * M(1) * M(2); - bubble_functions.push_back(b0 - - b0.value(centroid) * c_bubble); - auto b1 = 27 * M(0) * M(1) * M(3); - bubble_functions.push_back(b1 - - b1.value(centroid) * c_bubble); - auto b2 = 27 * M(0) * M(2) * M(3); - bubble_functions.push_back(b2 - - b2.value(centroid) * c_bubble); - auto b3 = 27 * M(1) * M(2) * M(3); - bubble_functions.push_back(b3 - - b3.value(centroid) * c_bubble); - - bubble_functions.push_back(c_bubble); - } - - // Extract out the support points for the extra bubble (both - // volume and face) functions: - const std::vector> support_points = - unit_support_points(degree); - const std::vector> bubble_support_points( - support_points.begin() + fe_p.n(), support_points.end()); - Assert(bubble_support_points.size() == bubble_functions.size(), - ExcInternalError()); - const unsigned int n_bubbles = bubble_support_points.size(); - - // Assemble the final basis: - std::vector> lump_polys; - for (unsigned int i = 0; i < fe_p.n(); ++i) - { - BarycentricPolynomial p = fe_p[i]; - - for (unsigned int j = 0; j < n_bubbles; ++j) - { - p = p - p.value(bubble_support_points[j]) * - bubble_functions[j]; - } - - lump_polys.push_back(p); - } - - for (auto &p : bubble_functions) - lump_polys.push_back(std::move(p)); - - // Sanity check: -#ifdef DEBUG - BarycentricPolynomial unity; - for (const auto &p : lump_polys) - unity = unity + p; - - Point test; - for (unsigned int d = 0; d < dim; ++d) - test[d] = 2.0; - Assert(std::abs(unity.value(test) - 1.0) < 1e-10, - ExcInternalError()); -#endif - - return BarycentricPolynomials(lump_polys); - } - default: - Assert(degree < 3, ExcNotImplemented()); - } - - Assert(degree < 3, ExcNotImplemented()); - // bogus return to placate compilers - return BarycentricPolynomials::get_fe_p_basis(degree); - } - - - - template - FiniteElementData - get_fe_data(const unsigned int degree) - { - // It's not efficient, but delegate computation of the degree of the - // finite element (which is different from the input argument) to the - // basis. - const auto polys = get_basis(degree); - return FiniteElementData(get_dpo_vector(degree), - ReferenceCells::get_simplex(), - 1, // n_components - polys.degree(), - FiniteElementData::H1); - } - } // namespace FE_P_BubblesImplementation - - - - template - FE_P_Bubbles::FE_P_Bubbles(const unsigned int degree) - : dealii::FE_Poly( - FE_P_BubblesImplementation::get_basis(degree), - FE_P_BubblesImplementation::get_fe_data(degree), - std::vector( - FE_P_BubblesImplementation::get_fe_data(degree).dofs_per_cell, - true), - std::vector( - FE_P_BubblesImplementation::get_fe_data(degree).dofs_per_cell, - std::vector(1, true))) - , approximation_degree(degree) - { - this->unit_support_points = - FE_P_BubblesImplementation::unit_support_points(degree); - - // TODO - // this->unit_face_support_points = - // unit_face_support_points_fe_poly(degree); - } - - - - template - std::string - FE_P_Bubbles::get_name() const - { - return "Simplex::FE_P_Bubbles<" + Utilities::dim_string(dim, spacedim) + - ">" + "(" + std::to_string(approximation_degree) + ")"; - } - - - - template - void - FE_P_Bubbles:: - convert_generalized_support_point_values_to_dof_values( - const std::vector> &support_point_values, - std::vector & nodal_values) const - { - AssertDimension(support_point_values.size(), - this->get_unit_support_points().size()); - AssertDimension(support_point_values.size(), nodal_values.size()); - AssertDimension(this->dofs_per_cell, nodal_values.size()); - - for (unsigned int i = 0; i < this->dofs_per_cell; ++i) - { - AssertDimension(support_point_values[i].size(), 1); - - nodal_values[i] = support_point_values[i](0); - } - } - - - - template - std::unique_ptr> - FE_P_Bubbles::clone() const - { - return std::make_unique>(*this); - } } // namespace Simplex // explicit instantiations -#include "fe_lib.inst" +#include "fe_simplex_p.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/source/simplex/grid_generator.inst.in b/source/fe/fe_simplex_p.inst.in similarity index 51% rename from source/simplex/grid_generator.inst.in rename to source/fe/fe_simplex_p.inst.in index b65ed758bb..0f29f676d1 100644 --- a/source/simplex/grid_generator.inst.in +++ b/source/fe/fe_simplex_p.inst.in @@ -15,19 +15,11 @@ -for (deal_II_dimension : DIMENSIONS) +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS) { - template void GridGenerator::subdivided_hyper_rectangle_with_simplices( - Triangulation & tria, - const std::vector &repetitions, - const Point & p1, - const Point & p2, - const bool colorize); - - template void GridGenerator::subdivided_hyper_cube_with_simplices( - Triangulation & tria, - const unsigned int repetitions, - const double p1, - const double p2, - const bool colorize); +#if deal_II_dimension <= deal_II_space_dimension + template class Simplex::FE_Poly; + template class Simplex::FE_P; + template class Simplex::FE_DGP; +#endif } diff --git a/source/fe/fe_simplex_p_bubbles.cc b/source/fe/fe_simplex_p_bubbles.cc new file mode 100644 index 0000000000..8c2ee43f53 --- /dev/null +++ b/source/fe/fe_simplex_p_bubbles.cc @@ -0,0 +1,390 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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 + +DEAL_II_NAMESPACE_OPEN + +namespace Simplex +{ + namespace + { + /** + * Set up a vector that contains the unit (reference) cell support points + * for FE_Poly and sufficiently similar elements. + */ + template + std::vector> + unit_support_points_fe_poly_bubbles(const unsigned int degree) + { + std::vector> unit_points; + + // Piecewise constants are a special case: use a support point at the + // centroid and only the centroid + if (degree == 0) + { + Point centroid; + std::fill(centroid.begin_raw(), + centroid.end_raw(), + 1.0 / double(dim + 1)); + unit_points.emplace_back(centroid); + return unit_points; + } + + if (dim == 1) + { + // We don't really have dim = 1 support for simplex elements yet, but + // its convenient for populating the face array + Assert(degree <= 2, ExcNotImplemented()); + if (degree >= 1) + { + unit_points.emplace_back(0.0); + unit_points.emplace_back(1.0); + + if (degree == 2) + unit_points.emplace_back(0.5); + } + } + else if (dim == 2) + { + Assert(degree <= 2, ExcNotImplemented()); + if (degree >= 1) + { + unit_points.emplace_back(0.0, 0.0); + unit_points.emplace_back(1.0, 0.0); + unit_points.emplace_back(0.0, 1.0); + + if (degree == 2) + { + unit_points.emplace_back(0.5, 0.0); + unit_points.emplace_back(0.5, 0.5); + unit_points.emplace_back(0.0, 0.5); + } + } + } + else if (dim == 3) + { + Assert(degree <= 2, ExcNotImplemented()); + if (degree >= 1) + { + unit_points.emplace_back(0.0, 0.0, 0.0); + unit_points.emplace_back(1.0, 0.0, 0.0); + unit_points.emplace_back(0.0, 1.0, 0.0); + unit_points.emplace_back(0.0, 0.0, 1.0); + + if (degree == 2) + { + unit_points.emplace_back(0.5, 0.0, 0.0); + unit_points.emplace_back(0.5, 0.5, 0.0); + unit_points.emplace_back(0.0, 0.5, 0.0); + unit_points.emplace_back(0.0, 0.0, 0.5); + unit_points.emplace_back(0.5, 0.0, 0.5); + unit_points.emplace_back(0.0, 0.5, 0.5); + } + } + } + else + { + Assert(false, ExcNotImplemented()); + } + + return unit_points; + } + } // namespace + + namespace FE_P_BubblesImplementation + { + template + std::vector + get_dpo_vector(const unsigned int degree) + { + std::vector dpo(dim + 1); + if (degree == 0) + { + dpo[dim] = 1; // single interior dof + } + else + { + Assert(degree == 1 || degree == 2, ExcNotImplemented()); + dpo[0] = 1; // vertex dofs + + if (degree == 2) + { + dpo[1] = 1; // line dofs + + if (dim > 1) + dpo[dim] = 1; // the internal bubble function + if (dim == 3) + dpo[dim - 1] = 1; // face bubble functions + } + } + + return dpo; + } + + + + template + std::vector> + unit_support_points(const unsigned int degree) + { + Assert(degree < 3, ExcNotImplemented()); + std::vector> points = + unit_support_points_fe_poly_bubbles(degree); + + Point centroid; + std::fill(centroid.begin_raw(), + centroid.end_raw(), + 1.0 / double(dim + 1)); + + switch (dim) + { + case 1: + // nothing more to do + return points; + case 2: + { + if (degree == 2) + points.push_back(centroid); + return points; + } + case 3: + { + if (degree == 2) + { + const double q13 = 1.0 / 3.0; + points.emplace_back(q13, q13, 0.0); + points.emplace_back(q13, 0.0, q13); + points.emplace_back(0.0, q13, q13); + points.emplace_back(q13, q13, q13); + points.push_back(centroid); + } + return points; + } + default: + Assert(false, ExcNotImplemented()); + } + return points; + } + + + + template + BarycentricPolynomials + get_basis(const unsigned int degree) + { + Point centroid; + std::fill(centroid.begin_raw(), + centroid.end_raw(), + 1.0 / double(dim + 1)); + + auto M = [](const unsigned int d) { + return BarycentricPolynomial::monomial(d); + }; + + switch (degree) + { + // we don't need to add bubbles to P0 or P1 + case 0: + case 1: + return BarycentricPolynomials::get_fe_p_basis(degree); + case 2: + { + const auto fe_p = + BarycentricPolynomials::get_fe_p_basis(degree); + // no further work is needed in 1D + if (dim == 1) + return fe_p; + + // in 2D and 3D we add a centroid bubble function + auto c_bubble = BarycentricPolynomial() + 1; + for (unsigned int d = 0; d < dim + 1; ++d) + c_bubble = c_bubble * M(d); + c_bubble = c_bubble / c_bubble.value(centroid); + + std::vector> bubble_functions; + if (dim == 2) + { + bubble_functions.push_back(c_bubble); + } + else if (dim == 3) + { + // need 'face bubble' functions in addition to the centroid. + // Furthermore we need to subtract them off from the other + // functions so that we end up with an interpolatory basis + auto b0 = 27 * M(0) * M(1) * M(2); + bubble_functions.push_back(b0 - + b0.value(centroid) * c_bubble); + auto b1 = 27 * M(0) * M(1) * M(3); + bubble_functions.push_back(b1 - + b1.value(centroid) * c_bubble); + auto b2 = 27 * M(0) * M(2) * M(3); + bubble_functions.push_back(b2 - + b2.value(centroid) * c_bubble); + auto b3 = 27 * M(1) * M(2) * M(3); + bubble_functions.push_back(b3 - + b3.value(centroid) * c_bubble); + + bubble_functions.push_back(c_bubble); + } + + // Extract out the support points for the extra bubble (both + // volume and face) functions: + const std::vector> support_points = + unit_support_points(degree); + const std::vector> bubble_support_points( + support_points.begin() + fe_p.n(), support_points.end()); + Assert(bubble_support_points.size() == bubble_functions.size(), + ExcInternalError()); + const unsigned int n_bubbles = bubble_support_points.size(); + + // Assemble the final basis: + std::vector> lump_polys; + for (unsigned int i = 0; i < fe_p.n(); ++i) + { + BarycentricPolynomial p = fe_p[i]; + + for (unsigned int j = 0; j < n_bubbles; ++j) + { + p = p - p.value(bubble_support_points[j]) * + bubble_functions[j]; + } + + lump_polys.push_back(p); + } + + for (auto &p : bubble_functions) + lump_polys.push_back(std::move(p)); + + // Sanity check: +#ifdef DEBUG + BarycentricPolynomial unity; + for (const auto &p : lump_polys) + unity = unity + p; + + Point test; + for (unsigned int d = 0; d < dim; ++d) + test[d] = 2.0; + Assert(std::abs(unity.value(test) - 1.0) < 1e-10, + ExcInternalError()); +#endif + + return BarycentricPolynomials(lump_polys); + } + default: + Assert(degree < 3, ExcNotImplemented()); + } + + Assert(degree < 3, ExcNotImplemented()); + // bogus return to placate compilers + return BarycentricPolynomials::get_fe_p_basis(degree); + } + + + + template + FiniteElementData + get_fe_data(const unsigned int degree) + { + // It's not efficient, but delegate computation of the degree of the + // finite element (which is different from the input argument) to the + // basis. + const auto polys = get_basis(degree); + return FiniteElementData(get_dpo_vector(degree), + ReferenceCells::get_simplex(), + 1, // n_components + polys.degree(), + FiniteElementData::H1); + } + } // namespace FE_P_BubblesImplementation + + + + template + FE_P_Bubbles::FE_P_Bubbles(const unsigned int degree) + : dealii::FE_Poly( + FE_P_BubblesImplementation::get_basis(degree), + FE_P_BubblesImplementation::get_fe_data(degree), + std::vector( + FE_P_BubblesImplementation::get_fe_data(degree).dofs_per_cell, + true), + std::vector( + FE_P_BubblesImplementation::get_fe_data(degree).dofs_per_cell, + std::vector(1, true))) + , approximation_degree(degree) + { + this->unit_support_points = + FE_P_BubblesImplementation::unit_support_points(degree); + + // TODO + // this->unit_face_support_points = + // unit_face_support_points_fe_poly(degree); + } + + + + template + std::string + FE_P_Bubbles::get_name() const + { + return "Simplex::FE_P_Bubbles<" + Utilities::dim_string(dim, spacedim) + + ">" + "(" + std::to_string(approximation_degree) + ")"; + } + + + + template + void + FE_P_Bubbles:: + convert_generalized_support_point_values_to_dof_values( + const std::vector> &support_point_values, + std::vector & nodal_values) const + { + AssertDimension(support_point_values.size(), + this->get_unit_support_points().size()); + AssertDimension(support_point_values.size(), nodal_values.size()); + AssertDimension(this->dofs_per_cell, nodal_values.size()); + + for (unsigned int i = 0; i < this->dofs_per_cell; ++i) + { + AssertDimension(support_point_values[i].size(), 1); + + nodal_values[i] = support_point_values[i](0); + } + } + + + + template + std::unique_ptr> + FE_P_Bubbles::clone() const + { + return std::make_unique>(*this); + } +} // namespace Simplex + +// explicit instantiations +#include "fe_simplex_p_bubbles.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/fe/fe_simplex_p_bubbles.inst.in b/source/fe/fe_simplex_p_bubbles.inst.in new file mode 100644 index 0000000000..ca82953a46 --- /dev/null +++ b/source/fe/fe_simplex_p_bubbles.inst.in @@ -0,0 +1,24 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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. +// +// --------------------------------------------------------------------- + + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + template class Simplex::FE_P_Bubbles; +#endif + } diff --git a/source/fe/fe_wedge_p.cc b/source/fe/fe_wedge_p.cc new file mode 100644 index 0000000000..1a1d1fd0ca --- /dev/null +++ b/source/fe/fe_wedge_p.cc @@ -0,0 +1,325 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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 + +DEAL_II_NAMESPACE_OPEN + +namespace Simplex +{ + namespace + { + /** + * Helper function to set up the dpo vector of FE_WedgeP for a given @p degree. + */ + internal::GenericDoFsPerObject + get_dpo_vector_fe_wedge_p(const unsigned int degree) + { + internal::GenericDoFsPerObject dpo; + + if (degree == 1) + { + dpo.dofs_per_object_exclusive = {{1}, {0}, {0, 0, 0, 0, 0}, {0}}; + dpo.dofs_per_object_inclusive = {{1}, {2}, {3, 3, 4, 4, 4}, {6}}; + dpo.object_index = {{}, {6}, {6}, {6}}; + dpo.first_object_index_on_face = {{}, + {3, 3, 4, 4, 4}, + {3, 3, 4, 4, 4}}; + } + else if (degree == 2) + { + dpo.dofs_per_object_exclusive = {{1}, {1}, {0, 0, 1, 1, 1}, {0}}; + dpo.dofs_per_object_inclusive = {{1}, {3}, {6, 6, 9, 9, 9}, {18}}; + dpo.object_index = {{}, {6}, {15, 15, 15, 16, 17}, {18}}; + dpo.first_object_index_on_face = {{}, + {3, 3, 4, 4, 4}, + {6, 6, 8, 8, 8}}; + } + else + { + Assert(false, ExcNotImplemented()); + } + + return dpo; + } + + /** + * Helper function to set up the dpo vector of FE_WedgeDGP for a given @p degree. + */ + internal::GenericDoFsPerObject + get_dpo_vector_fe_wedge_dgp(const unsigned int degree) + { + unsigned int n_dofs = 0; + + if (degree == 1) + n_dofs = 6; + else if (degree == 2) + n_dofs = 18; + else + Assert(false, ExcNotImplemented()); + + return internal::expand(3, {{0, 0, 0, n_dofs}}, ReferenceCells::Wedge); + } + } // namespace + + template + FE_Wedge::FE_Wedge( + const unsigned int degree, + const internal::GenericDoFsPerObject & dpos, + const typename FiniteElementData::Conformity conformity) + : dealii::FE_Poly( + Simplex::ScalarWedgePolynomial(degree), + FiniteElementData(dpos, + ReferenceCells::Wedge, + 1, + degree, + conformity), + std::vector( + FiniteElementData(dpos, ReferenceCells::Wedge, 1, degree) + .dofs_per_cell, + true), + std::vector( + FiniteElementData(dpos, ReferenceCells::Wedge, 1, degree) + .dofs_per_cell, + std::vector(1, true))) + { + AssertDimension(dim, 3); + + if (degree == 1) + { + this->unit_support_points.emplace_back(0.0, 0.0, 0.0); + this->unit_support_points.emplace_back(1.0, 0.0, 0.0); + this->unit_support_points.emplace_back(0.0, 1.0, 0.0); + this->unit_support_points.emplace_back(0.0, 0.0, 1.0); + this->unit_support_points.emplace_back(1.0, 0.0, 1.0); + this->unit_support_points.emplace_back(0.0, 1.0, 1.0); + } + } + + + + template + FE_WedgeP::FE_WedgeP(const unsigned int degree) + : FE_Wedge(degree, + get_dpo_vector_fe_wedge_p(degree), + FiniteElementData::H1) + {} + + + + template + std::unique_ptr> + FE_WedgeP::clone() const + { + return std::make_unique>(*this); + } + + + + template + std::string + FE_WedgeP::get_name() const + { + std::ostringstream namebuf; + namebuf << "FE_WedgeP<" << dim << ">(" << this->degree << ")"; + + return namebuf.str(); + } + + + + template + FiniteElementDomination::Domination + FE_WedgeP::compare_for_domination( + const FiniteElement &fe_other, + const unsigned int codim) const + { + Assert(codim <= dim, ExcImpossibleInDim(dim)); + + // vertex/line/face domination + // (if fe_other is derived from FE_DGP) + // ------------------------------------ + if (codim > 0) + if (dynamic_cast *>(&fe_other) != nullptr) + // there are no requirements between continuous and discontinuous + // elements + return FiniteElementDomination::no_requirements; + + + // vertex/line/face domination + // (if fe_other is not derived from FE_DGP) + // & cell domination + // ---------------------------------------- + if (const FE_WedgeP *fe_wp_other = + dynamic_cast *>(&fe_other)) + { + if (this->degree < fe_wp_other->degree) + return FiniteElementDomination::this_element_dominates; + else if (this->degree == fe_wp_other->degree) + return FiniteElementDomination::either_element_can_dominate; + else + return FiniteElementDomination::other_element_dominates; + } + else if (const FE_P *fe_p_other = + dynamic_cast *>(&fe_other)) + { + if (this->degree < fe_p_other->degree) + return FiniteElementDomination::this_element_dominates; + else if (this->degree == fe_p_other->degree) + return FiniteElementDomination::either_element_can_dominate; + else + return FiniteElementDomination::other_element_dominates; + } + else if (const FE_Q *fe_q_other = + dynamic_cast *>(&fe_other)) + { + if (this->degree < fe_q_other->degree) + return FiniteElementDomination::this_element_dominates; + else if (this->degree == fe_q_other->degree) + return FiniteElementDomination::either_element_can_dominate; + else + return FiniteElementDomination::other_element_dominates; + } + else if (const FE_Nothing *fe_nothing = + dynamic_cast *>(&fe_other)) + { + if (fe_nothing->is_dominating()) + return FiniteElementDomination::other_element_dominates; + else + // the FE_Nothing has no degrees of freedom and it is typically used + // in a context where we don't require any continuity along the + // interface + return FiniteElementDomination::no_requirements; + } + + Assert(false, ExcNotImplemented()); + return FiniteElementDomination::neither_element_dominates; + } + + + + template + std::vector> + FE_WedgeP::hp_vertex_dof_identities( + const FiniteElement &fe_other) const + { + (void)fe_other; + + Assert((dynamic_cast *>(&fe_other)) || + (dynamic_cast *>(&fe_other)), + ExcNotImplemented()); + + return {{0, 0}}; + } + + + + template + std::vector> + FE_WedgeP::hp_line_dof_identities( + const FiniteElement &fe_other) const + { + (void)fe_other; + + Assert((dynamic_cast *>(&fe_other)) || + (dynamic_cast *>(&fe_other)), + ExcNotImplemented()); + + std::vector> result; + + for (unsigned int i = 0; i < this->degree - 1; ++i) + result.emplace_back(i, i); + + return result; + } + + + + template + std::vector> + FE_WedgeP::hp_quad_dof_identities( + const FiniteElement &fe_other, + const unsigned int face_no) const + { + (void)fe_other; + + AssertIndexRange(face_no, 5); + + if (face_no < 2) + { + Assert((dynamic_cast *>(&fe_other)), + ExcNotImplemented()); + } + else + { + Assert((dynamic_cast *>(&fe_other)), + ExcNotImplemented()); + } + + std::vector> result; + + for (unsigned int i = 0; i < this->n_dofs_per_quad(face_no); ++i) + result.emplace_back(i, i); + + return result; + } + + + + template + FE_WedgeDGP::FE_WedgeDGP(const unsigned int degree) + : FE_Wedge(degree, + get_dpo_vector_fe_wedge_dgp(degree), + FiniteElementData::L2) + {} + + + + template + std::unique_ptr> + FE_WedgeDGP::clone() const + { + return std::make_unique>(*this); + } + + + + template + std::string + FE_WedgeDGP::get_name() const + { + std::ostringstream namebuf; + namebuf << "FE_WedgeDGP<" << dim << ">(" << this->degree << ")"; + + return namebuf.str(); + } +} // namespace Simplex + +// explicit instantiations +#include "fe_wedge_p.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/fe/fe_wedge_p.inst.in b/source/fe/fe_wedge_p.inst.in new file mode 100644 index 0000000000..9b6ee61535 --- /dev/null +++ b/source/fe/fe_wedge_p.inst.in @@ -0,0 +1,28 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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. +// +// --------------------------------------------------------------------- + + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + template class Simplex::FE_Wedge; + template class Simplex::FE_WedgeP; + template class Simplex::FE_WedgeDGP; +#endif + } diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index ce59fc0d2c..c845f2d454 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -8220,6 +8220,196 @@ namespace GridGenerator return surface_to_volume_mapping; } + + + template + void + subdivided_hyper_rectangle_with_simplices( + Triangulation & tria, + const std::vector &repetitions, + const Point & p1, + const Point & p2, + const bool colorize) + { +# ifndef DEAL_II_WITH_SIMPLEX_SUPPORT + Assert(false, ExcNeedsSimplexSupport()); +# endif + AssertDimension(dim, spacedim); + + AssertThrow(colorize == false, ExcNotImplemented()); + + std::vector> vertices; + std::vector> cells; + + if (dim == 2) + { + // determine cell sizes + const Point dx((p2[0] - p1[0]) / repetitions[0], + (p2[1] - p1[1]) / repetitions[1]); + + // create vertices + for (unsigned int j = 0; j <= repetitions[1]; ++j) + for (unsigned int i = 0; i <= repetitions[0]; ++i) + vertices.push_back( + Point(p1[0] + dx[0] * i, p1[1] + dx[1] * j)); + + // create cells + for (unsigned int j = 0; j < repetitions[1]; ++j) + for (unsigned int i = 0; i < repetitions[0]; ++i) + { + // create reference QUAD cell + std::array quad{{ + (j + 0) * (repetitions[0] + 1) + i + 0, // + (j + 0) * (repetitions[0] + 1) + i + 1, // + (j + 1) * (repetitions[0] + 1) + i + 0, // + (j + 1) * (repetitions[0] + 1) + i + 1 // + }}; // + + // TRI cell 0 + { + CellData tri; + tri.vertices = {quad[0], quad[1], quad[2]}; + cells.push_back(tri); + } + + // TRI cell 1 + { + CellData tri; + tri.vertices = {quad[3], quad[2], quad[1]}; + cells.push_back(tri); + } + } + } + else if (dim == 3) + { + // determine cell sizes + const Point dx((p2[0] - p1[0]) / repetitions[0], + (p2[1] - p1[1]) / repetitions[1], + (p2[2] - p1[2]) / repetitions[2]); + + // create vertices + for (unsigned int k = 0; k <= repetitions[2]; ++k) + for (unsigned int j = 0; j <= repetitions[1]; ++j) + for (unsigned int i = 0; i <= repetitions[0]; ++i) + vertices.push_back(Point(p1[0] + dx[0] * i, + p1[1] + dx[1] * j, + p1[2] + dx[2] * k)); + + // create cells + for (unsigned int k = 0; k < repetitions[2]; ++k) + for (unsigned int j = 0; j < repetitions[1]; ++j) + for (unsigned int i = 0; i < repetitions[0]; ++i) + { + // create reference HEX cell + std::array quad{ + {(k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 0, + (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 1, + (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 0, + (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 1, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 0, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 1, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 0, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 1}}; + + // TET cell 0 + { + CellData cell; + if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) + cell.vertices = {{quad[0], quad[1], quad[2], quad[4]}}; + else + cell.vertices = {{quad[0], quad[1], quad[3], quad[5]}}; + + cells.push_back(cell); + } + + // TET cell 1 + { + CellData cell; + if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) + cell.vertices = {{quad[2], quad[1], quad[3], quad[7]}}; + else + cell.vertices = {{quad[0], quad[3], quad[2], quad[6]}}; + cells.push_back(cell); + } + + // TET cell 2 + { + CellData cell; + if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) + cell.vertices = {{quad[1], quad[4], quad[5], quad[7]}}; + else + cell.vertices = {{quad[0], quad[4], quad[5], quad[6]}}; + cells.push_back(cell); + } + + // TET cell 3 + { + CellData cell; + if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) + cell.vertices = {{quad[2], quad[4], quad[7], quad[6]}}; + else + cell.vertices = {{quad[3], quad[5], quad[7], quad[6]}}; + cells.push_back(cell); + } + + // TET cell 4 + { + CellData cell; + if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) + cell.vertices = {{quad[1], quad[2], quad[4], quad[7]}}; + else + cell.vertices = {{quad[0], quad[3], quad[6], quad[5]}}; + cells.push_back(cell); + } + } + } + else + { + AssertThrow(colorize == false, ExcNotImplemented()); + } + + // actually create triangulation + tria.create_triangulation(vertices, cells, SubCellData()); + } + + + + template + void + subdivided_hyper_cube_with_simplices(Triangulation &tria, + const unsigned int repetitions, + const double p1, + const double p2, + const bool colorize) + { + if (dim == 2) + { + subdivided_hyper_rectangle_with_simplices( + tria, {{repetitions, repetitions}}, {p1, p1}, {p2, p2}, colorize); + } + else if (dim == 3) + { + subdivided_hyper_rectangle_with_simplices( + tria, + {{repetitions, repetitions, repetitions}}, + {p1, p1, p1}, + {p2, p2, p2}, + colorize); + } + else + { + AssertThrow(false, ExcNotImplemented()) + } + } + } // namespace GridGenerator // explicit instantiations diff --git a/source/grid/grid_generator.inst.in b/source/grid/grid_generator.inst.in index 17b7d02ad0..3f063e77d0 100644 --- a/source/grid/grid_generator.inst.in +++ b/source/grid/grid_generator.inst.in @@ -281,3 +281,20 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) #endif \} } + +for (deal_II_dimension : DIMENSIONS) + { + template void GridGenerator::subdivided_hyper_rectangle_with_simplices( + Triangulation & tria, + const std::vector &repetitions, + const Point & p1, + const Point & p2, + const bool colorize); + + template void GridGenerator::subdivided_hyper_cube_with_simplices( + Triangulation & tria, + const unsigned int repetitions, + const double p1, + const double p2, + const bool colorize); + } diff --git a/source/grid/reference_cell.cc b/source/grid/reference_cell.cc index 41b656fc6f..a9bc30f535 100644 --- a/source/grid/reference_cell.cc +++ b/source/grid/reference_cell.cc @@ -14,9 +14,14 @@ // --------------------------------------------------------------------- #include +#include #include #include +#include +#include +#include +#include #include #include #include @@ -25,10 +30,6 @@ #include #include -#include -#include -#include - #include DEAL_II_NAMESPACE_OPEN diff --git a/source/simplex/CMakeLists.txt b/source/simplex/CMakeLists.txt deleted file mode 100644 index 6814164168..0000000000 --- a/source/simplex/CMakeLists.txt +++ /dev/null @@ -1,48 +0,0 @@ -## --------------------------------------------------------------------- -## -## Copyright (C) 2020 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_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) - -SET(_unity_include_src - fe_lib.cc - grid_generator.cc - barycentric_polynomials.cc - polynomials.cc - quadrature_lib.cc - ) - -SET(_separate_src - ) - -# include all files in the unity file -SET(_n_includes_per_unity_file 20) - -SETUP_SOURCE_LIST("${_unity_include_src}" - "${_separate_src}" - ${_n_includes_per_unity_file} - _src - ) - -SET(_inst - fe_lib.inst.in - grid_generator.inst.in -) - -FILE(GLOB _header - ${CMAKE_SOURCE_DIR}/include/deal.II/simplex/*.h - ) - -DEAL_II_ADD_LIBRARY(obj_simplex OBJECT ${_src} ${_header} ${_inst}) -EXPAND_INSTANTIATIONS(obj_simplex "${_inst}") diff --git a/source/simplex/grid_generator.cc b/source/simplex/grid_generator.cc deleted file mode 100644 index aaed0379d3..0000000000 --- a/source/simplex/grid_generator.cc +++ /dev/null @@ -1,214 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2020 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 - -DEAL_II_NAMESPACE_OPEN - -namespace GridGenerator -{ - template - void - subdivided_hyper_rectangle_with_simplices( - Triangulation & tria, - const std::vector &repetitions, - const Point & p1, - const Point & p2, - const bool colorize) - { -#ifndef DEAL_II_WITH_SIMPLEX_SUPPORT - Assert(false, ExcNeedsSimplexSupport()); -#endif - AssertDimension(dim, spacedim); - - AssertThrow(colorize == false, ExcNotImplemented()); - - std::vector> vertices; - std::vector> cells; - - if (dim == 2) - { - // determine cell sizes - const Point dx((p2[0] - p1[0]) / repetitions[0], - (p2[1] - p1[1]) / repetitions[1]); - - // create vertices - for (unsigned int j = 0; j <= repetitions[1]; ++j) - for (unsigned int i = 0; i <= repetitions[0]; ++i) - vertices.push_back( - Point(p1[0] + dx[0] * i, p1[1] + dx[1] * j)); - - // create cells - for (unsigned int j = 0; j < repetitions[1]; ++j) - for (unsigned int i = 0; i < repetitions[0]; ++i) - { - // create reference QUAD cell - std::array quad{{ - (j + 0) * (repetitions[0] + 1) + i + 0, // - (j + 0) * (repetitions[0] + 1) + i + 1, // - (j + 1) * (repetitions[0] + 1) + i + 0, // - (j + 1) * (repetitions[0] + 1) + i + 1 // - }}; // - - // TRI cell 0 - { - CellData tri; - tri.vertices = {quad[0], quad[1], quad[2]}; - cells.push_back(tri); - } - - // TRI cell 1 - { - CellData tri; - tri.vertices = {quad[3], quad[2], quad[1]}; - cells.push_back(tri); - } - } - } - else if (dim == 3) - { - // determine cell sizes - const Point dx((p2[0] - p1[0]) / repetitions[0], - (p2[1] - p1[1]) / repetitions[1], - (p2[2] - p1[2]) / repetitions[2]); - - // create vertices - for (unsigned int k = 0; k <= repetitions[2]; ++k) - for (unsigned int j = 0; j <= repetitions[1]; ++j) - for (unsigned int i = 0; i <= repetitions[0]; ++i) - vertices.push_back(Point(p1[0] + dx[0] * i, - p1[1] + dx[1] * j, - p1[2] + dx[2] * k)); - - // create cells - for (unsigned int k = 0; k < repetitions[2]; ++k) - for (unsigned int j = 0; j < repetitions[1]; ++j) - for (unsigned int i = 0; i < repetitions[0]; ++i) - { - // create reference HEX cell - std::array quad{ - {(k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + - (j + 0) * (repetitions[0] + 1) + i + 0, - (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + - (j + 0) * (repetitions[0] + 1) + i + 1, - (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + - (j + 1) * (repetitions[0] + 1) + i + 0, - (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + - (j + 1) * (repetitions[0] + 1) + i + 1, - (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + - (j + 0) * (repetitions[0] + 1) + i + 0, - (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + - (j + 0) * (repetitions[0] + 1) + i + 1, - (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + - (j + 1) * (repetitions[0] + 1) + i + 0, - (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + - (j + 1) * (repetitions[0] + 1) + i + 1}}; - - // TET cell 0 - { - CellData cell; - if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) - cell.vertices = {{quad[0], quad[1], quad[2], quad[4]}}; - else - cell.vertices = {{quad[0], quad[1], quad[3], quad[5]}}; - - cells.push_back(cell); - } - - // TET cell 1 - { - CellData cell; - if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) - cell.vertices = {{quad[2], quad[1], quad[3], quad[7]}}; - else - cell.vertices = {{quad[0], quad[3], quad[2], quad[6]}}; - cells.push_back(cell); - } - - // TET cell 2 - { - CellData cell; - if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) - cell.vertices = {{quad[1], quad[4], quad[5], quad[7]}}; - else - cell.vertices = {{quad[0], quad[4], quad[5], quad[6]}}; - cells.push_back(cell); - } - - // TET cell 3 - { - CellData cell; - if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) - cell.vertices = {{quad[2], quad[4], quad[7], quad[6]}}; - else - cell.vertices = {{quad[3], quad[5], quad[7], quad[6]}}; - cells.push_back(cell); - } - - // TET cell 4 - { - CellData cell; - if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) - cell.vertices = {{quad[1], quad[2], quad[4], quad[7]}}; - else - cell.vertices = {{quad[0], quad[3], quad[6], quad[5]}}; - cells.push_back(cell); - } - } - } - else - { - AssertThrow(colorize == false, ExcNotImplemented()); - } - - // actually create triangulation - tria.create_triangulation(vertices, cells, SubCellData()); - } - - template - void - subdivided_hyper_cube_with_simplices(Triangulation &tria, - const unsigned int repetitions, - const double p1, - const double p2, - const bool colorize) - { - if (dim == 2) - { - subdivided_hyper_rectangle_with_simplices( - tria, {{repetitions, repetitions}}, {p1, p1}, {p2, p2}, colorize); - } - else if (dim == 3) - { - subdivided_hyper_rectangle_with_simplices( - tria, - {{repetitions, repetitions, repetitions}}, - {p1, p1, p1}, - {p2, p2, p2}, - colorize); - } - else - { - AssertThrow(false, ExcNotImplemented()) - } - } -} // namespace GridGenerator - -// explicit instantiations -#include "grid_generator.inst" - -DEAL_II_NAMESPACE_CLOSE diff --git a/source/simplex/quadrature_lib.cc b/source/simplex/quadrature_lib.cc deleted file mode 100644 index cbbd1c9ad6..0000000000 --- a/source/simplex/quadrature_lib.cc +++ /dev/null @@ -1,485 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2020 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 - - -DEAL_II_NAMESPACE_OPEN - -namespace Simplex -{ - template - QGauss::QGauss(const unsigned int n_points_1D) - : QSimplex(Quadrature()) - { - // fill quadrature points and quadrature weights - if (dim == 1) - { - const dealii::QGauss quad(n_points_1D); - - this->quadrature_points = quad.get_points(); - this->weights = quad.get_weights(); - } - else if (dim == 2) - { - if (n_points_1D == 1) - { - const double p = 1.0 / 3.0; - this->quadrature_points.emplace_back(p, p); - this->weights.emplace_back(0.5); - } - else if (n_points_1D == 2) - { - const double Q23 = 2.0 / 3.0; - const double Q16 = 1.0 / 6.0; - - this->quadrature_points.emplace_back(Q23, Q16); - this->quadrature_points.emplace_back(Q16, Q23); - this->quadrature_points.emplace_back(Q16, Q16); - this->weights.emplace_back(Q16); - this->weights.emplace_back(Q16); - this->weights.emplace_back(Q16); - } - else if (n_points_1D == 3) - { - const double q12 = 0.5; - - // clang-format off - this->quadrature_points.emplace_back(0.3333333333330, 0.3333333333330); - this->quadrature_points.emplace_back(0.7974269853530, 0.1012865073230); - this->quadrature_points.emplace_back(0.1012865073230, 0.7974269853530); - this->quadrature_points.emplace_back(0.1012865073230, 0.1012865073230); - this->quadrature_points.emplace_back(0.0597158717898, 0.4701420641050); - this->quadrature_points.emplace_back(0.4701420641050, 0.0597158717898); - this->quadrature_points.emplace_back(0.4701420641050, 0.4701420641050); - // clang-format on - - this->weights.emplace_back(q12 * 0.225); - this->weights.emplace_back(q12 * 0.125939180545); - this->weights.emplace_back(q12 * 0.125939180545); - this->weights.emplace_back(q12 * 0.125939180545); - this->weights.emplace_back(q12 * 0.132394152789); - this->weights.emplace_back(q12 * 0.132394152789); - this->weights.emplace_back(q12 * 0.132394152789); - } - else if (n_points_1D == 4) - { - Quadrature::operator=(QWitherdenVincent(n_points_1D)); - } - } - else if (dim == 3) - { - if (n_points_1D == 1) - { - const double Q14 = 1.0 / 4.0; - const double Q16 = 1.0 / 6.0; - - this->quadrature_points.emplace_back(Q14, Q14, Q14); - this->weights.emplace_back(Q16); - } - else if (n_points_1D == 2) - { - const double Q124 = 1.0 / 6.0 / 4.0; - - const double palpha = (5.0 + 3.0 * sqrt(5.0)) / 20.0; - const double pbeta = (5.0 - sqrt(5.0)) / 20.0; - this->quadrature_points.emplace_back(pbeta, pbeta, pbeta); - this->quadrature_points.emplace_back(palpha, pbeta, pbeta); - this->quadrature_points.emplace_back(pbeta, palpha, pbeta); - this->quadrature_points.emplace_back(pbeta, pbeta, palpha); - this->weights.emplace_back(Q124); - this->weights.emplace_back(Q124); - this->weights.emplace_back(Q124); - this->weights.emplace_back(Q124); - } - else if (n_points_1D == 3) - { - const double Q16 = 1.0 / 6.0; - - // clang-format off - this->quadrature_points.emplace_back(0.5684305841968444, 0.1438564719343852, 0.1438564719343852); - this->quadrature_points.emplace_back(0.1438564719343852, 0.1438564719343852, 0.1438564719343852); - this->quadrature_points.emplace_back(0.1438564719343852, 0.1438564719343852, 0.5684305841968444); - this->quadrature_points.emplace_back(0.1438564719343852, 0.5684305841968444, 0.1438564719343852); - this->quadrature_points.emplace_back(0.0000000000000000, 0.5000000000000000, 0.5000000000000000); - this->quadrature_points.emplace_back(0.5000000000000000, 0.0000000000000000, 0.5000000000000000); - this->quadrature_points.emplace_back(0.5000000000000000, 0.5000000000000000, 0.0000000000000000); - this->quadrature_points.emplace_back(0.5000000000000000, 0.0000000000000000, 0.0000000000000000); - this->quadrature_points.emplace_back(0.0000000000000000, 0.5000000000000000, 0.0000000000000000); - this->quadrature_points.emplace_back(0.0000000000000000, 0.0000000000000000, 0.5000000000000000); - // clang-format on - - this->weights.emplace_back(0.2177650698804054 * Q16); - this->weights.emplace_back(0.2177650698804054 * Q16); - this->weights.emplace_back(0.2177650698804054 * Q16); - this->weights.emplace_back(0.2177650698804054 * Q16); - this->weights.emplace_back(0.0214899534130631 * Q16); - this->weights.emplace_back(0.0214899534130631 * Q16); - this->weights.emplace_back(0.0214899534130631 * Q16); - this->weights.emplace_back(0.0214899534130631 * Q16); - this->weights.emplace_back(0.0214899534130631 * Q16); - this->weights.emplace_back(0.0214899534130631 * Q16); - } - else if (n_points_1D == 4) - { - Quadrature::operator=(QWitherdenVincent(n_points_1D)); - } - } - - AssertDimension(this->quadrature_points.size(), this->weights.size()); - Assert(this->quadrature_points.size() > 0, - ExcNotImplemented( - "Simplex::QGauss is currently only implemented for " - "n_points_1D = 1, 2, 3, and 4 while you are asking for " - "n_points_1D = " + - Utilities::to_string(n_points_1D))); - } - - namespace - { - template - std::vector> - all_permutations(const std::array &b_point) - { - std::vector> output; - - // We want all possible permutations of the barycentric coordinates. - // The easiest way to get all of them is to sort the input first and - // then use next_permutation to cycle through them all. - std::array temp = b_point; - std::sort(temp.begin(), temp.end()); - do - { - output.push_back(temp); - } - while (std::next_permutation(temp.begin(), temp.end())); - - return output; - } - } // namespace - - - - template - QWitherdenVincent::QWitherdenVincent(const unsigned int n_points_1D) - : QSimplex(Quadrature()) - { - Assert(1 <= dim && dim <= 3, ExcNotImplemented()); - // Just use Gauss in 1D: this is a high-order open rule so this is a - // reasonable equivalent for generic programming. - if (dim == 1) - { - Quadrature::operator=(dealii::QGauss(n_points_1D)); - return; - } - - std::array centroid; - std::fill(centroid.begin(), centroid.end(), 1.0 / (dim + 1.0)); - std::vector>> b_point_permutations; - std::vector b_weights; - - // We can simplify the implementation of these quadrature rules - // by quite a bit by exploiting symmetry - we do essentially the - // same thing for each barycentric coordinate, so we can express - // our quadrature rule as permutations of barycentric points - // instead of writing things out explicitly. - - // Apply a Barycentric permutation where one point is different. - auto process_point_1 = [&](const double a, const double w) { - const double b = 1.0 - dim * a; - std::array b_point; - std::fill(b_point.begin(), b_point.begin() + dim, a); - b_point[dim] = b; - - b_weights.push_back(w); - b_point_permutations.push_back(all_permutations(b_point)); - }; - - // Apply a Barycentric permutation where two points (in 3D) are different. - auto process_point_2 = [&](const double a, const double w) { - Assert(dim == 3, ExcInternalError()); - const double b = (1.0 - 2.0 * a) / 2.0; - std::array b_point; - std::fill(b_point.begin(), b_point.begin() + dim - 1, a); - b_point[dim - 1] = b; - b_point[dim] = b; - - b_weights.push_back(w); - b_point_permutations.push_back(all_permutations(b_point)); - }; - - // Apply a Barycentric permutation where three (or four) points - // are different (since there are two inputs). - auto process_point_3 = [&](const double a, const double b, const double w) { - const double c = 1.0 - (dim - 1.0) * a - b; - std::array b_point; - std::fill(b_point.begin(), b_point.begin() + dim - 1, a); - b_point[dim - 1] = b; - b_point[dim] = c; - - b_weights.push_back(w); - b_point_permutations.push_back(all_permutations(b_point)); - }; - - if (n_points_1D == 1) - { - b_point_permutations.push_back({centroid}); - b_weights.push_back(1.0); - } - else if (n_points_1D == 2) - { - // This is WV-4 in 2D and WV-3 in 3D - if (dim == 2) - { - process_point_1(9.1576213509770743e-02, 1.0995174365532187e-01); - process_point_1(4.4594849091596489e-01, 2.2338158967801147e-01); - } - else if (dim == 3) - { - process_point_1(3.281633025163817e-01, 1.362178425370874e-01); - process_point_1(1.080472498984286e-01, 1.137821574629126e-01); - } - } - else if (n_points_1D == 3) - { - // This is the WV-5 rule in both 2D and 3D - if (dim == 2) - { - b_weights.push_back(0.225); - b_point_permutations.push_back({centroid}); - - process_point_1(1.0128650732345634e-01, 1.2593918054482714e-01); - process_point_1(4.7014206410511511e-01, 1.3239415278850619e-01); - } - else if (dim == 3) - { - process_point_1(3.108859192633006e-01, 1.126879257180159e-01); - process_point_1(9.273525031089125e-02, 7.349304311636196e-02); - - process_point_2(4.550370412564964e-02, 4.254602077708147e-02); - } - } - else if (n_points_1D == 4) - { - // This is the WV-7 rule in both 2D and 3D - if (dim == 2) - { - process_point_1(3.3730648554587850e-02, 1.6545050110792131e-02); - process_point_1(4.7430969250471822e-01, 7.7086646185986069e-02); - process_point_1(2.4157738259540357e-01, 1.2794417123015558e-01); - process_point_3(4.7036644652595216e-02, - 1.9868331479735168e-01, - 5.5878732903199779e-02); - } - else if (dim == 3) - { - b_point_permutations.push_back({centroid}); - b_weights.push_back(9.548528946413085e-02); - - process_point_1(3.157011497782028e-01, 4.232958120996703e-02); - process_point_2(5.048982259839635e-02, 3.189692783285758e-02); - - process_point_3(1.888338310260010e-01, - 5.751716375870000e-01, - 3.720713072833462e-02); - process_point_3(2.126547254148314e-02, - 8.108302410985486e-01, - 8.110770829903342e-03); - } - } - else if (n_points_1D == 5) - { - // This is the WV-9 rule in both 2D and 3D - if (dim == 2) - { - b_point_permutations.push_back({centroid}); - b_weights.push_back(9.7135796282798836e-02); - - process_point_1(4.4729513394452691e-02, 2.5577675658698031e-02); - process_point_1(4.8968251919873762e-01, 3.1334700227139071e-02); - process_point_1(4.3708959149293664e-01, 7.7827541004774278e-02); - process_point_1(1.8820353561903275e-01, 7.9647738927210249e-02); - - process_point_3(3.6838412054736258e-02, - 2.2196298916076568e-01, - 4.3283539377289376e-02); - } - else if (dim == 3) - { - b_point_permutations.push_back({centroid}); - b_weights.push_back(5.801054891248025e-02); - - process_point_1(6.198169755222693e-10, 6.431928175925639e-05); - process_point_1(1.607745353952616e-01, 2.317333846242546e-02); - process_point_1(3.222765218214210e-01, 2.956291233542929e-02); - process_point_1(4.510891834541358e-02, 8.063979979616182e-03); - - process_point_2(1.122965460043761e-01, 3.813408010370246e-02); - - process_point_3(4.588714487524592e-01, - 2.554579233041310e-03, - 8.384422198298552e-03); - process_point_3(3.377587068533860e-02, - 7.183503264420745e-01, - 1.023455935274533e-02); - process_point_3(1.836413698099279e-01, - 3.441591057817528e-02, - 2.052491596798814e-02); - } - } - else if (n_points_1D == 6) - { - // There is no WV-11 rule in 3D yet - if (dim == 2) - { - b_point_permutations.push_back({centroid}); - b_weights.push_back(8.5761179732224219e-02); - - process_point_1(2.8485417614371900e-02, 1.0431870512894697e-02); - process_point_1(4.9589190096589092e-01, 1.6606273054585369e-02); - process_point_1(1.0263548271224643e-01, 3.8630759237019321e-02); - process_point_1(4.3846592676435220e-01, 6.7316154079468296e-02); - process_point_1(2.1021995670317828e-01, 7.0515684111716576e-02); - - process_point_3(7.3254276860644785e-03, - 1.4932478865208237e-01, - 1.0290289572953278e-02); - process_point_3(4.6010500165429957e-02, - 2.8958112563770588e-01, - 4.0332476640500554e-02); - } - else if (dim == 3) - { - Assert(false, ExcNotImplemented()); - } - } - else - { - Assert(false, ExcNotImplemented()); - } - - Assert(b_point_permutations.size() == b_weights.size(), ExcInternalError()); - for (unsigned int permutation_n = 0; permutation_n < b_weights.size(); - ++permutation_n) - { - for (const std::array &b_point : - b_point_permutations[permutation_n]) - { - const double volume = (dim == 2 ? 1.0 / 2.0 : 1.0 / 6.0); - this->weights.emplace_back(volume * b_weights[permutation_n]); - Point c_point; - std::copy(b_point.begin(), - b_point.begin() + dim, - c_point.begin_raw()); - this->quadrature_points.emplace_back(c_point); - } - } - } - - - - template - QGaussWedge::QGaussWedge(const unsigned int n_points) - : Quadrature() - { - AssertDimension(dim, 3); - - Simplex::QGauss<2> quad_tri(n_points); - QGauss<1> quad_line(n_points); - - for (unsigned int i = 0; i < quad_line.size(); ++i) - for (unsigned int j = 0; j < quad_tri.size(); ++j) - { - this->quadrature_points.emplace_back(quad_tri.point(j)[0], - quad_tri.point(j)[1], - quad_line.point(i)[0]); - this->weights.emplace_back(quad_tri.weight(j) * quad_line.weight(i)); - } - - AssertDimension(this->quadrature_points.size(), this->weights.size()); - Assert(this->quadrature_points.size() > 0, - ExcMessage("No valid quadrature points!")); - } - - - - template - QGaussPyramid::QGaussPyramid(const unsigned int n_points_1D) - : Quadrature() - { - AssertDimension(dim, 3); - - if (n_points_1D == 1) - { - const double Q14 = 1.0 / 4.0; - const double Q43 = 4.0 / 3.0; - - this->quadrature_points.emplace_back(0, 0, Q14); - this->weights.emplace_back(Q43); - } - else if (n_points_1D == 2) - { - // clang-format off - this->quadrature_points.emplace_back(-0.26318405556971, -0.26318405556971, 0.54415184401122); - this->quadrature_points.emplace_back(-0.50661630334979, -0.50661630334979, 0.12251482265544); - this->quadrature_points.emplace_back(-0.26318405556971, +0.26318405556971, 0.54415184401122); - this->quadrature_points.emplace_back(-0.50661630334979, +0.50661630334979, 0.12251482265544); - this->quadrature_points.emplace_back(+0.26318405556971, -0.26318405556971, 0.54415184401122); - this->quadrature_points.emplace_back(+0.50661630334979, -0.50661630334979, 0.12251482265544); - this->quadrature_points.emplace_back(+0.26318405556971, +0.26318405556971, 0.54415184401122); - this->quadrature_points.emplace_back(+0.50661630334979, +0.50661630334979, 0.12251482265544); - // clang-format on - - this->weights.emplace_back(0.10078588207983); - this->weights.emplace_back(0.23254745125351); - this->weights.emplace_back(0.10078588207983); - this->weights.emplace_back(0.23254745125351); - this->weights.emplace_back(0.10078588207983); - this->weights.emplace_back(0.23254745125351); - this->weights.emplace_back(0.10078588207983); - this->weights.emplace_back(0.23254745125351); - } - - AssertDimension(this->quadrature_points.size(), this->weights.size()); - Assert(this->quadrature_points.size() > 0, - ExcMessage("No valid quadrature points!")); - } - -} // namespace Simplex - - -template class Simplex::QGauss<1>; -template class Simplex::QGauss<2>; -template class Simplex::QGauss<3>; -template class Simplex::QGaussWedge<1>; -template class Simplex::QGaussWedge<2>; -template class Simplex::QGaussWedge<3>; -template class Simplex::QGaussPyramid<1>; -template class Simplex::QGaussPyramid<2>; -template class Simplex::QGaussPyramid<3>; - -template class Simplex::QWitherdenVincent<1>; -template class Simplex::QWitherdenVincent<2>; -template class Simplex::QWitherdenVincent<3>; - -DEAL_II_NAMESPACE_CLOSE diff --git a/tests/matrix_free/compute_diagonal_05.cc b/tests/matrix_free/compute_diagonal_05.cc index 79a5ec2fbe..5daa7c3d05 100644 --- a/tests/matrix_free/compute_diagonal_05.cc +++ b/tests/matrix_free/compute_diagonal_05.cc @@ -17,11 +17,12 @@ // Test MatrixFreeTools::compute_diagonal() for a Laplace operator on a simplex // mesh. +#include + +#include #include -#include -#include -#include +#include #include "compute_diagonal_util.h" diff --git a/tests/multigrid-global-coarsening/multigrid_util.h b/tests/multigrid-global-coarsening/multigrid_util.h index d52957f3bb..bc59afd3cf 100644 --- a/tests/multigrid-global-coarsening/multigrid_util.h +++ b/tests/multigrid-global-coarsening/multigrid_util.h @@ -27,6 +27,7 @@ #include #include +#include #include #include @@ -54,10 +55,6 @@ #include #include -#include -#include -#include - #include "../tests.h" template diff --git a/tests/simplex/barycentric_01.cc b/tests/simplex/barycentric_01.cc index 6bd73fbff9..277d7220be 100644 --- a/tests/simplex/barycentric_01.cc +++ b/tests/simplex/barycentric_01.cc @@ -16,11 +16,13 @@ // Test Simplex::BarycentricPolynomial and Simplex::BarycentricPolynomials. #include +#include #include -#include -#include -#include +#include +#include +#include +#include #include "../tests.h" diff --git a/tests/simplex/cell_measure_01.cc b/tests/simplex/cell_measure_01.cc index a21e5837c8..12488df21f 100644 --- a/tests/simplex/cell_measure_01.cc +++ b/tests/simplex/cell_measure_01.cc @@ -17,13 +17,15 @@ // Test TriaAccessor::measure() and TriaAccessor::diameter(). +#include +#include +#include +#include #include #include #include -#include - #include "../tests.h" template diff --git a/tests/simplex/compute_point_locations_01.cc b/tests/simplex/compute_point_locations_01.cc index 66dc56cb39..7f01defc34 100644 --- a/tests/simplex/compute_point_locations_01.cc +++ b/tests/simplex/compute_point_locations_01.cc @@ -11,17 +11,20 @@ // // --------------------------------------------------------------------- +#include + +#include +#include +#include #include #include +#include #include +#include #include #include -#include -#include -#include - #include "../tests.h" // Create a simplex mesh in the unit cube. Check points distribute to each cell diff --git a/tests/simplex/compute_projection_matrices_01.cc b/tests/simplex/compute_projection_matrices_01.cc index 02d85f17b5..0eb96bc768 100644 --- a/tests/simplex/compute_projection_matrices_01.cc +++ b/tests/simplex/compute_projection_matrices_01.cc @@ -18,19 +18,23 @@ // (and indirectly FETools::compute_embedding_matrices() for simplices). +#include + +#include +#include +#include #include #include +#include #include +#include + #include #include #include -#include -#include -#include - #include "../tests.h" using namespace dealii; diff --git a/tests/simplex/conv_hex_to_simplex.cc b/tests/simplex/conv_hex_to_simplex.cc index 10dab721d2..975c8a8e55 100644 --- a/tests/simplex/conv_hex_to_simplex.cc +++ b/tests/simplex/conv_hex_to_simplex.cc @@ -27,8 +27,6 @@ #include #include -#include - #include #include diff --git a/tests/simplex/data_out_write_hdf5_01.cc b/tests/simplex/data_out_write_hdf5_01.cc index 2c9b1b93c7..bcae95fde2 100644 --- a/tests/simplex/data_out_write_hdf5_01.cc +++ b/tests/simplex/data_out_write_hdf5_01.cc @@ -21,9 +21,14 @@ #include +#include +#include +#include #include +#include #include +#include #include #include #include @@ -31,9 +36,6 @@ #include #include -#include -#include - #include "../tests.h" using namespace dealii; diff --git a/tests/simplex/data_out_write_hdf5_02.cc b/tests/simplex/data_out_write_hdf5_02.cc index 9b17353106..038c15220b 100644 --- a/tests/simplex/data_out_write_hdf5_02.cc +++ b/tests/simplex/data_out_write_hdf5_02.cc @@ -24,9 +24,14 @@ #include #include +#include +#include +#include #include +#include #include +#include #include #include #include @@ -34,9 +39,6 @@ #include #include -#include -#include - #include "../tests.h" using namespace dealii; diff --git a/tests/simplex/data_out_write_vtk_01.cc b/tests/simplex/data_out_write_vtk_01.cc index 0886b0fc42..b2d6b36480 100644 --- a/tests/simplex/data_out_write_vtk_01.cc +++ b/tests/simplex/data_out_write_vtk_01.cc @@ -17,11 +17,18 @@ // Test DataOut::write_vtk() for simplex meshes. +#include + #include +#include +#include +#include #include +#include #include +#include #include #include #include @@ -29,10 +36,6 @@ #include #include -#include -#include -#include - #include "../tests.h" using namespace dealii; diff --git a/tests/simplex/data_out_write_vtk_02.cc b/tests/simplex/data_out_write_vtk_02.cc index 2e9ddbd96c..d88279770d 100644 --- a/tests/simplex/data_out_write_vtk_02.cc +++ b/tests/simplex/data_out_write_vtk_02.cc @@ -17,12 +17,19 @@ // Test DataOut::write_vtk() for simplex meshes. +#include + #include +#include #include +#include +#include #include +#include #include +#include #include #include #include @@ -30,10 +37,6 @@ #include #include -#include -#include -#include - #include "../tests.h" #include "./simplex_grids.h" diff --git a/tests/simplex/data_out_write_vtu_01.cc b/tests/simplex/data_out_write_vtu_01.cc index aa197b2b20..a58cbf2307 100644 --- a/tests/simplex/data_out_write_vtu_01.cc +++ b/tests/simplex/data_out_write_vtu_01.cc @@ -19,9 +19,14 @@ #include +#include +#include +#include #include +#include #include +#include #include #include #include @@ -29,9 +34,6 @@ #include #include -#include -#include - #include "../tests.h" using namespace dealii; diff --git a/tests/simplex/extract_boundary_dofs.cc b/tests/simplex/extract_boundary_dofs.cc index 711f944524..9329aca91d 100644 --- a/tests/simplex/extract_boundary_dofs.cc +++ b/tests/simplex/extract_boundary_dofs.cc @@ -22,8 +22,12 @@ #include #include -#include -#include +#include +#include +#include +#include + +#include #include "../tests.h" diff --git a/tests/simplex/fe_lib_01.cc b/tests/simplex/fe_lib_01.cc index 9c46e8cbac..50460a2f09 100644 --- a/tests/simplex/fe_lib_01.cc +++ b/tests/simplex/fe_lib_01.cc @@ -17,7 +17,10 @@ // Test n_dofs_per-methods of Simplex::FE_P and Simplex::FE_DGP. -#include +#include +#include +#include +#include #include "../tests.h" diff --git a/tests/simplex/fe_lib_02.cc b/tests/simplex/fe_lib_02.cc index f04b0bd660..bb96eef5f3 100644 --- a/tests/simplex/fe_lib_02.cc +++ b/tests/simplex/fe_lib_02.cc @@ -17,8 +17,12 @@ // Evaluate Simplex::FE_P and Simplex::FE_DGP at quadrature points. -#include -#include +#include + +#include +#include +#include +#include #include "../tests.h" diff --git a/tests/simplex/fe_p_bubbles_01.cc b/tests/simplex/fe_p_bubbles_01.cc index 4c90162db7..c79de2e615 100644 --- a/tests/simplex/fe_p_bubbles_01.cc +++ b/tests/simplex/fe_p_bubbles_01.cc @@ -17,8 +17,12 @@ // Evaluate Simplex::FE_P_Bubbles. -#include -#include +#include + +#include +#include +#include +#include #include "../tests.h" diff --git a/tests/simplex/fe_p_bubbles_02.cc b/tests/simplex/fe_p_bubbles_02.cc index 09e869623f..ca12e66851 100644 --- a/tests/simplex/fe_p_bubbles_02.cc +++ b/tests/simplex/fe_p_bubbles_02.cc @@ -17,10 +17,15 @@ // convergence rate. #include +#include #include +#include +#include +#include #include +#include #include #include @@ -30,9 +35,6 @@ #include #include -#include -#include - #include "../tests.h" using namespace dealii; diff --git a/tests/simplex/get_projection_matrix_01.cc b/tests/simplex/get_projection_matrix_01.cc index 095fee6f5e..327b17589d 100644 --- a/tests/simplex/get_projection_matrix_01.cc +++ b/tests/simplex/get_projection_matrix_01.cc @@ -17,19 +17,23 @@ // Test FETools::get_projection_matrix for simplices. +#include + +#include +#include +#include #include #include +#include #include +#include + #include #include #include -#include -#include -#include - #include "../tests.h" using namespace dealii; diff --git a/tests/simplex/hanging_nodes_01.cc b/tests/simplex/hanging_nodes_01.cc index ac96c5eb37..cd10837641 100644 --- a/tests/simplex/hanging_nodes_01.cc +++ b/tests/simplex/hanging_nodes_01.cc @@ -31,15 +31,17 @@ #include #include +#include +#include +#include +#include + #include #include #include #include -#include -#include - #include "../tests.h" #include "hanging_nodes.h" diff --git a/tests/simplex/hanging_nodes_02.cc b/tests/simplex/hanging_nodes_02.cc index 838092d0f8..3fce2d9bf6 100644 --- a/tests/simplex/hanging_nodes_02.cc +++ b/tests/simplex/hanging_nodes_02.cc @@ -31,6 +31,11 @@ #include #include +#include +#include +#include +#include + #include #include #include @@ -39,9 +44,6 @@ #include -#include -#include - #include "../tests.h" #include "hanging_nodes.h" diff --git a/tests/simplex/hanging_nodes_03.cc b/tests/simplex/hanging_nodes_03.cc index c0e802b554..120aecb5a5 100644 --- a/tests/simplex/hanging_nodes_03.cc +++ b/tests/simplex/hanging_nodes_03.cc @@ -31,6 +31,11 @@ #include #include +#include +#include +#include +#include + #include #include #include @@ -39,9 +44,6 @@ #include -#include -#include - #include "../tests.h" #include "hanging_nodes.h" diff --git a/tests/simplex/hanging_nodes_hybrid_01.cc b/tests/simplex/hanging_nodes_hybrid_01.cc index 6a98030110..27368b8f35 100644 --- a/tests/simplex/hanging_nodes_hybrid_01.cc +++ b/tests/simplex/hanging_nodes_hybrid_01.cc @@ -31,12 +31,14 @@ */ +#include #include +#include +#include +#include #include -#include - #include #include "../tests.h" diff --git a/tests/simplex/hanging_nodes_hybrid_02.cc b/tests/simplex/hanging_nodes_hybrid_02.cc index cf7dc743d0..7e2bd1d6f1 100644 --- a/tests/simplex/hanging_nodes_hybrid_02.cc +++ b/tests/simplex/hanging_nodes_hybrid_02.cc @@ -31,12 +31,14 @@ */ +#include #include +#include +#include +#include #include -#include - #include #include "../tests.h" diff --git a/tests/simplex/hanging_nodes_hybrid_03.cc b/tests/simplex/hanging_nodes_hybrid_03.cc index 63440a7801..4c56033805 100644 --- a/tests/simplex/hanging_nodes_hybrid_03.cc +++ b/tests/simplex/hanging_nodes_hybrid_03.cc @@ -18,12 +18,14 @@ // verify hanging node constraints on locally hp-refined hybrid meshes +#include #include +#include +#include +#include #include -#include - #include #include "../tests.h" diff --git a/tests/simplex/mapping_fe_01.cc b/tests/simplex/mapping_fe_01.cc index 6bcb8ce79f..dc23a7bc0c 100644 --- a/tests/simplex/mapping_fe_01.cc +++ b/tests/simplex/mapping_fe_01.cc @@ -16,9 +16,15 @@ // Distribute Simplex::FE_Wedge on a DoFHandler. +#include + #include +#include +#include +#include #include +#include #include #include @@ -26,9 +32,6 @@ #include -#include -#include - #include "../tests.h" using namespace dealii; diff --git a/tests/simplex/mapping_fe_fields_01.cc b/tests/simplex/mapping_fe_fields_01.cc index 96577c03fe..4bc83fa77b 100644 --- a/tests/simplex/mapping_fe_fields_01.cc +++ b/tests/simplex/mapping_fe_fields_01.cc @@ -17,22 +17,25 @@ // Test MappingFEField and VectorTools::get_position_vector() for simplex // meshes. +#include + #include +#include +#include +#include #include #include +#include #include #include +#include #include #include #include -#include -#include -#include - #include "../tests.h" using namespace dealii; diff --git a/tests/simplex/mapping_fe_fields_02.cc b/tests/simplex/mapping_fe_fields_02.cc index 9580d402a7..3e9051d4d4 100644 --- a/tests/simplex/mapping_fe_fields_02.cc +++ b/tests/simplex/mapping_fe_fields_02.cc @@ -17,22 +17,25 @@ // Like mapping_fe_fields_01 but for deformed meshes tested for linear and // quadratic mapping. +#include + #include +#include +#include +#include #include #include +#include #include #include +#include #include #include #include -#include -#include -#include - #include "../tests.h" using namespace dealii; diff --git a/tests/simplex/mapping_transformations_01.cc b/tests/simplex/mapping_transformations_01.cc index 830fa8a4e1..970dce5548 100644 --- a/tests/simplex/mapping_transformations_01.cc +++ b/tests/simplex/mapping_transformations_01.cc @@ -18,6 +18,10 @@ #include +#include +#include +#include +#include #include #include @@ -25,8 +29,6 @@ #include #include -#include - #include #include "../tests.h" diff --git a/tests/simplex/matrix_free_01.cc b/tests/simplex/matrix_free_01.cc index daeeb53115..1e607dca78 100644 --- a/tests/simplex/matrix_free_01.cc +++ b/tests/simplex/matrix_free_01.cc @@ -18,10 +18,16 @@ // continuous elements and compare results between matrix-free and matrix-based // implementations. +#include + #include #include +#include #include +#include +#include +#include #include #include @@ -39,10 +45,6 @@ #include #include -#include -#include -#include - #include "../tests.h" #include "./simplex_grids.h" diff --git a/tests/simplex/matrix_free_02.cc b/tests/simplex/matrix_free_02.cc index 9f128749ec..d73069d0a0 100644 --- a/tests/simplex/matrix_free_02.cc +++ b/tests/simplex/matrix_free_02.cc @@ -17,12 +17,18 @@ // Same as matrix_free_01 but testing mixed meshes (and also pure simplex and // hypercube mesh as special case of mixed meshs). +#include + #include #include #include +#include #include +#include +#include +#include #include #include @@ -41,10 +47,6 @@ #include #include -#include -#include -#include - #include "../tests.h" #include "./simplex_grids.h" diff --git a/tests/simplex/matrix_free_03.cc b/tests/simplex/matrix_free_03.cc index 5eee6625c6..5e9e00f137 100644 --- a/tests/simplex/matrix_free_03.cc +++ b/tests/simplex/matrix_free_03.cc @@ -18,6 +18,8 @@ // continuous elements and compare results between matrix-free and matrix-based // implementations. +#include + #include #include @@ -25,7 +27,11 @@ #include #include +#include #include +#include +#include +#include #include #include @@ -47,10 +53,6 @@ #include #include -#include -#include -#include - #include "../tests.h" using namespace dealii; diff --git a/tests/simplex/matrix_free_04.cc b/tests/simplex/matrix_free_04.cc index 5d613014ae..52e750f4ef 100644 --- a/tests/simplex/matrix_free_04.cc +++ b/tests/simplex/matrix_free_04.cc @@ -16,6 +16,8 @@ // Solve Poisson problem problem on a mixed mesh with DG and MatrixFree. +#include + #include #include @@ -23,7 +25,11 @@ #include #include +#include #include +#include +#include +#include #include #include @@ -46,10 +52,6 @@ #include #include -#include -#include -#include - #include "../tests.h" #include "./simplex_grids.h" diff --git a/tests/simplex/matrix_free_shape_info_01.cc b/tests/simplex/matrix_free_shape_info_01.cc index cfe1180855..126faf65d3 100644 --- a/tests/simplex/matrix_free_shape_info_01.cc +++ b/tests/simplex/matrix_free_shape_info_01.cc @@ -17,21 +17,24 @@ // Test ShapeData for Simplex::FE_P and Simplex::QGauss +#include + #include +#include +#include +#include #include +#include #include +#include #include #include #include -#include -#include -#include - #include "../tests.h" using namespace dealii; diff --git a/tests/simplex/poisson_01.cc b/tests/simplex/poisson_01.cc index 56a2020e41..2828aa1882 100644 --- a/tests/simplex/poisson_01.cc +++ b/tests/simplex/poisson_01.cc @@ -19,6 +19,7 @@ #include #include +#include #include #include @@ -27,8 +28,12 @@ #include #include +#include #include +#include +#include #include +#include #include #include @@ -53,10 +58,6 @@ #include #include -#include -#include -#include - #include "../tests.h" #include "simplex_grids.h" diff --git a/tests/simplex/poisson_02.cc b/tests/simplex/poisson_02.cc index c117671b3d..f2906722bd 100644 --- a/tests/simplex/poisson_02.cc +++ b/tests/simplex/poisson_02.cc @@ -27,8 +27,12 @@ #include #include +#include #include +#include +#include #include +#include #include #include @@ -51,10 +55,6 @@ #include #include -#include -#include -#include - #include #include diff --git a/tests/simplex/polynomials_01.cc b/tests/simplex/polynomials_01.cc index c501094a14..96ffeb06c0 100644 --- a/tests/simplex/polynomials_01.cc +++ b/tests/simplex/polynomials_01.cc @@ -18,10 +18,9 @@ // quadrature rule. +#include #include -#include - #include "../tests.h" using namespace dealii; diff --git a/tests/simplex/q_projection_01.cc b/tests/simplex/q_projection_01.cc index 7d6f6a1065..60f4cf89a5 100644 --- a/tests/simplex/q_projection_01.cc +++ b/tests/simplex/q_projection_01.cc @@ -18,8 +18,7 @@ #include - -#include +#include #include "../tests.h" diff --git a/tests/simplex/q_witherden_vincent_01.cc b/tests/simplex/q_witherden_vincent_01.cc index 67010b0a31..b521142c98 100644 --- a/tests/simplex/q_witherden_vincent_01.cc +++ b/tests/simplex/q_witherden_vincent_01.cc @@ -14,8 +14,7 @@ // --------------------------------------------------------------------- #include - -#include +#include #include "../tests.h" diff --git a/tests/simplex/quadrature_lib_01.cc b/tests/simplex/quadrature_lib_01.cc index fa05a42fd0..a3e417b5a7 100644 --- a/tests/simplex/quadrature_lib_01.cc +++ b/tests/simplex/quadrature_lib_01.cc @@ -17,7 +17,7 @@ // Test Simplex::QGauss: output its quadrature points and weights. -#include +#include #include "../tests.h" diff --git a/tests/simplex/reference_cell_kind_01.cc b/tests/simplex/reference_cell_kind_01.cc index 0d3e7fe611..7885dc9b34 100644 --- a/tests/simplex/reference_cell_kind_01.cc +++ b/tests/simplex/reference_cell_kind_01.cc @@ -17,9 +17,9 @@ // Test ReferenceCell::Kind::faces_for_given_vertex(). -#include +#include -#include +#include #include "../tests.h" diff --git a/tests/simplex/refinement_01.cc b/tests/simplex/refinement_01.cc index 5346fe9a19..aa52255570 100644 --- a/tests/simplex/refinement_01.cc +++ b/tests/simplex/refinement_01.cc @@ -17,11 +17,10 @@ // Test refinement of simplex mesh. +#include #include #include -#include - #include "../tests.h" #include "./simplex_grids.h" diff --git a/tests/simplex/refinement_02.cc b/tests/simplex/refinement_02.cc index d74bfe7339..984c499cd4 100644 --- a/tests/simplex/refinement_02.cc +++ b/tests/simplex/refinement_02.cc @@ -17,11 +17,10 @@ // Test the correct setup of neighbors during refinement of simplex mesh. +#include #include #include -#include - #include "../tests.h" #include "./simplex_grids.h" diff --git a/tests/simplex/refinement_03.cc b/tests/simplex/refinement_03.cc index 5236adfd48..b72cd17947 100644 --- a/tests/simplex/refinement_03.cc +++ b/tests/simplex/refinement_03.cc @@ -18,10 +18,9 @@ // Verify CellAccessor::neighbor_child_on_subface for a triangle mesh. +#include #include -#include - #include "../tests.h" diff --git a/tests/simplex/step-01.cc b/tests/simplex/step-01.cc index b9a41e2418..5ce2637c7b 100644 --- a/tests/simplex/step-01.cc +++ b/tests/simplex/step-01.cc @@ -27,8 +27,6 @@ #include #include -#include - #include #include #include diff --git a/tests/simplex/step-02.cc b/tests/simplex/step-02.cc index 32cf8295d4..7798f2e45c 100644 --- a/tests/simplex/step-02.cc +++ b/tests/simplex/step-02.cc @@ -26,6 +26,11 @@ #include #include +#include +#include +#include +#include + #include #include #include @@ -34,8 +39,6 @@ #include #include -#include - #include #include "../tests.h" diff --git a/tests/simplex/step-03.cc b/tests/simplex/step-03.cc index ae8a03e29d..6924433f29 100644 --- a/tests/simplex/step-03.cc +++ b/tests/simplex/step-03.cc @@ -29,8 +29,12 @@ #include #include +#include #include +#include +#include #include +#include #include #include @@ -51,10 +55,6 @@ #include #include -#include -#include -#include - #include #include diff --git a/tests/simplex/step-04.cc b/tests/simplex/step-04.cc index 892c1efb04..60d0b14681 100644 --- a/tests/simplex/step-04.cc +++ b/tests/simplex/step-04.cc @@ -35,8 +35,12 @@ #include #include +#include #include +#include +#include #include +#include #include #include @@ -55,10 +59,6 @@ #include #include -#include -#include -#include - #include #include diff --git a/tests/simplex/step-06.cc b/tests/simplex/step-06.cc index d355951591..d13a24d4d9 100644 --- a/tests/simplex/step-06.cc +++ b/tests/simplex/step-06.cc @@ -30,10 +30,13 @@ #ifdef USE_SIMPLEX +# include # include -# include -# include +# include +# include +# include +# include #else # include diff --git a/tests/simplex/step-07.cc b/tests/simplex/step-07.cc index 4cbbd6a9fa..2ffa48a28c 100644 --- a/tests/simplex/step-07.cc +++ b/tests/simplex/step-07.cc @@ -35,8 +35,12 @@ #include #include +#include #include +#include +#include #include +#include #include #include @@ -58,9 +62,6 @@ #include #include -#include -#include - #include #include #include diff --git a/tests/simplex/step-08.cc b/tests/simplex/step-08.cc index eececbf1b4..483a891611 100644 --- a/tests/simplex/step-08.cc +++ b/tests/simplex/step-08.cc @@ -31,9 +31,13 @@ #include #include +#include #include +#include +#include #include #include +#include #include #include @@ -55,10 +59,6 @@ #include #include -#include -#include -#include - #include #include diff --git a/tests/simplex/step-12.cc b/tests/simplex/step-12.cc index 94d64a0856..5e8936a356 100644 --- a/tests/simplex/step-12.cc +++ b/tests/simplex/step-12.cc @@ -31,7 +31,11 @@ #include #include +#include +#include +#include #include +#include #include #include @@ -54,10 +58,6 @@ #include #include -#include -#include -#include - #include #include diff --git a/tests/simplex/step-12a.cc b/tests/simplex/step-12a.cc index 29484a3aa2..40df59beef 100644 --- a/tests/simplex/step-12a.cc +++ b/tests/simplex/step-12a.cc @@ -49,16 +49,20 @@ // Finally, the new include file for using the mesh_loop from the MeshWorker // framework +#include + +#include +#include +#include +#include #include +#include + #include #include -#include -#include -#include - #include #include diff --git a/tests/simplex/step-18.cc b/tests/simplex/step-18.cc index f2459a1000..25cd897699 100644 --- a/tests/simplex/step-18.cc +++ b/tests/simplex/step-18.cc @@ -70,11 +70,15 @@ #include "../tests.h" // simplex +#include + +#include +#include +#include +#include #include -#include -#include -#include +#include //#define HEX diff --git a/tests/simplex/step-20.cc b/tests/simplex/step-20.cc index 34e482fe7c..c16551e69f 100644 --- a/tests/simplex/step-20.cc +++ b/tests/simplex/step-20.cc @@ -35,8 +35,12 @@ #include #include +#include +#include +#include #include #include +#include #include #include @@ -57,10 +61,6 @@ #include #include -#include -#include -#include - #include #include diff --git a/tests/simplex/step-23.cc b/tests/simplex/step-23.cc index 1a8dee444d..52af8814de 100644 --- a/tests/simplex/step-23.cc +++ b/tests/simplex/step-23.cc @@ -50,11 +50,15 @@ #include "../tests.h" // simplex +#include + +#include +#include +#include +#include #include -#include -#include -#include +#include // #define HEX diff --git a/tests/simplex/step-38.cc b/tests/simplex/step-38.cc index cda26794bb..e47803bce1 100644 --- a/tests/simplex/step-38.cc +++ b/tests/simplex/step-38.cc @@ -46,9 +46,14 @@ #include // simplex -#include -#include -#include +#include + +#include +#include +#include +#include + +#include #include #include diff --git a/tests/simplex/step-55.cc b/tests/simplex/step-55.cc index 76561d44bc..c657613c0c 100644 --- a/tests/simplex/step-55.cc +++ b/tests/simplex/step-55.cc @@ -86,12 +86,15 @@ namespace LA #ifdef HEX # include #else +# include + # include +# include +# include +# include +# include # include - -# include -# include #endif diff --git a/tests/simplex/step-68.cc b/tests/simplex/step-68.cc index 8bfbc11d63..861b224e97 100644 --- a/tests/simplex/step-68.cc +++ b/tests/simplex/step-68.cc @@ -30,6 +30,7 @@ #include #include #include +#include #include #include @@ -40,8 +41,12 @@ #include #include +#include #include +#include +#include #include +#include #include #include @@ -59,10 +64,6 @@ #include #include -#include -#include -#include - #include #include diff --git a/tests/simplex/step-74.cc b/tests/simplex/step-74.cc index 7c9ce0553c..340cd3c42a 100644 --- a/tests/simplex/step-74.cc +++ b/tests/simplex/step-74.cc @@ -32,7 +32,11 @@ #include #include +#include +#include +#include #include +#include #include #include @@ -56,9 +60,6 @@ #include #include -#include -#include - #include "../tests.h" #define USE_SIMPLEX diff --git a/tests/simplex/unit_tangential_vectors_01.cc b/tests/simplex/unit_tangential_vectors_01.cc index 61fec47d1d..d72e26a48a 100644 --- a/tests/simplex/unit_tangential_vectors_01.cc +++ b/tests/simplex/unit_tangential_vectors_01.cc @@ -18,9 +18,9 @@ // for all reference-cell types. -#include +#include -#include +#include #include "../tests.h" diff --git a/tests/simplex/variable_face_quadratures_03.cc b/tests/simplex/variable_face_quadratures_03.cc index 95855659f9..6560f1515a 100644 --- a/tests/simplex/variable_face_quadratures_03.cc +++ b/tests/simplex/variable_face_quadratures_03.cc @@ -18,19 +18,22 @@ #include +#include +#include +#include +#include #include #include +#include #include +#include + #include #include -#include -#include -#include - #include "../tests.h" using namespace dealii; diff --git a/tests/simplex/wedge_01.cc b/tests/simplex/wedge_01.cc index 5629248bee..62ecb8ddfb 100644 --- a/tests/simplex/wedge_01.cc +++ b/tests/simplex/wedge_01.cc @@ -18,11 +18,14 @@ #include +#include +#include +#include +#include + #include #include -#include - #include "../tests.h" using namespace dealii; -- 2.39.5