From: Peter Munch Date: Sat, 5 Dec 2020 12:00:33 +0000 (+0100) Subject: Implement ScalarPyramidPolynomial and FE_PyramidP X-Git-Tag: v9.3.0-rc1~781^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=953765f7af66e5634452301bee731591d014b493;p=dealii.git Implement ScalarPyramidPolynomial and FE_PyramidP --- diff --git a/include/deal.II/simplex/fe_lib.h b/include/deal.II/simplex/fe_lib.h index e7c45a8ed9..030b89d526 100644 --- a/include/deal.II/simplex/fe_lib.h +++ b/include/deal.II/simplex/fe_lib.h @@ -222,6 +222,66 @@ namespace Simplex }; + /** + * Implementation of a scalar Lagrange finite element on a pyramid that yields + * the finite element space of continuous, piecewise polynomials of + * degree p. + * + * @ingroup simplex + */ + template + class FE_PyramidP : public dealii::FE_Poly + { + 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; + }; + + } // namespace Simplex DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/simplex/polynomials.h b/include/deal.II/simplex/polynomials.h index a1c484aaa5..051877bdc7 100644 --- a/include/deal.II/simplex/polynomials.h +++ b/include/deal.II/simplex/polynomials.h @@ -266,6 +266,119 @@ 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() + * + * @note Not implemented yet. + */ + 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 functions template template @@ -303,6 +416,25 @@ 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/source/simplex/fe_lib.cc b/source/simplex/fe_lib.cc index abd6e59042..9312656e3c 100644 --- a/source/simplex/fe_lib.cc +++ b/source/simplex/fe_lib.cc @@ -115,6 +115,32 @@ 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_pyramid(const unsigned int degree) + { + internal::GenericDoFsPerObject dpo; + + // all dofs are internal + 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; + } } // namespace @@ -559,6 +585,153 @@ namespace Simplex + template + FE_PyramidP::FE_PyramidP(const unsigned int degree) + : dealii::FE_Poly( + Simplex::ScalarPyramidPolynomial(degree), + FiniteElementData(get_dpo_vector_fe_pyramid(degree), + ReferenceCell::Type::Pyramid, + 1, + degree, + FiniteElementData::H1), + std::vector(FiniteElementData(get_dpo_vector_fe_pyramid( + degree), + ReferenceCell::Type::Pyramid, + 1, + degree) + .dofs_per_cell, + true), + std::vector( + FiniteElementData(get_dpo_vector_fe_pyramid(degree), + ReferenceCell::Type::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 + 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 + { + (void)fe_other; + (void)codim; + + Assert((dynamic_cast *>(&fe_other)) || + (dynamic_cast *>(&fe_other)), + ExcNotImplemented()); + + return FiniteElementDomination::this_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; + } + + + } // namespace Simplex // explicit instantiations diff --git a/source/simplex/fe_lib.inst.in b/source/simplex/fe_lib.inst.in index 160d876f95..36b7279cd3 100644 --- a/source/simplex/fe_lib.inst.in +++ b/source/simplex/fe_lib.inst.in @@ -23,5 +23,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS) template class Simplex::FE_DGP; template class Simplex::FE_WedgeP; + template class Simplex::FE_PyramidP; #endif } diff --git a/source/simplex/polynomials.cc b/source/simplex/polynomials.cc index bae98e5f95..dc9ec63a43 100644 --- a/source/simplex/polynomials.cc +++ b/source/simplex/polynomials.cc @@ -52,6 +52,21 @@ namespace Simplex return 0; } + unsigned int + compute_n_polynomials_pyramid(const unsigned int dim, + const unsigned int degree) + { + if (dim == 3) + { + if (degree == 1) + return 5; + } + + Assert(false, ExcNotImplemented()); + + return 0; + } + unsigned int compute_n_polynomials_wedge(const unsigned int dim, const unsigned int degree) @@ -711,12 +726,255 @@ namespace Simplex + template + ScalarPyramidPolynomial::ScalarPyramidPolynomial( + const unsigned int degree) + : ScalarPolynomialsBase(degree, + compute_n_polynomials_pyramid(dim, degree)) + {} + + + template + double + ScalarPyramidPolynomial::compute_value(const unsigned int i, + const Point & p) const + { + AssertDimension(dim, 3); + AssertIndexRange(this->degree(), 2); + + const double Q14 = 0.25; + double ration; + + const double r = p[0]; + const double s = p[1]; + const double t = p[2]; + + if (fabs(t - 1.0) > 1.0e-14) + { + ration = (r * s * t) / (1.0 - t); + } + else + { + ration = 0.0; + } + + if (i == 0) + return Q14 * ((1.0 - r) * (1.0 - s) - t + ration); + if (i == 1) + return Q14 * ((1.0 + r) * (1.0 - s) - t - ration); + if (i == 2) + return Q14 * ((1.0 - r) * (1.0 + s) - t - ration); + if (i == 3) + return Q14 * ((1.0 + r) * (1.0 + s) - t + ration); + else + return t; + } + + + + template + Tensor<1, dim> + ScalarPyramidPolynomial::compute_grad(const unsigned int i, + const Point & p) const + { + AssertDimension(dim, 3); + AssertIndexRange(this->degree(), 4); + + Tensor<1, dim> grad; + + if (this->degree() == 1) + { + const double Q14 = 0.25; + + const double r = p[0]; + const double s = p[1]; + const double t = p[2]; + + double rationdr; + double rationds; + double rationdt; + + if (fabs(t - 1.0) > 1.0e-14) + { + rationdr = s * t / (1.0 - t); + rationds = r * t / (1.0 - t); + rationdt = r * s / ((1.0 - t) * (1.0 - t)); + } + else + { + rationdr = 1.0; + rationds = 1.0; + rationdt = 1.0; + } + + + if (i == 0) + { + grad[0] = Q14 * (-1.0 * (1.0 - s) + rationdr); + grad[1] = Q14 * (-1.0 * (1.0 - r) + rationds); + grad[2] = Q14 * (rationdt - 1.0); + } + else if (i == 1) + { + grad[0] = Q14 * (1.0 * (1.0 - s) - rationdr); + grad[1] = Q14 * (-1.0 * (1.0 + r) - rationds); + grad[2] = Q14 * (-1.0 * rationdt - 1.0); + } + else if (i == 2) + { + grad[0] = Q14 * (-1.0 * (1.0 + s) - rationdr); + grad[1] = Q14 * (1.0 * (1.0 - r) - rationds); + grad[2] = Q14 * (-1.0 * rationdt - 1.0); + } + else if (i == 3) + { + grad[0] = Q14 * (1.0 * (1.0 + s) + rationdr); + grad[1] = Q14 * (1.0 * (1.0 + r) + rationds); + grad[2] = Q14 * (rationdt - 1.0); + } + else if (i == 4) + { + grad[0] = 0.0; + grad[1] = 0.0; + grad[2] = 1.0; + } + else + { + Assert(false, ExcNotImplemented()); + } + } + + return grad; + } + + + + template + Tensor<2, dim> + ScalarPyramidPolynomial::compute_grad_grad(const unsigned int i, + const Point & p) const + { + (void)i; + (void)p; + + Assert(false, ExcNotImplemented()); + return Tensor<2, dim>(); + } + + + + template + void + ScalarPyramidPolynomial::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> + ScalarPyramidPolynomial::compute_1st_derivative( + const unsigned int i, + const Point & p) const + { + return compute_grad(i, p); + } + + + + template + Tensor<2, dim> + ScalarPyramidPolynomial::compute_2nd_derivative( + const unsigned int i, + const Point & p) const + { + (void)i; + (void)p; + + Assert(false, ExcNotImplemented()); + + return {}; + } + + + + template + Tensor<3, dim> + ScalarPyramidPolynomial::compute_3rd_derivative( + const unsigned int i, + const Point & p) const + { + (void)i; + (void)p; + + Assert(false, ExcNotImplemented()); + + return {}; + } + + + + template + Tensor<4, dim> + ScalarPyramidPolynomial::compute_4th_derivative( + const unsigned int i, + const Point & p) const + { + (void)i; + (void)p; + + Assert(false, ExcNotImplemented()); + + return {}; + } + + + + template + std::string + ScalarPyramidPolynomial::name() const + { + return "ScalarPyramidPolynomial"; + } + + + + template + std::unique_ptr> + ScalarPyramidPolynomial::clone() const + { + return std::make_unique>(*this); + } + + + template class ScalarPolynomial<1>; template class ScalarPolynomial<2>; template class ScalarPolynomial<3>; template class ScalarWedgePolynomial<1>; template class ScalarWedgePolynomial<2>; template class ScalarWedgePolynomial<3>; + template class ScalarPyramidPolynomial<1>; + template class ScalarPyramidPolynomial<2>; + template class ScalarPyramidPolynomial<3>; } // namespace Simplex