]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Clean-ups, annotations, etc.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 8 Mar 2001 10:35:32 +0000 (10:35 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 8 Mar 2001 10:35:32 +0000 (10:35 +0000)
git-svn-id: https://svn.dealii.org/trunk@4170 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/include/grid/tria_boundary.h
deal.II/deal.II/source/fe/fe_dgq.cc
deal.II/deal.II/source/fe/fe_q.cc
deal.II/deal.II/source/fe/mapping.cc
deal.II/deal.II/source/fe/mapping_cartesian.cc
deal.II/deal.II/source/fe/mapping_q.cc
deal.II/deal.II/source/fe/mapping_q1.cc

index fed9d47ca019f9419857d58aaa030951685e49b3..dfc0513c518adcaa4061fb171cd2e50f1f19e848 100644 (file)
@@ -111,6 +111,7 @@ class Mapping : public Subscriptor
                                      * Virtual destructor.
                                      */
     virtual ~Mapping ();
+    
                                     /**
                                      * Prepare internal data
                                      * structures and fill in values
@@ -129,6 +130,7 @@ class Mapping : public Subscriptor
     virtual InternalDataBase*
     get_face_data (const UpdateFlags flags,
                   const Quadrature<dim-1>& quadrature) const = 0;
+    
                                     /**
                                      * Prepare internal data
                                      * structure for transformation
index be6c2e2374161c2b3cc7bcd8da0dead0149097ed..5b5d82bc133b67b343c35679435bdec3b573d29b 100644 (file)
@@ -142,6 +142,7 @@ class MappingCartesian : public Mapping<dim>
                                      * @p{p_real} on the real cell
                                      * @p{cell} and returns @p{p_real}.
                                      */
+//TODO: document meaning of mdata argument    
     virtual Point<dim> transform_unit_to_real_cell (
       const typename Triangulation<dim>::cell_iterator cell,
       const Point<dim> &p,
@@ -209,14 +210,14 @@ class MappingCartesian : public Mapping<dim>
                                         /**
                                          * Constructor.
                                          */
-       InternalData (const Quadrature<dim>quadrature);
+       InternalData (const Quadrature<dim> &quadrature);
 
                                         /**
                                          * Length of the cell in
                                          * different coordinate
                                          * directions.
                                          */
-       std::vector<double> length;
+       Tensor<1,dim> length;
 
                                         /**
                                          * Vector of all quadrature
@@ -236,7 +237,7 @@ class MappingCartesian : public Mapping<dim>
         std::vector<std::vector<Tensor<1,dim> > > unit_tangentials;
        
                                         /**
-                                         * Auxuliary vectors for internal use.
+                                         * Auxiliary vectors for internal use.
                                          */
         std::vector<std::vector<Tensor<1,dim> > > aux;
     };
@@ -252,6 +253,14 @@ class MappingCartesian : public Mapping<dim>
                       InternalData& data,
                       std::vector<Point<dim> > &quadrature_points,
                       std::vector<Point<dim> >& normal_vectors) const;
+
+  private:
+                                    /**
+                                     * Value to indicate that a given
+                                     * face or subface number is
+                                     * invalid.
+                                     */
+    static const unsigned int invalid_face_number = static_cast<unsigned int>(-1);    
 };
 
 
index 5c33194d566d7e5f4754ce4c22fcf36b31bec1df..4938a56974916aa55dfcb20342a1ab7447a3dbaf 100644 (file)
@@ -25,7 +25,11 @@ class LagrangeEquidistant;
 
 
 //TODO: fill_fe_face_values should exist in a version doing all faces
-// to save initialization time.
+//TODO: to save initialization time.
+
+//TODO: doc: laplace_on_quad_vector, compute_laplace_on_quad, etc all
+//TODO: reference to each other, there is no description what they actually
+//TODO: do or are good for
 
 /**
  * Mapping class that uses Qp-mappings on boundary AND on inner
@@ -269,7 +273,7 @@ class MappingQ : public MappingQ1<dim>
                                      *
                                      * For @p{degree==2} this function
                                      * sets the
-                                     * @p{laplace_on_quad_vector} to
+                                     * @p{laplace_on_hex_vector} to
                                      * the hardcoded data. For
                                      * @p{degree>2} this vector is
                                      * computed.
@@ -358,6 +362,7 @@ class MappingQ : public MappingQ1<dim>
                                      * function. For @p{dim=1} and 2 this
                                      * function is empty.
                                      */
+//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;
     
index 46450cfae4c0ff118ac77b63cc493bbe251cc4c4..304fe5c7a7c9d2cd422cba077fbbc52afe1a12c6 100644 (file)
@@ -147,6 +147,7 @@ class MappingQ1 : public Mapping<dim>
                                      * @p{p_real} on the real cell
                                      * @p{cell} and returns @p{p_real}.
                                      */
