]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Minor clean-ups. Prepare Mapping_Q for possible creation of a Mapping_C1 class.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 13 Mar 2001 12:22:00 +0000 (12:22 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 13 Mar 2001 12:22:00 +0000 (12:22 +0000)
git-svn-id: https://svn.dealii.org/trunk@4195 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/mapping.h
deal.II/deal.II/include/fe/mapping_cartesian.h
deal.II/deal.II/include/fe/mapping_q.h
deal.II/deal.II/include/fe/mapping_q1.h
deal.II/deal.II/source/fe/mapping_q.cc
deal.II/deal.II/source/fe/mapping_q1.cc

index dfc0513c518adcaa4061fb171cd2e50f1f19e848..28c7afbd5003a5cfee067709e2dd64a2bc7b7fd1 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
+//    Copyright (C) 2000, 2001 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -51,6 +51,7 @@ class Mapping : public Subscriptor
                                    * Class for internal data of finite
                                    * element and mapping objects.
                                    */
+//TODO: can we make the following class protected?    
     class InternalDataBase: public Subscriptor
     {
       public:
@@ -111,7 +112,8 @@ class Mapping : public Subscriptor
                                      * Virtual destructor.
                                      */
     virtual ~Mapping ();
-    
+
+//TODO: why make the following functions public? they are only helpful for fevalues and maybe the finite elements?
                                     /**
                                      * Prepare internal data
                                      * structures and fill in values
@@ -246,7 +248,7 @@ class Mapping : public Subscriptor
     virtual void transform_contravariant (std::vector<Tensor<1,dim> >       &dst,
                                          const std::vector<Tensor<1,dim> > &src,
                                          const InternalDataBase& internal,
-                                     const unsigned int src_offset) const = 0;
+                                         const unsigned int src_offset) const = 0;
     
                                     /**
                                      * Tranform a field of covariant vectors.
index 5b5d82bc133b67b343c35679435bdec3b573d29b..e32c05fcf4376039756df14bb359a0c702e25027 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
+//    Copyright (C) 2001 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
index 4938a56974916aa55dfcb20342a1ab7447a3dbaf..88e6ab49ce85b6cfb9a79c56fa88c159516229dc 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
+//    Copyright (C) 2000, 2001 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -57,6 +57,7 @@ class MappingQ : public MappingQ1<dim>
                                      */
     ~MappingQ ();
 
+//TODO: why make the following functions public? they are only helpful for fevalues and maybe the finite elements?
                                     /**
                                      * Implementation of the interface in
                                      * @ref{Mapping}.
@@ -178,7 +179,7 @@ class MappingQ : public MappingQ1<dim>
                                         /**
                                          * Constructor.
                                          */
-       InternalData(unsigned int n_shape_functions);
+       InternalData (const unsigned int n_shape_functions);
        
                                         /**
                                          * Unit normal vectors. Used
@@ -292,11 +293,18 @@ class MappingQ : public MappingQ1<dim>
 
                                     /**
                                      * Takes a
-                                     * @p{set_laplace_on_hex(quad)_vector}
+                                     * @p{laplace_on_hex(quad)_vector}
                                      * and applies it to the vector
                                      * @p{a} to compute the inner
-                                     * support points. They are
-                                     * appended to the vector @p{a}.
+                                     * support points as a linear
+                                     * combination of the exterior
+                                     * points.
+                                     *
+                                     * The vector @p{a} initially
+                                     * containts the locations of the
+                                     * @p{n_outer} points, the
+                                     * @p{n_inner} computed inner
+                                     * points are appended.
                                      */
     void apply_laplace_vector(const std::vector<std::vector<double> > &lvs,
                              std::vector<Point<dim> > &a) const;
@@ -319,9 +327,9 @@ class MappingQ : public MappingQ1<dim>
                                      * solution of a Laplace equation
                                      * with the position of the outer
                                      * support points as boundary
-                                     * values. The outer support
-                                     * points are all support points
-                                     * except of the inner ones.
+                                     * values, in order to make the
+                                     * transformation as smooth as
+                                     * possible.
                                      */
     void compute_support_points_laplace(
       const typename Triangulation<dim>::cell_iterator &cell,
@@ -340,31 +348,59 @@ class MappingQ : public MappingQ1<dim>
       std::vector<Point<dim> > &a) const;
     
                                     /**
-                                     * For @p{dim=2,3}. Adds (appends) the
-                                     * support points of all lines to
-                                     * the vector a.
+                                     * For @p{dim=2,3}. Append
+                                     * (appends) the support points
+                                     * of all shape functions located
+                                     * on bounding lines to the
+                                     * vector @p{a}. Points located
+                                     * on the line but on vertices
+                                     * are not included.
                                      *
                                      * Needed by the
                                      * @p{compute_support_points_simple(laplace)}
                                      * functions. For @p{dim=1} this
                                      * function is empty.
+                                     *
+                                     * This function is made virtual
+                                     * in order to allow derived
+                                     * classes to choose shape
+                                     * function support points
+                                     * differently than the present
+                                     * class, which chooses the
+                                     * points as interpolation points
+                                     * on the boundary.
                                      */
-    void add_line_support_points (const Triangulation<dim>::cell_iterator &cell,
-                                 std::vector<Point<dim> > &a) const;
+    virtual void
+    add_line_support_points (const Triangulation<dim>::cell_iterator &cell,
+                            std::vector<Point<dim> > &a) const;
 
                                     /**
-                                     * For @p{dim=3}. Adds (appends) the
-                                     * support points of all faces (quads in 3d) to
-                                     * the vector a.
+                                     * For @p{dim=3}. Append the
+                                     * support points of all shape
+                                     * functions located on bounding
+                                     * faces (quads in 3d) to the
+                                     * vector @p{a}. Points located
+                                     * on the line but on vertices
+                                     * are not included.
                                      *
                                      * Needed by the
                                      * @p{compute_support_points_laplace}
-                                     * function. For @p{dim=1} and 2 this
-                                     * function is empty.
+                                     * function. For @p{dim=1} and 2
+                                     * this function is empty.
+                                     *
+                                     * This function is made virtual
+                                     * in order to allow derived
+                                     * classes to choose shape
+                                     * function support points
+                                     * differently than the present
+                                     * class, which chooses the
+                                     * points as interpolation points
+                                     * on the boundary.
                                      */
 //TODO: rename function to add_quad_support_points, to unify notation    
-    void add_face_support_points(const typename Triangulation<dim>::cell_iterator &cell,
-                                std::vector<Point<dim> > &a) const;
+    virtual void
+    add_face_support_points(const typename Triangulation<dim>::cell_iterator &cell,
+                           std::vector<Point<dim> > &a) const;
     
                                     /**
                                      * For @p{dim=2} and 3. Simple
@@ -375,6 +411,7 @@ class MappingQ : public MappingQ1<dim>
                                      * Needed by the
                                      * @p{compute_support_points_simple}
                                      */
+//TODO: remove this function altogether?    
     void fill_quad_support_points_simple (const Triangulation<dim>::cell_iterator &cell,
                                          std::vector<Point<dim> > &a) const;
     
@@ -445,7 +482,7 @@ class MappingQ : public MappingQ1<dim>
                                      * Number of the Qp tensor
                                      * product shape functions.
                                      */
-    unsigned int n_shape_functions;
+    const unsigned int n_shape_functions;
 
                                     /**
                                      * Mapping from lexicographic to
index 304fe5c7a7c9d2cd422cba077fbbc52afe1a12c6..bd52d90ab5bd12e9ecb81d6a43f620d221d3b033 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
+//    Copyright (C) 2000, 2001 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -36,6 +36,7 @@ template <int dim>
 class MappingQ1 : public Mapping<dim>
 {
   public:
+//TODO: why make the following functions public? they are only helpful for fevalues and maybe the finite elements?
                                     /**
                                      * Implementation of the interface in
                                      * @ref{Mapping}.
@@ -213,9 +214,10 @@ class MappingQ1 : public Mapping<dim>
     {
       public:
                                         /**
-                                         * Constructor.
+                                         * Constructor. Pass the
+                                         * number of shape functions.
                                          */
-       InternalData(unsigned int n_shape_functions);
+       InternalData(const unsigned int n_shape_functions);
 
                                         /**
                                          * Shape function at quadrature
@@ -224,31 +226,31 @@ class MappingQ1 : public Mapping<dim>
                                          * vertices must be reordered
                                          * to obtain transformation.
                                          */
-       double shape (unsigned int qpoint,
-                     unsigned int shape_nr) const;
+       double shape (const unsigned int qpoint,
+                     const unsigned int shape_nr) const;
        
                                         /**
                                          * Shape function at quadrature
                                          * point. See above.
                                          */
-       double &shape (unsigned int qpoint,
-                      unsigned int shape_nr);
+       double &shape (const unsigned int qpoint,
+                      const unsigned int shape_nr);
        
                                         /**
                                          * Gradient of shape function
                                          * in quadrature point. See
                                          * above.
                                          */
-       Tensor<1,dim> derivative (unsigned int qpoint,
-                                 unsigned int shape_nr) const;
+       Tensor<1,dim> derivative (const unsigned int qpoint,
+                                 const unsigned int shape_nr) const;
 
                                         /**
                                          * Gradient of shape function
                                          * in quadrature point. See
                                          * above.
                                          */
-       Tensor<1,dim> &derivative (unsigned int qpoint,
-                                  unsigned int shape_nr);
+       Tensor<1,dim> &derivative (const unsigned int qpoint,
+                                  const unsigned int shape_nr);
        
                                         /**
                                          * Values of shape
@@ -299,11 +301,6 @@ class MappingQ1 : public Mapping<dim>
                                          */
         std::vector<std::vector<Tensor<1,dim> > > aux;
 
-                                        /**
-                                         * Number of shape functions.
-                                         */
-       unsigned int n_shape_functions;
-
                                         /**
                                          * Stores the support points of
                                          * the mapping shape functions on
@@ -316,7 +313,7 @@ class MappingQ1 : public Mapping<dim>
                                          * @p{mapping_support_points} are
                                          * stored.
                                          */
-       DoFHandler<dim>::cell_iterator cell_of_current_support_points;
+       typename DoFHandler<dim>::cell_iterator cell_of_current_support_points;
        
                                         /**
                                          * Default value of this flag
@@ -326,6 +323,21 @@ class MappingQ1 : public Mapping<dim>
                                          * @p{false}.
                                          */
        bool is_mapping_q1_data;
+
+                                        /**
+                                         * Number of shape
+                                         * functions. If this is a Q1
+                                         * mapping, then it is simply
+                                         * the number of vertices per
+                                         * cell. However, since also
+                                         * derived classes use this
+                                         * class (e.g. the
+                                         * @ref{Mapping_Q} class),
+                                         * the number of shape
+                                         * functions may also be
+                                         * different.
+                                         */
+       unsigned int n_shape_functions;
     };
     
                                     /**
@@ -424,9 +436,11 @@ class MappingQ1 : public Mapping<dim>
       std::vector<Point<dim> > &a) const;
 
                                     /**
-                                     *Number of shape functions
+                                     * Number of shape functions. Is
+                                     * simply the number of vertices
+                                     * per cell for the Q1 mapping.
                                      */
-    static const unsigned int n_shape_functions = 1 << dim;
+    static const unsigned int n_shape_functions = GeometryInfo<dim>::vertices_per_cell;
 };
 
 
