From 0882fc182ab2099f404df28c8b35b4556d61eb20 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 15 Sep 2015 16:01:17 -0500 Subject: [PATCH] Move functions into MappingQGeneric. --- include/deal.II/fe/mapping_q1.h | 20 -- include/deal.II/fe/mapping_q_generic.h | 20 ++ source/fe/mapping_q1.cc | 362 +----------------------- source/fe/mapping_q_generic.cc | 373 ++++++++++++++++++++++++- 4 files changed, 394 insertions(+), 381 deletions(-) diff --git a/include/deal.II/fe/mapping_q1.h b/include/deal.II/fe/mapping_q1.h index 6f5850ec00..8951b5af9a 100644 --- a/include/deal.II/fe/mapping_q1.h +++ b/include/deal.II/fe/mapping_q1.h @@ -105,26 +105,6 @@ protected: */ MappingQ1 (const unsigned int degree); - /** - * Transforms the point @p p on the real cell to the corresponding point on - * the unit cell @p cell by a Newton iteration. - * - * Takes a reference to an @p InternalData that is assumed to be previously - * created by the @p get_data function with @p UpdateFlags including @p - * update_transformation_values and @p update_transformation_gradients and a - * one point Quadrature 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 initial_p_unit. - * - * @p mdata will be changed by this function. - */ - Point - transform_real_to_unit_cell_internal (const typename Triangulation::cell_iterator &cell, - const Point &p, - const Point &initial_p_unit, - InternalData &mdata) const; - /** * Computes the support points of the mapping. For @p MappingQ1 these are * the vertices, as obtained by calling Mapping::get_vertices(). diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h index 46fa2b4b85..0a5838766c 100644 --- a/include/deal.II/fe/mapping_q_generic.h +++ b/include/deal.II/fe/mapping_q_generic.h @@ -535,6 +535,26 @@ protected: compute_mapping_support_points (const typename Triangulation::cell_iterator &cell, std::vector > &a) const; + /** + * Transforms the point @p p on the real cell to the corresponding point on + * the unit cell @p cell by a Newton iteration. + * + * Takes a reference to an @p InternalData that is assumed to be previously + * created by the @p get_data function with @p UpdateFlags including @p + * update_transformation_values and @p update_transformation_gradients and a + * one point Quadrature 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 initial_p_unit. + * + * @p mdata will be changed by this function. + */ + Point + transform_real_to_unit_cell_internal (const typename Triangulation::cell_iterator &cell, + const Point &p, + const Point &initial_p_unit, + InternalData &mdata) const; + /** * For dim=2,3. Append the support points of all shape * functions located on bounding lines of the given cell to the diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc index 81181053a2..67b7ffc9f5 100644 --- a/source/fe/mapping_q1.cc +++ b/source/fe/mapping_q1.cc @@ -478,370 +478,12 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it // return the result. note that this // statement may throw an exception, which // we simply pass up to the caller - return transform_real_to_unit_cell_internal(cell, p, initial_p_unit, - *mdata); + return this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit, + *mdata); } } -namespace -{ - /** - * Using the relative weights of the shape functions evaluated at - * one point on the reference cell (and stored in data.shape_values - * and accessed via data.shape(0,i)) and the locations of mapping - * support points (stored in data.mapping_support_points), compute - * the mapped location of that point in real space. - */ - template - Point - compute_mapped_location_of_point (const typename MappingQ1::InternalData &data) - { - AssertDimension (data.shape_values.size(), - data.mapping_support_points.size()); - - // use now the InternalData to compute the point in real space. - Point p_real; - for (unsigned int i=0; i - 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 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; - Assert(n_shapes!=0, ExcInternalError()); - Assert(mdata.shape_derivatives.size()==n_shapes, ExcInternalError()); - Assert(mdata.shape_second_derivatives.size()==n_shapes, ExcInternalError()); - - std::vector > &points=mdata.mapping_support_points; - Assert(points.size()==n_shapes, ExcInternalError()); - - Point p_minus_F; - - Tensor<1,spacedim> DF[dim]; - Tensor<1,spacedim> D2F[dim][dim]; - - Point p_unit = initial_p_unit; - Point f; - Tensor<2,dim> df; - - // 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,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; jdiameter(); - const unsigned int loop_limit = 10; - - unsigned int loop=0; - - while (f.norm()>eps && loop++ d; - Tensor<2,dim> df_1; - - df_1 = invert(df); - contract (d, df_1, static_cast&>(f)); - p_unit -= d; - - for (unsigned int j=0; j > (1, p_unit)); - - for (unsigned int k=0; k &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:: -transform_real_to_unit_cell_internal -(const typename Triangulation::cell_iterator &cell, - const Point &p, - const Point &initial_p_unit, - InternalData &mdata) const -{ - // 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); -} - diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index adaac27fe7..5429fc3b4f 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -30,7 +30,7 @@ #include #include #include -#include +#include #include #include @@ -1049,6 +1049,377 @@ transform_unit_to_real_cell (const typename Triangulation::cell_it } +// In the code below, GCC tries to instantiate MappingQGeneric<3,4> when +// seeing which of the overloaded versions of +// do_transform_real_to_unit_cell_internal() to call. This leads to bad +// error messages and, generally, nothing very good. Avoid this by ensuring +// that this class exists, but does not have an inner InternalData +// type, thereby ruling out the codim-1 version of the function +// below when doing overload resolution. +template <> +class MappingQGeneric<3,4> +{}; + +namespace +{ + /** + * Using the relative weights of the shape functions evaluated at + * one point on the reference cell (and stored in data.shape_values + * and accessed via data.shape(0,i)) and the locations of mapping + * support points (stored in data.mapping_support_points), compute + * the mapped location of that point in real space. + */ + template + Point + compute_mapped_location_of_point (const typename MappingQGeneric::InternalData &data) + { + AssertDimension (data.shape_values.size(), + data.mapping_support_points.size()); + + // use now the InternalData to compute the point in real space. + Point p_real; + for (unsigned int i=0; i + 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 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 MappingQGeneric::InternalData &mdata) + { + const unsigned int spacedim = dim+1; + + const unsigned int n_shapes=mdata.shape_values.size(); + (void)n_shapes; + Assert(n_shapes!=0, ExcInternalError()); + Assert(mdata.shape_derivatives.size()==n_shapes, ExcInternalError()); + Assert(mdata.shape_second_derivatives.size()==n_shapes, ExcInternalError()); + + std::vector > &points=mdata.mapping_support_points; + Assert(points.size()==n_shapes, ExcInternalError()); + + Point p_minus_F; + + Tensor<1,spacedim> DF[dim]; + Tensor<1,spacedim> D2F[dim][dim]; + + Point p_unit = initial_p_unit; + Point f; + Tensor<2,dim> df; + + // 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,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; jdiameter(); + const unsigned int loop_limit = 10; + + unsigned int loop=0; + + while (f.norm()>eps && loop++ d; + Tensor<2,dim> df_1; + + df_1 = invert(df); + contract (d, df_1, static_cast&>(f)); + p_unit -= d; + + for (unsigned int j=0; j > (1, p_unit)); + + for (unsigned int k=0; k &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 + */ + Point<1> + do_transform_real_to_unit_cell_internal + (const typename Triangulation<1,3>::cell_iterator &, + const Point<3> &, + const Point<1> &, + MappingQGeneric<1,3>::InternalData &) + { + Assert (false, ExcNotImplemented()); + return Point<1>(); + } + +} + + + +template +Point +MappingQGeneric:: +transform_real_to_unit_cell_internal +(const typename Triangulation::cell_iterator &cell, + const Point &p, + const Point &initial_p_unit, + InternalData &mdata) const +{ + // 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 UpdateFlags -- 2.39.5