};
+ /**
+ * 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 <int dim, int spacedim = dim>
+ class FE_PyramidP : public dealii::FE_Poly<dim, spacedim>
+ {
+ public:
+ /**
+ * Constructor.
+ */
+ FE_PyramidP(const unsigned int degree);
+
+ /**
+ * @copydoc dealii::FiniteElement::clone()
+ */
+ std::unique_ptr<FiniteElement<dim, spacedim>>
+ clone() const override;
+
+ /**
+ * Return a string that uniquely identifies a finite element. This class
+ * returns <tt>Simplex::FE_PyramidP<dim>(degree)</tt>, 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<dim, spacedim> &fe_other,
+ const unsigned int codim) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::hp_vertex_dof_identities()
+ */
+ std::vector<std::pair<unsigned int, unsigned int>>
+ hp_vertex_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::hp_line_dof_identities()
+ */
+ std::vector<std::pair<unsigned int, unsigned int>>
+ hp_line_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::hp_quad_dof_identities()
+ */
+ std::vector<std::pair<unsigned int, unsigned int>>
+ hp_quad_dof_identities(const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int face_no = 0) const override;
+ };
+
+
} // namespace Simplex
DEAL_II_NAMESPACE_CLOSE
+ /**
+ * Polynomials defined on pyramid entities. This class is basis of
+ * Simplex::FE_PyramidP.
+ */
+ template <int dim>
+ class ScalarPyramidPolynomial : public ScalarPolynomialsBase<dim>
+ {
+ 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<dim> & unit_point,
+ std::vector<double> & values,
+ std::vector<Tensor<1, dim>> &grads,
+ std::vector<Tensor<2, dim>> &grad_grads,
+ std::vector<Tensor<3, dim>> &third_derivatives,
+ std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_value()
+ */
+ double
+ compute_value(const unsigned int i, const Point<dim> &p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_derivative()
+ *
+ * @note Currently, only implemented for first derivative.
+ */
+ template <int order>
+ Tensor<order, dim>
+ compute_derivative(const unsigned int i, const Point<dim> &p) const;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_1st_derivative()
+ */
+ Tensor<1, dim>
+ compute_1st_derivative(const unsigned int i,
+ const Point<dim> & p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_2nd_derivative()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<2, dim>
+ compute_2nd_derivative(const unsigned int i,
+ const Point<dim> & p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_3rd_derivative()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<3, dim>
+ compute_3rd_derivative(const unsigned int i,
+ const Point<dim> & p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_4th_derivative()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<4, dim>
+ compute_4th_derivative(const unsigned int i,
+ const Point<dim> & p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_grad()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<1, dim>
+ compute_grad(const unsigned int i, const Point<dim> &p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_grad_grad()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<2, dim>
+ compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::name()
+ */
+ std::string
+ name() const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::clone()
+ */
+ virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
+ clone() const override;
+ };
+
+
+
// template functions
template <int dim>
template <int order>
return der;
}
+
+
+ template <int dim>
+ template <int order>
+ Tensor<order, dim>
+ ScalarPyramidPolynomial<dim>::compute_derivative(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ Tensor<order, dim> 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
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
+ template <int dim, int spacedim>
+ FE_PyramidP<dim, spacedim>::FE_PyramidP(const unsigned int degree)
+ : dealii::FE_Poly<dim, spacedim>(
+ Simplex::ScalarPyramidPolynomial<dim>(degree),
+ FiniteElementData<dim>(get_dpo_vector_fe_pyramid(degree),
+ ReferenceCell::Type::Pyramid,
+ 1,
+ degree,
+ FiniteElementData<dim>::H1),
+ std::vector<bool>(FiniteElementData<dim>(get_dpo_vector_fe_pyramid(
+ degree),
+ ReferenceCell::Type::Pyramid,
+ 1,
+ degree)
+ .dofs_per_cell,
+ true),
+ std::vector<ComponentMask>(
+ FiniteElementData<dim>(get_dpo_vector_fe_pyramid(degree),
+ ReferenceCell::Type::Pyramid,
+ 1,
+ degree)
+ .dofs_per_cell,
+ std::vector<bool>(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 <int dim, int spacedim>
+ std::unique_ptr<FiniteElement<dim, spacedim>>
+ FE_PyramidP<dim, spacedim>::clone() const
+ {
+ return std::make_unique<FE_PyramidP<dim, spacedim>>(*this);
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::string
+ FE_PyramidP<dim, spacedim>::get_name() const
+ {
+ std::ostringstream namebuf;
+ namebuf << "FE_PyramidP<" << dim << ">(" << this->degree << ")";
+
+ return namebuf.str();
+ }
+
+
+
+ template <int dim, int spacedim>
+ FiniteElementDomination::Domination
+ FE_PyramidP<dim, spacedim>::compare_for_domination(
+ const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int codim) const
+ {
+ (void)fe_other;
+ (void)codim;
+
+ Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)) ||
+ (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
+ ExcNotImplemented());
+
+ return FiniteElementDomination::this_element_dominates;
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::vector<std::pair<unsigned int, unsigned int>>
+ FE_PyramidP<dim, spacedim>::hp_vertex_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const
+ {
+ (void)fe_other;
+
+ Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)) ||
+ (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
+ ExcNotImplemented());
+
+ return {{0, 0}};
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::vector<std::pair<unsigned int, unsigned int>>
+ FE_PyramidP<dim, spacedim>::hp_line_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const
+ {
+ (void)fe_other;
+
+ Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)) ||
+ (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
+ ExcNotImplemented());
+
+ std::vector<std::pair<unsigned int, unsigned int>> result;
+
+ for (unsigned int i = 0; i < this->degree - 1; ++i)
+ result.emplace_back(i, i);
+
+ return result;
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::vector<std::pair<unsigned int, unsigned int>>
+ FE_PyramidP<dim, spacedim>::hp_quad_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int face_no) const
+ {
+ (void)fe_other;
+
+
+ AssertIndexRange(face_no, 5);
+
+ if (face_no == 0)
+ {
+ Assert((dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
+ ExcNotImplemented());
+ }
+ else
+ {
+ Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)),
+ ExcNotImplemented());
+ }
+
+ std::vector<std::pair<unsigned int, unsigned int>> 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
template class Simplex::FE_DGP<deal_II_dimension, deal_II_space_dimension>;
template class Simplex::FE_WedgeP<deal_II_dimension,
deal_II_space_dimension>;
+ template class Simplex::FE_PyramidP<deal_II_dimension,
+ deal_II_space_dimension>;
#endif
}
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)
+ template <int dim>
+ ScalarPyramidPolynomial<dim>::ScalarPyramidPolynomial(
+ const unsigned int degree)
+ : ScalarPolynomialsBase<dim>(degree,
+ compute_n_polynomials_pyramid(dim, degree))
+ {}
+
+
+ template <int dim>
+ double
+ ScalarPyramidPolynomial<dim>::compute_value(const unsigned int i,
+ const Point<dim> & 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 <int dim>
+ Tensor<1, dim>
+ ScalarPyramidPolynomial<dim>::compute_grad(const unsigned int i,
+ const Point<dim> & 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 <int dim>
+ Tensor<2, dim>
+ ScalarPyramidPolynomial<dim>::compute_grad_grad(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ (void)i;
+ (void)p;
+
+ Assert(false, ExcNotImplemented());
+ return Tensor<2, dim>();
+ }
+
+
+
+ template <int dim>
+ void
+ ScalarPyramidPolynomial<dim>::evaluate(
+ const Point<dim> & unit_point,
+ std::vector<double> & values,
+ std::vector<Tensor<1, dim>> &grads,
+ std::vector<Tensor<2, dim>> &grad_grads,
+ std::vector<Tensor<3, dim>> &third_derivatives,
+ std::vector<Tensor<4, dim>> &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 <int dim>
+ Tensor<1, dim>
+ ScalarPyramidPolynomial<dim>::compute_1st_derivative(
+ const unsigned int i,
+ const Point<dim> & p) const
+ {
+ return compute_grad(i, p);
+ }
+
+
+
+ template <int dim>
+ Tensor<2, dim>
+ ScalarPyramidPolynomial<dim>::compute_2nd_derivative(
+ const unsigned int i,
+ const Point<dim> & p) const
+ {
+ (void)i;
+ (void)p;
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+
+
+ template <int dim>
+ Tensor<3, dim>
+ ScalarPyramidPolynomial<dim>::compute_3rd_derivative(
+ const unsigned int i,
+ const Point<dim> & p) const
+ {
+ (void)i;
+ (void)p;
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+
+
+ template <int dim>
+ Tensor<4, dim>
+ ScalarPyramidPolynomial<dim>::compute_4th_derivative(
+ const unsigned int i,
+ const Point<dim> & p) const
+ {
+ (void)i;
+ (void)p;
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+
+
+ template <int dim>
+ std::string
+ ScalarPyramidPolynomial<dim>::name() const
+ {
+ return "ScalarPyramidPolynomial";
+ }
+
+
+
+ template <int dim>
+ std::unique_ptr<ScalarPolynomialsBase<dim>>
+ ScalarPyramidPolynomial<dim>::clone() const
+ {
+ return std::make_unique<ScalarPyramidPolynomial<dim>>(*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
projects / dealii.git / commitdiff