From 71359c746ff8d082678df71e0719a5c453e487d5 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@math.tamu.edu>
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<dim,spacedim>::
 transform_unit_to_real_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                              const Point<dim> &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<dim> point_quadrature(p);
-  std_cxx11::unique_ptr<InternalData> 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<dim>
+  tensor_pols (Polynomials::generate_complete_Lagrange_basis(line_support_points.get_points()));
+  Assert (tensor_pols.n() == Utilities::fixed_power<dim>(polynomial_degree+1),
+          ExcInternalError());
+
+  // then also construct the mapping from lexicographic to the Qp shape function numbering
+  const std::vector<unsigned int>
+  renumber (FETools::
+            lexicographic_to_hierarchic_numbering (
+              FiniteElementData<dim> (get_dpo_vector<dim>(polynomial_degree), 1,
+                                      polynomial_degree)));
+
+  std::vector<Point<spacedim> > support_points (tensor_pols.n());
+  compute_mapping_support_points(cell, support_points);
+
+  Point<spacedim> mapped_point;
+  for (unsigned int i=0; i<tensor_pols.n(); ++i)
+    mapped_point += support_points[renumber[i]] * tensor_pols.compute_value (i, p);
+
+  return mapped_point;
 }
 
 
-- 
2.39.5