From: Luca Heltai Date: Thu, 3 Mar 2016 18:52:56 +0000 (+0100) Subject: Fixing Jacobians. X-Git-Tag: v8.5.0-rc1~1131^2~20 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=acb7d0730c7a1a61fde154bf4283624f18a5a49d;p=dealii.git Fixing Jacobians. --- diff --git a/include/deal.II/fe/mapping_manifold.h b/include/deal.II/fe/mapping_manifold.h index 52371c38e0..4a16e5866f 100644 --- a/include/deal.II/fe/mapping_manifold.h +++ b/include/deal.II/fe/mapping_manifold.h @@ -205,17 +205,30 @@ public: void compute_manifold_quadrature_weights(const Quadrature &quadrature); + /** + * Store vertices internally. + */ + void + store_vertices(const typename Triangulation::cell_iterator &cell) const; + /** * Return an estimate (in bytes) or the memory consumption of this object. */ virtual std::size_t memory_consumption () const; /** - * Store the current cell. + * The current cell vertices. + * + * Computed each. + */ + mutable std::vector > vertices; + + /** + * The current cell. * * Computed each. */ - typename Triangulation::cell_iterator current_cell; + mutable typename Triangulation::cell_iterator cell; /** * The actual quadrature on the reference cell. @@ -331,11 +344,16 @@ public: // */ // mutable typename Triangulation::cell_iterator cell_of_current_support_points; - // /** - // * The determinant of the Jacobian in each quadrature point. Filled if - // * #update_volume_elements. - // */ - // mutable std::vector volume_elements; + /** + * The determinant of the Jacobian in each quadrature point. Filled if + * #update_volume_elements. + */ + mutable std::vector volume_elements; + + /** + * A Q1 Finite element, to compute weights. + */ + const FE_Q fe_q; }; @@ -399,63 +417,6 @@ protected: */ const FE_Q fe_q; - /** - * A table of weights by which we multiply the locations of the support - * points on the perimeter of a quad to get the location of interior support - * points. - * - * Sizes: support_point_weights_on_quad.size()= number of inner - * unit_support_points support_point_weights_on_quad[i].size()= number of - * outer unit_support_points, i.e. unit_support_points on the boundary of - * the quad - * - * For the definition of this vector see equation (8) of the `mapping' - * report. - */ - Table<2,double> support_point_weights_on_quad; - - /** - * A table of weights by which we multiply the locations of the support - * points on the perimeter of a hex to get the location of interior support - * points. - * - * For the definition of this vector see equation (8) of the `mapping' - * report. - */ - Table<2,double> support_point_weights_on_hex; - - /** - * Return the locations of support points for the mapping. For example, for - * $Q_1$ mappings these are the vertices, and for higher order polynomial - * mappings they are the vertices plus interior points on edges, faces, and - * the cell interior that are placed in consultation with the Manifold - * description of the domain and its boundary. However, other classes may - * override this function differently. In particular, the MappingQ1Eulerian - * class does exactly this by not computing the support points from the - * geometry of the current cell but instead evaluating an externally given - * displacement field in addition to the geometry of the cell. - * - * The default implementation of this function is appropriate for most - * cases. It takes the locations of support points on the boundary of the - * cell from the underlying manifold. Interior support points (ie. support - * points in quads for 2d, in hexes for 3d) are then computed using the - * solution of a Laplace equation with the position of the outer support - * points as boundary values, in order to make the transformation as smooth - * as possible. - * - * The function works its way from the vertices (which it takes from the - * given cell) via the support points on the line (for which it calls the - * add_line_support_points() function) and the support points on the quad - * faces (in 3d, for which it calls the add_quad_support_points() function). - * It then adds interior support points that are either computed by - * interpolation from the surrounding points using weights computed by - * solving a Laplace equation, or if dim > - // compute_mapping_support_points (const typename Triangulation::cell_iterator &cell) const; - /** * Transforms the point @p p on the real cell to the corresponding point on * the unit cell @p cell by a Newton iteration. @@ -464,13 +425,6 @@ protected: // transform_real_to_unit_cell_internal (const typename Triangulation::cell_iterator &cell, // const Point &p, // const Point &initial_p_unit) const; - - /** - * Make MappingQ a friend since it needs to call the fill_fe_values() - * functions on its MappingManifold(1) sub-object. - */ - template friend class MappingQ; - }; @@ -481,17 +435,16 @@ protected: #ifndef DOXYGEN -// template -// inline -// const Tensor<1,dim> & -// MappingManifold::InternalData::derivative (const unsigned int qpoint, -// const unsigned int shape_nr) const -// { -// Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(), -// ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, -// shape_derivatives.size())); -// return shape_derivatives [qpoint*n_shape_functions + shape_nr]; -// } +template +inline +void +MappingManifold::InternalData::store_vertices (const typename Triangulation::cell_iterator &cell) const +{ + vertices.resize(GeometryInfo::vertices_per_cell); + for (unsigned int i=0; i::vertices_per_cell; ++i) + vertices[i] = cell->vertex(i); + this->cell = cell; +} diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc index d6216c6779..dfa7449188 100644 --- a/source/fe/mapping_manifold.cc +++ b/source/fe/mapping_manifold.cc @@ -41,7 +41,8 @@ DEAL_II_NAMESPACE_OPEN template -MappingManifold::InternalData::InternalData () +MappingManifold::InternalData::InternalData () : + fe_q(1) {} @@ -81,7 +82,7 @@ initialize (const UpdateFlags update_flags, // see if we need the (transformation) shape function values // and/or gradients and resize the necessary arrays - if (this->update_each & update_quadrature_points) + if (this->update_each & (update_quadrature_points | update_contravariant_transformation) ) compute_manifold_quadrature_weights(q); if (this->update_each & update_covariant_transformation) @@ -178,8 +179,7 @@ initialize_face (const UpdateFlags update_flags, template -MappingManifold::MappingManifold () - : +MappingManifold::MappingManifold () : fe_q(1) // support_point_weights_on_quad (compute_support_point_weights_on_quad(this->polynomial_degree)), // support_point_weights_on_hex (compute_support_point_weights_on_hex(this->polynomial_degree)), @@ -189,8 +189,7 @@ MappingManifold::MappingManifold () template -MappingManifold::MappingManifold (const MappingManifold &mapping) - : +MappingManifold::MappingManifold (const MappingManifold &mapping) : fe_q(1) {} @@ -359,28 +358,23 @@ namespace internal */ template void - maybe_compute_q_points (const typename dealii::Triangulation::cell_iterator &cell, - const typename QProjector::DataSetDescriptor data_set, + maybe_compute_q_points (const typename QProjector::DataSetDescriptor data_set, const typename dealii::MappingManifold::InternalData &data, std::vector > &quadrature_points) { const UpdateFlags update_flags = data.update_each; + AssertDimension(data.vertices.size(), GeometryInfo::vertices_per_cell); if (update_flags & update_quadrature_points) { - AssertDimension(quadrature_points.size(), data.cell_manifold_quadrature_weights.size()); - std::vector > vertices; - for (unsigned int v=0; v::vertices_per_cell; ++v) - vertices.push_back(cell->vertex(v)); - for (unsigned int point=0; pointget_manifold(). - get_new_point(Quadrature(vertices, + quadrature_points[point] = data.cell->get_manifold(). + get_new_point(Quadrature(data.vertices, data.cell_manifold_quadrature_weights[point])); } } @@ -394,8 +388,7 @@ namespace internal */ template void - maybe_update_Jacobians (const typename dealii::Triangulation::cell_iterator &cell, - const typename dealii::QProjector::DataSetDescriptor data_set, + maybe_update_Jacobians (const typename dealii::QProjector::DataSetDescriptor data_set, const typename dealii::MappingManifold::InternalData &data) { const UpdateFlags update_flags = data.update_each; @@ -407,13 +400,22 @@ namespace internal std::fill(data.contravariant.begin(), data.contravariant.end(), DerivativeForm<1,dim,spacedim>()); + // Cache of weights used to compute points on the reference cell + std::vector weights(GeometryInfo::vertices_per_cell); + AssertDimension(GeometryInfo::vertices_per_cell, + data.vertices.size()); for (unsigned int point=0; point &p = data.quad.point(point+data_set); + // And get its image on the manifold: + const Point P = data.cell->get_manifold(). + get_new_point(Quadrature(data.vertices, + data.cell_manifold_quadrature_weights[point+data_set])); + // Always get the maximum length from the point to the // boundary of the reference element, to compute the // tangent vectors from the Manifold object @@ -429,7 +431,17 @@ namespace internal // which is positive, if the coordinate is < .5, double L = ai > .5 ? -ai: 1-ai; - data.contravariant[point][i] = cell->get_manifold().get_tangent_vector(p, np)/L; + // Get the weights to compute the np point in real space + for (unsigned int j=0; j::vertices_per_cell; ++j) + weights[j] = data.fe_q.shape_value(j, np); + + Point NP=data.cell->get_manifold(). + get_new_point(Quadrature(data.vertices, weights)); + + Tensor<1,spacedim> T = data.cell->get_manifold().get_tangent_vector(P, NP); + + for (unsigned int d=0; d::cell_iterator &cell, const unsigned int n_q_points=quadrature.size(); - internal::maybe_compute_q_points (cell, - QProjector::DataSetDescriptor::cell (), + data.store_vertices(cell); + + internal::maybe_compute_q_points (QProjector::DataSetDescriptor::cell (), data, output_data.quadrature_points); -// internal::maybe_update_Jacobians (cell_similarity, -// QProjector::DataSetDescriptor::cell (), -// data); + internal::maybe_update_Jacobians (QProjector::DataSetDescriptor::cell (), + data); // const UpdateFlags update_flags = data.update_each; // const std::vector &weights=quadrature.get_weights();