Check orientation when caching for required families

Author: nmnobre

include/fe/fe.h

@@ -790,10 +790,23 @@ class FE : public FEGenericBase<typename FEOutputType<T>::type>
                                       const unsigned vdim = 1);
 
   /**
-   * An array of the node locations on the last
-   * element we computed on
+   * Vectors holding the node locations, edge and
+   * face orientations of the last element we cached.
    */
   std::vector<Point> cached_nodes;
+  std::vector<bool> cached_edges, cached_faces;
+
+  /**
+   * Repopulate the element cache with the node locations,
+   * edge and face orientations of the element \p elem.
+   */
+  void cache(const Elem * elem);
+
+  /**
+   * Check if the node locations, edge and face orientations
+   * held in the element cache match those of element \p elem.
+   */
+  bool match(const Elem * elem);
 
   /**
    * The last side and last edge we did a reinit on

include/fe/fe_interface.h

@@ -845,6 +845,11 @@ class FEInterface
    */
   static FEFieldType field_type (const FEFamily & fe_family);
 
+  /**
+   * \returns Whether the element's shape functions are orientation-dependent.
+   */
+  static bool orientation_dependent (const FEFamily & fe_family);
+
   /**
    * \returns The number of components of a vector-valued element.
    * Scalar-valued elements return 1.

src/fe/fe.C

@@ -42,8 +42,9 @@ namespace {
 namespace libMesh
 {
 // ------------------------------------------------------------
-// Whether we cache the node locations on the last element we computed on
-// to try to avoid calling init_shape_functions and compute_shape_functions
+// Whether we cache the node locations, edge and face orientations of the last
+// element we computed on as needed to avoid calling init_shape_functions and
+// compute_shape_functions
 static const bool * caching = nullptr;
 
 class CachingSetup: public Singleton::Setup
@@ -148,6 +149,50 @@ void FE<Dim,T>::dofs_on_edge(const Elem * const elem,
 
 
 
+template <unsigned int Dim, FEFamily T>
+void FE<Dim,T>::cache(const Elem * elem)
+{
+  cached_nodes.resize(elem->n_nodes());
+  for (auto n : elem->node_index_range())
+    cached_nodes[n] = elem->point(n);
+
+  if (FEInterface::orientation_dependent(T))
+    {
+      cached_edges.resize(elem->n_edges());
+      for (auto n : elem->edge_index_range())
+        cached_edges[n] = elem->positive_edge_orientation(n);
+
+      cached_faces.resize(elem->n_faces());
+      for (auto n : elem->face_index_range())
+        cached_faces[n] = elem->positive_face_orientation(n);
+    }
+}
+
+
+
+template <unsigned int Dim, FEFamily T>
+bool FE<Dim,T>::match(const Elem * elem)
+{
+  bool m = cached_nodes.size() == elem->n_nodes();
+  for (unsigned n = 1; m && n < elem->n_nodes(); n++)
+    m = (elem->point(n) - elem->point(0)).relative_fuzzy_equals(cached_nodes[n] - cached_nodes[0]);
+
+  if (FEInterface::orientation_dependent(T))
+    {
+      m &= cached_edges.size() == elem->n_edges();
+      for (unsigned n = 0; m && n < elem->n_edges(); n++)
+        m = elem->positive_edge_orientation(n) == cached_edges[n];
+
+      m &= cached_faces.size() == elem->n_faces();
+      for (unsigned n = 0; m && n < elem->n_faces(); n++)
+        m = elem->positive_face_orientation(n) == cached_faces[n];
+    }
+
+  return m;
+}
+
+
+
 template <unsigned int Dim, FEFamily T>
 void FE<Dim,T>::reinit(const Elem * elem,
                        const std::vector<Point> * const pts,
@@ -162,7 +207,7 @@ void FE<Dim,T>::reinit(const Elem * elem,
 
   // Try to avoid calling init_shape_functions
   // even when shapes_need_reinit
-  bool cached_nodes_still_fit = false;
+  bool cached_elem_still_fits = false;
 
   // Most of the hard work happens when we have an actual element
   if (elem)
@@ -215,6 +260,7 @@ void FE<Dim,T>::reinit(const Elem * elem,
               this->_elem_type = elem->type();
               this->_elem_p_level = elem->p_level();
               this->_p_level = this->_add_p_level_in_reinit * elem->p_level();
+
               // Initialize the shape functions
               this->_fe_map->template init_reference_to_physical_map<Dim>
                 (this->qrule->get_points(), elem);
@@ -223,27 +269,22 @@ void FE<Dim,T>::reinit(const Elem * elem,
           else
             {
               this->_elem = elem;
-              // Check if cached elem geometry is sufficiently similar
-              cached_nodes_still_fit = cached_nodes.size() == elem->n_nodes();
-              for (unsigned n = 1; cached_nodes_still_fit && n < elem->n_nodes(); n++)
-                cached_nodes_still_fit = (elem->point(n) - elem->point(0)).relative_fuzzy_equals
-                                         (cached_nodes[n] - cached_nodes[0]);
+
+              // Check if cached element's nodes, edge and face orientations still fit
+              cached_elem_still_fits = this->match(elem);
+
               // Initialize the shape functions if needed
-              if (this->shapes_need_reinit() && !cached_nodes_still_fit)
+              if (this->shapes_need_reinit() && !cached_elem_still_fits)
                 {
                   this->_fe_map->template init_reference_to_physical_map<Dim>
                     (this->qrule->get_points(), elem);
                   this->init_shape_functions (this->qrule->get_points(), elem);
                 }
             }
 
-          // Replace cached nodes if no longer fitting
-          if (this->shapes_need_reinit() && !cached_nodes_still_fit && *caching)
-            {
-              cached_nodes.resize(elem->n_nodes());
-              for (auto n : elem->node_index_range())
-                cached_nodes[n] = elem->point(n);
-            }
+          // Replace cached nodes, edge and face orientations if no longer fitting
+          if (this->shapes_need_reinit() && !cached_elem_still_fits && *caching)
+            this->cache(elem);
 
           // The shape functions correspond to the qrule
           this->shapes_on_quadrature = true;
@@ -300,7 +341,7 @@ void FE<Dim,T>::reinit(const Elem * elem,
 
   // Compute the shape functions and the derivatives at all of the
   // quadrature points.
-  if (!cached_nodes_still_fit)
+  if (!cached_elem_still_fits)
     {
       if (pts != nullptr)
         this->compute_shape_functions (elem,*pts);

src/fe/fe_interface.C

@@ -2655,6 +2655,25 @@ FEFieldType FEInterface::field_type (const FEFamily & fe_family)
     }
 }
 
+bool FEInterface::orientation_dependent (const FEFamily & fe_family)
+{
+  switch (fe_family)
+    {
+    case HIERARCHIC:
+    case L2_HIERARCHIC:
+    case HIERARCHIC_VEC:
+    case L2_HIERARCHIC_VEC:
+    case BERNSTEIN:
+    case RATIONAL_BERNSTEIN:
+    case SZABAB:
+    case NEDELEC_ONE:
+    case RAVIART_THOMAS:
+      return true;
+    default:
+      return false;
+    }
+}
+
 unsigned int FEInterface::n_vec_dim (const MeshBase & mesh,
                                      const FEType & fe_type)
 {