From: Graham Harper Date: Tue, 3 Jul 2018 21:17:38 +0000 (+0000) Subject: Add the implementation of the Bernardi-Raugel element. X-Git-Tag: v9.1.0-rc1~846^2~20 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5d56330ab09723e677e36f28767b5023f4ee7e25;p=dealii.git Add the implementation of the Bernardi-Raugel element. --- diff --git a/include/deal.II/fe/fe_bernardi_raugel.h b/include/deal.II/fe/fe_bernardi_raugel.h new file mode 100644 index 0000000000..cce14f340b --- /dev/null +++ b/include/deal.II/fe/fe_bernardi_raugel.h @@ -0,0 +1,103 @@ +// --------------------------------------------------------------------- +// +// 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_fe_bernardi_raugel_h +#define dealii_fe_bernardi_raugel_h + +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +/** + * The Bernardi-Raugel element. + * + *

Degrees of freedom

+ * The BR1 element has dim degrees of freedom on each node and 1 on each + * face. The shape functions are ordered by the $(Q_1)^d$ shape functions supported on each + * vertex, increasing according to vertex ordering on the element in GeometryInfo, then the + * bubble functions follow in the ordering given in PolynomialsBernardiRaugel. + * + * This element only has 1 degree (degree $p=1$) because it yields an LBB stable pair BR1-P0 + * for Stokes problems which is lower degree than the Taylor-Hood element. The pair is + * sometimes referred to as an enriched P1-P0 element or a reduced P2-P0 element. + * + * This element does not support hanging nodes or multigrid in the current implementation. + * + */ +template +class FE_BernardiRaugel + : public FE_PolyTensor, dim> +{ +public: + /** + * Constructor for the Bernardi-Raugel element of degree @p p. The only + * supported degree is 1. + * + * @arg p: The degree of the element $p=1$ for $BR_1$. + */ + FE_BernardiRaugel(const unsigned int p = 1); + + /** + * Return a string that uniquely identifies a finite element. This class + * returns FE_BR(degree), with @p dim and @p degree replaced + * by appropriate values. + */ + virtual std::string + get_name() const; + + virtual std::unique_ptr> + clone() const; + + // documentation inherited from the base class + virtual void + convert_generalized_support_point_values_to_dof_values( + const std::vector> &support_point_values, + std::vector & nodal_values) const; + +private: + /** + * Only for internal use. Its full name is @p get_dofs_per_object_vector + * function and it creates the @p dofs_per_object vector that is needed + * within the constructor to be passed to the constructor of @p + * FiniteElementData. + */ + static std::vector + get_dpo_vector(); + + /** + * Initialize the FiniteElement::generalized_support_points and + * FiniteElement::generalized_face_support_points fields. Called from + * the constructor. See the + * @ref GlossGeneralizedSupport "glossary entry on generalized support points" + * for more information. + */ + void + initialize_support_points(); +}; + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/fe/fe_bernardi_raugel.cc b/source/fe/fe_bernardi_raugel.cc new file mode 100644 index 0000000000..bcc27fb2bd --- /dev/null +++ b/source/fe/fe_bernardi_raugel.cc @@ -0,0 +1,183 @@ +// --------------------------------------------------------------------- +// +// 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 +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + + +DEAL_II_NAMESPACE_OPEN + +template +FE_BernardiRaugel::FE_BernardiRaugel(const unsigned int p) + : FE_PolyTensor, dim>( + p, + FiniteElementData(get_dpo_vector(), + dim, + 2, + FiniteElementData::Hdiv), + std::vector(PolynomialsBernardiRaugel::compute_n_pols(p), + true), + std::vector(PolynomialsBernardiRaugel::compute_n_pols( + p), + std::vector(dim, true))) +{ + Assert(dim == 2 || dim == 3, ExcImpossibleInDim(dim)); + Assert(p == 1, ExcMessage("Only BR1 elements are available")); + + // const unsigned int n_dofs = this->dofs_per_cell; + + this->mapping_type = mapping_none; + // These must be done first, since + // they change the evaluation of + // basis functions + + // Set up the generalized support + // points + initialize_support_points(); +} + + + +template +std::string +FE_BernardiRaugel::get_name() const +{ + std::ostringstream namebuf; + namebuf << "FE_BR<" << dim << ">(" << 1 << ")"; + + return namebuf.str(); +} + + +template +std::unique_ptr> +FE_BernardiRaugel::clone() const +{ + return std_cxx14::make_unique>(*this); +} + + + +template +void +FE_BernardiRaugel::convert_generalized_support_point_values_to_dof_values( + const std::vector> &support_point_values, + std::vector & nodal_values) const +{ + Assert(support_point_values.size() == this->generalized_support_points.size(), + ExcDimensionMismatch(support_point_values.size(), + this->generalized_support_points.size())); + AssertDimension(support_point_values[0].size(), dim); + Assert(nodal_values.size() == this->dofs_per_cell, + ExcDimensionMismatch(nodal_values.size(), this->dofs_per_cell)); + + static 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); + } + + for (unsigned int i = 0; i < dim * GeometryInfo::vertices_per_cell; ++i) + nodal_values[i] = support_point_values[i][i % dim]; + + // Compute the support points values for the bubble functions + for (unsigned int i = dim * GeometryInfo::vertices_per_cell; + i < dim * GeometryInfo::vertices_per_cell + + GeometryInfo::faces_per_cell; + ++i) + { + nodal_values[i] = 0; + for (unsigned int j = 0; j < dim; ++j) + nodal_values[i] += + support_point_values[i][j] * + normals[i - dim * GeometryInfo::vertices_per_cell][j]; + } +} + + +template +std::vector +FE_BernardiRaugel::get_dpo_vector() +{ + // compute the number of unknowns per cell interior/face/edge + // + // there are dim degrees of freedom per node and there + // is 1 degree of freedom per edge in 2D (face in 3D) + std::vector dpo(dim + 1, 0u); + dpo[0] = dim; + dpo[dim - 1] = 1u; + + return dpo; +} + + +template +void +FE_BernardiRaugel::initialize_support_points() +{ + // The support points for our shape functions are the nodes and + // the face midpoints, for a total of #vertices + #faces points + this->generalized_support_points.resize(this->dofs_per_cell); + + // We need dim copies of each vertex for the first dim*vertices_per_cell + // generalized support points + for (unsigned int i = 0; i < dim * GeometryInfo::vertices_per_cell; ++i) + this->generalized_support_points[i] = + GeometryInfo::unit_cell_vertex(i / dim); + + // The remaining 2*dim points are the edge midpoints + for (unsigned int i = 0; i < dim; ++i) + { + for (unsigned int j = 0; j < 2; ++j) + { + Point p; + p[0] = 0.5; + p[1] = 0.5; + if (dim == 3) + p[2] = 0.5; + p[i] = j; + + unsigned int k = + dim * GeometryInfo::vertices_per_cell + i * 2 + j; + this->generalized_support_points[k] = p; + } + } +} + +template class FE_BernardiRaugel<1>; +template class FE_BernardiRaugel<2>; +template class FE_BernardiRaugel<3>; + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/fe/fe_poly_tensor.cc b/source/fe/fe_poly_tensor.cc index 98da782066..ae9506d861 100644 --- a/source/fe/fe_poly_tensor.cc +++ b/source/fe/fe_poly_tensor.cc @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include diff --git a/source/fe/fe_poly_tensor.inst.in b/source/fe/fe_poly_tensor.inst.in index b073e9c8e5..7501ece4f9 100644 --- a/source/fe/fe_poly_tensor.inst.in +++ b/source/fe/fe_poly_tensor.inst.in @@ -23,6 +23,8 @@ for (deal_II_dimension : DIMENSIONS) deal_II_dimension>; template class FE_PolyTensor, deal_II_dimension>; + template class FE_PolyTensor, + deal_II_dimension>; template class FE_PolyTensor, deal_II_dimension>; template class FE_PolyTensor,