index 269673debcd37d3fd359758fcf0f167a3f1b2040..136e07ddcc608b331542f9d37c806359b9b4645c 100644 (file)
@@ -85,7 +85,7 @@ MappingQ<dim>::MappingQ (const unsigned int p):
                n_outer((dim==2) ? 4+4*(degree-1)
                        :8+12*(degree-1)+6*(degree-1)*(degree-1)),
                tensor_pols(0),
-               n_shape_functions(0),
+               n_shape_functions(power(degree+1,dim)),
                renumber(0),
 //TODO: why have two ways to compute? if they both work, choose one and remove the other    
                alternative_normals_computation(false),
@@ -101,7 +101,8 @@ MappingQ<dim>::MappingQ (const unsigned int p):
     v.push_back(LagrangeEquidistant(degree,i));
 
   tensor_pols = new TensorProductPolynomials<dim> (v);
-  n_shape_functions=tensor_pols->n_tensor_product_polynomials();
+  Assert (n_shape_functions==tensor_pols->n_tensor_product_polynomials(),
+         ExcInternalError());
   Assert(n_inner+n_outer==n_shape_functions, ExcInternalError());
   
                                   // build the renumbering of the
@@ -189,6 +190,7 @@ MappingQ<dim>::compute_shapes_virtual (const std::vector<Point<dim> > &unit_poin
 }
 
 