+//TODO: document meaning of mdata argument    
     virtual Point<dim> transform_unit_to_real_cell (
       const typename Triangulation<dim>::cell_iterator cell,
       const Point<dim> &p,
index 4503cd8030320b0a260dd20a4afcfaa5a277f503..80ab2385645a2d628e4d8daf3243e4f51e458de1 100644 (file)
@@ -121,9 +121,11 @@ class Boundary : public Subscriptor
     get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const;
 
                                     /**
-                                     * Return intermediate points
-                                     * on the boundary line.
-                                     *
+                                     * Return equally spaced
+                                     * intermediate points on a
+                                     * boundary line.
+//TODO: this function is also called for inner lines!?
+                                      *
                                      * The number of points requested
                                      * is given by the size of the
                                      * vector @p{points}. It is the
@@ -145,8 +147,10 @@ class Boundary : public Subscriptor
                                     typename std::vector<Point<dim> > &points) const;
     
                                     /**
-                                     * Return intermediate points
-                                     * on the boundary quad.
+                                     * Return equally spaced
+                                     * intermediate points on a
+                                     * boundary quad.
+//TODO: this function is also called for inner quads!?
                                      *
                                      * The number of points requested
                                      * is given by the size of the
@@ -235,8 +239,9 @@ class StraightBoundary : public Boundary<dim> {
                                      * p{n+1} partitions of equal
                                      * lengths.
                                      *
-                                     * Refer to the general documentation of
-                                     * this class and the documentation of the
+                                     * Refer to the general
+                                     * documentation of this class
+                                     * and the documentation of the
                                      * base class.
                                      */
     virtual void
@@ -250,8 +255,9 @@ class StraightBoundary : public Boundary<dim> {
                                      * @p{(m+1)(m+1)} subquads of equal
                                      * size.
                                      *
-                                     * Refer to the general documentation of
-                                     * this class and the documentation of the
+                                     * Refer to the general
+                                     * documentation of this class
+                                     * and the documentation of the
                                      * base class.
                                      */
     virtual void
index 670f8a122d7a9fe8c14deba3764b1cd4c103319c..e9d98c7a3ff529af00e1ae0fec1770132425e02b 100644 (file)
@@ -60,7 +60,9 @@ FE_DGQ<dim>::FE_DGQ (unsigned int degree)
       poly = new TensorProductPolynomials<dim> (v);
     }
 
-                                  // generate permutation/rotation index sets to generate some matrices from others
+                                  // generate permutation/rotation
+                                  // index sets to generate some
+                                  // matrices from others
   std::vector<unsigned int> right;
   std::vector<unsigned int> top;
   rotate_indices (right, 'Z');
index 449145ec4ab705fae2c1a0a7030ae1d6f02977e3..90b17638d29e817b196368e6a89e6fd6db58c69a 100644 (file)
@@ -49,6 +49,7 @@ FE_Q<dim>::FE_Q (const unsigned int degree)
 
                                   // do some internal book-keeping
   build_renumbering (*this, degree, renumber);
+//TODO: externalize this to a proper template function  
 #if deal_II_dimension > 1
   build_face_renumbering (FiniteElementData<dim-1>(FE_Q<dim-1>::get_dpo_vector(degree),1),
                          degree,
index 8452ebf17e3e624f60bbea0452496decab574f98..0a8fd3918d46165b68998051f570b8dc762a48b5 100644 (file)
@@ -18,7 +18,7 @@
 #if (deal_II_dimension == 1)
 
 template<>
-const unsigned int Mapping<deal_II_dimension>::normal_directions[2] =
+const unsigned int Mapping<1>::normal_directions[2] =
 {
   1, 0
 };
@@ -28,7 +28,7 @@ const unsigned int Mapping<deal_II_dimension>::normal_directions[2] =
 #if (deal_II_dimension == 2)
 
 template<>
-const unsigned int Mapping<deal_II_dimension>::normal_directions[4] =
+const unsigned int Mapping<2>::normal_directions[4] =
 {
   2, 0, 3, 1
 };
@@ -38,7 +38,7 @@ const unsigned int Mapping<deal_II_dimension>::normal_directions[4] =
 #if (deal_II_dimension == 3)
 
 template<>
-const unsigned int Mapping<deal_II_dimension>::normal_directions[6] =
+const unsigned int Mapping<3>::normal_directions[6] =
 {
   3, 2, 5, 0, 4, 1
 };
@@ -51,12 +51,7 @@ Mapping<dim>::~Mapping ()
 {}
 
 
-/*------------------------------ InternalData ------------------------------*/
-
-
-template <int dim>
-Mapping<dim>::InternalDataBase::~InternalDataBase ()
-{}
+/*------------------------------ InternalDataBase ------------------------------*/
 
 
 template <int dim>
@@ -69,6 +64,11 @@ Mapping<dim>::InternalDataBase::InternalDataBase ():
 
 
 
+template <int dim>
+Mapping<dim>::InternalDataBase::~InternalDataBase ()
+{}
+
+
 /*------------------------------ InternalData ------------------------------*/
 
 
index 49fad8d038af6300cf6e2c860d68829061b8911d..805a8a4515da6791b03d924b6c9943516c9b50ce 100644 (file)
 
 
 
+template <int dim>
+const unsigned int MappingCartesian<dim>::invalid_face_number;
+
+
+
 template<int dim>
 MappingCartesian<dim>::InternalData::InternalData (const Quadrature<dim>& q)
                :
-               length (dim, 0.),
                quadrature_points (q.get_points ())
 {}
 
@@ -42,6 +46,8 @@ MappingCartesian<dim>::update_once (const UpdateFlags) const
   return update_default;
 }
 
+
+
 template <int dim>
 UpdateFlags
 MappingCartesian<dim>::update_each (const UpdateFlags in) const
@@ -54,10 +60,11 @@ MappingCartesian<dim>::update_each (const UpdateFlags in) const
 }
 
 
+
 template <int dim>
 Mapping<dim>::InternalDataBase*
