From: hartmann Date: Wed, 28 Mar 2001 15:03:15 +0000 (+0000) Subject: Fix a bug in the MappingQ::transform_unit_to_real(real_to_unit)_cell functions. This... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c9c30ff407245cd3bcde1a8aa0497668a7cfa162;p=dealii-svn.git Fix a bug in the MappingQ::transform_unit_to_real(real_to_unit)_cell functions. This needed some reorganization that includes giving MappingQ its own implementation of these virtual functions, introducing a new MappingQ1::transform_real_to_unit_internal and removing the MappingQ1::get_cell_data function. git-svn-id: https://svn.dealii.org/trunk@4316 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/mapping_q.h b/deal.II/deal.II/include/fe/mapping_q.h index 7b897388fa..d894a80537 100644 --- a/deal.II/deal.II/include/fe/mapping_q.h +++ b/deal.II/deal.II/include/fe/mapping_q.h @@ -50,6 +50,30 @@ class MappingQ : public MappingQ1 * Destructor. */ virtual ~MappingQ (); + + /** + * Transforms the point @p{p} on + * the unit cell to the point + * @p{p_real} on the real cell + * @p{cell} and returns @p{p_real}. + */ + virtual Point transform_unit_to_real_cell ( + const typename Triangulation::cell_iterator cell, + const Point &p) const; + + /** + * Transforms the point @p{p} on + * the real cell to the point + * @p{p_unit} on the unit cell + * @p{cell} and returns @p{p_unit}. + * + * Uses Newton iteration and the + * @p{transform_unit_to_real_cell} + * function. + */ + virtual Point transform_real_to_unit_cell ( + const typename Triangulation::cell_iterator cell, + const Point &p) const; /** * Implementation of the interface in diff --git a/deal.II/deal.II/include/fe/mapping_q1.h b/deal.II/deal.II/include/fe/mapping_q1.h index 49e59156e2..7d15d8310a 100644 --- a/deal.II/deal.II/include/fe/mapping_q1.h +++ b/deal.II/deal.II/include/fe/mapping_q1.h @@ -388,15 +388,6 @@ class MappingQ1 : public Mapping virtual void compute_shapes_virtual (const std::vector > &unit_points, InternalData &data) const; - /** - * Mapping between tensor product - * ordering and real ordering of - * the vertices. - */ - static const unsigned int vertex_mapping[1< * @p{transform_unit_to_real_cell} * and multiply (through the * Newton iteration) by - * @p{transform_real_to_unit_cell}. + * @p{transform_real_to_unit_cell_internal}. * * Takes a reference to an * @p{InternalData} that must @@ -423,31 +414,55 @@ class MappingQ1 : public Mapping * computations of the mapping * support points. */ - virtual Point transform_unit_to_real_cell_internal (const InternalData &m_data) const; - + Point transform_unit_to_real_cell_internal (const InternalData &mdata) const; + /** - * Returns an @p{InternalData} - * whose data vectors are resized - * corresponding to the - * @p{update_flags} and a - * one-point - * quadrature. Furthermore the - * @p{InternalData} stores the - * mapping support points of the - * given @p{cell}. + * Transforms the point @p{p} on + * the real cell to the point + * @p{p_unit} on the unit cell + * @p{cell} by a Newton + * iteration. * - * This function is called by - * @p{transform_unit_to_real_cell} - * and by - * @p{transform_real_to_unit_cell}. - * The resulting @p{InternalData} - * is given to the - * @p{transform_unit_to_real_internal} - * function. + * 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 + * including the given point + * @p{p_unit}. Hence this + * function assumes that + * @p{mdata} already includes the + * transformation shape values + * and gradients computed at + * @p{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. */ - InternalData* get_cell_data(const typename Triangulation::cell_iterator cell, - const UpdateFlags update_flags) const; + void transform_real_to_unit_cell_internal (const typename Triangulation::cell_iterator cell, + const Point &p, + InternalData &mdata, + Point &p_unit) const; + /** + * Mapping between tensor product + * ordering and real ordering of + * the vertices. + */ + static const unsigned int vertex_mapping[1<::compute_shapes_virtual (const std::vector > &unit_points, MappingQ1::InternalData &data) const { - const unsigned int n_points=unit_points.size(); std::vector values; std::vector > grads; @@ -1208,6 +1207,76 @@ MappingQ::transform_contravariant (std::vector > &dst, } +template +Point MappingQ::transform_unit_to_real_cell ( + const typename Triangulation::cell_iterator cell, + const Point &p) const +{ + // Use the get_data function to + // create an InternalData with data + // vectors of the right size and + // transformation shape values + // already computed at point p. + const Quadrature point_quadrature(p); + InternalData *mdata=dynamic_cast ( + get_data(update_transformation_values, point_quadrature)); + Assert(mdata!=0, ExcInternalError()); + + mdata->use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells + || cell->has_boundary_lines()); + + MappingQ1::InternalData *p_data=0; + if (mdata->use_mapping_q1_on_current_cell) + p_data=&mdata->mapping_q1_data; + else + p_data=mdata; + + compute_mapping_support_points(cell, p_data->mapping_support_points); + + return transform_unit_to_real_cell_internal(*p_data); +} + +template +Point MappingQ::transform_real_to_unit_cell ( + const typename Triangulation::cell_iterator cell, + const Point &p) const +{ + // Let the start value of the + // newton iteration be the center + // of the unit cell + Point p_unit; + for (unsigned int i=0; i point_quadrature(p_unit); + InternalData *mdata=dynamic_cast ( + get_data(update_transformation_values | update_transformation_gradients, + point_quadrature)); + Assert(mdata!=0, ExcInternalError()); + + mdata->use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells + || cell->has_boundary_lines()); + + MappingQ1::InternalData *p_data=0; + if (mdata->use_mapping_q1_on_current_cell) + p_data=&mdata->mapping_q1_data; + else + p_data=mdata; + + // perform the newton iteration. + transform_real_to_unit_cell_internal(cell, p, *p_data, p_unit); + + delete mdata; + return p_unit; +} + + // explicit instantiation template class MappingQ; diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc index c4c3118356..85730cf02e 100644 --- a/deal.II/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/deal.II/source/fe/mapping_q1.cc @@ -845,15 +845,19 @@ Point MappingQ1::transform_unit_to_real_cell ( const typename Triangulation::cell_iterator cell, const Point &p) const { - // Use the get_cell_data function - // to create an InternalData with - // data vectors already of the - // right size and with mapping - // support points set. - InternalData *mdata=get_cell_data(cell, update_transformation_values); + // Use the get_data function to + // create an InternalData with data + // vectors of the right size and + // transformation shape values + // already computed at point p. + const Quadrature point_quadrature(p); + InternalData *mdata=dynamic_cast ( + get_data(update_transformation_values, point_quadrature)); Assert(mdata!=0, ExcInternalError()); - compute_shapes(std::vector > (1, p), *mdata); + // compute the mapping support + // points + compute_mapping_support_points(cell, mdata->mapping_support_points); return transform_unit_to_real_cell_internal(*mdata); } @@ -876,40 +880,68 @@ Point MappingQ1::transform_unit_to_real_cell_internal ( } - template Point MappingQ1::transform_real_to_unit_cell ( const typename Triangulation::cell_iterator cell, const Point &p) const { - // Use the get_cell_data function - // to create an InternalData with - // data vectors already of the - // right size and with mapping - // support points set. - InternalData *mdata=get_cell_data(cell, - update_transformation_values - | update_transformation_gradients); + // Let the start value of the + // newton iteration be the center + // of the unit cell + Point p_unit; + for (unsigned int i=0; i point_quadrature(p_unit); + InternalData *mdata=dynamic_cast ( + get_data(update_transformation_values | update_transformation_gradients, + point_quadrature)); Assert(mdata!=0, ExcInternalError()); - std::vector > &points=mdata->mapping_support_points; + + // perform the newton iteration. + transform_real_to_unit_cell_internal(cell, p, *mdata, p_unit); + + delete mdata; + return p_unit; +} + + + +template +void MappingQ1::transform_real_to_unit_cell_internal ( + const typename Triangulation::cell_iterator cell, + const Point &p, + InternalData &mdata, + Point &p_unit) const +{ + const unsigned int n_shapes=mdata.shape_values.size(); + Assert(n_shapes!=0, ExcInternalError()); + Assert(mdata.shape_derivatives.size()==n_shapes, ExcInternalError()); + + std::vector > &points=mdata.mapping_support_points; + compute_mapping_support_points(cell, points); + Assert(points.size()==n_shapes, ExcInternalError()); // Newton iteration to solve // f(x)=p(x)-p=0 // x_{n+1}=x_n-[f'(x)]^{-1}f(x) - // Let the start value be the - // center of the unit cell - // (p_unit stands for x) - Point p_unit; - for (unsigned int i=0; i > (1, p_unit), *mdata); + // The shape values and derivatives + // of the mapping at this point are + // previously computed. + // f(x) - Point p_real(transform_unit_to_real_cell_internal(*mdata)); + Point p_real(transform_unit_to_real_cell_internal(mdata)); Point f = p_real-p; const double eps=1e-15*cell->diameter(); @@ -918,9 +950,9 @@ Point MappingQ1::transform_real_to_unit_cell ( { // f'(x) Tensor<2,dim> df; - for (unsigned int k=0; kn_shape_functions; ++k) + for (unsigned int k=0; k &grad_transform=mdata->derivative(0,k); + const Tensor<1,dim> &grad_transform=mdata.derivative(0,k); const Point &point=points[k]; for (unsigned int i=0; i MappingQ1::transform_real_to_unit_cell ( p_unit -= d; // shape values and derivatives // at new p_unit point - compute_shapes(std::vector > (1, p_unit), *mdata); + compute_shapes(std::vector > (1, p_unit), mdata); // f(x) - p_real=transform_unit_to_real_cell_internal(*mdata); + p_real=transform_unit_to_real_cell_internal(mdata); f = p_real-p; } - - delete mdata; - return p_unit; -} - - -template -MappingQ1::InternalData* -MappingQ1::get_cell_data(const typename Triangulation::cell_iterator cell, - const UpdateFlags update_flags) const -{ - static const Point dummy_p; - static const Quadrature dummy_quadrature(dummy_p); - - InternalData *mdata=dynamic_cast ( - get_data(update_flags, dummy_quadrature)); - Assert(mdata!=0, ExcInternalError()); - - // compute the mapping support - // points - std::vector > &points=mdata->mapping_support_points; - compute_mapping_support_points(cell, points); - - return mdata; } - - template