+
 template <int dim>
 UpdateFlags
 MappingQ<dim>::update_each (const UpdateFlags in) const
@@ -260,10 +262,10 @@ MappingQ<dim>::get_subface_data (const UpdateFlags update_flags,
 
 template <int dim>
 void
-MappingQ<dim>::compute_face_data (UpdateFlags update_flags,
-                                 const Quadrature<dim>q,
-                                 const unsigned int n_original_q_points,
-                                 MappingQ1<dim>::InternalDatamapping_q1_data) const
+MappingQ<dim>::compute_face_data (UpdateFlags                   update_flags,
+                                 const Quadrature<dim>        &q,
+                                 const unsigned int            n_original_q_points,
+                                 MappingQ1<dim>::InternalData &mapping_q1_data) const
 {
                                   // convert data object to internal
                                   // data for this class. fails with
@@ -328,9 +330,17 @@ MappingQ<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
                                   // possible
   InternalData &data = dynamic_cast<InternalData&> (mapping_data);
 
-  data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells
-                                       || cell->has_boundary_lines());
-  
+                                  // check whether this cell needs
+                                  // the full mapping or can be
+                                  // treated by a reduced Q1 mapping,
+                                  // e.g. if the cell is in the
+                                  // interior of the domain
+  data.use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells
+                                         || cell->has_boundary_lines());
+
+                                  // depending on this result, use
+                                  // this or the other data object
+                                  // for the mapping
   if (data.use_mapping_q1_on_current_cell)
     MappingQ1<dim>::fill_fe_values(cell, q, data.mapping_q1_data,
                                   quadrature_points, JxW_values);
