--- /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_fe_bernardi_raugel_h
+#define dealii_fe_bernardi_raugel_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/geometry_info.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/polynomials_bernardi_raugel.h>
+#include <deal.II/base/table.h>
+#include <deal.II/base/tensor_product_polynomials.h>
+
+#include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_poly_tensor.h>
+
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * The Bernardi-Raugel element.
+ *
+ * <h3>Degrees of freedom</h3>
+ * The BR1 element has <i>dim</i> 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 <int dim>
+class FE_BernardiRaugel
+ : public FE_PolyTensor<PolynomialsBernardiRaugel<dim>, 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 <tt>FE_BR<dim>(degree)</tt>, with @p dim and @p degree replaced
+ * by appropriate values.
+ */
+ virtual std::string
+ get_name() const;
+
+ virtual std::unique_ptr<FiniteElement<dim, dim>>
+ clone() const;
+
+ // documentation inherited from the base class
+ virtual void
+ convert_generalized_support_point_values_to_dof_values(
+ const std::vector<Vector<double>> &support_point_values,
+ std::vector<double> & 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<unsigned int>
+ get_dpo_vector();
+
+ /**
+ * Initialize the FiniteElement<dim>::generalized_support_points and
+ * FiniteElement<dim>::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
--- /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_p.h>
+#include <deal.II/base/qprojector.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/std_cxx14/memory.h>
+#include <deal.II/base/table.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+
+#include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_bernardi_raugel.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <iostream>
+#include <sstream>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+template <int dim>
+FE_BernardiRaugel<dim>::FE_BernardiRaugel(const unsigned int p)
+ : FE_PolyTensor<PolynomialsBernardiRaugel<dim>, dim>(
+ p,
+ FiniteElementData<dim>(get_dpo_vector(),
+ dim,
+ 2,
+ FiniteElementData<dim>::Hdiv),
+ std::vector<bool>(PolynomialsBernardiRaugel<dim>::compute_n_pols(p),
+ true),
+ std::vector<ComponentMask>(PolynomialsBernardiRaugel<dim>::compute_n_pols(
+ p),
+ std::vector<bool>(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 <int dim>
+std::string
+FE_BernardiRaugel<dim>::get_name() const
+{
+ std::ostringstream namebuf;
+ namebuf << "FE_BR<" << dim << ">(" << 1 << ")";
+
+ return namebuf.str();
+}
+
+
+template <int dim>
+std::unique_ptr<FiniteElement<dim, dim>>
+FE_BernardiRaugel<dim>::clone() const
+{
+ return std_cxx14::make_unique<FE_BernardiRaugel<dim>>(*this);
+}
+
+
+
+template <int dim>
+void
+FE_BernardiRaugel<dim>::convert_generalized_support_point_values_to_dof_values(
+ const std::vector<Vector<double>> &support_point_values,
+ std::vector<double> & 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<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);
+ }
+
+ for (unsigned int i = 0; i < dim * GeometryInfo<dim>::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<dim>::vertices_per_cell;
+ i < dim * GeometryInfo<dim>::vertices_per_cell +
+ GeometryInfo<dim>::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<dim>::vertices_per_cell][j];
+ }
+}
+
+
+template <int dim>
+std::vector<unsigned int>
+FE_BernardiRaugel<dim>::get_dpo_vector()
+{
+ // compute the number of unknowns per cell interior/face/edge
+ //
+ // there are <tt>dim</tt> degrees of freedom per node and there
+ // is 1 degree of freedom per edge in 2D (face in 3D)
+ std::vector<unsigned int> dpo(dim + 1, 0u);
+ dpo[0] = dim;
+ dpo[dim - 1] = 1u;
+
+ return dpo;
+}
+
+
+template <int dim>
+void
+FE_BernardiRaugel<dim>::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<dim>::vertices_per_cell; ++i)
+ this->generalized_support_points[i] =
+ GeometryInfo<dim>::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<dim> p;
+ p[0] = 0.5;
+ p[1] = 0.5;
+ if (dim == 3)
+ p[2] = 0.5;
+ p[i] = j;
+
+ unsigned int k =
+ dim * GeometryInfo<dim>::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
projects / dealii.git / commitdiff