From 71359c746ff8d082678df71e0719a5c453e487d5 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 9 Sep 2015 21:57:24 -0500 Subject: [PATCH] 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. --- source/fe/mapping_q_generic.cc | 36 +++++++++++++++++++++------------- 1 file changed, 22 insertions(+), 14 deletions(-) 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