From a13893459431d81de0a5b4fc6f49699351fa70e2 Mon Sep 17 00:00:00 2001
From: David Wells <drwells@vt.edu>
Date: Fri, 17 Jul 2015 23:54:21 -0400
Subject: [PATCH] Use a separate function to extract vertices.

---
 include/deal.II/fe/mapping.h             | 11 +++++
 include/deal.II/fe/mapping_q1_eulerian.h | 20 ++++----
 source/fe/mapping.cc                     | 15 ++++++
 source/fe/mapping_q1.cc                  |  5 +-
 source/fe/mapping_q1_eulerian.cc         | 60 +++++++++---------------
 5 files changed, 61 insertions(+), 50 deletions(-)

diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h
index 6fc969415f..144e769594 100644
--- a/include/deal.II/fe/mapping.h
+++ b/include/deal.II/fe/mapping.h
@@ -19,6 +19,7 @@
 
 #include <deal.II/base/config.h>
 #include <deal.II/base/derivative_form.h>
+#include <deal.II/base/std_cxx11/array.h>
 #include <deal.II/base/vector_slice.h>
 #include <deal.II/grid/tria.h>
 #include <deal.II/fe/fe_update_flags.h>
@@ -230,6 +231,16 @@ public:
    */
   virtual ~Mapping ();
 
+  /**
+   * Return the mapped vertices of a cell. These values are not equal to the
+   * vertex coordinates stored by the triangulation for MappingQEulerian and
+   * MappingQ1Eulerian.
+   */
+  virtual
+  std_cxx11::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+  get_vertices (
+    const typename Triangulation<dim,spacedim>::cell_iterator &cell) const;
+
   /**
    * Transforms the point @p p on the unit cell to the point @p p_real on the
    * real cell @p cell and returns @p p_real.
diff --git a/include/deal.II/fe/mapping_q1_eulerian.h b/include/deal.II/fe/mapping_q1_eulerian.h
index 8b31afdc44..06fdb8c7be 100644
--- a/include/deal.II/fe/mapping_q1_eulerian.h
+++ b/include/deal.II/fe/mapping_q1_eulerian.h
@@ -17,6 +17,7 @@
 #define dealii__mapping_q1_eulerian_h
 
 #include <deal.II/base/config.h>
+#include <deal.II/base/std_cxx11/array.h>
 #include <deal.II/base/smartpointer.h>
 #include <deal.II/fe/mapping_q1.h>
 
@@ -99,6 +100,14 @@ public:
   MappingQ1Eulerian (const VECTOR  &euler_transform_vectors,
                      const DoFHandler<dim,spacedim> &shiftmap_dof_handler);
 
+  /**
+   * Return the mapped vertices of the cell.
+   */
+  virtual
+  std_cxx11::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+  get_vertices
+  (const typename Triangulation<dim,spacedim>::cell_iterator &cell) const;
+
   /**
    * Return a pointer to a copy of the present object. The caller of this copy
    * then assumes ownership of it.
@@ -146,17 +155,6 @@ protected:
    * Pointer to the DoFHandler to which the mapping vector is associated.
    */
   SmartPointer<const DoFHandler<dim,spacedim>,MappingQ1Eulerian<dim,VECTOR,spacedim> > shiftmap_dof_handler;
-
-
-private:
-  /**
-   * Computes the support points of the mapping. For @p MappingQ1Eulerian
-   * these are the vertices.
-   */
-  virtual void compute_mapping_support_points(
-    const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-    std::vector<Point<spacedim> > &a) const;
-
 };
 
 /*@}*/
diff --git a/source/fe/mapping.cc b/source/fe/mapping.cc
index af4a2c4359..6160a21dba 100644
--- a/source/fe/mapping.cc
+++ b/source/fe/mapping.cc
@@ -25,6 +25,21 @@ Mapping<dim, spacedim>::~Mapping ()
 {}
 
 
+template <int dim, int spacedim>
+std_cxx11::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+Mapping<dim, spacedim>::get_vertices (
+  const typename Triangulation<dim,spacedim>::cell_iterator &cell) const
+{
+  std_cxx11::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertices;
+  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
+    {
+      vertices[i] = cell->vertex(i);
+    }
+  return vertices;
+}
+
+
+
 /*------------------------------ InternalDataBase ------------------------------*/
 
 
diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc
index 638d46d588..210f01965f 100644
--- a/source/fe/mapping_q1.cc
+++ b/source/fe/mapping_q1.cc
@@ -18,6 +18,7 @@
 #include <deal.II/base/quadrature.h>
 #include <deal.II/base/qprojector.h>
 #include <deal.II/base/memory_consumption.h>
