--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2018 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_bernardi_raugel_h
+#define dealii_polynomials_bernardi_raugel_h
+
+#include <deal.II/base/geometry_info.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/polynomial_space.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/tensor_product_polynomials.h>
+
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+/**
+ * This class implements the Bernardi-Raugel polynomials similarly to the
+ * description in the <i>Mathematics of Computation</i> paper from 1985 by
+ * Christine Bernardi and Geneviève Raugel.
+ *
+ * The Bernardi-Raugel polynomials are originally defined as an enrichment
+ * of the $(P_1)^d$ elements on simplicial meshes for Stokes problems by the
+ * addition of bubble functions, yielding a locking-free finite element which
+ * is a subset of $(P_2)^d$ elements. This implementation is an enrichment of
+ * $(Q_1)^d$ elements which is a subset of $(Q_2)^d$ elements for
+ * quadrilateral and hexahedral meshes.
+ *
+ * The $BR_1$ bubble functions are defined to have magnitude 1 at the center
+ * of face $e_i$ and direction $\mathbf{n}_i$ normal to face $e_i$, and
+ * magnitude 0 on all other vertices and faces. Ordering is consistent with
+ * the face numbering in GeometryInfo. The vector $\mathbf{n}_i$ points in
+ * the positive axis direction and not necessarily normal to the element for
+ * consistent orientation across edges.
+ *<dl>
+ * <dt> 2D bubble functions (in order)
+ * <dd> $x=0$ edge: $\mathbf{p}_1 = \mathbf{n}_1 (1-x)(y)(1-y)$
+ *
+ * $x=1$ edge: $\mathbf{p}_2 = \mathbf{n}_2 (x)(y)(1-y)$
+ *
+ * $y=0$ edge: $\mathbf{p}_3 = \mathbf{n}_3 (x)(1-x)(1-y)$
+ *
+ * $y=1$ edge: $\mathbf{p}_4 = \mathbf{n}_4 (x)(1-x)(y)$
+ *
+ * <dt> 3D bubble functions (in order)
+ * <dd> $x=0$ edge: $\mathbf{p}_1 = \mathbf{n}_1 (1-x)(y)(1-y)(z)(1-z)$
+ *
+ * $x=1$ edge: $\mathbf{p}_2 = \mathbf{n}_2 (x)(y)(1-y)(z)(1-z)$
+ *
+ * $y=0$ edge: $\mathbf{p}_3 = \mathbf{n}_3 (x)(1-x)(1-y)(z)(1-z)$
+ *
+ * $y=1$ edge: $\mathbf{p}_4 = \mathbf{n}_4 (x)(1-x)(y)(z)(1-z)$
+ *
+ * $z=0$ edge: $\mathbf{p}_5 = \mathbf{n}_5 (x)(1-x)(y)(1-y)(1-z)$
+ *
+ * $z=1$ edge: $\mathbf{p}_6 = \mathbf{n}_6 (x)(1-x)(y)(1-y)(z)$
+ *
+ *</dl>
+ *
+ * Then the $BR_1(E)$ polynomials are defined on quadrilaterals and hexahedra
+ * by $BR_1(E) = Q_1(E) \oplus \mbox{span}\{\mathbf{p}_i, i=1,...,2d\}$.
+ *
+ *
+ * @ingroup Polynomials
+ * @author Graham Harper
+ * @date 2018
+ */
+template <int dim>
+class PolynomialsBernardiRaugel
+{
+public:
+ /**
+ * Constructor. Creates all basis functions for Bernardi-Raugel polynomials
+ * of given degree.
+ *
+ * @arg k: the degree of the Bernardi-Raugel-space, which is currently
+ * limited to the case <tt>k=1</tt>
+ */
+ PolynomialsBernardiRaugel(const unsigned int k);
+
+ /**
+ * Return the number of Bernardi-Raugel polynomials.
+ */
+ unsigned int
+ n() const;
+
+
+ /**
+ * Return the degree of Bernardi-Raugel polynomials.
+ * Since the bubble functions are quadratic in at least one variable,
+ * the degree of the Bernardi-Raugel polynomials is two.
+ */
+ unsigned int
+ degree() const;
+
+ /**
+ * Return the name of the space, which is <tt>BernardiRaugel</tt>.
+ */
+ std::string
+ name() const;
+
+ /**
+ * Compute the value and derivatives of each Bernardi-
+ * Raugel polynomial at @p unit_point.
+ *
+ * The size of the vectors must either be zero or equal <tt>n()</tt>. In
+ * the first case, the function will not compute these values.
+ *
+ * If you need values or derivatives of all tensor product polynomials then
+ * use this function, rather than using any of the <tt>compute_value</tt>,
+ * <tt>compute_grad</tt> or <tt>compute_grad_grad</tt> functions, see below,
+ * in a loop over all tensor product polynomials.
+ */
+ void
+ compute(const Point<dim> & unit_point,
+ std::vector<Tensor<1, dim>> &values,
+ std::vector<Tensor<2, dim>> &grads,
+ std::vector<Tensor<3, dim>> &grad_grads,
+ std::vector<Tensor<4, dim>> &third_derivatives,
+ std::vector<Tensor<5, dim>> &fourth_derivatives) const;
+ /**
+ * Return the number of polynomials in the space <tt>BR(degree)</tt> without
+ * requiring to build an object of PolynomialsBernardiRaugel. This is
+ * required by the FiniteElement classes.
+ */
+ static unsigned int
+ compute_n_pols(unsigned int degree);
+
+
+private:
+ /**
+ * The degree of this object given to the constructor (must be 1).
+ */
+ const unsigned int my_degree;
+
+ /**
+ * The number of Bernardi-Raugel polynomials.
+ */
+ const unsigned int n_pols;
+
+ /**
+ * An object representing the polynomial space of Q
+ * functions which forms the <tt>BR</tt> polynomials through
+ * outer products of these with the corresponding unit ijk
+ * vectors.
+ */
+ const AnisotropicPolynomials<dim> polynomial_space_Q;
+ /**
+ * An object representing the polynomial space of bubble
+ * functions which forms the <tt>BR</tt> polynomials through
+ * outer products of these with the corresponding normals.
+ */
+ const AnisotropicPolynomials<dim> polynomial_space_bubble;
+
+ /**
+ * A static member function that creates the polynomial space we use to
+ * initialize the #polynomial_space_Q member variable.
+ */
+ static std::vector<std::vector<Polynomials::Polynomial<double>>>
+ create_polynomials_Q();
+
+ /**
+ * A static member function that creates the polynomial space we use to
+ * initialize the #polynomial_space_bubble member variable.
+ */
+ static std::vector<std::vector<Polynomials::Polynomial<double>>>
+ create_polynomials_bubble();
+};
+
+
+template <int dim>
+inline unsigned int
+PolynomialsBernardiRaugel<dim>::n() const
+{
+ return n_pols;
+}
+
+template <int dim>
+inline unsigned int
+PolynomialsBernardiRaugel<dim>::degree() const
+{
+ return 2;
+}
+
+template <int dim>
+inline std::string
+PolynomialsBernardiRaugel<dim>::name() const
+{
+ return "BernardiRaugel";
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2018 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 <deal.II/base/polynomials_bernardi_raugel.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+template <int dim>
+PolynomialsBernardiRaugel<dim>::PolynomialsBernardiRaugel(const unsigned int k)
+ : my_degree(k)
+ , n_pols(compute_n_pols(k))
+ , polynomial_space_Q(create_polynomials_Q())
+ , polynomial_space_bubble(create_polynomials_bubble())
+{}
+
+
+template <int dim>
+std::vector<std::vector<Polynomials::Polynomial<double>>>
+PolynomialsBernardiRaugel<dim>::create_polynomials_bubble()
+{
+ std::vector<std::vector<Polynomials::Polynomial<double>>> pols;
+ std::vector<Polynomials::Polynomial<double>> bubble_shapes;
+ bubble_shapes.push_back(Polynomials::LagrangeEquidistant(1, 0));
+ bubble_shapes.push_back(Polynomials::LagrangeEquidistant(1, 1));
+ bubble_shapes.push_back(Polynomials::LagrangeEquidistant(2, 1));
+
+ for (unsigned int d = 0; d < dim; ++d)
+ pols.push_back(bubble_shapes);
+ // In 2D, the only q_ij polynomials we will use are 31,32,13,23
+ // where ij corresponds to index (i-1)+3*(j-1) (2,5,6,7)
+
+ // In 3D, the only q_ijk polynomials we will use are 331,332,313,323,133,233
+ // where ijk corresponds to index (i-1)+3*(j-1)+9*(k-1) (8,17,20,23,24,25)
+ return pols;
+}
+
+
+template <int dim>
+std::vector<std::vector<Polynomials::Polynomial<double>>>
+PolynomialsBernardiRaugel<dim>::create_polynomials_Q()
+{
+ std::vector<std::vector<Polynomials::Polynomial<double>>> pols;
+ std::vector<Polynomials::Polynomial<double>> Q_shapes;
+ Q_shapes.push_back(Polynomials::LagrangeEquidistant(1, 0));
+ Q_shapes.push_back(Polynomials::LagrangeEquidistant(1, 1));
+ for (unsigned int d = 0; d < dim; ++d)
+ pols.push_back(Q_shapes);
+
+ return pols;
+}
+
+
+template <int dim>
+void
+PolynomialsBernardiRaugel<dim>::compute(
+ const Point<dim> & unit_point,
+ std::vector<Tensor<1, dim>> &values,
+ std::vector<Tensor<2, dim>> &grads,
+ std::vector<Tensor<3, dim>> &grad_grads,
+ std::vector<Tensor<4, dim>> &third_derivatives,
+ std::vector<Tensor<5, dim>> &fourth_derivatives) const
+{
+ Assert(values.size() == n_pols || values.size() == 0,
+ ExcDimensionMismatch(values.size(), n_pols));
+ Assert(grads.size() == n_pols || grads.size() == 0,
+ ExcDimensionMismatch(grads.size(), n_pols));
+ Assert(grad_grads.size() == n_pols || grad_grads.size() == 0,
+ ExcDimensionMismatch(grad_grads.size(), n_pols));
+ Assert(third_derivatives.size() == n_pols || third_derivatives.size() == 0,
+ ExcDimensionMismatch(third_derivatives.size(), n_pols));
+ Assert(fourth_derivatives.size() == n_pols || fourth_derivatives.size() == 0,
+ ExcDimensionMismatch(fourth_derivatives.size(), n_pols));
+
+ std::vector<double> Q_values;
+ std::vector<Tensor<1, dim>> Q_grads;
+ std::vector<Tensor<2, dim>> Q_grad_grads;
+ std::vector<Tensor<3, dim>> Q_third_derivatives;
+ std::vector<Tensor<4, dim>> Q_fourth_derivatives;
+ std::vector<double> bubble_values;
+ std::vector<Tensor<1, dim>> bubble_grads;
+ std::vector<Tensor<2, dim>> bubble_grad_grads;
+ std::vector<Tensor<3, dim>> bubble_third_derivatives;
+ std::vector<Tensor<4, dim>> bubble_fourth_derivatives;
+
+ int n_bubbles = std::pow(3, dim); // size for create_polynomials_bubble
+ int n_q = 1 << dim; // size for create_polynomials_q
+
+ // don't resize if the provided vector has 0 length
+ Q_values.resize((values.size() == 0) ? 0 : n_q);
+ Q_grads.resize((grads.size() == 0) ? 0 : n_q);
+ Q_grad_grads.resize((grad_grads.size() == 0) ? 0 : n_q);
+ Q_third_derivatives.resize((third_derivatives.size() == 0) ? 0 : n_q);
+ Q_fourth_derivatives.resize((fourth_derivatives.size() == 0) ? 0 : n_q);
+ bubble_values.resize((values.size() == 0) ? 0 : n_bubbles);
+ bubble_grads.resize((grads.size() == 0) ? 0 : n_bubbles);
+ bubble_grad_grads.resize((grad_grads.size() == 0) ? 0 : n_bubbles);
+ bubble_third_derivatives.resize((third_derivatives.size() == 0) ? 0 :
+ n_bubbles);
+ bubble_fourth_derivatives.resize(
+ (fourth_derivatives.size() == 0) ? 0 : n_bubbles);
+
+ // 1 normal vector per face, ordering consistent with GeometryInfo
+ // Normal vectors point in the +x, +y, and +z directions for
+ // consistent orientation across edges
+ std::vector<Tensor<1, dim>> normals;
+ for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+ {
+ Tensor<1, dim> normal;
+ normal[i / 2] = 1;
+ normals.push_back(normal);
+ }
+
+ // dim standard basis vectors for R^dim, usual ordering
+ std::vector<Tensor<1, dim>> units;
+ for (unsigned int i = 0; i < dim; ++i)
+ {
+ Tensor<1, dim> unit;
+ unit[i] = 1;
+ units.push_back(unit);
+ }
+
+ // set indices for the anisotropic polynomials to find
+ // them after polynomial_space_bubble.compute is called
+ std::vector<int> aniso_indices;
+ if (dim == 2)
+ {
+ aniso_indices.push_back(6);
+ aniso_indices.push_back(7);
+ aniso_indices.push_back(2);
+ aniso_indices.push_back(5);
+ }
+ else if (dim == 3)
+ {
+ aniso_indices.push_back(24);
+ aniso_indices.push_back(25);
+ aniso_indices.push_back(20);
+ aniso_indices.push_back(23);
+ aniso_indices.push_back(8);
+ aniso_indices.push_back(17);
+ }
+
+ polynomial_space_bubble.compute(unit_point,
+ bubble_values,
+ bubble_grads,
+ bubble_grad_grads,
+ bubble_third_derivatives,
+ bubble_fourth_derivatives);
+ polynomial_space_Q.compute(unit_point,
+ Q_values,
+ Q_grads,
+ Q_grad_grads,
+ Q_third_derivatives,
+ Q_fourth_derivatives);
+
+ // first dim*vertices_per_cell functions are Q_1^2 functions
+ for (unsigned int i = 0; i < dim * GeometryInfo<dim>::vertices_per_cell; ++i)
+ {
+ if (values.size() != 0)
+ {
+ values[i] = units[i % dim] * Q_values[i / dim];
+ }
+ if (grads.size() != 0)
+ {
+ grads[i] = outer_product(units[i % dim], Q_grads[i / dim]);
+ }
+ if (grad_grads.size() != 0)
+ {
+ grad_grads[i] = outer_product(units[i % dim], Q_grad_grads[i / dim]);
+ }
+ if (third_derivatives.size() != 0)
+ {
+ third_derivatives[i] =
+ outer_product(units[i % dim], Q_third_derivatives[i / dim]);
+ }
+ if (fourth_derivatives.size() != 0)
+ {
+ fourth_derivatives[i] =
+ outer_product(units[i % dim], Q_fourth_derivatives[i / dim]);
+ }
+ }
+
+ // last faces_per_cell functions are bubble functions
+ for (unsigned int i = dim * GeometryInfo<dim>::vertices_per_cell;
+ i < dim * GeometryInfo<dim>::vertices_per_cell +
+ GeometryInfo<dim>::faces_per_cell;
+ ++i)
+ {
+ unsigned int j =
+ i -
+ dim *
+ GeometryInfo<dim>::vertices_per_cell; // ranges 0 to faces_per_cell-1
+ if (values.size() != 0)
+ {
+ values[i] = normals[j] * bubble_values[aniso_indices[j]];
+ }
+ if (grads.size() != 0)
+ {
+ grads[i] = outer_product(normals[j], bubble_grads[aniso_indices[j]]);
+ }
+ if (grad_grads.size() != 0)
+ {
+ grad_grads[i] =
+ outer_product(normals[j], bubble_grad_grads[aniso_indices[j]]);
+ }
+ if (third_derivatives.size() != 0)
+ {
+ third_derivatives[i] =
+ outer_product(normals[j],
+ bubble_third_derivatives[aniso_indices[j]]);
+ }
+ if (fourth_derivatives.size() != 0)
+ {
+ fourth_derivatives[i] =
+ outer_product(normals[j],
+ bubble_fourth_derivatives[aniso_indices[j]]);
+ }
+ }
+}
+
+template <int dim>
+unsigned int
+PolynomialsBernardiRaugel<dim>::compute_n_pols(unsigned int k)
+{
+ (void)k;
+ Assert(k == 1, ExcNotImplemented());
+ if (dim == 2 || dim == 3)
+ return dim * GeometryInfo<dim>::vertices_per_cell +
+ GeometryInfo<dim>::faces_per_cell;
+ // 2*4+4=12 polynomials in 2D and 3*8+6=30 polynomials in 3D
+
+ Assert(false, ExcNotImplemented());
+ return 0;
+}
+
+template class PolynomialsBernardiRaugel<1>; // to prevent errors
+template class PolynomialsBernardiRaugel<2>;
+template class PolynomialsBernardiRaugel<3>;
+
+
+DEAL_II_NAMESPACE_CLOSE
projects / dealii.git / commitdiff