-MappingCartesian<dim>::get_data (const UpdateFlags flags,
-                         const Quadrature<dim>& q) const
+MappingCartesian<dim>::get_data (const UpdateFlags      flags,
+                                const Quadrature<dim> &q) const
 {
   Assert (flags & update_normal_vectors == 0, ExcNotImplemented());
   InternalData* data = new InternalData (q);
@@ -68,7 +75,7 @@ MappingCartesian<dim>::get_data (const UpdateFlags flags,
 template <int dim>
 Mapping<dim>::InternalDataBase*
 MappingCartesian<dim>::get_face_data (const UpdateFlags,
-                              const Quadrature<dim-1>& quadrature) const
+                                     const Quadrature<dim-1>& quadrature) const
 {
   QProjector<dim> q (quadrature, false);
   InternalData* data = new InternalData (q);
@@ -80,7 +87,7 @@ MappingCartesian<dim>::get_face_data (const UpdateFlags,
 template <int dim>
 Mapping<dim>::InternalDataBase*
 MappingCartesian<dim>::get_subface_data (const UpdateFlags,
-                                 const Quadrature<dim-1>& quadrature) const
+                                        const Quadrature<dim-1> &quadrature) const
 {
   QProjector<dim> q (quadrature, true);
   InternalData* data = new InternalData (q);
@@ -93,45 +100,55 @@ MappingCartesian<dim>::get_subface_data (const UpdateFlags,
 template <int dim>
 void
 MappingCartesian<dim>::compute_fill (const typename DoFHandler<dim>::cell_iterator &cell,
-                                    const unsigned int face_no,
-                                    const unsigned int sub_no,
-                                    InternalDatadata,
+                                    const unsigned int        face_no,
+                                    const unsigned int        sub_no,
+                                    InternalData             &data,
                                     std::vector<Point<dim> > &quadrature_points,
-                                    std::vector<Point<dim> >normal_vectors) const
+                                    std::vector<Point<dim> > &normal_vectors) const
 {
-  const UpdateFlags update_flags(data.current_update_flags());
+//TODO: does it make sense to query the update flags in data if we did not set any in the get_*_data functions?
+  UpdateFlags update_flags(data.current_update_flags());
 
   const unsigned int npts = quadrature_points.size ();
   unsigned int offset = 0;
-  bool onface = false;
   
-  if (face_no == static_cast<unsigned int> (-1))
+  if (face_no != invalid_face_number)
     {
-      Assert (sub_no == static_cast<unsigned int> (-1), ExcInternalError());
-    } else {
-      onface = true;
                                       // Add 1 on both sides of
                                       // assertion to avoid compiler
                                       // warning about testing
                                       // unsigned int < 0 in 1d.
       Assert (face_no+1 < GeometryInfo<dim>::faces_per_cell+1,
              ExcIndexRange (face_no, 0, GeometryInfo<dim>::faces_per_cell));
-      if (sub_no == static_cast<unsigned int> (-1))
+
+      if (sub_no == invalid_face_number)
+                                        // called from FEFaceValues
        offset = face_no * quadrature_points.size();
       else
        {
+                                          // called from FESubfaceValue
          Assert (sub_no+1 < GeometryInfo<dim>::subfaces_per_face+1,
                  ExcIndexRange (sub_no, 0, GeometryInfo<dim>::subfaces_per_face));
          offset = (face_no * GeometryInfo<dim>::subfaces_per_face + sub_no)
                   * quadrature_points.size();
        }
     }
+  else
+                                    // invalid face number, so
+                                    // subface should be invalid as
+                                    // well
+    Assert (sub_no == invalid_face_number, ExcInternalError());
+
+                                  // let @p{start} be the origin of a
+                                  // local coordinate system. it is
+                                  // chosen as the (lower) left
+                                  // vertex
+  const Point<dim> start = cell->vertex(0);
   
-                                  // Compute start point and sizes along axes.
-                                  // Strange vertex numbering makes this complicated again.
-  
-  Point<dim> start = cell->vertex (0);
-
+                                  // Compute start point and sizes
+                                  // along axes.  Strange vertex
+                                  // numbering makes this complicated
+                                  // again.
   switch (dim)
     {
       case 1:
@@ -150,23 +167,34 @@ MappingCartesian<dim>::compute_fill (const typename DoFHandler<dim>::cell_iterat
        Assert(false, ExcNotImplemented());
     }
   
-  
+
+                                  // transform quadrature point. this
+                                  // is obtained simply by scaling
+                                  // unit coordinates with lengths in
+                                  // each direction
   if (update_flags & update_q_points)
     {
       Assert (quadrature_points.size() == npts,
              ExcDimensionMismatch(quadrature_points.size(), npts));
-      for (unsigned int i=0;i<npts;++i)
+      for (unsigned int i=0; i<npts; ++i)
        {
-         Point<dim> p = start;
-         for (unsigned int d=0;d<dim;++d)
-           p(d) += data.length[d]*data.quadrature_points[i+offset](d);
-         quadrature_points[i] = p;
+         quadrature_points[i] = start;
+         for (unsigned int d=0; d<dim; ++d)
+           quadrature_points[i](d) += data.length[d] *
+                                      data.quadrature_points[i+offset](d);
        }
     }
 
-  
+
+                                  // compute normal vectors. since
+                                  // cells are aligned to coordinate
+                                  // axes, they are simply vectors
+                                  // with exactly one entry equal to
+                                  // 1 or -1
   if (update_flags & update_normal_vectors)
     {
+      Assert (normal_vectors.size() == npts,
+             ExcDimensionMismatch(normal_vectors.size(), npts));
       Point<dim> n;
       switch (100*dim+face_no)
        {
@@ -205,8 +233,9 @@ MappingCartesian<dim>::compute_fill (const typename DoFHandler<dim>::cell_iterat
          default:
                Assert (false, ExcInternalError());
        }
-      for (unsigned int i=0;i<npts;++i)
-       normal_vectors[i] = n;
+                                      // furthermore, all normal
+                                      // vectors on a face are equal
+      std::fill_n (normal_vectors.begin(), normal_vectors.end(), n);
     }
   
   data.first_cell = false;
@@ -221,17 +250,25 @@ MappingCartesian<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator& cel
                                std::vector<Point<dim> >& quadrature_points,
                                std::vector<double>& JxW_values) const
 {
-  InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
-  Assert(data_ptr!=0, ExcInternalError());
-  InternalData &data=*data_ptr;
+                                  // convert data object to internal
+                                  // data for this class. fails with
+                                  // an exception if that is not
+                                  // possible
+  InternalData &data = dynamic_cast<InternalData&> (mapping_data);
 
   typename std::vector<Point<dim> > dummy;
   
-  compute_fill (cell, static_cast<unsigned int> (-1), static_cast<unsigned int> (-1),
+  compute_fill (cell, invalid_face_number, invalid_face_number,
                data,
                quadrature_points,
                dummy);
 
+                                  // compute Jacobian
+                                  // determinant. all values are
+                                  // equal and are the product of the
+                                  // local lengths in each coordinate
+                                  // direction
+//TODO: does it make sense to query the update flags in data if we did not set any in the get_*_data functions?
   if (data.current_update_flags() & update_JxW_values)
     {
       double J = data.length[0];
@@ -247,34 +284,42 @@ MappingCartesian<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator& cel
 template <int dim>
 void
 MappingCartesian<dim>::fill_fe_face_values (const typename DoFHandler<dim>::cell_iterator &cell,
-                                    const unsigned int       face_no,
-                                    const Quadrature<dim-1> &q,
-                                    typename Mapping<dim>::InternalDataBase &mapping_data,
-                                    std::vector<Point<dim> >     &quadrature_points,
-                                    std::vector<double>          &JxW_values,
-                                    std::vector<Tensor<1,dim> >  &boundary_forms,
-                                    std::vector<Point<dim> >     &normal_vectors) const
+                                           const unsigned int            face_no,
+                                           const Quadrature<dim-1>      &q,
+                                           typename Mapping<dim>::InternalDataBase &mapping_data,
+                                           std::vector<Point<dim> >     &quadrature_points,
+                                           std::vector<double>          &JxW_values,
+                                           std::vector<Tensor<1,dim> >  &boundary_forms,
+                                           std::vector<Point<dim> >     &normal_vectors) const
 {
-  InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
-  Assert(data_ptr!=0, ExcInternalError());
-  InternalData &data=*data_ptr;
+                                  // convert data object to internal
+                                  // data for this class. fails with
+                                  // an exception if that is not
+                                  // possible
+  InternalData &data = dynamic_cast<InternalData&> (mapping_data);
 
-  compute_fill (cell, face_no, static_cast<unsigned int> (-1),
+  compute_fill (cell, face_no, invalid_face_number,
                data,
                quadrature_points,
                normal_vectors);
 
+                                  // first compute Jacobian
+                                  // determinant, which is simply the
+                                  // product of the local lengths
+                                  // since the jacobian is diagonal
   double J = 1.;
   for (unsigned int d=0;d<dim;++d)
     if (d != (normal_directions[face_no]/2))
       J *= data.length[d];
   
+//TODO: does it make sense to query the update flags in data if we did not set any in the get_*_data functions?
   if (data.current_update_flags() & update_JxW_values)
     {
       for (unsigned int i=0; i<JxW_values.size();++i)
        JxW_values[i] = J * q.weight(i);
     }
 
+//TODO: does it make sense to query the update flags in data if we did not set any in the get_*_data functions?
   if (data.current_update_flags() & update_boundary_forms)
     {
       for (unsigned int i=0; i<boundary_forms.size();++i)
@@ -295,26 +340,34 @@ MappingCartesian<dim>::fill_fe_subface_values (const typename DoFHandler<dim>::c
                                        std::vector<Tensor<1,dim> >  &boundary_forms,
                                        std::vector<Point<dim> >     &normal_vectors) const
 {
-  InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
-  Assert(data_ptr!=0, ExcInternalError());
-  InternalData &data=*data_ptr;
+                                  // convert data object to internal
+                                  // data for this class. fails with
+                                  // an exception if that is not
+                                  // possible
+  InternalData &data = dynamic_cast<InternalData&> (mapping_data);
 
   compute_fill (cell, face_no, sub_no,
                data,
                quadrature_points,
                normal_vectors);
 
+                                  // first compute Jacobian
+                                  // determinant, which is simply the
+                                  // product of the local lengths
+                                  // since the jacobian is diagonal
   double J = 1.;
   for (unsigned int d=0;d<dim;++d)
     if (d != (normal_directions[face_no]/2))
       J *= data.length[d];
   
+//TODO: does it make sense to query the update flags in data if we did not set any in the get_*_data functions?
   if (data.current_update_flags() & update_JxW_values)
     {
       for (unsigned int i=0; i<JxW_values.size();++i)
        JxW_values[i] = J * q.weight(i) / GeometryInfo<dim>::subfaces_per_face;
     }
 
+//TODO: does it make sense to query the update flags in data if we did not set any in the get_*_data functions?
   if (data.current_update_flags() & update_boundary_forms)
     {
       for (unsigned int i=0; i<boundary_forms.size();++i)
@@ -328,13 +381,13 @@ MappingCartesian<dim>::fill_fe_subface_values (const typename DoFHandler<dim>::c
 template <>
 void
 MappingCartesian<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator &,
-                                  const unsigned,
-                                  const Quadrature<0>&,
-                                  Mapping<1>::InternalDataBase&,
-                                  std::vector<Point<1> >&,
-                                  std::vector<double>&,
-                                  std::vector<Tensor<1,1> >&,
-                                  std::vector<Point<1> >&) const
+                                         const unsigned,
+                                         const Quadrature<0>&,
+                                         Mapping<1>::InternalDataBase&,
+                                         std::vector<Point<1> >&,
+                                         std::vector<double>&,
+                                         std::vector<Tensor<1,1> >&,
+                                         std::vector<Point<1> >&) const
 {
   Assert(false, ExcNotImplemented());
 }
@@ -343,14 +396,14 @@ MappingCartesian<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator &,
 template <>
 void
 MappingCartesian<1>::fill_fe_subface_values (const DoFHandler<1>::cell_iterator &,
-                                     const unsigned,
-                                     const unsigned,
-                                     const Quadrature<0>&,
-                                     Mapping<1>::InternalDataBase&,
-                                     std::vector<Point<1> >&,
-                                     std::vector<double>&,
-                                     std::vector<Tensor<1,1> >&,
-                                     std::vector<Point<1> >&) const
+                                            const unsigned,
+                                            const unsigned,
+                                            const Quadrature<0>&,
+                                            Mapping<1>::InternalDataBase&,
+                                            std::vector<Point<1> >&,
+                                            std::vector<double>&,
+                                            std::vector<Tensor<1,1> >&,
+                                            std::vector<Point<1> >&) const
 {
   Assert(false, ExcNotImplemented());
 }
@@ -407,6 +460,8 @@ Point<dim> MappingCartesian<dim>::transform_unit_to_real_cell (
   const Point<dim> &p,
   const typename Mapping<dim>::InternalDataBase *const m_data) const
 {
+//TODO: wouldn't it be simpler in this function to compute local lengths ourselved, rather than relying on an InternalDataBase object?
+  
                                   // If m_data!=0 use this
                                   // InternalData.
                                   //
@@ -432,12 +487,19 @@ Point<dim> MappingCartesian<dim>::transform_unit_to_real_cell (
     mdata = dynamic_cast<const InternalData *> (m_data);
   Assert(mdata!=0, ExcInternalError());
 
+//TODO: don't we need to reinit() mdata on the present cell or initialize the mdata->length array in some other way?
+  
                                   // use now the InternalData, that
                                   // mdata is pointing to, to compute
-                                  // the point in real space.
+                                  // the point in real space. it is
+                                  // simply a scaled version of the
+                                  // unit cell point shifted by the
+                                  // (lower) left corner of the cell
   Point<dim> p_real = cell->vertex(0);
-  for (unsigned int d=0;d<dim;++d)
+  for (unsigned int d=0; d<dim; ++d)
     p_real(d) += mdata->length[d]*p(d);
+
+//TODO: memory leak: if initially mdata==0, then an object is created but never deleted. see mapping_q1 for a better implementation  
   return p_real;
 }
 
@@ -448,7 +510,7 @@ Point<dim> MappingCartesian<dim>::transform_real_to_unit_cell (
   const typename Triangulation<dim>::cell_iterator cell,
   const Point<dim> &p) const
 {
-  const Point<dim>& start =  cell->vertex (0);
+  const Point<dim> &start = cell->vertex(0);
   Point<dim> real = p;
   real -= start;
 
@@ -479,53 +541,64 @@ template <typename tensor_>
 inline
 void
 MappingCartesian<dim>::contravariant_transformation (std::vector<tensor_>       &dst,
-                                             const std::vector<tensor_> &src,
-                                             const Mapping<dim>::InternalDataBase &mapping_data,
-                                             const unsigned int src_offset) const
+                                                    const std::vector<tensor_> &src,
+                                                    const Mapping<dim>::InternalDataBase &mapping_data,
+                                                    const unsigned int src_offset) const
 {
   Assert(tensor_::dimension==dim && tensor_::rank==1, ExcInvalidData());
-  const InternalData* data_ptr = dynamic_cast<const InternalData *> (&mapping_data);
-  Assert(data_ptr!=0, ExcInternalError());
-  const InternalData &data=*data_ptr;
+
+                                  // convert data object to internal
+                                  // data for this class. fails with
+                                  // an exception if that is not
+                                  // possible
+  const InternalData &data = dynamic_cast<const InternalData&> (mapping_data);
 
   Assert (data.update_flags & update_contravariant_transformation,
          typename FEValuesBase<dim>::ExcAccessToUninitializedField());
 
-  typename std::vector<tensor_>::const_iterator vec = src.begin() + src_offset;
-  typename std::vector<tensor_>::iterator result = dst.begin();
-  typename std::vector<tensor_>::const_iterator end = dst.end();
+  typename std::vector<tensor_>::const_iterator vec    = src.begin() + src_offset;
+  typename std::vector<tensor_>::iterator       result = dst.begin();
+  typename std::vector<tensor_>::const_iterator end    = dst.end();
   
+                                  // simply scale by Jacobian
+                                  // (which is diagonal here)
   while (result!=end)
     {
-      for (unsigned int d=0;d<dim;++d)
+      for (unsigned int d=0; d<dim; ++d)
        (*result)[d] = (*vec)[d]*data.length[d];
-      vec++;
-      result++;
+      ++vec;
+      ++result;
     }
 }
 
 
+
 template <int dim>
 template <typename tensor_>
 inline
 void
 MappingCartesian<dim>::covariant_transformation (std::vector<tensor_>       &dst,
-                                         const std::vector<tensor_> &src,
-                                         const Mapping<dim>::InternalDataBase &mapping_data,
-                                         const unsigned int src_offset) const
+                                                const std::vector<tensor_> &src,
+                                                const Mapping<dim>::InternalDataBase &mapping_data,
+                                                const unsigned int src_offset) const
 {
   Assert(tensor_::dimension==dim && tensor_::rank==1, ExcInvalidData());
-  const InternalData *data_ptr = dynamic_cast<const InternalData *> (&mapping_data);
-  Assert(data_ptr!=0, ExcInternalError());
-  const InternalData &data=*data_ptr;
+
+                                  // convert data object to internal
+                                  // data for this class. fails with
+                                  // an exception if that is not
+                                  // possible
+  const InternalData &data = dynamic_cast<const InternalData&> (mapping_data);
 
   Assert (data.update_flags & update_covariant_transformation,
          typename FEValuesBase<dim>::ExcAccessToUninitializedField());
   
-  typename std::vector<tensor_>::const_iterator vec = src.begin() + src_offset;
-  typename std::vector<tensor_>::iterator result = dst.begin();
-  typename std::vector<tensor_>::const_iterator end = dst.end();
-  
+  typename std::vector<tensor_>::const_iterator vec    = src.begin() + src_offset;
+  typename std::vector<tensor_>::iterator       result = dst.begin();
+  typename std::vector<tensor_>::const_iterator end    = dst.end();
+
+                                  // simply scale by inverse Jacobian
+                                  // (which is diagonal here)
   while (result!=end)
     {
       for (unsigned int d=0;d<dim;++d)
@@ -538,6 +611,7 @@ MappingCartesian<dim>::covariant_transformation (std::vector<tensor_>       &dst
 
 
 //----------------------------------------------------------------------//
+// explicit instantiations
 
 template class MappingCartesian<deal_II_dimension>;
 
index f42fe23efb3e9036108eca72c042856b4fecfb38..269673debcd37d3fd359758fcf0f167a3f1b2040 100644 (file)
 
 
 template<int dim>
-MappingQ<dim>::InternalData::InternalData (unsigned int n_shape_functions):
+MappingQ<dim>::InternalData::InternalData (const unsigned int n_shape_functions):
                MappingQ1<dim>::InternalData(n_shape_functions),
-                               use_mapping_q1_on_current_cell(false),
-                               mapping_q1_data(1 << dim)
+//TODO: in 1d, use_mapping_q1_on_current_cell is always false, but should be true.
+               use_mapping_q1_on_current_cell(false),
+               mapping_q1_data(1 << dim)
 {
   is_mapping_q1_data=false;
 }
@@ -38,8 +39,10 @@ MappingQ<dim>::InternalData::InternalData (unsigned int n_shape_functions):
 
 #if deal_II_dimension == 1
 
+// in 1d, it is irrelevant which polynomial degree to use, since all
+// cells are scaled linearly
 template<>
-MappingQ<1>::MappingQ (unsigned int):
+MappingQ<1>::MappingQ (const unsigned int):
                laplace_on_quad_vector(0),
                laplace_on_hex_vector(0),
                degree(1),
@@ -94,7 +97,7 @@ MappingQ<dim>::MappingQ (const unsigned int p):
                                   // functions for the Qp mapping of
                                   // cells at the boundary.
   std::vector<LagrangeEquidistant> v;
-  for (unsigned int i=0;i<=degree;++i)
+  for (unsigned int i=0; i<=degree; ++i)
     v.push_back(LagrangeEquidistant(degree,i));
 
   tensor_pols = new TensorProductPolynomials<dim> (v);
@@ -107,15 +110,16 @@ MappingQ<dim>::MappingQ (const unsigned int p):
   renumber.resize(n_shape_functions,0);
   std::vector<unsigned int> dpo(dim+1, 1);
   for (unsigned int i=1; i<dpo.size(); ++i)
-    dpo[i]=dpo[i-1]*(degree-1);
+    dpo[i] = dpo[i-1]*(degree-1);
   FiniteElementData<dim> fe_data(dpo, 1);
   FE_Q<dim>::build_renumbering (fe_data, p, renumber);
 
                                   // build laplace_on_quad_vector
   if (degree>1)
     {
-      set_laplace_on_quad_vector(laplace_on_quad_vector);
-      if (dim==3)
+      if (dim >= 2)
+       set_laplace_on_quad_vector(laplace_on_quad_vector);
+      if (dim >= 3)
        set_laplace_on_hex_vector(laplace_on_hex_vector);
     }
 }
@@ -134,7 +138,7 @@ MappingQ<dim>::~MappingQ ()
 template<>
 void
 MappingQ<1>::compute_shapes_virtual (const std::vector<Point<1> > &unit_points,
-                                    MappingQ1<1>::InternalData &data) const
+                                    MappingQ1<1>::InternalData   &data) const
 {
   MappingQ1<1>::compute_shapes_virtual(unit_points, data);
 }
@@ -199,6 +203,7 @@ MappingQ<dim>::update_each (const UpdateFlags in) const
 }
 
 
+
 template <int dim>
 Mapping<dim>::InternalDataBase*
 MappingQ<dim>::get_data (const UpdateFlags update_flags,
@@ -260,13 +265,16 @@ MappingQ<dim>::compute_face_data (UpdateFlags update_flags,
                                  const unsigned int n_original_q_points,
                                  MappingQ1<dim>::InternalData& mapping_q1_data) const
 {
-  InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_q1_data);
-  Assert(data_ptr!=0, ExcInternalError());
-  InternalData &data=*data_ptr;
+                                  // convert data object to internal
+                                  // data for this class. fails with
+                                  // an exception if that is not
+                                  // possible
+  InternalData &data = dynamic_cast<InternalData&> (mapping_q1_data);
   
   MappingQ1<dim>::compute_face_data(update_flags, q,
                                    n_original_q_points, data);
-  
+
+//TODO: externalize this to a proper template function. or: scrap the whole function for 1d since it doesn't make much sense anyway?
 #if (deal_II_dimension>1)
   if ((data.update_flags & update_normal_vectors)
       && alternative_normals_computation)
@@ -310,13 +318,15 @@ template <int dim>
 void
 MappingQ<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
                               const Quadrature<dim>                &q,
-                              Mapping<dim>::InternalDataBase      &mapping_data,
-                              std::vector<Point<dim> >                  &quadrature_points,
-                              std::vector<double>                       &JxW_values) const
+                              Mapping<dim>::InternalDataBase       &mapping_data,
+                              std::vector<Point<dim> >             &quadrature_points,
+                              std::vector<double>                  &JxW_values) const
 {
-  InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
-  Assert(data_ptr!=0, ExcInternalError());
-  InternalData &data=*data_ptr;
+                                  // convert data object to internal
+                                  // data for this class. fails with
+                                  // an exception if that is not
+                                  // 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());
@@ -341,10 +351,13 @@ MappingQ<dim>::fill_fe_face_values (const typename DoFHandler<dim>::cell_iterato
                                    std::vector<Tensor<1,dim> >  &exterior_forms,
                                    std::vector<Point<dim> >     &normal_vectors) const
 {
-  InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
-  Assert(data_ptr!=0, ExcInternalError());
-  InternalData &data=*data_ptr;
+                                  // convert data object to internal
+                                  // data for this class. fails with
+                                  // an exception if that is not
+                                  // possible
+  InternalData &data = dynamic_cast<InternalData&> (mapping_data);
   
+//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());
 
@@ -378,10 +391,13 @@ MappingQ<dim>::fill_fe_subface_values (const typename DoFHandler<dim>::cell_iter
                                       std::vector<Tensor<1,dim> >  &exterior_forms,
                                       std::vector<Point<dim> >     &normal_vectors) const
 {
-  InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
-  Assert(data_ptr!=0, ExcInternalError());
-  InternalData &data=*data_ptr;
+                                  // convert data object to internal
+                                  // data for this class. fails with
+                                  // an exception if that is not
+                                  // possible
+  InternalData &data = dynamic_cast<InternalData&> (mapping_data);
 
+//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());
 
@@ -422,7 +438,10 @@ MappingQ<dim>::set_laplace_on_quad_vector(std::vector<std::vector<double> > &loq
   Assert(degree>1, ExcInternalError());
   const unsigned int n_inner_2d=(degree-1)*(degree-1);
   const unsigned int n_outer_2d=4+4*(degree-1);
-  
+
+                                  // first check whether we have
+                                  // precomputed the values for some
+                                  // polynomial degree
   double const *loqv_ptr=0;
   if (degree==2)
     {
@@ -432,7 +451,7 @@ MappingQ<dim>::set_laplace_on_quad_vector(std::vector<std::vector<double> > &loq
     }
   else if (degree==3)
     {
-      static double loqv3[4*12]
+      static const double loqv3[4*12]
        ={80/1053., 1/81., 11/1053., 1/81., 25/117., 44/351.,
          7/117., 16/351., 7/117., 16/351., 25/117., 44/351.,
          1/81., 80/1053., 1/81., 11/1053., 44/351., 25/117.,
@@ -444,29 +463,33 @@ MappingQ<dim>::set_laplace_on_quad_vector(std::vector<std::vector<double> > &loq
       
       loqv_ptr=&loqv3[0];
     }
-  
+
   if (loqv_ptr!=0)
     {
+                                      // precomputed. copy values to
+                                      // the loqvs array
       loqvs.resize(n_inner_2d);
       for (unsigned int unit_point=0; unit_point<n_inner_2d; ++unit_point)
        {
          loqvs[unit_point].resize(n_outer_2d, 0);
-         std::vector<double> &loqv=loqvs[unit_point];
          for (unsigned int k=0; k<n_outer_2d; ++k)
-           loqv[k]=loqv_ptr[unit_point*n_outer_2d+k];
+           loqvs[unit_point][k]=loqv_ptr[unit_point*n_outer_2d+k];
        }
     }
   else
     {
+                                  // not precomputed, then do so now
       if (dim==2)
        compute_laplace_vector(loqvs);
+//TODO: what if dim==3?
     }
 
-  
+                                  // the sum of weights of the points
+                                  // at the outer rim should be
+                                  // one. check this
   for (unsigned int unit_point=0; unit_point<loqvs.size(); ++unit_point)
-    Assert(fabs(accumulate(
-      loqvs[unit_point].begin(),
-      loqvs[unit_point].end(),0.)-1)<1e-13,
+    Assert(fabs(std::accumulate(loqvs[unit_point].begin(),
+                               loqvs[unit_point].end(),0.)-1)<1e-13,
           ExcInternalError());
   
                                   // TEST output
@@ -490,6 +513,9 @@ MappingQ<3>::set_laplace_on_hex_vector(std::vector<std::vector<double> > &lohvs)
 {
   Assert(degree>1, ExcInternalError());
 
+                                  // first check whether we have
+                                  // precomputed the values for some
+                                  // polynomial degree
   double const *lohv_ptr=0;
   if (degree==2)
     {
@@ -504,22 +530,26 @@ MappingQ<3>::set_laplace_on_hex_vector(std::vector<std::vector<double> > &lohvs)
   
   if (lohv_ptr!=0)
     {
+                                      // precomputed. copy values to
+                                      // the lohvs array
       lohvs.resize(n_inner);
       for (unsigned int unit_point=0; unit_point<n_inner; ++unit_point)
        {
          lohvs[unit_point].resize(n_outer, 0);
-         std::vector<double> &loqv=lohvs[unit_point];
          for (unsigned int k=0; k<n_outer; ++k)
-           loqv[k]=lohv_ptr[unit_point*n_outer+k];
+           lohvs[unit_point][k]=lohv_ptr[unit_point*n_outer+k];
        }
     }
   else
+                                  // not precomputed, then do so now
     compute_laplace_vector(lohvs);
     
+                                  // the sum of weights of the points
+                                  // at the outer rim should be
+                                  // one. check this
   for (unsigned int unit_point=0; unit_point<n_inner; ++unit_point)
-    Assert(fabs(accumulate(
-      lohvs[unit_point].begin(),
-      lohvs[unit_point].end(),0.)-1)<1e-13,
+    Assert(fabs(std::accumulate(lohvs[unit_point].begin(),
+                               lohvs[unit_point].end(),0.) - 1)<1e-13,
           ExcInternalError());
   
                                   // TEST output
@@ -568,12 +598,11 @@ MappingQ<dim>::compute_laplace_vector(std::vector<std::vector<double> > &lvs) co
                                   // compute the shape
                                   // gradients at the quadrature
                                   // points on the unit cell
-  QGauss4<dim> quadrature;
+  const QGauss4<dim> quadrature;
   const unsigned int n_q_points=quadrature.n_quadrature_points;
   
   InternalData quadrature_data(n_shape_functions);
-  quadrature_data.shape_derivatives
-    .resize(n_shape_functions * n_q_points);
+  quadrature_data.shape_derivatives.resize(n_shape_functions * n_q_points);
   compute_shapes(quadrature.get_points(), quadrature_data);
   
                                   // Compute the stiffness matrix of
@@ -706,10 +735,10 @@ void
 MappingQ<dim>::compute_support_points_laplace(const typename Triangulation<dim>::cell_iterator &cell,
                                              std::vector<Point<dim> > &a) const
 {
-  a.resize(0);
+  a.resize(GeometryInfo<dim>::vertices_per_cell);
                                   // the vertices first
   for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-    a.push_back(cell->vertex(i));
+    a[i] = cell->vertex(i);
   
   if (degree>1)
     {
@@ -742,6 +771,7 @@ MappingQ<dim>::compute_support_points_simple(const typename Triangulation<dim>::
                                             std::vector<Point<dim> > &a) const
 {
                                   // the vertices first
+// TODO: check for size of 'a' first?
   for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
     a.push_back(cell->vertex(i));
   
@@ -755,64 +785,67 @@ MappingQ<dim>::compute_support_points_simple(const typename Triangulation<dim>::
                                       // (for dim=3)
       fill_quad_support_points_simple (cell, a);
       
-                                      // then the points in cell    
-      Point<dim> middle;
-      compute_midpoint(a, middle);
-      
-      if (degree==2)
-       a.push_back(middle);
-      else if (degree==3)
-                                        // The four points in the
-                                        // cell are located at the
-                                        // midpoint between the
-                                        // middle point and the 4
-                                        // vertices
-
-                                        // TODO: better position of
-                                        // points: transform them by
-                                        // a Q2 transformation.
-       for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-         a.push_back(middle*2./3.+cell->vertex(vertex_mapping[i])/3.);
-      else if (degree==4)
-       {
-         Assert(a.size()==16, ExcInternalError());
-         a.insert(a.end(), 9, Point<dim>());
+                                      // then the points in cell. for
+                                      // this we need the midpoint of
+                                      // the points already in @p{a}
+      const Point<dim> middle = std::accumulate(a.begin(), a.end(), Point<dim>())
+                               / a.size();
 
-         const unsigned int inner_map[8]=
-         { 0, 1, 2, 5, 8, 7, 6, 3 };
-         
-         
-                                          // The nine points in the
-                                          // cell are located at the
-                                          // midpoint between the
-                                          // middle point and (the 4
-                                          // vertices and the face
-                                          // midpoints)
-
-         a[16+4]=middle;
-         for (unsigned int i=0, j=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-           {
-             a[16+inner_map[j++]]=(middle+cell->vertex(i))/2.;
-             a[16+inner_map[j++]]=(middle+(cell->vertex(i)+cell->vertex((i+1)%4))/2.)/2.;
-           }
-       }
-      else
-       Assert(false, ExcNotImplemented());
-    }
+      switch (degree)
+       {
+         case 2:
+         {
+           a.push_back(middle);
+           break;
+         };
+
+         case 3:
+         {
+                                            // The four points in the
+                                            // cell are located at
+                                            // the midpoint between
+                                            // the middle point and
+                                            // the 4 vertices
+           
+                                            // TODO: better position of
+                                            // points: transform them by
+                                            // a Q2 transformation.
+           for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+             a.push_back(middle*2./3.+cell->vertex(vertex_mapping[i])/3.);
+           break;
+         };
+
+         case 4:
+         {
+           Assert(a.size()==16, ExcInternalError());
+           a.insert(a.end(), 9, Point<dim>());
+           
+           const unsigned int inner_map[8]=
+           { 0, 1, 2, 5, 8, 7, 6, 3 };
+           
+           
+                                            // The nine points in the
+                                            // cell are located at the
+                                            // midpoint between the
+                                            // middle point and (the 4
+                                            // vertices and the face
+                                            // midpoints)
+           
+           a[16+4]=middle;
+           for (unsigned int i=0, j=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+             {
+               a[16+inner_map[j++]]=(middle+cell->vertex(i))/2.;
+               a[16+inner_map[j++]]=(middle+(cell->vertex(i)+cell->vertex((i+1)%4))/2.)/2.;
+             }
+           break;
+         };
+
+         default:
+               Assert(false, ExcNotImplemented());
+       };
+    };
+  
   Assert(a.size()==n_shape_functions, ExcInternalError());
-
-}
-
-
-
-template <int dim> static
-void compute_midpoint(const std::vector<Point<dim> > &points, Point<dim> &m)
-{
-  typename std::vector<Point<dim> >::const_iterator p_iter=points.begin();
-  const typename std::vector<Point<dim> >::const_iterator p_end=points.end();
-  for (; p_iter!=p_end; ++p_iter)
-    m+=*p_iter;
-  m/=points.size();
 }
 
 
@@ -822,8 +855,11 @@ void compute_midpoint(const std::vector<Point<dim> > &points, Point<dim> &m)
 template <>
 void
 MappingQ<1>::add_line_support_points (const Triangulation<1>::cell_iterator &,
-                                  std::vector<Point<1> > &) const
-{}
+                                     std::vector<Point<1> > &) const
+{
+                                  // there are no points on bounding
+                                  // lines which are to be added
+}
 
 #endif
 
@@ -831,30 +867,37 @@ MappingQ<1>::add_line_support_points (const Triangulation<1>::cell_iterator &,
 template <int dim>
 void
 MappingQ<dim>::add_line_support_points (const Triangulation<dim>::cell_iterator &cell,
-                                    std::vector<Point<dim> > &a) const
+                                       std::vector<Point<dim> > &a) const
 {
-  const Boundary<dim> *boundary;
+  const Boundary<dim> *boundary = 0;
 
   std::vector<Point<dim> > line_points;
   if (degree>2)
     line_points.resize(degree-1);
-  
-  Triangulation<dim>::line_iterator line;
+
+                                  // 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)
     {
-      line = cell->line(line_no);
+      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)
        a.push_back(boundary->get_new_point_on_line(line));
       else
        {
          boundary->get_intermediate_points_on_line (line, line_points);
-         for (unsigned int i=0; i<line_points.size(); ++i)
-           a.push_back(line_points[i]);
+         a.insert (a.end(), line_points.begin(), line_points.end());
        } 
     }
 }
@@ -864,15 +907,16 @@ MappingQ<dim>::add_line_support_points (const Triangulation<dim>::cell_iterator
 
 #if deal_II_dimension==3
 
+//TODO: rename function to add_quad_support_points, to unify notation
 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,
-                    lines_per_face=GeometryInfo<3>::lines_per_face,
-                 vertices_per_cell=GeometryInfo<3>::vertices_per_cell;
+  const unsigned int faces_per_cell    = GeometryInfo<3>::faces_per_cell,
+                    vertices_per_face = GeometryInfo<3>::vertices_per_face,
+                    lines_per_face    = GeometryInfo<3>::lines_per_face,
+                    vertices_per_cell = GeometryInfo<3>::vertices_per_cell;
   
   static const unsigned int face_vertex_to_cell_vertex
     [faces_per_cell][vertices_per_face]={{0,1,2,3},
@@ -891,10 +935,11 @@ MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell
                                      {8,7,11,3}};
   
   
-  Triangulation<3>::face_iterator face;
+                                  // loop over all faces and collect points on them
   for (unsigned int face_no=0; face_no<faces_per_cell; ++face_no)
     {
-      face=cell->face(face_no);
+      const Triangulation<3>::face_iterator face=cell->face(face_no);
+      
       for (unsigned int i=0; i<vertices_per_face; ++i)
        Assert(face->vertex_index(i)==
               cell->vertex_index(face_vertex_to_cell_vertex[face_no][i]),
@@ -904,30 +949,41 @@ MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell
        Assert(face->line(i)==
               cell->line(face_line_to_cell_line[face_no][i]),
               ExcInternalError());
-      
+
+                                      // if face at boundary, then
+                                      // ask boundary object to
+                                      // return intermediate points
+                                      // on it
       if (face->at_boundary())
        {
-         std::vector<Point<3> > quad_points;
-         quad_points.resize((degree-1)*(degree-1));
+         std::vector<Point<3> > quad_points ((degree-1)*(degree-1));
 
          face->get_triangulation().get_boundary(face->boundary_indicator())
            .get_intermediate_points_on_quad (face, quad_points);
-         
-         for (unsigned int i=0; i<quad_points.size(); ++i)
-           a.push_back(quad_points[i]);
+         a.insert (a.end(), quad_points.begin(), quad_points.end());
        }
       else
        {
+                                          // face is not at boundary,
+                                          // but maybe some of its
+                                          // lines are. count them
          unsigned int lines_at_boundary=0;
          for (unsigned int i=0; i<lines_per_face; ++i)
            if (face->line(i)->at_boundary())
              ++lines_at_boundary;
          
          Assert(lines_at_boundary<lines_per_face, ExcInternalError());
-         
+
+                                          // if at least one of the
+                                          // lines bounding this quad
+                                          // is at the boundary, then
+                                          // collect points
+                                          // separately
          if (lines_at_boundary>0)
            {
                                               // sort the points into b
+//TODO: this is not thread-safe!!! b might be used for objects with
+//TODO: different degrees at the same time!
              static std::vector<Point<3> > b;
              b.resize(4*degree);
              Assert(4*degree==vertices_per_face+lines_per_face*(degree-1),
@@ -950,13 +1006,15 @@ MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell
            }
          else
            {
-             std::vector<Point<3> > quad_points;
-             quad_points.resize((degree-1)*(degree-1));
+                                              // face is entirely in
+                                              // the interior. get
+                                              // intermediate points
+                                              // from a straight
+                                              // boundary object
+             std::vector<Point<3> > quad_points ((degree-1)*(degree-1));
              
              straight_boundary.get_intermediate_points_on_quad (face, quad_points);
-             
-             for (unsigned int i=0; i<quad_points.size(); ++i)
-               a.push_back(quad_points[i]);
+             a.insert (a.end(), quad_points.begin(), quad_points.end());
            }
        }
     }
@@ -980,26 +1038,24 @@ MappingQ<dim>::add_face_support_points(const typename Triangulation<dim>::cell_i
 template <>
 void
 MappingQ<3>::fill_quad_support_points_simple (const Triangulation<3>::cell_iterator &cell,
-                                  std::vector<Point<3> > &a) const
+                                             std::vector<Point<3> > &a) const
 {
-  const Boundary<3> *boundary;
+  const Boundary<3> *boundary = 0;
 
   std::vector<Point<3> > quad_points;
   Assert(degree>1, ExcInternalError());
   quad_points.resize((degree-1)*(degree-1));
   
-  Triangulation<3>::quad_iterator quad;
   for (unsigned int quad_no=0; quad_no<GeometryInfo<3>::quads_per_cell; ++quad_no)
     {
-      quad = cell->face(quad_no);
+      const Triangulation<3>::quad_iterator quad = cell->face(quad_no);
       if (quad->at_boundary())
        boundary=&quad->get_triangulation().get_boundary(quad->boundary_indicator());
       else
        boundary=&straight_boundary;
 
       boundary->get_intermediate_points_on_quad (quad, quad_points);
-      for (unsigned int i=0; i<quad_points.size(); ++i)
-       a.push_back(quad_points[i]);
+      a.insert (a.end(), quad_points.begin(), quad_points.end());
     }
 }
 
@@ -1180,5 +1236,5 @@ MappingQ<dim>::transform_contravariant (std::vector<Point<dim> >       &dst,
 
 
 
-
+// explicit instantiation
 template class MappingQ<deal_II_dimension>;
index 79b70e7c79827c76b52d3a0cd77a75c1763a6f16..29713e2389a22f257d8789dd297a699b0a6116a9 100644 (file)
@@ -864,6 +864,7 @@ Point<dim> MappingQ1<dim>::transform_unit_to_real_cell (
       mdata=mdata_local;
     }
   else
+//TODO: can we assume that m_data->mapping_support_points is initialized    
     mdata = dynamic_cast<const InternalData *> (m_data);
   Assert(mdata!=0, ExcInternalError());
 

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.