From 40c37b5b3a82945abbbc2c72c6fff9a745815998 Mon Sep 17 00:00:00 2001 From: Graham Harper Date: Tue, 3 Jul 2018 21:21:14 +0000 Subject: [PATCH] Implement the Bernardi-Raugel polynomials. --- .../base/polynomials_bernardi_raugel.h | 210 +++++++++++++++ source/base/polynomials_bernardi_raugel.cc | 254 ++++++++++++++++++ 2 files changed, 464 insertions(+) create mode 100644 include/deal.II/base/polynomials_bernardi_raugel.h create mode 100644 source/base/polynomials_bernardi_raugel.cc diff --git a/include/deal.II/base/polynomials_bernardi_raugel.h b/include/deal.II/base/polynomials_bernardi_raugel.h new file mode 100644 index 0000000000..eff77c3847 --- /dev/null +++ b/include/deal.II/base/polynomials_bernardi_raugel.h @@ -0,0 +1,210 @@ +// --------------------------------------------------------------------- +// +// 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 +#include +#include +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + + +/** + * This class implements the Bernardi-Raugel polynomials similarly to the + * description in the Mathematics of Computation 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. + *
+ *
2D bubble functions (in order) + *
$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)$ + * + *
3D bubble functions (in order) + *
$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)$ + * + *
+ * + * 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 +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 k=1 + */ + 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 BernardiRaugel. + */ + 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 n(). 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 compute_value, + * compute_grad or compute_grad_grad functions, see below, + * in a loop over all tensor product polynomials. + */ + void + compute(const Point & unit_point, + std::vector> &values, + std::vector> &grads, + std::vector> &grad_grads, + std::vector> &third_derivatives, + std::vector> &fourth_derivatives) const; + /** + * Return the number of polynomials in the space BR(degree) 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 BR polynomials through + * outer products of these with the corresponding unit ijk + * vectors. + */ + const AnisotropicPolynomials polynomial_space_Q; + /** + * An object representing the polynomial space of bubble + * functions which forms the BR polynomials through + * outer products of these with the corresponding normals. + */ + const AnisotropicPolynomials 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>> + 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>> + create_polynomials_bubble(); +}; + + +template +inline unsigned int +PolynomialsBernardiRaugel::n() const +{ + return n_pols; +} + +template +inline unsigned int +PolynomialsBernardiRaugel::degree() const +{ + return 2; +} + +template +inline std::string +PolynomialsBernardiRaugel::name() const +{ + return "BernardiRaugel"; +} + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/base/polynomials_bernardi_raugel.cc b/source/base/polynomials_bernardi_raugel.cc new file mode 100644 index 0000000000..c94ff1fb4f --- /dev/null +++ b/source/base/polynomials_bernardi_raugel.cc @@ -0,0 +1,254 @@ +// --------------------------------------------------------------------- +// +// 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_NAMESPACE_OPEN + + +template +PolynomialsBernardiRaugel::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 +std::vector>> +PolynomialsBernardiRaugel::create_polynomials_bubble() +{ + std::vector>> pols; + std::vector> 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 +std::vector>> +PolynomialsBernardiRaugel::create_polynomials_Q() +{ + std::vector>> pols; + std::vector> 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 +void +PolynomialsBernardiRaugel::compute( + const Point & unit_point, + std::vector> &values, + std::vector> &grads, + std::vector> &grad_grads, + std::vector> &third_derivatives, + std::vector> &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 Q_values; + std::vector> Q_grads; + std::vector> Q_grad_grads; + std::vector> Q_third_derivatives; + std::vector> Q_fourth_derivatives; + std::vector bubble_values; + std::vector> bubble_grads; + std::vector> bubble_grad_grads; + std::vector> bubble_third_derivatives; + std::vector> 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> normals; + for (unsigned int i = 0; i < GeometryInfo::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> 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 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::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::vertices_per_cell; + i < dim * GeometryInfo::vertices_per_cell + + GeometryInfo::faces_per_cell; + ++i) + { + unsigned int j = + i - + dim * + GeometryInfo::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 +unsigned int +PolynomialsBernardiRaugel::compute_n_pols(unsigned int k) +{ + (void)k; + Assert(k == 1, ExcNotImplemented()); + if (dim == 2 || dim == 3) + return dim * GeometryInfo::vertices_per_cell + + GeometryInfo::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 -- 2.39.5