From c35f04452455a8002b55f3b9787b676d96e8f167 Mon Sep 17 00:00:00 2001
From: Luca Heltai <luca.heltai@sissa.it>
Date: Tue, 24 Apr 2018 13:27:56 +0200
Subject: [PATCH] Added implementation of MappingFEField::get_vertices

---
 include/deal.II/fe/mapping_fe_field.h        | 28 ++++++--
 source/fe/mapping_fe_field.cc                | 71 ++++++++++++++++++--
 tests/grid/grid_tools_min_max_diameter_02.cc |  3 +
 3 files changed, 90 insertions(+), 12 deletions(-)

diff --git a/include/deal.II/fe/mapping_fe_field.h b/include/deal.II/fe/mapping_fe_field.h
index c7235ef21e..f27680e7d6 100644
--- a/include/deal.II/fe/mapping_fe_field.h
+++ b/include/deal.II/fe/mapping_fe_field.h
@@ -24,6 +24,7 @@
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/fe/mapping.h>
 #include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_values.h>
 
 #include <array>
 
@@ -125,8 +126,6 @@ public:
   virtual
   std::unique_ptr<Mapping<dim,spacedim>> clone () const override;
 
-
-
   /**
    * See the documentation of Mapping::preserves_vertex_locations()
    * for the purpose of this function. The implementation in this
@@ -135,6 +134,16 @@ public:
   virtual
   bool preserves_vertex_locations () const override;
 
+  /**
+   * Return the mapped vertices of a cell.
+   *
+   * This mapping ignores the vertices of the Triangulation it is associated to,
+   * and constructs the position of the vertices according to the @p euler_vector
+   * that was passed at construction time.
+   */
+  virtual
+  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+  get_vertices (const typename Triangulation<dim,spacedim>::cell_iterator &cell) const override;
 
   /**
    * @name Mapping points between reference and real cells
@@ -206,7 +215,6 @@ public:
    * @}
    */
 
-
   /**
    * Return the degree of the mapping, i.e. the value which was passed to the
    * constructor.
@@ -574,6 +582,17 @@ private:
    */
   std::vector<unsigned int> fe_to_real;
 
+  /**
+   * FEValues object used to query the given finite element field at the
+   * support points in the reference configuration.
+   */
+  mutable FEValues<dim,spacedim> fe_values;
+
+  /**
+   * A variable to guard access to the fe_values variable.
+   */
+  mutable Threads::Mutex fe_values_mutex;
+
   void
   compute_data (const UpdateFlags      update_flags,
                 const Quadrature<dim>  &q,
@@ -599,9 +618,6 @@ private:
 
 #ifndef DOXYGEN
 
-
-
-
 #endif // DOXYGEN
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/source/fe/mapping_fe_field.cc b/source/fe/mapping_fe_field.cc
index 0d339c5fad..3009f8cc8f 100644
--- a/source/fe/mapping_fe_field.cc
+++ b/source/fe/mapping_fe_field.cc
@@ -202,6 +202,27 @@ MappingFEField<dim,spacedim,DoFHandlerType,VectorType>::InternalData::fourth_der
 
 
 
+namespace
+{
+  // Construct a quadrature formula containing the vertices of the reference
+  // cell in dimension dim (with invalid weights)
+  template<int dim>
+  Quadrature<dim> &get_vertex_quadrature()
+  {
+    static Quadrature<dim> quad;
+    if (quad.size() == 0)
+      {
+        std::vector<Point<dim> > points(GeometryInfo<dim>::vertices_per_cell);
+        for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+          points[i] = GeometryInfo<dim>::unit_cell_vertex(i);
+        quad = Quadrature<dim>(points);
+      }
+    return quad;
+  }
+}
+
+
+
 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
 MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::MappingFEField
 (const DoFHandlerType            &euler_dof_handler,
@@ -212,7 +233,8 @@ MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::MappingFEField
   euler_dof_handler(&euler_dof_handler),
   fe_mask(mask.size() ? mask :
           ComponentMask(euler_dof_handler.get_fe().get_nonzero_components(0).size(), true)),
-  fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
+  fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int),
+  fe_values(this->euler_dof_handler->get_fe(), get_vertex_quadrature<dim>(), update_values)
 {
   unsigned int size = 0;
   for (unsigned int i=0; i<fe_mask.size(); ++i)
@@ -231,15 +253,16 @@ MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::MappingFEField
   euler_vector(mapping.euler_vector),
   euler_dof_handler(mapping.euler_dof_handler),
   fe_mask(mapping.fe_mask),
-  fe_to_real(mapping.fe_to_real)
+  fe_to_real(mapping.fe_to_real),
+  fe_values(mapping.euler_dof_handler->get_fe(), get_vertex_quadrature<dim>(), update_values)
 {}
 
 
 
-template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
+template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
 inline
 const double &
-MappingFEField<dim,spacedim,DoFHandlerType,VectorType>::InternalData::shape
+MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::InternalData::shape
 (const unsigned int qpoint,
  const unsigned int shape_nr) const
 {
@@ -251,15 +274,51 @@ MappingFEField<dim,spacedim,DoFHandlerType,VectorType>::InternalData::shape
 
 
 
-template <int dim, int spacedim, typename DoFHandlerType, typename VectorType>
+template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
 bool
-MappingFEField<dim,spacedim,DoFHandlerType,VectorType>::preserves_vertex_locations () const
+MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::preserves_vertex_locations () const
 {
   return false;
 }
 
 
 
+template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
+std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::
+get_vertices
+(const typename Triangulation<dim,spacedim>::cell_iterator &cell) const
+{
+  // we transform our tria iterator into a dof iterator so we can access
+  // data not associated with triangulations
+  const typename DoFHandler<dim,spacedim>::cell_iterator dof_cell(*cell,
+      euler_dof_handler);
+
+  Assert (dof_cell->active() == true, ExcInactiveCell());
+  AssertDimension(GeometryInfo<dim>::vertices_per_cell, fe_values.n_quadrature_points);
+  AssertDimension(fe_to_real.size(), euler_dof_handler->get_fe().n_components());
+
+  std::vector<Vector<typename VectorType::value_type> >
+  values(fe_values.n_quadrature_points,
+         Vector<typename VectorType::value_type> (euler_dof_handler->get_fe().n_components()));
+
+  {
+    Threads::Mutex::ScopedLock lock(fe_values_mutex);
+    fe_values.reinit(dof_cell);
+    fe_values.get_function_values(*euler_vector, values);
+  }
+  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertices;
+
+  for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+    for (unsigned int j=0; j<fe_to_real.size(); ++j)
+      if (fe_to_real[j] != numbers::invalid_unsigned_int)
+        vertices[i][fe_to_real[j]] = values[i][j];
+
+  return vertices;
+}
+
+
+
 template <int dim, int spacedim, typename VectorType, typename DoFHandlerType>
 void
 MappingFEField<dim,spacedim,VectorType,DoFHandlerType>::
diff --git a/tests/grid/grid_tools_min_max_diameter_02.cc b/tests/grid/grid_tools_min_max_diameter_02.cc
index 8dd263bad8..f653eae67b 100644
--- a/tests/grid/grid_tools_min_max_diameter_02.cc
+++ b/tests/grid/grid_tools_min_max_diameter_02.cc
@@ -13,6 +13,9 @@
 //
 // ---------------------------------------------------------------------
 
+// test GridTools::minimal_cell_diameter and GridTools::maximal_cell_diameter with a
+// MappingFEField
+
 #include "../tests.h"
 #include <deal.II/base/quadrature_lib.h>
 #include <deal.II/lac/vector.h>
-- 
2.39.5