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 751fec4621abe5fb05d3ae07fde63b9232f9dc3d..7b3cd001ed241af2ec3bfa52265744d31807d426 100644 (file)
--- 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;
 }