1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2003 - 2019 by the deal.II authors 4 // 5 // This file is part of the deal.II library. 6 // 7 // The deal.II library is free software; you can use it, redistribute 8 // it, and/or modify it under the terms of the GNU Lesser General 9 // Public License as published by the Free Software Foundation; either 10 // version 2.1 of the License, or (at your option) any later version. 11 // The full text of the license can be found in the file LICENSE.md at 12 // the top level directory of deal.II. 13 // 14 // --------------------------------------------------------------------- 15 16 #ifndef dealii_fe_abf_h 17 #define dealii_fe_abf_h 18 19 #include <deal.II/base/config.h> 20 21 #include <deal.II/base/geometry_info.h> 22 #include <deal.II/base/polynomial.h> 23 #include <deal.II/base/polynomials_abf.h> 24 #include <deal.II/base/table.h> 25 #include <deal.II/base/tensor_product_polynomials.h> 26 27 #include <deal.II/fe/fe.h> 28 #include <deal.II/fe/fe_poly_tensor.h> 29 30 #include <vector> 31 32 DEAL_II_NAMESPACE_OPEN 33 34 35 /*!@addtogroup fe */ 36 /*@{*/ 37 38 /** 39 * Implementation of Arnold-Boffi-Falk (ABF) elements, conforming with the 40 * space H<sup>div</sup>. These elements generate vector fields with normal 41 * components continuous between mesh cells. 42 * 43 * These elements are based on an article from Arnold, Boffi and Falk: 44 * Quadrilateral H(div) finite elements, SIAM J. Numer. Anal. Vol.42, No.6, 45 * pp.2429-2451 46 * 47 * In this article, the authors demonstrate that the usual RT elements and 48 * also BDM and other proposed finite dimensional subspaces of H(div) do not 49 * work properly on arbitrary FE grids. I.e. the convergence rates deteriorate 50 * on these meshes. As a solution the authors propose the ABF elements, which 51 * are implemented in this module. 52 * 53 * This class is not implemented for the codimension one case (<tt>spacedim != 54 * dim</tt>). 55 * 56 * @todo Even if this element is implemented for two and three space 57 * dimensions, the definition of the node values relies on consistently 58 * oriented faces in 3D. Therefore, care should be taken on complicated 59 * meshes. 60 * 61 * <h3>Interpolation</h3> 62 * 63 * The 64 * @ref GlossInterpolation "interpolation" 65 * operators associated with the RT element are constructed such that 66 * interpolation and computing the divergence are commuting operations. We 67 * require this from interpolating arbitrary functions as well as the 68 * #restriction matrices. It can be achieved by two interpolation schemes, 69 * the simplified one in FE_RaviartThomasNodal and the original one here: 70 * 71 * <h4>Node values on edges/faces</h4> 72 * 73 * On edges or faces, the 74 * @ref GlossNodes "node values" 75 * are the moments of the normal component of the interpolated function with 76 * respect to the traces of the RT polynomials. Since the normal trace of the 77 * RT space of degree <i>k</i> on an edge/face is the space 78 * <i>Q<sub>k</sub></i>, the moments are taken with respect to this space. 79 * 80 * <h4>Interior node values</h4> 81 * 82 * Higher order RT spaces have interior nodes. These are moments taken with 83 * respect to the gradient of functions in <i>Q<sub>k</sub></i> on the cell 84 * (this space is the matching space for RT<sub>k</sub> in a mixed 85 * formulation). 86 * 87 * <h4>Generalized support points</h4> 88 * 89 * The node values above rely on integrals, which will be computed by 90 * quadrature rules themselves. The generalized support points are a set of 91 * points such that this quadrature can be performed with sufficient accuracy. 92 * The points needed are those of QGauss<sub>k+1</sub> on each face as well as 93 * QGauss<sub>k</sub> in the interior of the cell (or none for 94 * RT<sub>0</sub>). See the 95 * @ref GlossGeneralizedSupport "glossary entry on generalized support points" 96 * for more information. 97 */ 98 template <int dim> 99 class FE_ABF : public FE_PolyTensor<dim> 100 { 101 public: 102 /** 103 * Constructor for the ABF element of degree @p p. 104 */ 105 FE_ABF(const unsigned int p); 106 107 /** 108 * Return a string that uniquely identifies a finite element. This class 109 * returns <tt>FE_ABF<dim>(degree)</tt>, with @p dim and @p degree replaced 110 * by appropriate values. 111 */ 112 virtual std::string 113 get_name() const override; 114 115 /** 116 * This function returns @p true, if the shape function @p shape_index has 117 * non-zero function values somewhere on the face @p face_index. 118 * 119 * Right now, this is only implemented for RT0 in 1D. Otherwise, returns 120 * always @p true. 121 */ 122 virtual bool 123 has_support_on_face(const unsigned int shape_index, 124 const unsigned int face_index) const override; 125 126 // documentation inherited from the base class 127 virtual void 128 convert_generalized_support_point_values_to_dof_values( 129 const std::vector<Vector<double>> &support_point_values, 130 std::vector<double> & nodal_values) const override; 131 132 virtual std::size_t 133 memory_consumption() const override; 134 135 virtual std::unique_ptr<FiniteElement<dim, dim>> 136 clone() const override; 137 138 private: 139 /** 140 * The order of the ABF element. The lowest order elements are usually 141 * referred to as RT0, even though their shape functions are piecewise 142 * quadratics. 143 */ 144 const unsigned int rt_order; 145 146 /** 147 * Only for internal use. Its full name is @p get_dofs_per_object_vector 148 * function and it creates the @p dofs_per_object vector that is needed 149 * within the constructor to be passed to the constructor of @p 150 * FiniteElementData. 151 */ 152 static std::vector<unsigned int> 153 get_dpo_vector(const unsigned int degree); 154 155 /** 156 * Initialize the @p generalized_support_points field of the FiniteElement 157 * class and fill the tables with interpolation weights (#boundary_weights 158 * and #interior_weights). Called from the constructor. 159 * 160 * See the 161 * @ref GlossGeneralizedSupport "glossary entry on generalized support points" 162 * for more information. 163 */ 164 void 165 initialize_support_points(const unsigned int rt_degree); 166 167 /** 168 * Initialize the interpolation from functions on refined mesh cells onto 169 * the father cell. According to the philosophy of the Raviart-Thomas 170 * element, this restriction operator preserves the divergence of a function 171 * weakly. 172 */ 173 void 174 initialize_restriction(); 175 176 /** 177 * Fields of cell-independent data. 178 * 179 * For information about the general purpose of this class, see the 180 * documentation of the base class. 181 */ 182 class InternalData : public FiniteElement<dim>::InternalDataBase 183 { 184 public: 185 /** 186 * Array with shape function values in quadrature points. There is one row 187 * for each shape function, containing values for each quadrature point. 188 * Since the shape functions are vector-valued (with as many components as 189 * there are space dimensions), the value is a tensor. 190 * 191 * In this array, we store the values of the shape function in the 192 * quadrature points on the unit cell. The transformation to the real 193 * space cell is then simply done by multiplication with the Jacobian of 194 * the mapping. 195 */ 196 std::vector<std::vector<Tensor<1, dim>>> shape_values; 197 198 /** 199 * Array with shape function gradients in quadrature points. There is one 200 * row for each shape function, containing values for each quadrature 201 * point. 202 * 203 * We store the gradients in the quadrature points on the unit cell. We 204 * then only have to apply the transformation (which is a matrix-vector 205 * multiplication) when visiting an actual cell. 206 */ 207 std::vector<std::vector<Tensor<2, dim>>> shape_gradients; 208 }; 209 210 /** 211 * These are the factors multiplied to a function in the 212 * #generalized_face_support_points when computing the integration. They are 213 * organized such that there is one row for each generalized face support 214 * point and one column for each degree of freedom on the face. 215 */ 216 Table<2, double> boundary_weights; 217 /** 218 * Precomputed factors for interpolation of interior degrees of freedom. The 219 * rationale for this Table is the same as for #boundary_weights. Only, this 220 * table has a third coordinate for the space direction of the component 221 * evaluated. 222 */ 223 Table<3, double> interior_weights; 224 225 226 227 /** 228 * These are the factors multiplied to a function in the 229 * #generalized_face_support_points when computing the integration. They are 230 * organized such that there is one row for each generalized face support 231 * point and one column for each degree of freedom on the face. 232 */ 233 Table<2, double> boundary_weights_abf; 234 /** 235 * Precomputed factors for interpolation of interior degrees of freedom. The 236 * rationale for this Table is the same as for #boundary_weights. Only, this 237 * table has a third coordinate for the space direction of the component 238 * evaluated. 239 */ 240 Table<3, double> interior_weights_abf; 241 242 243 // Allow access from other dimensions. 244 template <int dim1> 245 friend class FE_ABF; 246 }; 247 248 249 250 /*@}*/ 251 252 253 DEAL_II_NAMESPACE_CLOSE 254 255 #endif 256