+#include <deal.II/base/std_cxx11/array.h>
 #include <deal.II/base/std_cxx11/unique_ptr.h>
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/grid/tria.h>
@@ -725,10 +726,12 @@ MappingQ1<dim,spacedim>::compute_mapping_support_points(
   const typename Triangulation<dim,spacedim>::cell_iterator &cell,
   std::vector<Point<spacedim> > &a) const
 {
+  std_cxx11::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+  vertices = this->get_vertices(cell);
   a.resize(GeometryInfo<dim>::vertices_per_cell);
 
   for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-    a[i] = cell->vertex(i);
+    a[i] = vertices[i];
 }
 
 
diff --git a/source/fe/mapping_q1_eulerian.cc b/source/fe/mapping_q1_eulerian.cc
index 0903bd44da..77101747d1 100644
--- a/source/fe/mapping_q1_eulerian.cc
+++ b/source/fe/mapping_q1_eulerian.cc
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 2001 - 2014 by the deal.II authors
+// Copyright (C) 2001 - 2015 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -13,6 +13,7 @@
 //
 // ---------------------------------------------------------------------
 
+#include <deal.II/base/std_cxx11/array.h>
 #include <deal.II/fe/mapping_q1_eulerian.h>
 #include <deal.II/lac/vector.h>
 #include <deal.II/lac/petsc_vector.h>
@@ -36,67 +37,50 @@ MappingQ1Eulerian (const EulerVectorType  &euler_transform_vectors,
 
 
 template <int dim, class EulerVectorType, int spacedim>
-void
+std_cxx11::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
 MappingQ1Eulerian<dim, EulerVectorType, spacedim>::
-compute_mapping_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-                               std::vector<Point<spacedim> > &a) const
+get_vertices
+(const typename Triangulation<dim,spacedim>::cell_iterator &cell) const
 {
+  std_cxx11::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertices;
+  // The assertions can not be in the constructor, since this would
+  // require to call dof_handler.distribute_dofs(fe) *before* the mapping
+  // object is constructed, which is not necessarily what we want.
 
-  // The assertions can not be in the
-  // constructor, since this would
-  // require to call
-  // dof_handler.distribute_dofs(fe)
-  // *before* the mapping object is
-  // constructed, which is not
-  // necessarily what we want.
-
-//TODO: Only one of these two assertions should be relevant
+  //TODO: Only one of these two assertions should be relevant
   AssertDimension (spacedim, shiftmap_dof_handler->get_fe().n_dofs_per_vertex());
-  AssertDimension(shiftmap_dof_handler->get_fe().n_components(), spacedim);
+  AssertDimension (shiftmap_dof_handler->get_fe().n_components(), spacedim);
 
   AssertDimension (shiftmap_dof_handler->n_dofs(), euler_transform_vectors->size());
 
-  // cast the
-  // Triangulation<dim>::cell_iterator
-  // into a
-  // DoFHandler<dim>::cell_iterator
-  // which is necessary for access to
+  // cast the Triangulation<dim>::cell_iterator into a
+  // DoFHandler<dim>::cell_iterator which is necessary for access to
   // DoFCellAccessor::get_dof_values()
   typename DoFHandler<dim,spacedim>::cell_iterator dof_cell (*cell, shiftmap_dof_handler);
 
-  // We require the cell to be active
-  // since we can only then get nodal
+  // We require the cell to be active since we can only then get nodal
   // values for the shifts
   Assert (dof_cell->active() == true, ExcInactiveCell());
 
-  // for Q1 elements, the number of
-  // support points should equal the
-  // number of vertices
-  a.resize(GeometryInfo<dim>::vertices_per_cell);
-
-  // now get the values of the shift
-  // vectors at the vertices
+  // now get the values of the shift vectors at the vertices
   Vector<double> mapping_values (shiftmap_dof_handler->get_fe().dofs_per_cell);
   dof_cell->get_dof_values (*euler_transform_vectors, mapping_values);
 
-
   for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
     {
       Point<spacedim> shift_vector;
 
-      // pick out the value of the
-      // shift vector at the present
-      // vertex. since vertex dofs
-      // are always numbered first,
-      // we can access them easily
+      // pick out the value of the shift vector at the present
+      // vertex. since vertex dofs are always numbered first, we can
+      // access them easily
       for (unsigned int j=0; j<spacedim; ++j)
         shift_vector[j] = mapping_values(i*spacedim+j);
 
-      // compute new support point by
-      // old (reference) value and
-      // added shift
-      a[i] = cell->vertex(i) + shift_vector;
+      // compute new support point by old (reference) value and added
+      // shift
+      vertices[i] = cell->vertex(i) + shift_vector;
     }
+  return vertices;
 }
 
 
-- 
2.39.5