*/
namespace Simplex
{
- /**
- * Polynomials defined on dim-dimensional simplex entities. This class is
- * basis of Simplex::FE_P.
- *
- * @ingroup simplex
- */
- template <int dim>
- class ScalarPolynomial : 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 (degree=1) and quadratic polynomials
- * (degree=2) are implemented.
- */
- ScalarPolynomial(const unsigned int degree);
-
- /**
- * @copydoc ScalarPolynomialsBase::evaluate()
- *
- * @note Currently, only the vectors @p values, @p grads, and @p grad_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 and second 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()
- */
- 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;
- };
-
-
-
/**
* Polynomials defined on wedge entities. This class is basis of
* Simplex::FE_WedgeP.
*
* The polynomials are created via a tensor product of a
- * Simplex::ScalarPolynomial<2>(degree) and a
- * Simplex::ScalarPolynomial<1>(degree), however, are re-numerated to better
- * match the definition of FiniteElement.
+ * Simplex::BarycentricPolynomials<2>::get_fe_p_basis(degree) and a
+ * Simplex::BarycentricPolynomials<1>::get_fe_p_basis(degree), however, are
+ * re-numerated to better match the definition of FiniteElement.
*/
template <int dim>
class ScalarWedgePolynomial : public ScalarPolynomialsBase<dim>
- // template functions
- template <int dim>
- template <int order>
- Tensor<order, dim>
- ScalarPolynomial<dim>::compute_derivative(const unsigned int i,
- const Point<dim> & p) const
- {
- Tensor<order, dim> derivative;
-
- if (order == 1)
- {
- Tensor<1, dim> &derivative_1 =
- *reinterpret_cast<Tensor<1, dim> *>(&derivative);
-
- const auto grad = compute_grad(i, p);
- for (unsigned int i = 0; i < dim; ++i)
- derivative_1[i] = grad[i];
- }
- else if (order == 2)
- {
- Tensor<2, dim> &derivative_2 =
- *reinterpret_cast<Tensor<2, dim> *>(&derivative);
-
- const auto grad_grad = compute_grad_grad(i, p);
-
- for (unsigned int i = 0; i < dim; ++i)
- for (unsigned int j = 0; j < dim; ++j)
- derivative_2[i][j] = grad_grad[i][j];
- }
- else
- {
- Assert(false, ExcNotImplemented());
- }
-
- return derivative;
- }
-
-
-
template <int dim>
template <int order>
Tensor<order, dim>
{
namespace
{
- unsigned int
- compute_n_polynomials(const unsigned int dim, const unsigned int degree)
- {
- if (dim == 1)
- {
- if (degree == 1)
- return 2;
- if (degree == 2)
- return 3;
- }
- else if (dim == 2)
- {
- if (degree == 1)
- return 3;
- if (degree == 2)
- return 6;
- }
- else if (dim == 3)
- {
- if (degree == 1)
- return 4;
- if (degree == 2)
- return 10;
- }
-
- Assert(false, ExcNotImplemented());
-
- return 0;
- }
-
unsigned int
compute_n_polynomials_pyramid(const unsigned int dim,
const unsigned int degree)
- template <int dim>
- ScalarPolynomial<dim>::ScalarPolynomial(const unsigned int degree)
- : ScalarPolynomialsBase<dim>(degree, compute_n_polynomials(dim, degree))
- {}
-
-
-
- template <int dim>
- double
- ScalarPolynomial<dim>::compute_value(const unsigned int i,
- const Point<dim> & p) const
- {
- if (dim == 1)
- {
- if (this->degree() == 1)
- {
- if (i == 0)
- return 1.0 - p[0];
- else if (i == 1)
- return p[0];
- }
- else if (this->degree() == 2)
- {
- if (i == 0)
- return 2.0 * p[0] * p[0] - 3.0 * p[0] + 1;
- else if (i == 1)
- return 2.0 * p[0] * p[0] - p[0];
- else if (i == 2)
- return -4.0 * p[0] * p[0] + 4.0 * p[0];
- }
- }
- else if (dim == 2)
- {
- if (this->degree() == 1)
- {
- if (i == 0)
- return 1.0 - p[0] - p[1];
- else if (i == 1)
- return p[0];
- else if (i == 2)
- return p[1];
- }
- else if (this->degree() == 2)
- {
- const double t1 = 1.0 - p[0] - p[1];
- const double t2 = p[0];
- const double t3 = p[1];
-
- if (i == 0)
- return t1 * (2.0 * t1 - 1.0);
- else if (i == 1)
- return t2 * (2.0 * t2 - 1.0);
- else if (i == 2)
- return t3 * (2.0 * t3 - 1.0);
- else if (i == 3)
- return 4.0 * t2 * t1;
- else if (i == 4)
- return 4.0 * t2 * t3;
- else if (i == 5)
- return 4.0 * t3 * t1;
- }
- }
- else if (dim == 3)
- {
- if (this->degree() == 1)
- {
- if (i == 0)
- return 1.0 - p[0] - p[1] - p[2];
- else if (i == 1)
- return p[0];
- else if (i == 2)
- return p[1];
- else if (i == 3)
- return p[2];
- }
- else if (this->degree() == 2)
- {
- const double r = p[0];
- const double s = p[1];
- const double t = p[2];
- const double u = 1.0 - p[0] - p[1] - p[2];
- if (i == 0)
- return u * (2.0 * u - 1.0);
- else if (i == 1)
- return r * (2.0 * r - 1.0);
- else if (i == 2)
- return s * (2.0 * s - 1.0);
- else if (i == 3)
- return t * (2.0 * t - 1.0);
- else if (i == 4)
- return 4.0 * r * u;
- else if (i == 5)
- return 4.0 * r * s;
- else if (i == 6)
- return 4.0 * s * u;
- else if (i == 7)
- return 4.0 * t * u;
- else if (i == 8)
- return 4.0 * r * t;
- else if (i == 9)
- return 4.0 * s * t;
- }
- }
-
- Assert(false, ExcNotImplemented());
-
- return 0;
- }
-
-
-
- template <int dim>
- Tensor<1, dim>
- ScalarPolynomial<dim>::compute_grad(const unsigned int i,
- const Point<dim> & p) const
- {
- Tensor<1, dim> grad;
-
- if (dim == 1)
- {
- if (this->degree() == 1)
- {
- if (i == 0)
- grad[0] = -1.0;
- else if (i == 1)
- grad[0] = 1.0;
- }
- else if (this->degree() == 2)
- {
- if (i == 0)
- grad[0] = 4.0 * p[0] - 3.0;
- else if (i == 1)
- grad[0] = 4.0 * p[0] - 1.0;
- else if (i == 2)
- grad[0] = -8.0 * p[0] + 4.0;
- }
- else
- {
- Assert(false, ExcNotImplemented());
- }
- }
- else if (dim == 2)
- {
- if (this->degree() == 1)
- {
- if (i == 0)
- {
- grad[0] = -1.0;
- grad[1] = -1.0;
- }
- else if (i == 1)
- {
- grad[0] = +1.0;
- grad[1] = +0.0;
- }
- else if (i == 2)
- {
- grad[0] = +0.0;
- grad[1] = +1.0;
- }
- else
- {
- Assert(false, ExcNotImplemented());
- }
- }
- else if (this->degree() == 2)
- {
- if (i == 0)
- {
- grad[0] = -3.0 + 4.0 * (p[0] + p[1]);
- grad[1] = -3.0 + 4.0 * (p[0] + p[1]);
- }
- else if (i == 1)
- {
- grad[0] = 4.0 * p[0] - 1.0;
- grad[1] = 0.0;
- }
- else if (i == 2)
- {
- grad[0] = 0.0;
- grad[1] = 4.0 * p[1] - 1.0;
- }
- else if (i == 3)
- {
- grad[0] = 4.0 * (1.0 - 2.0 * p[0] - p[1]);
- grad[1] = -4.0 * p[0];
- }
- else if (i == 4)
- {
- grad[0] = 4.0 * p[1];
- grad[1] = 4.0 * p[0];
- }
- else if (i == 5)
- {
- grad[0] = -4.0 * p[1];
- grad[1] = 4.0 * (1.0 - p[0] - 2.0 * p[1]);
- }
- else
- {
- Assert(false, ExcNotImplemented());
- }
- }
- else
- {
- Assert(false, ExcNotImplemented());
- }
- }
- else if (dim == 3)
- {
- if (this->degree() == 1)
- {
- if (i == 0)
- {
- grad[0] = -1.0;
- grad[1] = -1.0;
- grad[2] = -1.0;
- }
- else if (i == 1)
- {
- grad[0] = +1.0;
- grad[1] = +0.0;
- grad[2] = +0.0;
- }
- else if (i == 2)
- {
- grad[0] = +0.0;
- grad[1] = +1.0;
- grad[2] = +0.0;
- }
- else if (i == 3)
- {
- grad[0] = +0.0;
- grad[1] = +0.0;
- grad[2] = +1.0;
- }
- }
- else if (this->degree() == 2)
- {
- const double r = p[0];
- const double s = p[1];
- const double t = p[2];
- const double u = 1.0 - p[0] - p[1] - p[2];
-
- if (i == 0)
- {
- grad[0] = -4.0 * u + 1.;
- grad[1] = grad[0];
- grad[2] = grad[0];
- }
- else if (i == 1)
- {
- grad[0] = +4.0 * r - 1.;
- grad[1] = +0.0;
- grad[2] = +0.0;
- }
- else if (i == 2)
- {
- grad[0] = +0.0;
- grad[1] = +4.0 * s - 1.;
- grad[2] = +0.0;
- }
- else if (i == 3)
- {
- grad[0] = +0.0;
- grad[1] = +0.0;
- grad[2] = +4.0 * t - 1.;
- }
- else if (i == 4)
- {
- grad[0] = +4.0 * (u - r);
- grad[1] = -4.0 * r;
- grad[2] = -4.0 * r;
- }
- else if (i == 5)
- {
- grad[0] = +4.0 * s;
- grad[1] = +4.0 * r;
- grad[2] = +0.0;
- }
- else if (i == 6)
- {
- grad[0] = -4.0 * s;
- grad[1] = +4.0 * (u - s);
- grad[2] = -4.0 * s;
- }
- else if (i == 7)
- {
- grad[0] = -4.0 * t;
- grad[1] = -4.0 * t;
- grad[2] = +4.0 * (u - t);
- }
- else if (i == 8)
- {
- grad[0] = +4.0 * t;
- grad[1] = +0.0;
- grad[2] = +4.0 * r;
- }
- else if (i == 9)
- {
- grad[0] = +0.0;
- grad[1] = +4.0 * t;
- grad[2] = +4.0 * s;
- }
- }
- else
- {
- Assert(false, ExcNotImplemented());
- }
- }
- else
- {
- Assert(false, ExcNotImplemented());
- }
-
- return grad;
- }
-
-
-
- template <int dim>
- Tensor<2, dim>
- ScalarPolynomial<dim>::compute_grad_grad(const unsigned int i,
- const Point<dim> & p) const
- {
- (void)p;
-
- Tensor<2, dim> result;
-
- if (this->degree() < 2)
- return result;
-
- if (dim == 1)
- {
- if (i == 0)
- {
- result[0][0] = 4.0;
- }
- else if (i == 1)
- {
- result[0][0] = 4.0;
- }
- else if (i == 2)
- {
- result[0][0] = -8.0;
- }
- else
- {
- Assert(false, ExcNotImplemented());
- }
- }
- else if (dim == 2)
- {
- if (i == 0)
- {
- for (unsigned j = 0; j < dim; ++j)
- for (unsigned k = 0; k < dim; ++k)
- result[j][k] = 4.0;
- }
- else if (i == 1)
- {
- result[0][0] = 4.0;
- }
- else if (i == 2)
- {
- result[1][1] = 4.0;
- }
- else if (i == 3)
- {
- result[0][0] = -8.0;
- result[0][1] = result[1][0] = -4.0;
- }
- else if (i == 4)
- {
- result[0][1] = result[1][0] = 4.0;
- }
- else if (i == 5)
- {
- result[1][1] = -8.0;
- result[0][1] = result[1][0] = -4.0;
- }
- else
- {
- Assert(false, ExcNotImplemented());
- }
- }
- else if (dim == 3)
- {
- if (i == 0)
- {
- for (unsigned j = 0; j < dim; ++j)
- for (unsigned k = 0; k < dim; ++k)
- result[j][k] = 4.0;
- }
- else if (i == 1)
- {
- result[0][0] = 4.0;
- }
- else if (i == 2)
- {
- result[1][1] = 4.0;
- }
- else if (i == 3)
- {
- result[2][2] = 4.0;
- }
- else if (i == 4)
- {
- result[0][0] = -8.0;
- result[0][1] = result[0][2] = result[1][0] = result[2][0] = -4.0;
- }
- else if (i == 5)
- {
- result[0][1] = result[1][0] = 4.0;
- }
- else if (i == 6)
- {
- result[1][1] = -8.0;
- result[0][1] = result[1][0] = result[1][2] = result[2][1] = -4.0;
- }
- else if (i == 7)
- {
- result[2][2] = -8.0;
- result[2][0] = result[2][1] = result[0][2] = result[1][2] = -4.0;
- }
- else if (i == 8)
- {
- result[0][2] = result[2][0] = 4.0;
- }
- else if (i == 9)
- {
- result[1][2] = result[2][1] = 4.0;
- }
- else
- {
- Assert(false, ExcNotImplemented());
- }
- }
- else
- {
- Assert(false, ExcNotImplemented());
- }
-
- return result;
- }
-
-
-
- template <int dim>
- void
- ScalarPolynomial<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);
-
- if (grad_grads.size() == this->n())
- for (unsigned int i = 0; i < this->n(); i++)
- grad_grads[i] = compute_grad_grad(i, unit_point);
- }
-
-
-
- template <int dim>
- Tensor<1, dim>
- ScalarPolynomial<dim>::compute_1st_derivative(const unsigned int i,
- const Point<dim> & p) const
- {
- return compute_grad(i, p);
- }
-
-
-
- template <int dim>
- Tensor<2, dim>
- ScalarPolynomial<dim>::compute_2nd_derivative(const unsigned int i,
- const Point<dim> & p) const
- {
- return compute_grad_grad(i, p);
- }
-
-
-
- template <int dim>
- Tensor<3, dim>
- ScalarPolynomial<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>
- ScalarPolynomial<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
- ScalarPolynomial<dim>::name() const
- {
- return "Simplex";
- }
-
-
-
- template <int dim>
- std::unique_ptr<ScalarPolynomialsBase<dim>>
- ScalarPolynomial<dim>::clone() const
- {
- return std::make_unique<ScalarPolynomial<dim>>(*this);
- }
-
-
-
template <int dim>
ScalarWedgePolynomial<dim>::ScalarWedgePolynomial(const unsigned int degree)
: ScalarPolynomialsBase<dim>(degree,
- template class ScalarPolynomial<1>;
- template class ScalarPolynomial<2>;
- template class ScalarPolynomial<3>;
template class ScalarWedgePolynomial<1>;
template class ScalarWedgePolynomial<2>;
template class ScalarWedgePolynomial<3>;
projects / dealii.git / commitdiff