Clean up code: clearly separate input arguments and return values, rather than using...

git-svn-id: https://svn.dealii.org/trunk@25580 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/include/deal.II/fe/mapping_q1.h b/deal.II/include/deal.II/fe/mapping_q1.h
index befd17bb2df84f36bfefa02ff3bc398fd15da32e..ce7d36d57422a5487fe792e0806e3d5fd7ae1adf 100644 (file)
--- a/deal.II/include/deal.II/fe/mapping_q1.h
+++ b/deal.II/include/deal.II/fe/mapping_q1.h
@@ -492,7 +492,7 @@ class MappingQ1 : public Mapping<dim,spacedim>
                                       *
                                       * This function is called by
                                       * @p transform_unit_to_real_cell
-                                      * and multiply (through the
+                                      * and multiple times (through the
                                       * Newton iteration) by
                                       * @p transform_real_to_unit_cell_internal.
                                       *
@@ -510,12 +510,13 @@ class MappingQ1 : public Mapping<dim,spacedim>
                                       * computations of the mapping
                                       * support points.
                                       */
-    Point<spacedim> transform_unit_to_real_cell_internal (const InternalData &mdata) const;
+    Point<spacedim>
+    transform_unit_to_real_cell_internal (const InternalData &mdata) const;
 
                                      /**
                                       * Transforms the point @p p on
-                                      * the real cell to the point
-                                      * @p p_unit on the unit cell
+                                      * the real cell to the corresponding
+                                      * point on the unit cell
                                       * @p cell by a Newton
                                       * iteration.
                                       *
@@ -529,34 +530,24 @@ class MappingQ1 : public Mapping<dim,spacedim>
                                       * and
                                       * @p update_transformation_gradients
                                       * and a one point Quadrature
-                                      * including the given point
-                                      * @p p_unit.  Hence this
+                                      * that includes the given
+                                      * initial guess for the
+                                      * transformation
+                                      * @p initial_p_unit.  Hence this
                                       * function assumes that
                                       * @p mdata already includes the
                                       * transformation shape values
                                       * and gradients computed at
-                                      * @p p_unit.
+                                      * @p initial_p_unit.
                                       *
-                                      * These assumptions should be
-                                      * fulfilled by the calling
-                                      * function. That is up to now
-                                      * only the function
-                                      * @p transform_real_to_unit_cell
-                                      * and its overloaded versions.
                                       * @p mdata will be changed by
                                       * this function.
                                       */
-    void transform_real_to_unit_cell_internal (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-                                               const Point<spacedim> &p,
-                                               InternalData &mdata,
-                                               Point<dim> &p_unit) const;
-
-
-
-
-
-
-
+    Point<dim>
+    transform_real_to_unit_cell_internal (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                                          const Point<spacedim> &p,
+                                          const Point<dim> &initial_p_unit,
+                                          InternalData &mdata) const;
 
                                        /**
                                         * Always returns @p true because
@@ -581,11 +572,12 @@ class MappingQ1 : public Mapping<dim,spacedim>
   protected:
                                      /* Trick to templatize transform_real_to_unit_cell<dim, dim+1> */
     template<int dim_>
-    void transform_real_to_unit_cell_internal_codim1
+    Point<dim_>
+    transform_real_to_unit_cell_internal_codim1
     (const typename Triangulation<dim_,dim_+1>::cell_iterator &cell,
      const Point<dim_+1> &p,
-     InternalData        &mdata,
-     Point<dim_>         &p_unit) const;
+     const Point<dim_>         &initial_p_unit,
+     InternalData        &mdata) const;
 
 /**
    Compute an initial guess to pass to the Newton method in
@@ -752,32 +744,34 @@ class MappingQ1 : public Mapping<dim,spacedim>
 };
 
 
+// explicit specializations
+
 template<>
-void
+Point<2>
 MappingQ1<2,3>::
 transform_real_to_unit_cell_internal
 (const Triangulation<2,3>::cell_iterator &cell,
- const Point<3>                            &p,
- InternalData                                     &mdata,
- Point<2>                                       &p_unit) const;
+ const Point<3> &p,
+ const Point<2> &initial_p_unit,
+ InternalData    &mdata) const;
 
-/* Only used in mapping Q if degree > 1 */
 template<>
-void
+Point<1>
 MappingQ1<1,2>::
-transform_real_to_unit_cell_internal (const Triangulation<1,2>::cell_iterator &cell,
-                                      const Point<2> &p,
-                                      InternalData   &mdata,
-                                      Point<1>       &p_unit) const;
+transform_real_to_unit_cell_internal
+(const Triangulation<1,2>::cell_iterator &cell,
+ const Point<2> &p,
+ const Point<1> &initial_p_unit,
+ InternalData    &mdata) const;
 
-/* Only used in mapping Q if degree > 1 */
 template<>
-void
+Point<1>
 MappingQ1<1,3>::