@@ -340,6 +350,7 @@ MappingQ<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
 }
 
 
+
 template <int dim>
 void
 MappingQ<dim>::fill_fe_face_values (const typename DoFHandler<dim>::cell_iterator &cell,
@@ -357,6 +368,11 @@ MappingQ<dim>::fill_fe_face_values (const typename DoFHandler<dim>::cell_iterato
                                   // possible
   InternalData &data = dynamic_cast<InternalData&> (mapping_data);
   
+                                  // check whether this cell needs
+                                  // the full mapping or can be
+                                  // treated by a reduced Q1 mapping,
+                                  // e.g. if the cell is in the
+                                  // interior of the domain
 //TODO: shouldn't we ask whether the face is at the boundary, rather than the cell?  
   data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells
                                        || cell->has_boundary_lines());
@@ -364,6 +380,9 @@ MappingQ<dim>::fill_fe_face_values (const typename DoFHandler<dim>::cell_iterato
   const unsigned int npts=q.n_quadrature_points;
   const unsigned int offset=face_no*npts;
 
+                                  // depending on this result, use
+                                  // this or the other data object
+                                  // for the mapping
   if (data.use_mapping_q1_on_current_cell)
     MappingQ1<dim>::compute_fill_face (cell, face_no, false,
                                       npts, offset, q.get_weights(),
@@ -397,6 +416,11 @@ MappingQ<dim>::fill_fe_subface_values (const typename DoFHandler<dim>::cell_iter
                                   // possible
   InternalData &data = dynamic_cast<InternalData&> (mapping_data);
 
+                                  // check whether this cell needs
+                                  // the full mapping or can be
+                                  // treated by a reduced Q1 mapping,
+                                  // e.g. if the cell is in the
+                                  // interior of the domain
 //TODO: shouldn't we ask whether the face is at the boundary, rather than the cell?  
   data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells
                                        || cell->has_boundary_lines());
@@ -405,6 +429,9 @@ MappingQ<dim>::fill_fe_subface_values (const typename DoFHandler<dim>::cell_iter
   const unsigned int offset=
     (face_no*GeometryInfo<dim>::subfaces_per_face + sub_no)*npts;
 
+                                  // depending on this result, use
+                                  // this or the other data object
+                                  // for the mapping
   if (data.use_mapping_q1_on_current_cell)
     MappingQ1<dim>::compute_fill_face (cell, face_no, true,
                                       npts, offset, q.get_weights(),
@@ -698,13 +725,18 @@ MappingQ<dim>::apply_laplace_vector(const std::vector<std::vector<double> > &lvs
   const unsigned int n_outer_apply=lvs[0].size();
   Assert(a.size()==n_outer_apply, ExcDimensionMismatch(a.size(), n_outer_apply));
 
+                                  // compute each inner point as
+                                  // linear combination of the outer
+                                  // points. the weights are given by
+                                  // the lvs entries, the outer
+                                  // points are the first (existing)
+                                  // elements of a
   for (unsigned int unit_point=0; unit_point<n_inner_apply; ++unit_point)
     {
-      const std::vector<double> &lv=lvs[unit_point];
-      Assert(lv.size()==n_outer_apply, ExcInternalError());
+      Assert(lvs[unit_point].size()==n_outer_apply, ExcInternalError());
       Point<dim> p;
       for (unsigned int k=0; k<n_outer_apply; ++k)
-       p+=lv[k]*a[k];
+       p+=lvs[unit_point][k]*a[k];
 
       a.push_back(p);
     }
@@ -717,10 +749,21 @@ MappingQ<dim>::compute_mapping_support_points(
   const typename Triangulation<dim>::cell_iterator &cell,
   std::vector<Point<dim> > &a) const
 {
+                                  // if this is a cell for which we
+                                  // want to compute the full
+                                  // mapping, then get them from the
+                                  // following function
   if (use_mapping_q_on_all_cells || cell->has_boundary_lines())
     compute_support_points_laplace(cell, a);
 //  compute_support_points_simple(cell, a);
+//TODO: can we delete the previous line?  
+
   else
+                                    // otherwise: use a Q1 mapping
+                                    // for which the mapping shape
+                                    // function support points are
+                                    // simply the vertices of the
+                                    // cell
     {
       a.resize(GeometryInfo<dim>::vertices_per_cell);
       
@@ -735,36 +778,43 @@ void
 MappingQ<dim>::compute_support_points_laplace(const typename Triangulation<dim>::cell_iterator &cell,
                                              std::vector<Point<dim> > &a) const
 {
+                                  // in any case, we need the
+                                  // vertices first
   a.resize(GeometryInfo<dim>::vertices_per_cell);
-                                  // the vertices first
   for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
     a[i] = cell->vertex(i);
   
   if (degree>1)
-    {
-      if (dim==1)
-       {
-         Assert(false, ExcNotImplemented());
-       }
-      else
-       {
-                                          // then the points on lines
-                                          // (for dim=2,3)
-         add_line_support_points (cell, a);
-         
-         if (dim==2)
-           apply_laplace_vector(laplace_on_quad_vector,a);
-         else if (dim==3)
-           {
-             add_face_support_points(cell, a);
+    switch (dim)
+      {
+       case 2:
+                                              // in 2d, add the
+                                              // points on the four
+                                              // bounding lines to
+                                              // the exterior (outer)
+                                              // points
+             add_line_support_points (cell, a);
+             apply_laplace_vector (laplace_on_quad_vector,a);
+             break;
+
+       case 3:
+                                              // in 3d also add the
+                                              // points located on
+                                              // the boundary faces
+             add_line_support_points (cell, a);
+             add_face_support_points (cell, a);
+              apply_laplace_vector (laplace_on_hex_vector, a);
+             break;
              
-              apply_laplace_vector(laplace_on_hex_vector, a);
-           }
-       }
-    }    
+       default 1:
+             Assert(false, ExcNotImplemented());
+             break;
+      };
 }
 
 
+
+//TODO: remove the following function altogether
 template <int dim>
 void
 MappingQ<dim>::compute_support_points_simple(const typename Triangulation<dim>::cell_iterator &cell,
@@ -869,38 +919,41 @@ void
 MappingQ<dim>::add_line_support_points (const Triangulation<dim>::cell_iterator &cell,
                                        std::vector<Point<dim> > &a) const
 {
-  const Boundary<dim> *boundary = 0;
-
-  std::vector<Point<dim> > line_points;
-  if (degree>2)
-    line_points.resize(degree-1);
-
-                                  // loop over each of the lines, and
-                                  // if it is at the boundary, then
-                                  // first get the boundary
-                                  // description and second compute
-                                  // the points on it
-  for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
+                                  // if we only need the midpoint,
+                                  // then ask for it.
+  if (degree==2)
     {
-      const typename Triangulation<dim>::line_iterator line = cell->line(line_no);
-      if (line->at_boundary())
-       boundary=&line->get_triangulation().get_boundary(line->boundary_indicator());
-      else
-       boundary=&straight_boundary;
-
-                                      // if we only need the
-                                      // midpoint, then ask for
-                                      // it. otherwise call the more
-                                      // complicated functions
-      if (degree==2)
+      for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
        a.push_back(boundary->get_new_point_on_line(line));
-      else
+    }
+  else
+                                    // otherwise call the more
+                                    // complicated functions and ask
+                                    // for inner points from the
+                                    // boundary description
+    {
+      std::vector<Point<dim> > line_points (degree-1);
+      
+                                      // loop over each of the lines,
+                                      // and if it is at the
+                                      // boundary, then first get the
+                                      // boundary description and
+                                      // second compute the points on
+                                      // it
+      for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
        {
+         const typename Triangulation<dim>::line_iterator line = cell->line(line_no);
+         
+         const Boundary<dim> * const boundary
+           = (line->at_boundary() ?
+              &line->get_triangulation().get_boundary(line->boundary_indicator()) :
+              &straight_boundary);
+         
          boundary->get_intermediate_points_on_line (line, line_points);
          a.insert (a.end(), line_points.begin(), line_points.end());
-       } 
+       }
     }
-}
+};
 
 
 
@@ -911,7 +964,7 @@ MappingQ<dim>::add_line_support_points (const Triangulation<dim>::cell_iterator
 template<>
 void
 MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell,
-                                    std::vector<Point<3> > &a) const
+                                    std::vector<Point<3> >                &a) const
 {
   const unsigned int faces_per_cell    = GeometryInfo<3>::faces_per_cell,
                     vertices_per_face = GeometryInfo<3>::vertices_per_face,
@@ -956,6 +1009,7 @@ MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell
                                       // on it
       if (face->at_boundary())
        {
+//TODO: move this variable out of the inner loop         
          std::vector<Point<3> > quad_points ((degree-1)*(degree-1));
 
          face->get_triangulation().get_boundary(face->boundary_indicator())
@@ -1011,6 +1065,7 @@ MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell
                                               // intermediate points
                                               // from a straight
                                               // boundary object
+//TODO: move this variable out of the loop           
              std::vector<Point<3> > quad_points ((degree-1)*(degree-1));
              
              straight_boundary.get_intermediate_points_on_quad (face, quad_points);
index 29713e2389a22f257d8789dd297a699b0a6116a9..19a896fa1c7164e65f6581d4ddfd190b09afdff3 100644 (file)
 
 
 
+template <int dim>
+const unsigned int MappingQ1<dim>::n_shape_functions;
+
+
+
 template<int dim>
-MappingQ1<dim>::InternalData::InternalData (unsigned int n_shape_functions):
-               n_shape_functions(n_shape_functions),
-               is_mapping_q1_data(true)
+MappingQ1<dim>::InternalData::InternalData (const unsigned int n_shape_functions):
+               is_mapping_q1_data(true),
+               n_shape_functions (n_shape_functions)
 {}
 
 
 template<int dim> inline
 double
-MappingQ1<dim>::InternalData::shape (unsigned int qpoint,
-                                    unsigned int shape_nr) const
+MappingQ1<dim>::InternalData::shape (const unsigned int qpoint,
+                                    const unsigned int shape_nr) const
 {
   Assert(qpoint*n_shape_functions + shape_nr < shape_values.size(),
         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, shape_values.size()));
@@ -300,10 +305,10 @@ MappingQ1<dim>::update_each (const UpdateFlags in) const
 
 template <int dim>
 void
-MappingQ1<dim>::compute_data (const UpdateFlags update_flags,
-                             const Quadrature<dim>q,
-                             const unsigned int n_original_q_points,
-                             InternalDatadata) const
+MappingQ1<dim>::compute_data (const UpdateFlags      update_flags,
+                             const Quadrature<dim> &q,
+                             const unsigned int     n_original_q_points,
+                             InternalData          &data) const
 {
   const unsigned int npts = q.n_quadrature_points;
 
@@ -313,8 +318,6 @@ MappingQ1<dim>::compute_data (const UpdateFlags update_flags,
 
   const UpdateFlags flags(data.update_flags);
   
-  //  cerr << "Data: " << hex << flags << dec << endl;
-
   if (flags & update_transformation_values)
     data.shape_values.resize(data.n_shape_functions * npts);
 
@@ -411,8 +414,8 @@ MappingQ1<dim>::compute_face_data (UpdateFlags update_flags,
 
 template <int dim>
 Mapping<dim>::InternalDataBase*
-MappingQ1<dim>::get_face_data (const UpdateFlags update_flags,
-                              const Quadrature<dim-1>quadrature) const
+MappingQ1<dim>::get_face_data (const UpdateFlags        update_flags,
+                              const Quadrature<dim-1> &quadrature) const
 {
   InternalData* data = new InternalData(n_shape_functions);
   QProjector<dim> q (quadrature, false);
@@ -1022,7 +1025,7 @@ MappingQ1<dim>::covariant_transformation (std::vector<tensor_>       &dst,
   typename std::vector<tensor_>::const_iterator vec = src.begin() + src_offset;
   typename std::vector<Tensor<2,dim> >::const_iterator tensor = data.covariant.begin();
   typename std::vector<tensor_>::iterator result = dst.begin();
-  typename std::vector<tensor_>::const_iterator end = dst.end();
+  const typename std::vector<tensor_>::const_iterator end = dst.end();
   
   while (result!=end)
     {

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.