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_raviart_thomas_h 17 #define dealii_fe_raviart_thomas_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_raviart_thomas.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 /*!@addtogroup fe */ 35 /*@{*/ 36 37 /** 38 * Implementation of Raviart-Thomas (RT) elements. The Raviart-Thomas space 39 * is designed to solve problems in which the solution only lives in the 40 * space 41 * $H^\text{div}=\{ {\mathbf u} \in L_2: \text{div}\, {\mathbf u} \in L_2\}$, 42 * rather than in the more commonly used space 43 * $H^1=\{ u \in L_2: \nabla u \in L_2\}$. In other words, the solution must 44 * be a vector field whose divergence is square integrable, but for which the 45 * gradient may not be square integrable. The typical application for this 46 * space (and these elements) is to the mixed formulation of the Laplace 47 * equation and related situations, see for example step-20. The defining 48 * characteristic of functions in $H^\text{div}$ is that they are in 49 * general discontinuous -- but that if you draw a line in 2d (or a 50 * surface in 3d), then the <i>normal</i> component of the vector 51 * field must be continuous across the line (or surface) even though 52 * the tangential component may not be. As a consequence, the 53 * Raviart-Thomas element is constructed in such a way that (i) it is 54 * @ref vector_valued "vector-valued", (ii) the shape functions are 55 * discontinuous, but (iii) the normal component of the vector field 56 * represented by each shape function is continuous across the faces 57 * of cells. 58 * 59 * Other properties of the Raviart-Thomas element are that (i) it is 60 * @ref GlossPrimitive "not a primitive element"; (ii) the shape functions 61 * are defined so that certain integrals over the faces are either zero 62 * or one, rather than the common case of certain point values being 63 * either zero or one. (There is, however, the FE_RaviartThomasNodal 64 * element that uses point values.) 65 * 66 * We follow the commonly used -- though confusing -- definition of the "degree" 67 * of RT elements. Specifically, the "degree" of the element denotes 68 * the polynomial degree of the <i>largest complete polynomial subspace</i> 69 * contained in the finite element space, even if the space may contain shape 70 * functions of higher polynomial degree. The lowest order element is 71 * consequently FE_RaviartThomas(0), i.e., the Raviart-Thomas element "of 72 * degree zero", even though the functions of this space are in general 73 * polynomials of degree one in each variable. This choice of "degree" 74 * implies that the approximation order of the function itself is 75 * <i>degree+1</i>, as with usual polynomial spaces. The numbering so chosen 76 * implies the sequence 77 * @f[ 78 * Q_{k+1} 79 * \stackrel{\text{grad}}{\rightarrow} 80 * \text{Nedelec}_k 81 * \stackrel{\text{curl}}{\rightarrow} 82 * \text{RaviartThomas}_k 83 * \stackrel{\text{div}}{\rightarrow} 84 * DGQ_{k} 85 * @f] 86 * 87 * This class is not implemented for the codimension one case (<tt>spacedim != 88 * dim</tt>). 89 * 90 * 91 * <h3>Interpolation</h3> 92 * 93 * The 94 * @ref GlossInterpolation "interpolation" 95 * operators associated with the RT element are constructed such that 96 * interpolation and computing the divergence are commuting operations. We 97 * require this from interpolating arbitrary functions as well as the 98 * #restriction matrices. It can be achieved by two interpolation schemes, 99 * the simplified one in FE_RaviartThomasNodal and the original one here: 100 * 101 * <h4>Node values on edges/faces</h4> 102 * 103 * On edges or faces, the 104 * @ref GlossNodes "node values" 105 * are the moments of the normal component of the interpolated function with 106 * respect to the traces of the RT polynomials. Since the normal trace of the 107 * RT space of degree <i>k</i> on an edge/face is the space 108 * <i>Q<sub>k</sub></i>, the moments are taken with respect to this space. 109 * 110 * <h4>Interior node values</h4> 111 * 112 * Higher order RT spaces have interior nodes. These are moments taken with 113 * respect to the gradient of functions in <i>Q<sub>k</sub></i> on the cell 114 * (this space is the matching space for RT<sub>k</sub> in a mixed 115 * formulation). 116 * 117 * <h4>Generalized support points</h4> 118 * 119 * The node values above rely on integrals, which will be computed by 120 * quadrature rules themselves. The generalized support points are a set of 121 * points such that this quadrature can be performed with sufficient accuracy. 122 * The points needed are those of QGauss<sub>k+1</sub> on each face as well as 123 * QGauss<sub>k+1</sub> in the interior of the cell (or none for 124 * RT<sub>0</sub>). 125 */ 126 template <int dim> 127 class FE_RaviartThomas : public FE_PolyTensor<dim> 128 { 129 public: 130 /** 131 * Constructor for the Raviart-Thomas element of degree @p p. 132 */ 133 FE_RaviartThomas(const unsigned int p); 134 135 /** 136 * Return a string that uniquely identifies a finite element. This class 137 * returns <tt>FE_RaviartThomas<dim>(degree)</tt>, with @p dim and @p degree 138 * replaced by appropriate values. 139 */ 140 virtual std::string 141 get_name() const override; 142 143 // documentation inherited from the base class 144 virtual std::unique_ptr<FiniteElement<dim, dim>> 145 clone() const override; 146 147 /** 148 * This function returns @p true, if the shape function @p shape_index has 149 * non-zero function values somewhere on the face @p face_index. 150 * 151 * Right now, this is only implemented for RT0 in 1D. Otherwise, returns 152 * always @p true. 153 */ 154 virtual bool 155 has_support_on_face(const unsigned int shape_index, 156 const unsigned int face_index) const override; 157 158 // documentation inherited from the base class 159 virtual void 160 convert_generalized_support_point_values_to_dof_values( 161 const std::vector<Vector<double>> &support_point_values, 162 std::vector<double> & nodal_values) const override; 163 164 /** 165 * Return a list of constant modes of the element. This method is currently 166 * not correctly implemented because it returns ones for all components. 167 */ 168 virtual std::pair<Table<2, bool>, std::vector<unsigned int>> 169 get_constant_modes() const override; 170 171 virtual std::size_t 172 memory_consumption() const override; 173 174 private: 175 /** 176 * Only for internal use. Its full name is @p get_dofs_per_object_vector 177 * function and it creates the @p dofs_per_object vector that is needed 178 * within the constructor to be passed to the constructor of @p 179 * FiniteElementData. 180 */ 181 static std::vector<unsigned int> 182 get_dpo_vector(const unsigned int degree); 183 184 /** 185 * Initialize the @p generalized_support_points field of the FiniteElement 186 * class and fill the tables with interpolation weights (#boundary_weights 187 * and #interior_weights). Called from the constructor. 188 */ 189 void 190 initialize_support_points(const unsigned int rt_degree); 191 192 /** 193 * Initialize the interpolation from functions on refined mesh cells onto 194 * the father cell. According to the philosophy of the Raviart-Thomas 195 * element, this restriction operator preserves the divergence of a function 196 * weakly. 197 */ 198 void 199 initialize_restriction(); 200 201 /** 202 * These are the factors multiplied to a function in the 203 * #generalized_face_support_points when computing the integration. They are 204 * organized such that there is one row for each generalized face support 205 * point and one column for each degree of freedom on the face. 206 * 207 * See the 208 * @ref GlossGeneralizedSupport "glossary entry on generalized support points" 209 * for more information. 210 */ 211 Table<2, double> boundary_weights; 212 213 /** 214 * Precomputed factors for interpolation of interior degrees of freedom. The 215 * rationale for this Table is the same as for #boundary_weights. Only, this 216 * table has a third coordinate for the space direction of the component 217 * evaluated. 218 */ 219 Table<3, double> interior_weights; 220 221 // Allow access from other dimensions. 222 template <int dim1> 223 friend class FE_RaviartThomas; 224 }; 225 226 227 228 /** 229 * The Raviart-Thomas elements with node functionals defined as point values 230 * in Gauss points. 231 * 232 * <h3>Description of node values</h3> 233 * 234 * For this Raviart-Thomas element, the node values are not cell and face 235 * moments with respect to certain polynomials, but the values in quadrature 236 * points. Following the general scheme for numbering degrees of freedom, the 237 * node values on edges are first, edge by edge, according to the natural 238 * ordering of the edges of a cell. The interior degrees of freedom are last. 239 * 240 * For an RT-element of degree <i>k</i>, we choose <i>(k+1)<sup>d-1</sup></i> 241 * Gauss points on each face. These points are ordered lexicographically with 242 * respect to the orientation of the face. This way, the normal component 243 * which is in <i>Q<sub>k</sub></i> is uniquely determined. Furthermore, since 244 * this Gauss-formula is exact on <i>Q<sub>2k+1</sub></i>, these node values 245 * correspond to the exact integration of the moments of the RT-space. 246 * 247 * In the interior of the cells, the moments are with respect to an 248 * anisotropic <i>Q<sub>k</sub></i> space, where the test functions are one 249 * degree lower in the direction corresponding to the vector component under 250 * consideration. This is emulated by using an anisotropic Gauss formula for 251 * integration. 252 * 253 * @todo The current implementation is for Cartesian meshes only. You must use 254 * MappingCartesian. 255 * 256 * @todo Even if this element is implemented for two and three space 257 * dimensions, the definition of the node values relies on consistently 258 * oriented faces in 3D. Therefore, care should be taken on complicated 259 * meshes. 260 * 261 * @note The degree stored in the member variable 262 * FiniteElementData<dim>::degree is higher by one than the constructor 263 * argument! 264 */ 265 template <int dim> 266 class FE_RaviartThomasNodal : public FE_PolyTensor<dim> 267 { 268 public: 269 /** 270 * Constructor for the Raviart-Thomas element of degree @p p. 271 */ 272 FE_RaviartThomasNodal(const unsigned int p); 273 274 /** 275 * Return a string that uniquely identifies a finite element. This class 276 * returns <tt>FE_RaviartThomasNodal<dim>(degree)</tt>, with @p dim and @p 277 * degree replaced by appropriate values. 278 */ 279 virtual std::string 280 get_name() const override; 281 282 // documentation inherited from the base class 283 virtual std::unique_ptr<FiniteElement<dim, dim>> 284 clone() const override; 285 286 virtual void 287 convert_generalized_support_point_values_to_dof_values( 288 const std::vector<Vector<double>> &support_point_values, 289 std::vector<double> & nodal_values) const override; 290 291 virtual void 292 get_face_interpolation_matrix(const FiniteElement<dim> &source, 293 FullMatrix<double> & matrix, 294 const unsigned int face_no = 0) const override; 295 296 virtual void 297 get_subface_interpolation_matrix( 298 const FiniteElement<dim> &source, 299 const unsigned int subface, 300 FullMatrix<double> & matrix, 301 const unsigned int face_no = 0) const override; 302 virtual bool 303 hp_constraints_are_implemented() const override; 304 305 virtual std::vector<std::pair<unsigned int, unsigned int>> 306 hp_vertex_dof_identities(const FiniteElement<dim> &fe_other) const override; 307 308 virtual std::vector<std::pair<unsigned int, unsigned int>> 309 hp_line_dof_identities(const FiniteElement<dim> &fe_other) const override; 310 311 virtual std::vector<std::pair<unsigned int, unsigned int>> 312 hp_quad_dof_identities(const FiniteElement<dim> &fe_other, 313 const unsigned int face_no = 0) const override; 314 315 /** 316 * @copydoc FiniteElement::compare_for_domination() 317 */ 318 virtual FiniteElementDomination::Domination 319 compare_for_domination(const FiniteElement<dim> &fe_other, 320 const unsigned int codim = 0) const override final; 321 322 private: 323 /** 324 * Only for internal use. Its full name is @p get_dofs_per_object_vector 325 * function and it creates the @p dofs_per_object vector that is needed 326 * within the constructor to be passed to the constructor of @p 327 * FiniteElementData. 328 */ 329 static std::vector<unsigned int> 330 get_dpo_vector(const unsigned int degree); 331 332 /** 333 * Compute the vector used for the @p restriction_is_additive field passed 334 * to the base class's constructor. 335 */ 336 static std::vector<bool> 337 get_ria_vector(const unsigned int degree); 338 339 /** 340 * This function returns @p true, if the shape function @p shape_index has 341 * non-zero function values somewhere on the face @p face_index. 342 * 343 * Right now, this is only implemented for RT0 in 1D. Otherwise, returns 344 * always @p true. 345 */ 346 virtual bool 347 has_support_on_face(const unsigned int shape_index, 348 const unsigned int face_index) const override; 349 /** 350 * Initialize the FiniteElement<dim>::generalized_support_points and 351 * FiniteElement<dim>::generalized_face_support_points fields. Called from 352 * the constructor. 353 * 354 * See the 355 * @ref GlossGeneralizedSupport "glossary entry on generalized support points" 356 * for more information. 357 */ 358 void 359 initialize_support_points(const unsigned int rt_degree); 360 }; 361 362 363 /*@}*/ 364 365 /* -------------- declaration of explicit specializations ------------- */ 366 367 #ifndef DOXYGEN 368 369 template <> 370 void 371 FE_RaviartThomas<1>::initialize_restriction(); 372 373 #endif // DOXYGEN 374 375 DEAL_II_NAMESPACE_CLOSE 376 377 #endif 378