deal.II/include/deal.II/fe/mapping_q.h

   1 // ---------------------------------------------------------------------
   2 // $Id$
   3 //
   4 // Copyright (C) 2001 - 2013 by the deal.II authors
   5 //
   6 // This file is part of the deal.II library.
   7 //
   8 // The deal.II library is free software; you can use it, redistribute
   9 // it, and/or modify it under the terms of the GNU Lesser General
  10 // Public License as published by the Free Software Foundation; either
  11 // version 2.1 of the License, or (at your option) any later version.
  12 // The full text of the license can be found in the file LICENSE at
  13 // the top level of the deal.II distribution.
  14 //
  15 // ---------------------------------------------------------------------
  16
  17 #ifndef __deal2__mapping_q_h
  18 #define __deal2__mapping_q_h
  19
  20
  21 #include <deal.II/base/config.h>
  22 #include <deal.II/base/table.h>
  23 #include <deal.II/fe/mapping_q1.h>
  24 #include <deal.II/fe/fe_q.h>
  25
  26 DEAL_II_NAMESPACE_OPEN
  27
  28 template <int dim, typename POLY> class TensorProductPolynomials;
  29
  30
  31 /*!@addtogroup mapping */
  32 /*@{*/
  33
  34 /**
  35  * Mapping class that uses Qp-mappings on boundary cells. The mapping shape
  36  * functions make use of tensor product polynomials with unit cell support
  37  * points equal to the points of the Gauss-Lobatto quadrature formula. These
  38  * points give a well-conditioned interpolation also for very high orders and
  39  * are therefore preferred over equidistant support points.
  40  *
  41  * For more details about Qp-mappings, see the `mapping' report at
  42  * <tt>deal.II/doc/reports/mapping_q/index.html</tt> in the `Reports'
  43  * section of `Documentation'.
  44  *
  45  * For more information about the <tt>spacedim</tt> template parameter
  46  * check the documentation of FiniteElement or the one of
  47  * Triangulation.
  48  *
  49  * @note Since the boundary description is closely tied to the unit cell
  50  * support points, new boundary descriptions need to explicitly use the
  51  * Gauss-Lobatto points.
  52  *
  53  * @author Ralf Hartmann, 2000, 2001, 2005; Guido Kanschat 2000, 2001
  54  */
  55 template <int dim, int spacedim=dim>
  56 class MappingQ : public MappingQ1<dim,spacedim>
  57 {
  58 public:
  59   /**
  60    * Constructor.  @p p gives the degree of mapping polynomials on boundary
  61    * cells.
  62    *
  63    * The second argument determines whether the higher order mapping should
  64    * also be used on interior cells. If its value is <code>false</code> (the
  65    * default), the a lower-order mapping is used in the interior. This is
  66    * sufficient for most cases where higher order mappings are only used to
  67    * better approximate the boundary. In that case, cells bounded by straight
  68    * lines are acceptable in the interior. However, there are cases where one
  69    * would also like to use a higher order mapping in the interior. The
  70    * MappingQEulerian class is one such case.
  71    */
  72   MappingQ (const unsigned int p,
  73             const bool use_mapping_q_on_all_cells = false);
  74
  75   /**
  76    * Copy constructor. Performs a deep copy, i.e. duplicates what #tensor_pols
  77    * points to instead of simply copying the #tensor_pols pointer as done by a
  78    * default copy constructor.
  79    */
  80   MappingQ (const MappingQ<dim,spacedim> &mapping);
  81
  82   /**
  83    * Destructor.
  84    */
  85   virtual ~MappingQ ();
  86
  87   /**
  88    * Transforms the point @p p on the unit cell to the point @p p_real on the
  89    * real cell @p cell and returns @p p_real.
  90    */
  91   virtual Point<spacedim>
  92   transform_unit_to_real_cell (
  93     const typename Triangulation<dim,spacedim>::cell_iterator &cell,
  94     const Point<dim>                                 &p) const;
  95
  96   /**
  97    * Transforms the point @p p on the real cell to the point @p p_unit on the
  98    * unit cell @p cell and returns @p p_unit.
  99    *
 100    * Uses Newton iteration and the @p transform_unit_to_real_cell function.
 101    *
 102    * In the codimension one case, this function returns the normal projection
 103    * of the real point @p p on the curve or surface identified by the @p cell.
 104    *
 105    * @note Polynomial mappings from the reference (unit) cell coordinates to
 106    * the coordinate system of a real cell are not always invertible if the
 107    * point for which the inverse mapping is to be computed lies outside the
 108    * cell's boundaries.  In such cases, the current function may fail to
 109    * compute a point on the reference cell whose image under the mapping
 110    * equals the given point @p p.  If this is the case then this function
 111    * throws an exception of type Mapping::ExcTransformationFailed .  Whether
 112    * the given point @p p lies outside the cell can therefore be determined by
 113    * checking whether the return reference coordinates lie inside of outside
 114    * the reference cell (e.g., using GeometryInfo::is_inside_unit_cell) or
 115    * whether the exception mentioned above has been thrown.
 116    */
 117   virtual Point<dim>
 118   transform_real_to_unit_cell (
 119     const typename Triangulation<dim,spacedim>::cell_iterator &cell,
 120     const Point<spacedim>                            &p) const;
 121
 122   virtual void
 123   transform (const VectorSlice<const std::vector<Tensor<1,dim> > > input,
 124              VectorSlice<std::vector<Tensor<1,spacedim> > > output,
 125              const typename Mapping<dim,spacedim>::InternalDataBase &internal,
 126              const MappingType type) const;
 127
 128   virtual void
 129   transform (const VectorSlice<const std::vector<DerivativeForm<1, dim, spacedim> > >    input,
 130              VectorSlice<std::vector<Tensor<2,spacedim> > > output,
 131              const typename Mapping<dim,spacedim>::InternalDataBase &internal,
 132              const MappingType type) const;
 133
 134   virtual
 135   void
 136   transform (const VectorSlice<const std::vector<Tensor<2, dim> > >     input,
 137              VectorSlice<std::vector<Tensor<2,spacedim> > >             output,
 138              const typename Mapping<dim,spacedim>::InternalDataBase &internal,
 139              const MappingType type) const;
 140
 141   /**
 142    * Return the degree of the mapping, i.e. the value which was passed to the
 143    * constructor.
 144    */
 145   unsigned int get_degree () const;
 146
 147   /**
 148    * Return a pointer to a copy of the present object. The caller of this copy
 149    * then assumes ownership of it.
 150    */
 151   virtual
 152   Mapping<dim,spacedim> *clone () const;
 153
 154   /**
 155    * Storage for internal data of Q_degree transformation.
 156    */
 157   class InternalData : public MappingQ1<dim,spacedim>::InternalData
 158   {
 159   public:
 160     /**
 161      * Constructor.
 162      */
 163     InternalData (const unsigned int n_shape_functions);
 164
 165
 166     /**
 167      * Return an estimate (in bytes) or the memory consumption of this object.
 168      */
 169     virtual std::size_t memory_consumption () const;
 170
 171     /**
 172      * Unit normal vectors. Used for the alternative computation of the normal
 173      * vectors. See doc of the @p alternative_normals_computation flag.
 174      *
 175      * Filled (hardcoded) once in @p get_face_data.
 176      */
 177     std::vector<std::vector<Point<dim> > > unit_normals;
 178
 179     /**
 180      * Flag that is set by the <tt>fill_fe_[[sub]face]_values</tt> function.
 181      *
 182      * If this flag is @p true we are on an interior cell and the @p
 183      * mapping_q1_data is used.
 184      */
 185     bool use_mapping_q1_on_current_cell;
 186
 187     /**
 188      * On interior cells @p MappingQ1 is used.
 189      */
 190     typename MappingQ1<dim,spacedim>::InternalData mapping_q1_data;
 191   };
 192
 193 protected:
 194   /**
 195    * Implementation of the interface in Mapping.
 196    */
 197   virtual void
 198   fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
 199                   const Quadrature<dim>                                     &quadrature,
 200                   typename Mapping<dim,spacedim>::InternalDataBase          &mapping_data,
 201                   typename std::vector<Point<spacedim> >                    &quadrature_points,
 202                   std::vector<double>                                       &JxW_values,
 203                   std::vector<DerivativeForm<1,dim,spacedim> >       &jacobians,
 204                   std::vector<DerivativeForm<2,dim,spacedim> >       &jacobian_grads,
 205                   std::vector<DerivativeForm<1,spacedim,dim> >      &inverse_jacobians,
 206                   std::vector<Point<spacedim> >                             &cell_normal_vectors,
 207                   CellSimilarity::Similarity                           &cell_similarity) const ;
 208
 209   /**
 210    * Implementation of the interface in Mapping.
 211    */
 212   virtual void
 213   fill_fe_face_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
 214                        const unsigned int face_no,
 215                        const Quadrature<dim-1>& quadrature,
 216                        typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
 217                        typename std::vector<Point<spacedim> >        &quadrature_points,
 218                        std::vector<double>             &JxW_values,
 219                        typename std::vector<Tensor<1,spacedim> >        &exterior_form,
 220                        typename std::vector<Point<spacedim> >        &normal_vectors) const ;
 221
 222   /**
 223    * Implementation of the interface in Mapping.
 224    */
 225   virtual void
 226   fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
 227                           const unsigned int face_no,
 228                           const unsigned int sub_no,
 229                           const Quadrature<dim-1>& quadrature,
 230                           typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
 231                           typename std::vector<Point<spacedim> >        &quadrature_points,
 232                           std::vector<double>             &JxW_values,
 233                           typename std::vector<Tensor<1,spacedim> >        &exterior_form,
 234                           typename std::vector<Point<spacedim> >        &normal_vectors) const ;
 235
 236   /**
 237    * For <tt>dim=2,3</tt>. Append the support points of all shape functions
 238    * located on bounding lines to the vector @p a. Points located on the line
 239    * but not on vertices are not included.
 240    *
 241    * Needed by the @p compute_support_points_laplace function . For
 242    * <tt>dim=1</tt> this function is empty.
 243    *
 244    * This function is made virtual in order to allow derived classes to choose
 245    * shape function support points differently than the present class, which
 246    * chooses the points as interpolation points on the boundary.
 247    */
 248   virtual void
 249   add_line_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
 250                            std::vector<Point<spacedim> > &a) const;
 251
 252   /**
 253    * For <tt>dim=3</tt>. Append the support points of all shape functions
 254    * located on bounding faces (quads in 3d) to the vector @p a. Points
 255    * located on the quad but not on vertices are not included.
 256    *
 257    * Needed by the @p compute_support_points_laplace function. For
 258    * <tt>dim=1</tt> and <tt>dim=2</tt> this function is empty.
 259    *
 260    * This function is made virtual in order to allow derived classes to choose
 261    * shape function support points differently than the present class, which
 262    * chooses the points as interpolation points on the boundary.
 263    */
 264   virtual void
 265   add_quad_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
 266                           std::vector<Point<spacedim> > &a) const;
 267
 268 private:
 269
 270   virtual
 271   typename Mapping<dim,spacedim>::InternalDataBase *
 272   get_data (const UpdateFlags,
 273             const Quadrature<dim> &quadrature) const;
 274
 275   virtual
 276   typename Mapping<dim,spacedim>::InternalDataBase *
 277   get_face_data (const UpdateFlags flags,
 278                  const Quadrature<dim-1>& quadrature) const;
 279
 280   virtual
 281   typename Mapping<dim,spacedim>::InternalDataBase *
 282   get_subface_data (const UpdateFlags flags,
 283                     const Quadrature<dim-1>& quadrature) const;
 284
 285   /**
 286    * Compute shape values and/or derivatives.
 287    */
 288   virtual void
 289   compute_shapes_virtual (const std::vector<Point<dim> > &unit_points,
 290                           typename MappingQ1<dim,spacedim>::InternalData &data) const;
 291
 292   /**
 293    * This function is needed by the constructor of
 294    * <tt>MappingQ<dim,spacedim></tt> for <tt>dim=</tt> 2 and 3.
 295    *
 296    * For <tt>degree<4</tt> this function sets the @p laplace_on_quad_vector to
 297    * the hardcoded data. For <tt>degree>=4</tt> and MappingQ<2> this vector is
 298    * computed.
 299    *
 300    * For the definition of the @p laplace_on_quad_vector please refer to
 301    * equation (8) of the `mapping' report.
 302    */
 303   void
 304   set_laplace_on_quad_vector(Table<2,double> &loqvs) const;
 305
 306   /**
 307    * This function is needed by the constructor of <tt>MappingQ<3></tt>.
 308    *
 309    * For <tt>degree==2</tt> this function sets the @p laplace_on_hex_vector to
 310    * the hardcoded data. For <tt>degree>2</tt> this vector is computed.
 311    *
 312    * For the definition of the @p laplace_on_hex_vector please refer to
 313    * equation (8) of the `mapping' report.
 314    */
 315   void set_laplace_on_hex_vector(Table<2,double> &lohvs) const;
 316
 317   /**
 318    * Computes the <tt>laplace_on_quad(hex)_vector</tt>.
 319    *
 320    * Called by the <tt>set_laplace_on_quad(hex)_vector</tt> functions if the
 321    * data is not yet hardcoded.
 322    *
 323    * For the definition of the <tt>laplace_on_quad(hex)_vector</tt> please
 324    * refer to equation (8) of the `mapping' report.
 325    */
 326   void compute_laplace_vector(Table<2,double> &lvs) const;
 327
 328   /**
 329    * Takes a <tt>laplace_on_hex(quad)_vector</tt> and applies it to the vector
 330    * @p a to compute the inner support points as a linear combination of the
 331    * exterior points.
 332    *
 333    * The vector @p a initially contains the locations of the @p n_outer
 334    * points, the @p n_inner computed inner points are appended.
 335    *
 336    * See equation (7) of the `mapping' report.
 337    */
 338   void apply_laplace_vector(const Table<2,double>   &lvs,
 339                             std::vector<Point<spacedim> > &a) const;
 340
 341   /**
 342    * Computes the support points of the mapping.
 343    */
 344   virtual void compute_mapping_support_points(
 345     const typename Triangulation<dim,spacedim>::cell_iterator &cell,
 346     std::vector<Point<spacedim> > &a) const;
 347
 348   /**
 349    * Computes all support points of the mapping shape functions. The inner
 350    * support points (ie. support points in quads for 2d, in hexes for 3d) are
 351    * computed using the solution of a Laplace equation with the position of
 352    * the outer support points as boundary values, in order to make the
 353    * transformation as smooth as possible.
 354    */
 355   void compute_support_points_laplace(
 356     const typename Triangulation<dim,spacedim>::cell_iterator &cell,
 357     std::vector<Point<spacedim> > &a) const;
 358
 359   /**
 360    * Needed by the @p laplace_on_quad function (for <tt>dim==2</tt>). Filled
 361    * by the constructor.
 362    *
 363    * Sizes:
 364    * laplace_on_quad_vector.size()= number of inner unit_support_points
 365    * laplace_on_quad_vector[i].size()= number of outer unit_support_points,
 366    *   i.e.  unit_support_points on the boundary of the quad
 367    *
 368    * For the definition of this vector see equation (8) of the `mapping'
 369    * report.
 370    */
 371   Table<2,double> laplace_on_quad_vector;
 372
 373   /**
 374    * Needed by the @p laplace_on_hex function (for <tt>dim==3</tt>). Filled by
 375    * the constructor.
 376    *
 377    * For the definition of this vector see equation (8) of the `mapping'
 378    * report.
 379    */
 380   Table<2,double> laplace_on_hex_vector;
 381
 382   /**
 383    * Exception.
 384    */
 385   DeclException1 (ExcLaplaceVectorNotSet,
 386                   int,
 387                   << "laplace_vector not set for degree=" << arg1 << ".");
 388
 389   /**
 390    * Degree @p p of the polynomials used as shape functions for the Qp mapping
 391    * of cells at the boundary.
 392    */
 393   const unsigned int degree;
 394
 395   /**
 396    * Number of inner mapping shape functions.
 397    */
 398   const unsigned int n_inner;
 399
 400   /**
 401    * Number of mapping shape functions on the boundary.
 402    */
 403   const unsigned int n_outer;
 404
 405   /**
 406    * Pointer to the @p dim-dimensional tensor product polynomials used as
 407    * shape functions for the Qp mapping of cells at the boundary.
 408    */
 409   const TensorProductPolynomials<dim> *tensor_pols;
 410
 411   /**
 412    * Number of the Qp tensor product shape functions.
 413    */
 414   const unsigned int n_shape_functions;
 415
 416   /**
 417    * Mapping from lexicographic to to the Qp shape function numbering. Its
 418    * size is @p dofs_per_cell.
 419    */
 420   const std::vector<unsigned int> renumber;
 421
 422   /**
 423    * If this flag is set @p true then @p MappingQ is used on all cells, not
 424    * only on boundary cells.
 425    */
 426   const bool use_mapping_q_on_all_cells;
 427
 428   /**
 429    * An FE_Q object which is only needed in 3D, since it knows how to reorder
 430    * shape functions/DoFs on non-standard faces. This is used to reorder
 431    * support points in the same way. We could make this a pointer to prevent
 432    * construction in 1D and 2D, but since memory and time requirements are not
 433    * particularly high this seems unnecessary at the moment.
 434    */
 435   const FE_Q<dim> feq;
 436
 437   /**
 438    * Declare other MappingQ classes friends.
 439    */
 440   template <int,int> friend class MappingQ;
 441 };
 442
 443 /*@}*/
 444
 445 /* -------------- declaration of explicit specializations ------------- */
 446
 447 #ifndef DOXYGEN
 448
 449 template<> MappingQ<1>::MappingQ (const unsigned int,
 450                                   const bool);
 451 template<> MappingQ<1>::~MappingQ ();
 452
 453 template<>
 454 void MappingQ<1>::compute_shapes_virtual (const std::vector<Point<1> > &unit_points,
 455                                           MappingQ1<1>::InternalData   &data) const;
 456 template <>
 457 void MappingQ<1>::set_laplace_on_quad_vector(Table<2,double> &) const;
 458
 459 template <>
 460 void MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const;
 461
 462 template <>
 463 void MappingQ<1>::compute_laplace_vector(Table<2,double> &) const;
 464 template <>
 465 void MappingQ<1>::add_line_support_points (const Triangulation<1>::cell_iterator &,
 466                                            std::vector<Point<1> > &) const;
 467 template <>
 468 void MappingQ<1,2>::add_line_support_points (const Triangulation<1,2>::cell_iterator &,
 469                                              std::vector<Point<2> > &) const;
 470 template <>
 471 void MappingQ<1,3>::add_line_support_points (const Triangulation<1,3>::cell_iterator &,
 472                                              std::vector<Point<3> > &) const;
 473
 474 template<>
 475 void MappingQ<3>::add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
 476                                           std::vector<Point<3> >                &a) const;
 477
 478 #endif // DOXYGEN
 479
 480 DEAL_II_NAMESPACE_CLOSE
 481
 482 #endif