-transform_real_to_unit_cell_internal (const Triangulation<1,3>::cell_iterator &cell,
-                                      const Point<3> &p,
-                                      InternalData   &mdata,
-                                      Point<1>       &p_unit) const;
+transform_real_to_unit_cell_internal
+(const Triangulation<1,3>::cell_iterator &cell,
+ const Point<3> &p,
+ const Point<1> &initial_p_unit,
+ InternalData    &mdata) const;
 
 
 /**

diff --git a/deal.II/source/fe/mapping_q.cc b/deal.II/source/fe/mapping_q.cc
index 0beb8330d2042c4148d2466245e9f2329bf00345..1d12cbbd1673cac6d97234051c25b14aa241d72f 100644 (file)
--- a/deal.II/source/fe/mapping_q.cc
+++ b/deal.II/source/fe/mapping_q.cc
@@ -1447,7 +1447,7 @@ transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_it
 {
                                    // first a Newton iteration based on a Q1
                                    // mapping.
-  Point<dim> p_unit =
+  Point<dim> initial_p_unit =
     MappingQ1<dim,spacedim>::transform_real_to_unit_cell(cell, p);
 
                                    // then a Newton iteration based on the
@@ -1472,9 +1472,9 @@ transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_it
                                       // cell in hopes that this gives a
                                       // better starting point to the
                                       // following iteration
-      p_unit = GeometryInfo<dim>::project_to_unit_cell(p_unit);
+      initial_p_unit = GeometryInfo<dim>::project_to_unit_cell(initial_p_unit);
 
-      const Quadrature<dim> point_quadrature(p_unit);
+      const Quadrature<dim> point_quadrature(initial_p_unit);
 
       UpdateFlags update_flags = update_transformation_values|update_transformation_gradients;
       if (spacedim>dim)
@@ -1494,10 +1494,10 @@ transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_it
         mdata->mapping_support_points.resize(GeometryInfo<dim>::vertices_per_cell);
 
 
-      this->transform_real_to_unit_cell_internal(cell, p, *mdata, p_unit);
+      return this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit, *mdata);
     }
-
-  return p_unit;
+  else
+    return initial_p_unit;
 }
 
 

diff --git a/deal.II/source/fe/mapping_q1.cc b/deal.II/source/fe/mapping_q1.cc
index 87485d368c518939075e8d1e445910dee0890327..604be5579fb8c614b815335c29c0acc8f368df09 100644 (file)
--- a/deal.II/source/fe/mapping_q1.cc
+++ b/deal.II/source/fe/mapping_q1.cc
@@ -1566,7 +1566,6 @@ MappingQ1<dim,spacedim>::
 transform_real_to_unit_cell_initial_guess (const std::vector<Point<spacedim> > &vertex,
                                            const Point<spacedim>               &p) const
 {
-
   Point<dim> p_unit;
 
   FullMatrix<double>  KA(GeometryInfo<dim>::vertices_per_cell, dim);
@@ -1600,44 +1599,43 @@ transform_real_to_unit_cell_initial_guess (const std::vector<Point<spacedim> > &
 
   A_1.vmult(dest,b); //A^{-1}*b
 
-
   for (unsigned int i=0;i<dim;++i)
     p_unit[i]=dest(i);
 
-
   return p_unit;
-
 }
 
+
+
 template<int dim, int spacedim>
 Point<dim>
 MappingQ1<dim,spacedim>::
 transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                              const Point<spacedim>                            &p) const
 {
-
-
                                    // Find the initial value for the
                                    // Newton iteration by a normal projection
                                    // to the least square plane determined by
                                    // the vertices of the cell
   std::vector<Point<spacedim> > a;
   compute_mapping_support_points (cell,a);
-  Point<dim> p_unit =
+  const Point<dim> initial_p_unit =
     transform_real_to_unit_cell_initial_guess(a,p);
 
                                    // if dim==1 there is nothing
-                                   // else to do the initial value is the answer.
-  if (dim>1)
+                                   // else to do to the initial
+                                   // value, and it is the answer
+  if (dim == 1)
+    return initial_p_unit;
+  else
     {
-
                                    // Use the get_data function to
                                    // create an InternalData with data
                                    // vectors of the right size and
                                    // transformation shape values and
                                    // derivatives already computed at
-                                   // point p_unit.
-      const Quadrature<dim> point_quadrature(p_unit);
+                                   // point initial_p_unit.
+      const Quadrature<dim> point_quadrature(initial_p_unit);
 
       UpdateFlags update_flags = update_transformation_values| update_transformation_gradients;
       if (spacedim>dim)
@@ -1657,23 +1655,23 @@ transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_it
                                        // Ignore non vertex support points.
       mdata->mapping_support_points.resize(GeometryInfo<dim>::vertices_per_cell);
 
-                                       // perform the newton iteration.
-      transform_real_to_unit_cell_internal(cell, p, *mdata, p_unit);
+                                       // perform the Newton iteration and
+                                       // return the result
+      return transform_real_to_unit_cell_internal(cell, p, initial_p_unit,
+                                                  *mdata);
     }
-
-  return p_unit;
 }
 
 
 
 template<int dim, int spacedim>
-void
+Point<dim>
 MappingQ1<dim,spacedim>::
 transform_real_to_unit_cell_internal
 (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
  const Point<spacedim>                            &p,
- InternalData                                     &mdata,
- Point<dim>                                       &p_unit) const
+ const Point<dim>                                 &initial_p_unit,
+ InternalData                                     &mdata) const
 {
   const unsigned int n_shapes=mdata.shape_values.size();
   Assert(n_shapes!=0, ExcInternalError());
@@ -1699,6 +1697,8 @@ transform_real_to_unit_cell_internal
 
                                    // f(x)
 
+  Point<dim> p_unit = initial_p_unit;
+
   compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata);
   Point<spacedim> p_real(transform_unit_to_real_cell_internal(mdata));
   Point<spacedim> f = p_real-p;
@@ -1747,7 +1747,6 @@ transform_real_to_unit_cell_internal
                                        // f(x)
       p_real = transform_unit_to_real_cell_internal(mdata);
       f = p_real-p;
-
     }
 
                                    // Here we check that in the last
@@ -1759,6 +1758,8 @@ transform_real_to_unit_cell_internal
                                    // increased and tested, and thus
                                    // havereached the limit.
   AssertThrow(loop<loop_limit, ExcTransformationFailed());
+
+  return p_unit;
 }
 
 
@@ -1780,36 +1781,45 @@ transform_real_to_unit_cell_internal
  */
 
 template<>
