From e2bfc7759098bcef8a4b0453fc0c4886ebb5a047 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 15 Sep 2015 08:51:23 -0500 Subject: [PATCH] Move even more functions. This also allows to simplify the logic because we can use function overloading rather than specializations. --- include/deal.II/fe/mapping_q1.h | 31 -- source/fe/mapping_q1.cc | 499 +++++++++++++++----------------- 2 files changed, 238 insertions(+), 292 deletions(-) diff --git a/include/deal.II/fe/mapping_q1.h b/include/deal.II/fe/mapping_q1.h index 33c5b83023..6f5850ec00 100644 --- a/include/deal.II/fe/mapping_q1.h +++ b/include/deal.II/fe/mapping_q1.h @@ -143,37 +143,6 @@ protected: }; -#ifndef DOXYGEN -// explicit specializations - -template<> -Point<2> -MappingQ1<2,3>:: -transform_real_to_unit_cell_internal -(const Triangulation<2,3>::cell_iterator &cell, - const Point<3> &p, - const Point<2> &initial_p_unit, - InternalData &mdata) const; - -template<> -Point<1> -MappingQ1<1,2>:: -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; - -template<> -Point<1> -MappingQ1<1,3>:: -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; - -#endif /** * In order to avoid creation of static MappingQ1 objects at several places in diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc index b41ae6e789..81181053a2 100644 --- a/source/fe/mapping_q1.cc +++ b/source/fe/mapping_q1.cc @@ -509,17 +509,188 @@ namespace } - template - Point - transform_real_to_unit_cell_internal_codim1 - (const typename Triangulation::cell_iterator &cell, - const Point &p, - const Point &initial_p_unit, - typename MappingQ1::InternalData &mdata) + /** + * Implementation of transform_real_to_unit_cell for dim==spacedim + */ + template + Point + do_transform_real_to_unit_cell_internal + (const typename Triangulation::cell_iterator &cell, + const Point &p, + const Point &initial_p_unit, + typename MappingQGeneric::InternalData &mdata) { - const unsigned int spacedim1 = dim_+1; - const unsigned int dim1 = dim_; + const unsigned int spacedim = dim; + + const unsigned int n_shapes=mdata.shape_values.size(); + (void)n_shapes; + Assert(n_shapes!=0, ExcInternalError()); + AssertDimension (mdata.shape_derivatives.size(), n_shapes); + std::vector > &points=mdata.mapping_support_points; + AssertDimension (points.size(), n_shapes); + + + // Newton iteration to solve + // f(x)=p(x)-p=0 + // where we are looking for 'x' and p(x) is the forward transformation + // from unit to real cell. We solve this using a Newton iteration + // x_{n+1}=x_n-[f'(x)]^{-1}f(x) + // The start value is set to be the linear approximation to the cell + + // The shape values and derivatives of the mapping at this point are + // previously computed. + + Point p_unit = initial_p_unit; + + mdata.compute_shape_function_values(std::vector > (1, p_unit)); + + Point p_real = compute_mapped_location_of_point(mdata); + Tensor<1,spacedim> f = p_real-p; + + // early out if we already have our point + if (f.norm_square() < 1e-24 * cell->diameter() * cell->diameter()) + return p_unit; + + // we need to compare the position of the computed p(x) against the given + // point 'p'. We will terminate the iteration and return 'x' if they are + // less than eps apart. The question is how to choose eps -- or, put maybe + // more generally: in which norm we want these 'p' and 'p(x)' to be eps + // apart. + // + // the question is difficult since we may have to deal with very elongated + // cells where we may achieve 1e-12*h for the distance of these two points + // in the 'long' direction, but achieving this tolerance in the 'short' + // direction of the cell may not be possible + // + // what we do instead is then to terminate iterations if + // \| p(x) - p \|_A < eps + // where the A-norm is somehow induced by the transformation of the cell. + // in particular, we want to measure distances relative to the sizes of + // the cell in its principal directions. + // + // to define what exactly A should be, note that to first order we have + // the following (assuming that x* is the solution of the problem, i.e., + // p(x*)=p): + // p(x) - p = p(x) - p(x*) + // = -grad p(x) * (x*-x) + higher order terms + // This suggest to measure with a norm that corresponds to + // A = {[grad p(x]^T [grad p(x)]}^{-1} + // because then + // \| p(x) - p \|_A \approx \| x - x* \| + // Consequently, we will try to enforce that + // \| p(x) - p \|_A = \| f \| <= eps + // + // Note that using this norm is a bit dangerous since the norm changes + // in every iteration (A isn't fixed by depends on xk). However, if the + // cell is not too deformed (it may be stretched, but not twisted) then + // the mapping is almost linear and A is indeed constant or nearly so. + const double eps = 1.e-11; + const unsigned int newton_iteration_limit = 20; + + unsigned int newton_iteration = 0; + double last_f_weighted_norm; + do + { +#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL + std::cout << "Newton iteration " << newton_iteration << std::endl; +#endif + + // f'(x) + Tensor<2,spacedim> df; + for (unsigned int k=0; k &grad_transform=mdata.derivative(0,k); + const Point &point=points[k]; + + for (unsigned int i=0; i delta; + Tensor<2,spacedim> df_inverse = invert(df); + contract (delta, df_inverse, static_cast&>(f)); + +#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL + std::cout << " delta=" << delta << std::endl; +#endif + + // do a line search + double step_length = 1; + do + { + // update of p_unit. The spacedim-th component of transformed point + // is simply ignored in codimension one case. When this component is + // not zero, then we are projecting the point to the surface or + // curve identified by the cell. + Point p_unit_trial = p_unit; + for (unsigned int i=0; i > (1, p_unit_trial)); + + // f(x) + Point p_real_trial = compute_mapped_location_of_point(mdata); + const Tensor<1,spacedim> f_trial = p_real_trial-p; + +#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL + std::cout << " step_length=" << step_length << std::endl + << " ||f || =" << f.norm() << std::endl + << " ||f*|| =" << f_trial.norm() << std::endl + << " ||f*||_A =" << (df_inverse * f_trial).norm() << std::endl; +#endif + + // see if we are making progress with the current step length + // and if not, reduce it by a factor of two and try again + // + // strictly speaking, we should probably use the same norm as we use + // for the outer algorithm. in practice, line search is just a + // crutch to find a "reasonable" step length, and so using the l2 + // norm is probably just fine + if (f_trial.norm() < f.norm()) + { + p_real = p_real_trial; + p_unit = p_unit_trial; + f = f_trial; + break; + } + else if (step_length > 0.05) + step_length /= 2; + else + AssertThrow (false, + (typename Mapping::ExcTransformationFailed())); + } + while (true); + + ++newton_iteration; + if (newton_iteration > newton_iteration_limit) + AssertThrow (false, + (typename Mapping::ExcTransformationFailed())); + last_f_weighted_norm = (df_inverse * f).norm(); + } + while (last_f_weighted_norm > eps); + + return p_unit; + } + + + + /** + * Implementation of transform_real_to_unit_cell for dim==spacedim-1 + */ + template + Point + do_transform_real_to_unit_cell_internal + (const typename Triangulation::cell_iterator &cell, + const Point &p, + const Point &initial_p_unit, + typename MappingQ1::InternalData &mdata) + { + const unsigned int spacedim = dim+1; const unsigned int n_shapes=mdata.shape_values.size(); (void)n_shapes; @@ -527,31 +698,31 @@ namespace Assert(mdata.shape_derivatives.size()==n_shapes, ExcInternalError()); Assert(mdata.shape_second_derivatives.size()==n_shapes, ExcInternalError()); - std::vector > &points=mdata.mapping_support_points; + std::vector > &points=mdata.mapping_support_points; Assert(points.size()==n_shapes, ExcInternalError()); - Point p_minus_F; + Point p_minus_F; - Tensor<1,spacedim1> DF[dim1]; - Tensor<1,spacedim1> D2F[dim1][dim1]; + Tensor<1,spacedim> DF[dim]; + Tensor<1,spacedim> D2F[dim][dim]; - Point p_unit = initial_p_unit; - Point f; - Tensor<2,dim1> df; + Point p_unit = initial_p_unit; + Point f; + Tensor<2,dim> df; - //Evaluate first and second derivatives + // Evaluate first and second derivatives mdata.compute_shape_function_values(std::vector > (1, p_unit)); for (unsigned int k=0; k &grad_phi_k = mdata.derivative(0,k); - const Tensor<2,dim1> &hessian_k = mdata.second_derivative(0,k); - const Point &point_k = points[k]; + const Tensor<1,dim> &grad_phi_k = mdata.derivative(0,k); + const Tensor<2,dim> &hessian_k = mdata.second_derivative(0,k); + const Point &point_k = points[k]; - for (unsigned int j=0; j(mdata); - for (unsigned int j=0; jeps && loop++ d; - Tensor<2,dim1> df_1; + Tensor<1,dim> d; + Tensor<2,dim> df_1; df_1 = invert(df); - contract (d, df_1, static_cast&>(f)); + contract (d, df_1, static_cast&>(f)); p_unit -= d; - for (unsigned int j=0; j &grad_phi_k = mdata.derivative(0,k); - const Tensor<2,dim1> &hessian_k = mdata.second_derivative(0,k); - const Point &point_k = points[k]; + const Tensor<1,dim> &grad_phi_k = mdata.derivative(0,k); + const Tensor<2,dim> &hessian_k = mdata.second_derivative(0,k); + const Point &point_k = points[k]; - for (unsigned int j=0; j(mdata); - for (unsigned int j=0; j::ExcTransformationFailed())); return p_unit; } + + + + + + /** + * Implementation of transform_real_to_unit_cell for other values of + * dim, spacedim + */ + template + Point + do_transform_real_to_unit_cell_internal + (const typename Triangulation::cell_iterator &, + const Point &, + const Point &, + typename MappingQ1::InternalData &) + { + Assert (false, ExcNotImplemented()); + return Point(); + } + } + template Point MappingQ1:: @@ -648,226 +837,14 @@ transform_real_to_unit_cell_internal const Point &initial_p_unit, InternalData &mdata) const { - const unsigned int n_shapes=mdata.shape_values.size(); - (void)n_shapes; - Assert(n_shapes!=0, ExcInternalError()); - AssertDimension (mdata.shape_derivatives.size(), n_shapes); - - std::vector > &points=mdata.mapping_support_points; - AssertDimension (points.size(), n_shapes); - - - // Newton iteration to solve - // f(x)=p(x)-p=0 - // where we are looking for 'x' and p(x) is the forward transformation - // from unit to real cell. We solve this using a Newton iteration - // x_{n+1}=x_n-[f'(x)]^{-1}f(x) - // The start value is set to be the linear approximation to the cell - - // The shape values and derivatives of the mapping at this point are - // previously computed. - - Point p_unit = initial_p_unit; - - mdata.compute_shape_function_values(std::vector > (1, p_unit)); - - Point p_real = compute_mapped_location_of_point(mdata); - Tensor<1,spacedim> f = p_real-p; - - // early out if we already have our point - if (f.norm_square() < 1e-24 * cell->diameter() * cell->diameter()) - return p_unit; - - // we need to compare the position of the computed p(x) against the given - // point 'p'. We will terminate the iteration and return 'x' if they are - // less than eps apart. The question is how to choose eps -- or, put maybe - // more generally: in which norm we want these 'p' and 'p(x)' to be eps - // apart. - // - // the question is difficult since we may have to deal with very elongated - // cells where we may achieve 1e-12*h for the distance of these two points - // in the 'long' direction, but achieving this tolerance in the 'short' - // direction of the cell may not be possible - // - // what we do instead is then to terminate iterations if - // \| p(x) - p \|_A < eps - // where the A-norm is somehow induced by the transformation of the cell. - // in particular, we want to measure distances relative to the sizes of - // the cell in its principal directions. - // - // to define what exactly A should be, note that to first order we have - // the following (assuming that x* is the solution of the problem, i.e., - // p(x*)=p): - // p(x) - p = p(x) - p(x*) - // = -grad p(x) * (x*-x) + higher order terms - // This suggest to measure with a norm that corresponds to - // A = {[grad p(x]^T [grad p(x)]}^{-1} - // because then - // \| p(x) - p \|_A \approx \| x - x* \| - // Consequently, we will try to enforce that - // \| p(x) - p \|_A = \| f \| <= eps - // - // Note that using this norm is a bit dangerous since the norm changes - // in every iteration (A isn't fixed by depends on xk). However, if the - // cell is not too deformed (it may be stretched, but not twisted) then - // the mapping is almost linear and A is indeed constant or nearly so. - const double eps = 1.e-11; - const unsigned int newton_iteration_limit = 20; - - unsigned int newton_iteration = 0; - double last_f_weighted_norm; - do - { -#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL - std::cout << "Newton iteration " << newton_iteration << std::endl; -#endif - - // f'(x) - Tensor<2,spacedim> df; - for (unsigned int k=0; k &grad_transform=mdata.derivative(0,k); - const Point &point=points[k]; - - for (unsigned int i=0; i df_inverse = invert(df); - const Tensor<1, spacedim> delta = - df_inverse * static_cast &>(f); - -#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL - std::cout << " delta=" << delta << std::endl; -#endif - - // do a line search - double step_length = 1; - do - { - // update of p_unit. The spacedim-th component of transformed point - // is simply ignored in codimension one case. When this component is - // not zero, then we are projecting the point to the surface or - // curve identified by the cell. - Point p_unit_trial = p_unit; - for (unsigned int i=0; i > (1, p_unit_trial)); - - // f(x) - Point p_real_trial = compute_mapped_location_of_point(mdata); - const Tensor<1,spacedim> f_trial = p_real_trial-p; - -#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL - std::cout << " step_length=" << step_length << std::endl - << " ||f || =" << f.norm() << std::endl - << " ||f*|| =" << f_trial.norm() << std::endl - << " ||f*||_A =" << (df_inverse * f_trial).norm() << std::endl; -#endif - - // see if we are making progress with the current step length - // and if not, reduce it by a factor of two and try again - // - // strictly speaking, we should probably use the same norm as we use - // for the outer algorithm. in practice, line search is just a - // crutch to find a "reasonable" step length, and so using the l2 - // norm is probably just fine - if (f_trial.norm() < f.norm()) - { - p_real = p_real_trial; - p_unit = p_unit_trial; - f = f_trial; - break; - } - else if (step_length > 0.05) - step_length /= 2; - else - AssertThrow (false, - (typename Mapping::ExcTransformationFailed())); - } - while (true); - - ++newton_iteration; - if (newton_iteration > newton_iteration_limit) - AssertThrow (false, - (typename Mapping::ExcTransformationFailed())); - last_f_weighted_norm = (df_inverse * f).norm(); - } - while (last_f_weighted_norm > eps); - - return p_unit; -} - - - -/* - This function becomes a little tricky in dimension <2,3>. - There is a surface embedded in R^3 and we pass a point p in R^3, that - is most likely not lying on the surface. - We then ask, - what point in R^2 (hopefully in the unit cell) satisfies that - map(x) = p. - - An appropriate modification of this question is: - Find x in R^2 and alpha in R such that - - map(x) + alpha * normal(x) = p - - - */ - -template<> -Point<2> -MappingQ1<2,3>:: -transform_real_to_unit_cell_internal (const Triangulation<2,3>::cell_iterator &cell, - const Point<3> &p, - const Point<2> &initial_p_unit, - InternalData &mdata) const -{ - return - transform_real_to_unit_cell_internal_codim1<2,3>(cell, p, initial_p_unit, - mdata); + // dispatch to the various specializations for spacedim=dim, + // spacedim=dim+1, etc + return do_transform_real_to_unit_cell_internal (cell, p, initial_p_unit, mdata); } -template<> -Point<1> -MappingQ1<1,2>:: -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 -{ - return - transform_real_to_unit_cell_internal_codim1<1,2>(cell, p, initial_p_unit, - mdata); -} - - -template<> -Point<1> -MappingQ1<1,3>:: -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 -{ - Assert(false, ExcNotImplemented()); - return Point<1>(); -} - - - - - - template Mapping * MappingQ1::clone () const -- 2.39.5