1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2015 - 2020 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_nedelec_sz_h 17 #define dealii_fe_nedelec_sz_h 18 19 #include <deal.II/base/config.h> 20 21 #include <deal.II/base/derivative_form.h> 22 #include <deal.II/base/polynomials_integrated_legendre_sz.h> 23 #include <deal.II/base/qprojector.h> 24 25 #include <deal.II/fe/fe.h> 26 #include <deal.II/fe/fe_values.h> 27 #include <deal.II/fe/mapping.h> 28 29 DEAL_II_NAMESPACE_OPEN 30 31 /*!@addtogroup fe */ 32 /*@{*/ 33 34 /** 35 * This class represents an implementation of the 36 * H<sup>curl</sup>-conforming Nédélec element described in the 37 * PhD thesis of S. Zaglmayr, <b>High Order Finite Element Methods for 38 * Electromagnetic Field Computation</b>, Johannes Kepler Universität Linz, 39 * 2006. It its used in the same context as described at the top of the 40 * description for the FE_Nedelec class. 41 * 42 * This element overcomes the sign conflict issues present in 43 * traditional Nédélec elements that arise from the edge 44 * and face parameterizations used in the basis functions. Therefore, 45 * this element should provide consistent results for general 46 * quadrilateral and hexahedral elements for which the relative 47 * orientations of edges and faces as seen from all adjacent cells are 48 * often difficult to establish. 49 * 50 * The way this element addresses the sign conflict problem is to 51 * assign local edges and faces a globally defined orientation. The 52 * local edge orientation is always chosen such that the first vertex 53 * defining the edge is the one that has the highest global vertex 54 * numbering, with the second edge vertex being that which has the 55 * lowest global vertex numbering. 56 * 57 * Similarly, the face orientation is always chosen such that the first 58 * vertex is chosen to be that with the highest global vertex numbering of the 59 * four vertices making up the face. The third vertex is then chosen to be that 60 * which is geometrically opposite the first vertex, and the second and fourth 61 * vertices are decided such that the second has a higher global vertex 62 * numbering than the fourth. 63 * 64 * Note that this element does not support non-conforming meshes at this time. 65 * 66 * Further details on this element, including some benchmarking, can be 67 * in the paper R. Kynch, P. Ledger: <b>Resolving the sign conflict 68 * problem for hp-hexahedral Nédélec elements with application to 69 * eddy current problems</b>, Computers & Structures 181, 41-54, 2017 (see 70 * https://doi.org/10.1016/j.compstruc.2016.05.021). 71 */ 72 template <int dim, int spacedim = dim> 73 class FE_NedelecSZ : public FiniteElement<dim, dim> 74 { 75 public: 76 static_assert(dim == spacedim, 77 "FE_NedelecSZ is only implemented for dim==spacedim!"); 78 79 /** 80 * Constructor for the NedelecSZ element of given @p order. The maximal 81 * polynomial degree of the shape functions is `order+1` (in each variable; 82 * the total polynomial degree may be higher). If `order = 0`, the element is 83 * linear and has degrees of freedom only on the edges. If `order >= 1` the 84 * element has degrees of freedom on the edges, faces and volume. For example 85 * the 3D version of FE_NedelecSZ has 12 degrees of freedom for `order = 0` 86 * and 54 for `degree = 1`. It is important to have enough quadrature points 87 * in order to perform the quadrature with sufficient accuracy. 88 * For example [QGauss<dim>(order + 2)](@ref QGauss) can be used for the 89 * quadrature formula, where `order` is the order of FE_NedelecSZ. 90 */ 91 FE_NedelecSZ(const unsigned int order); 92 93 virtual UpdateFlags 94 requires_update_flags(const UpdateFlags update_flags) const override; 95 96 virtual std::string 97 get_name() const override; 98 99 virtual std::unique_ptr<FiniteElement<dim, dim>> 100 clone() const override; 101 102 /** 103 * This element is vector-valued so this function will 104 * throw an exception. 105 */ 106 virtual double 107 shape_value(const unsigned int i, const Point<dim> &p) const override; 108 109 /** 110 * Not implemented. 111 */ 112 virtual double 113 shape_value_component(const unsigned int i, 114 const Point<dim> & p, 115 const unsigned int component) const override; 116 117 /** 118 * This element is vector-valued so this function will 119 * throw an exception. 120 */ 121 virtual Tensor<1, dim> 122 shape_grad(const unsigned int i, const Point<dim> &p) const override; 123 124 /** 125 * Not implemented. 126 */ 127 virtual Tensor<1, dim> 128 shape_grad_component(const unsigned int i, 129 const Point<dim> & p, 130 const unsigned int component) const override; 131 132 /** 133 * This element is vector-valued so this function will 134 * throw an exception. 135 */ 136 virtual Tensor<2, dim> 137 shape_grad_grad(const unsigned int i, const Point<dim> &p) const override; 138 139 /** 140 * Not implemented. 141 */ 142 virtual Tensor<2, dim> 143 shape_grad_grad_component(const unsigned int i, 144 const Point<dim> & p, 145 const unsigned int component) const override; 146 147 protected: 148 /** 149 * The mapping kind to be used to map shape functions from the reference 150 * cell to the mesh cell. 151 */ 152 MappingKind mapping_kind; 153 154 virtual std::unique_ptr< 155 typename dealii::FiniteElement<dim, spacedim>::InternalDataBase> 156 get_data( 157 const UpdateFlags update_flags, 158 const Mapping<dim, spacedim> &mapping, 159 const Quadrature<dim> & quadrature, 160 dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim, 161 spacedim> 162 &output_data) const override; 163 164 /** 165 * Compute information about the shape functions on the cell denoted by the 166 * first argument. Note that this function must recompute the cell-dependent 167 * degrees of freedom, and so is not thread-safe at this time. 168 */ 169 virtual void 170 fill_fe_values( 171 const typename Triangulation<dim, dim>::cell_iterator &cell, 172 const CellSimilarity::Similarity cell_similarity, 173 const Quadrature<dim> & quadrature, 174 const Mapping<dim, dim> & mapping, 175 const typename Mapping<dim, dim>::InternalDataBase & mapping_internal, 176 const dealii::internal::FEValuesImplementation::MappingRelatedData<dim, dim> 177 & mapping_data, 178 const typename FiniteElement<dim, dim>::InternalDataBase &fedata, 179 dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim> 180 &data) const override; 181 182 /** 183 * Compute information about the shape functions on the cell and face denoted 184 * by the first two arguments. Note that this function must recompute the 185 * cell-dependent degrees of freedom, and so is not thread-safe at this time. 186 */ 187 virtual void 188 fill_fe_face_values( 189 const typename Triangulation<dim, dim>::cell_iterator &cell, 190 const unsigned int face_no, 191 const Quadrature<dim - 1> & quadrature, 192 const Mapping<dim, dim> & mapping, 193 const typename Mapping<dim, dim>::InternalDataBase & mapping_internal, 194 const dealii::internal::FEValuesImplementation::MappingRelatedData<dim, dim> 195 & mapping_data, 196 const typename FiniteElement<dim, dim>::InternalDataBase &fedata, 197 dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim> 198 &data) const override; 199 200 /** 201 * Not implemented. 202 */ 203 virtual void 204 fill_fe_subface_values( 205 const typename Triangulation<dim, dim>::cell_iterator &cell, 206 const unsigned int face_no, 207 const unsigned int sub_no, 208 const Quadrature<dim - 1> & quadrature, 209 const Mapping<dim, dim> & mapping, 210 const typename Mapping<dim, dim>::InternalDataBase & mapping_internal, 211 const dealii::internal::FEValuesImplementation::MappingRelatedData<dim, dim> 212 & mapping_data, 213 const typename FiniteElement<dim, dim>::InternalDataBase &fedata, 214 dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim> 215 &data) const override; 216 217 /** 218 * Derived Internal data which is used to store cell-independent data. 219 * Note that due to the nature of this element, a number of useful 220 * pre-computed quantities are stored for the computation of cell-dependent 221 * shape functions. 222 * 223 * The main quantities which are stored are associated with edge and face 224 * parameterizations. These are: 225 * <ul> 226 * <li> $\lambda_{i}$ - trilinear function, equal to one at the $i$-th vertex 227 * and zero at all other vertices.</li> 228 * <li> $\sigma_{i}$ - linear functional associated with the $i$-th vertex.</li> 229 * </ul> 230 * 231 * The definitions of these functionals, as well as the edge and face 232 * parameterizations and edge and face extension parameters, can be found on 233 * page 82 of Zaglmayr's thesis. The details of the definition of the 234 * globally-defined edge and face orientations can be found on page 67. 235 */ 236 class InternalData : public FiniteElement<dim, dim>::InternalDataBase 237 { 238 public: 239 /** 240 * Storage for shape functions on the reference element. We only pre-compute 241 * cell-based DoFs, as the edge- and face-based DoFs depend on the cell. 242 * 243 * Due to the cell-dependent DoFs, this variable is declared mutable. 244 */ 245 mutable std::vector<std::vector<Tensor<1, dim>>> shape_values; 246 247 /** 248 * Storage for shape function gradients on the reference element. We only 249 * pre-compute cell-based DoFs, as the edge- and face-based DoFs depend on 250 * the cell. 251 * 252 * Due to the cell-dependent DoFs, this variable is declared mutable. 253 */ 254 mutable std::vector<std::vector<DerivativeForm<1, dim, dim>>> shape_grads; 255 256 /** 257 * Storage for all possible edge parameterization between vertices. These 258 * are required in the computation of edge- and face-based DoFs, which are 259 * cell-dependent. 260 * 261 * The edge parameterization of an edge, E, starting at vertex i and ending 262 * at vertex $j$ is given by $\sigma_{E} = \sigma_{i} - \sigma{j}$. 263 * 264 * sigma_imj_values[q][i][j] stores the value of the edge parameterization 265 * connected by vertices $i$ and $j$ at the q-th quadrature point. 266 * 267 * Note that not all of the $i$ and $j$ combinations result in valid edges 268 * on the hexahedral cell, but they are computed in this fashion for use 269 * with non-standard edge and face orientations. 270 */ 271 std::vector<std::vector<std::vector<double>>> sigma_imj_values; 272 273 /** 274 * Storage for gradients of all possible edge parameterizations between 275 * vertices. These are required in the computation of edge- and face-based 276 * DoFs, which are cell-dependent. Note that the components of the gradient 277 * are constant. 278 * 279 * The edge parameterization of an edge, $E$, starting at vertex $i$ and 280 * ending at vertex $j$ is given by $\sigma_{E} = \sigma_{i} - \sigma{j}$. 281 * 282 * sigma_imj_grads[i][j][d] stores the gradient of the edge parameterization 283 * connected by vertices $i$ and $j$ in component $d$. 284 * 285 * Note that the gradient of the edge parameterization is constant on an 286 * edge, so we do not need to store it at every quadrature point. 287 */ 288 std::vector<std::vector<std::vector<double>>> sigma_imj_grads; 289 290 /** 291 * Storage for values of edge parameterizations at quadrature points. These 292 * are stored for the 12 edges such that the global vertex numbering would 293 * follow the order defined by the "standard" deal.II cell. 294 * 295 * edge_sigma_values[m][q] stores the edge parameterization value at the 296 * q-th quadrature point on edge m. 297 * 298 * These values change with the orientation of the edges of a physical cell 299 * and so must take the "sign" into account when used for computation. 300 */ 301 std::vector<std::vector<double>> edge_sigma_values; 302 303 /** 304 * Storage for gradients of edge parameterization at quadrature points. 305 * These are stored for the 12 edges such that the global vertex numbering 306 * would follow the order defined by the "standard" deal.II cell. 307 * 308 * edge_sigma_grads[m][d] stores the gradient of the edge parameterization 309 * for component d on edge m. 310 * 311 * These values change with the orientation of the edges of a physical cell 312 * and so must take the "sign" into account when used for computation. 313 */ 314 std::vector<std::vector<double>> edge_sigma_grads; 315 316 /** 317 * Storage for edge extension parameters at quadrature points. These are 318 * stored for the 12 edges such that the global vertex numbering would 319 * follow the order defined by the "standard" deal.II cell. 320 * 321 * The edge extension parameter of an edge, $E$, starting at vertex $i$ and 322 * ending at vertex $j$ is given by $\lambda_{E} = \lambda_{i} + 323 * \lambda_{j}$. 324 * 325 * Note that under this definition, the values of $\lambda_{E}$ do not 326 * change with the orientation of the edge. 327 * 328 * edge_lambda_values[m][q] stores the edge extension parameter value at 329 * the $q$-th quadrature point on edge $m$. 330 */ 331 std::vector<std::vector<double>> edge_lambda_values; 332 333 /** 334 * Storage for gradients of edge extension parameters in 2D. In this case 335 * they are constant. These are stored for the 12 edges such that the global 336 * vertex numbering* would follow the order defined by the "standard" 337 * deal.II cell. 338 * 339 * edge_lambda_grads_2d[m][d] stores the gradient of the edge extension 340 * parameter for component $d$ on edge $m$. 341 */ 342 std::vector<std::vector<double>> edge_lambda_grads_2d; 343 344 /** 345 * Storage for gradients of edge extension parameters in 3D. In this case 346 * they are non-constant. These are stored for the 12 edges such that the 347 * global vertex numbering* would follow the order defined by the 348 * "standard" deal.II cell. 349 * 350 * edge_lambda_grads_3d[m][q][d] stores the gradient of the edge extension 351 * parameter for component $d$ at the $q$-th quadrature point on edge m. 352 */ 353 std::vector<std::vector<std::vector<double>>> edge_lambda_grads_3d; 354 355 /** 356 * Storage for 2nd derivatives of edge extension parameters in 3D, which are 357 * constant across the cell. These are stored for the 12 edges such that the 358 * global vertex numbering* would follow the order defined by the 359 * "standard" deal.II cell. 360 * 361 * edge_lambda_gradgrads_3d[m][d1][d2] stores the 2nd derivatives of the 362 * edge extension parameters with respect to components d1 and d2 on edge 363 * $m$. 364 */ 365 std::vector<std::vector<std::vector<double>>> edge_lambda_gradgrads_3d; 366 367 /** 368 * Storage for the face extension parameters. These are stored for the 6 369 * faces such that the global vertex numbering would follow the order 370 * defined by the "standard" deal.II cell. 371 * 372 * The face extension parameter of a face, F, defined by the vertices 373 * v1, v2, v3, v4 is given by 374 * $\lambda_{F} = \lambda_{v1} + \lambda_{v2} + \lambda_{v3} + 375 * \lambda_{v4}$. 376 * 377 * Note that under this definition, the values of $\lambda_{F}$ do not 378 * change with the orientation of the face. 379 * 380 * face_lambda_values[m][q] stores the face extension parameter value at 381 * the $q$-th quadrature point on face $m$. 382 */ 383 std::vector<std::vector<double>> face_lambda_values; 384 385 /** 386 * Storage for gradients of face extension parameters. These are stored for 387 * the 6 faces such that the global vertex numbering would follow the order 388 * defined by the "standard" deal.II cell. 389 * 390 * face_lambda_grads[m][d] stores the gradient of the face extension 391 * parameters for component $d$ on face $m$. 392 */ 393 std::vector<std::vector<double>> face_lambda_grads; 394 }; 395 396 private: 397 /** 398 * Internal function to return a vector of "dofs per object" 399 * where the components of the returned vector refer to: 400 * 0 = vertex 401 * 1 = edge 402 * 2 = face (which is a cell in 2D) 403 * 3 = cell 404 */ 405 static std::vector<unsigned int> 406 get_dpo_vector(const unsigned int degree); 407 408 /** 409 * Internal storage for all required integrated Legendre polynomials. 410 */ 411 std::vector<Polynomials::Polynomial<double>> IntegratedLegendrePolynomials; 412 413 /** 414 * Internal function to populate the internal array of integrated Legendre 415 * polynomials. 416 */ 417 void 418 create_polynomials(const unsigned int degree); 419 420 /** 421 * Returns the number of DoFs in the basis set. 422 */ 423 unsigned int 424 compute_num_dofs(const unsigned int degree) const; 425 426 /** 427 * Populates cell-dependent edge-based shape functions on the given 428 * InternalData object. 429 */ 430 void 431 fill_edge_values(const typename Triangulation<dim, dim>::cell_iterator &cell, 432 const Quadrature<dim> &quadrature, 433 const InternalData & fedata) const; 434 435 /** 436 * Populates the cell-dependent face-based shape functions on the given 437 * InternalData object. 438 */ 439 void 440 fill_face_values(const typename Triangulation<dim, dim>::cell_iterator &cell, 441 const Quadrature<dim> &quadrature, 442 const InternalData & fedata) const; 443 }; 444 445 446 447 /*@}*/ 448 449 DEAL_II_NAMESPACE_CLOSE 450 451 #endif 452