-void MappingQ1<2,3>::
+Point<2>
+MappingQ1<2,3>::
 transform_real_to_unit_cell_internal (const Triangulation<2,3>::cell_iterator &cell,
                                       const Point<3> &p,
-                                      InternalData   &mdata,
-                                      Point<2>       &p_unit) const
+                                      const Point<2> &initial_p_unit,
+                                      InternalData   &mdata) const
 {
-  transform_real_to_unit_cell_internal_codim1(cell,p, mdata, p_unit);
+  return
+      transform_real_to_unit_cell_internal_codim1(cell, p, initial_p_unit,
+                                                  mdata);
 }
 
 
 
 
 template<>
-void MappingQ1<1,2>::
+Point<1>
+MappingQ1<1,2>::
 transform_real_to_unit_cell_internal (const Triangulation<1,2>::cell_iterator &cell,
                                       const Point<2> &p,
-                                      InternalData   &mdata,
-                                      Point<1>       &p_unit) const
+                                      const Point<1> &initial_p_unit,
+                                      InternalData   &mdata) const
 {
-  transform_real_to_unit_cell_internal_codim1(cell,p, mdata, p_unit);
+  return
+      transform_real_to_unit_cell_internal_codim1(cell, p, initial_p_unit,
+                                                  mdata);
 }
 
+
 template<>
-void MappingQ1<1,3>::
+Point<1>
+MappingQ1<1,3>::
 transform_real_to_unit_cell_internal (const Triangulation<1,3>::cell_iterator &/*cell*/,
                                       const Point<3> &/*p*/,
-                                      InternalData   &/*mdata*/,
-                                      Point<1>       &/*p_unit*/) const
+                                      const Point<1> &/*initial_p_unit*/,
+                                      InternalData   &/*mdata*/) const
 {
   Assert(false, ExcNotImplemented());
+  return Point<1>();
 }
 
 
@@ -1817,14 +1827,14 @@ transform_real_to_unit_cell_internal (const Triangulation<1,3>::cell_iterator &/
 
 template<int dim, int spacedim>
 template<int dim_>
-void MappingQ1<dim,spacedim>::
+Point<dim_>
+MappingQ1<dim,spacedim>::
 transform_real_to_unit_cell_internal_codim1
 (const typename Triangulation<dim_,dim_ + 1>::cell_iterator &cell,
  const Point<dim_ + 1>                                      &p,
- MappingQ1<dim,spacedim>::InternalData                      &mdata,
- Point<dim_ >                                               &p_unit) const
+ const Point<dim_ >                                         &initial_p_unit,
+ MappingQ1<dim,spacedim>::InternalData                      &mdata) const
 {
-
   const unsigned int spacedim1 = dim_+1;
   const unsigned int dim1 = dim_;
 
@@ -1842,6 +1852,7 @@ transform_real_to_unit_cell_internal_codim1
   Point<spacedim1>  DF[dim1];
   Point<spacedim1>  D2F[dim1][dim1];
 
+  Point<dim1> p_unit = initial_p_unit;
   Point<dim1> f;
   Tensor<2,dim1>  df;
 
@@ -1935,6 +1946,8 @@ transform_real_to_unit_cell_internal_codim1
                                    // increased and tested, and thus
                                    // have reached the limit.
   AssertThrow (loop<loop_limit, ExcTransformationFailed());
+
+  return p_unit;
 }