From 6ecb1e40cc1f47f0f3e4c0c759544403dd20994a Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 30 May 2012 21:45:09 +0000 Subject: [PATCH] Clean up code: clearly separate input arguments and return values, rather than using one function argument as both an initial hint and the place where a value is returned. git-svn-id: https://svn.dealii.org/trunk@25580 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/fe/mapping_q1.h | 80 +++++++++++------------- deal.II/source/fe/mapping_q.cc | 12 ++-- deal.II/source/fe/mapping_q1.cc | 83 ++++++++++++++----------- 3 files changed, 91 insertions(+), 84 deletions(-) diff --git a/deal.II/include/deal.II/fe/mapping_q1.h b/deal.II/include/deal.II/fe/mapping_q1.h index befd17bb2d..ce7d36d574 100644 --- 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 * * 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 * computations of the mapping * support points. */ - Point transform_unit_to_real_cell_internal (const InternalData &mdata) const; + Point + 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 * 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::cell_iterator &cell, - const Point &p, - InternalData &mdata, - Point &p_unit) const; - - - - - - - + Point + transform_real_to_unit_cell_internal (const typename Triangulation::cell_iterator &cell, + const Point &p, + const Point &initial_p_unit, + InternalData &mdata) const; /** * Always returns @p true because @@ -581,11 +572,12 @@ class MappingQ1 : public Mapping protected: /* Trick to templatize transform_real_to_unit_cell */ template - void transform_real_to_unit_cell_internal_codim1 + Point + transform_real_to_unit_cell_internal_codim1 (const typename Triangulation::cell_iterator &cell, const Point &p, - InternalData &mdata, - Point &p_unit) const; + const Point &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 }; +// 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 0beb8330d2..1d12cbbd16 100644 --- 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::cell_it { // first a Newton iteration based on a Q1 // mapping. - Point p_unit = + Point initial_p_unit = MappingQ1::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::cell_it // cell in hopes that this gives a // better starting point to the // following iteration - p_unit = GeometryInfo::project_to_unit_cell(p_unit); + initial_p_unit = GeometryInfo::project_to_unit_cell(initial_p_unit); - const Quadrature point_quadrature(p_unit); + const Quadrature 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::cell_it mdata->mapping_support_points.resize(GeometryInfo::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 87485d368c..604be5579f 100644 --- a/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/source/fe/mapping_q1.cc @@ -1566,7 +1566,6 @@ MappingQ1:: transform_real_to_unit_cell_initial_guess (const std::vector > &vertex, const Point &p) const { - Point p_unit; FullMatrix KA(GeometryInfo::vertices_per_cell, dim); @@ -1600,44 +1599,43 @@ transform_real_to_unit_cell_initial_guess (const std::vector > & A_1.vmult(dest,b); //A^{-1}*b - for (unsigned int i=0;i Point MappingQ1:: transform_real_to_unit_cell (const typename Triangulation::cell_iterator &cell, const Point &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 > a; compute_mapping_support_points (cell,a); - Point p_unit = + const Point 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 point_quadrature(p_unit); + // point initial_p_unit. + const Quadrature 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::cell_it // Ignore non vertex support points. mdata->mapping_support_points.resize(GeometryInfo::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 -void +Point MappingQ1:: transform_real_to_unit_cell_internal (const typename Triangulation::cell_iterator &cell, const Point &p, - InternalData &mdata, - Point &p_unit) const + const Point &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 p_unit = initial_p_unit; + compute_shapes(std::vector > (1, p_unit), mdata); Point p_real(transform_unit_to_real_cell_internal(mdata)); Point 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 -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 template -void MappingQ1:: +Point +MappingQ1:: transform_real_to_unit_cell_internal_codim1 (const typename Triangulation::cell_iterator &cell, const Point &p, - MappingQ1::InternalData &mdata, - Point &p_unit) const + const Point &initial_p_unit, + MappingQ1::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 DF[dim1]; Point D2F[dim1][dim1]; + Point p_unit = initial_p_unit; Point 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