From: Wolfgang Bangerth Date: Thu, 10 Sep 2015 02:57:24 +0000 (-0500) Subject: Avoid the use of an InternalData object in transform_unit_to_real_cell(). X-Git-Tag: v8.4.0-rc2~439^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=71359c746ff8d082678df71e0719a5c453e487d5;p=dealii.git Avoid the use of an InternalData object in transform_unit_to_real_cell(). This avoids a significant overhead in the computation of the forward map for arbitrary points. There is potential for more optimization if we made the TensorProductPolynomials object and the permutation map member variables of the class. But that's for a later patch. Fixes #1569. --- diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 751fec4621..7b3cd001ed 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -791,20 +791,28 @@ MappingQGeneric:: transform_unit_to_real_cell (const typename Triangulation::cell_iterator &cell, const Point &p) const { - //TODO: This function is inefficient -- it might as well evaluate - //the shape functions directly, rather than first building the - //InternalData object and then working with it - - const Quadrature point_quadrature(p); - std_cxx11::unique_ptr mdata (new InternalData(polynomial_degree)); - mdata->initialize (this->requires_update_flags(update_quadrature_points), - point_quadrature, - 1); - - // compute the mapping support - // points - compute_mapping_support_points(cell, mdata->mapping_support_points); - return transform_unit_to_real_cell_internal(*mdata); + // set up the polynomial space + const QGaussLobatto<1> line_support_points (polynomial_degree + 1); + const TensorProductPolynomials + tensor_pols (Polynomials::generate_complete_Lagrange_basis(line_support_points.get_points())); + Assert (tensor_pols.n() == Utilities::fixed_power(polynomial_degree+1), + ExcInternalError()); + + // then also construct the mapping from lexicographic to the Qp shape function numbering + const std::vector + renumber (FETools:: + lexicographic_to_hierarchic_numbering ( + FiniteElementData (get_dpo_vector(polynomial_degree), 1, + polynomial_degree))); + + std::vector > support_points (tensor_pols.n()); + compute_mapping_support_points(cell, support_points); + + Point mapped_point; + for (unsigned int i=0; i