StaticCondensation can also be a DofMap

.gitignore

@@ -239,3 +239,4 @@ installed
 # VSCode related things
 .cache
 .vscode
+compile_commands.json

examples/miscellaneous/miscellaneous_ex16/miscellaneous_ex16.C

@@ -171,7 +171,7 @@ main(int argc, char ** argv)
   libmesh_error_msg_if(!libMesh::relative_fuzzy_equals(*sys.solution, *sc_sys.solution, 1e-4),
                        "mismatching solution");
   libMesh::out << "Static condensation reduced problem size to "
-               << sc_sys.get_static_condensation().get_condensed_mat().m() << std::endl;
+               << sc_sys.get_static_condensation_system_matrix().get_condensed_mat().m() << std::endl;
 
 #if defined(LIBMESH_HAVE_VTK) && !defined(LIBMESH_ENABLE_PARMESH)
 
@@ -208,7 +208,7 @@ assemble_poisson(EquationSystems & es, const std::string & system_name)
   // Get a pointer to the StaticCondensation class if it exists
   StaticCondensation * sc = nullptr;
   if (system.has_static_condensation())
-    sc = &system.get_static_condensation();
+    sc = &system.get_static_condensation_system_matrix();
 
   // A reference to the  DofMap object for this system.  The  DofMap
   // object handles the index translation from node and element numbers

examples/vector_fe/vector_fe_ex9/vector_fe_ex9.C

@@ -52,6 +52,7 @@
 #include "libmesh/equation_systems.h"
 #include "libmesh/nonlinear_implicit_system.h"
 #include "libmesh/nonlinear_solver.h"
+#include "libmesh/static_condensation_dof_map.h"
 #include "libmesh/static_condensation.h"
 
 // The exact solution and error computation.
@@ -151,11 +152,13 @@ main(int argc, char ** argv)
   const FEType scalar_fe_type(FIRST, L2_LAGRANGE);
   const FEType lm_fe_type(FIRST, SIDE_HIERARCHIC);
 
-  StaticCondensation * sc = nullptr;
+  StaticCondensationDofMap * sc_dof_map = nullptr;
+  StaticCondensation * sc_system_matrix = nullptr;
   if (system.has_static_condensation())
     {
-      sc = &system.get_static_condensation();
-      sc->dont_condense_vars({p_num});
+      sc_dof_map = &system.get_static_condensation_dof_map();
+      sc_dof_map->dont_condense_vars({p_num});
+      sc_system_matrix = cast_ptr<StaticCondensation *>(system.matrix);
     }
 
   HDGProblem hdg(nu, cavity);
@@ -170,7 +173,7 @@ main(int argc, char ** argv)
   hdg.scalar_fe_face = FEBase::build(dimension, scalar_fe_type);
   hdg.lm_fe_face = FEBase::build(dimension, lm_fe_type);
   hdg.mms = mms;
-  hdg.sc = sc;
+  hdg.sc = sc_system_matrix;
 
   system.nonlinear_solver->residual_object = &hdg;
   system.nonlinear_solver->jacobian_object = &hdg;

include/Makefile.in

@@ -621,6 +621,7 @@ include_HEADERS = \
         base/auto_ptr.h \
         base/dirichlet_boundaries.h \
         base/dof_map.h \
+        base/dof_map_base.h \
         base/dof_object.h \
         base/factory.h \
         base/float128_shims.h \
@@ -653,7 +654,6 @@ include_HEADERS = \
         enums/enum_fe_family.h \
         enums/enum_inf_map_type.h \
         enums/enum_io_package.h \
-        enums/enum_matrix_build_type.h \
         enums/enum_norm_type.h \
         enums/enum_order.h \
         enums/enum_parallel_type.h \
@@ -888,6 +888,7 @@ include_HEADERS = \
         numerics/sparse_matrix.h \
         numerics/sparse_shell_matrix.h \
         numerics/static_condensation.h \
+        numerics/static_condensation_dof_map.h \
         numerics/static_condensation_preconditioner.h \
         numerics/sum_shell_matrix.h \
         numerics/tensor_shell_matrix.h \

include/base/dof_map.h

@@ -38,6 +38,7 @@
 #include "libmesh/libmesh_logging.h"
 #include "libmesh/enum_elem_type.h"
 #include "libmesh/mesh_subdivision_support.h"
+#include "libmesh/dof_map_base.h"
 
 // TIMPI includes
 #include "timpi/parallel_implementation.h"
@@ -68,7 +69,7 @@ class PeriodicBoundaryBase;
 class PeriodicBoundaries;
 class System;
 class NonlinearImplicitSystem;
-class StaticCondensation;
+class StaticCondensationDofMap;
 template <typename T> class DenseVectorBase;
 template <typename T> class DenseVector;
 template <typename T> class DenseMatrix;
@@ -175,8 +176,8 @@ class NodeConstraints : public std::map<const Node *,
  * \date 2002-2007
  * \brief Manages the degrees of freedom (DOFs) in a simulation.
  */
-class DofMap : public ReferenceCountedObject<DofMap>,
-               public ParallelObject
+class DofMap : public DofMapBase,
+               public ReferenceCountedObject<DofMap>
 {
 public:
 
@@ -594,10 +595,7 @@ class DofMap : public ReferenceCountedObject<DofMap>,
    */
   const VariableGroup & variable_group (const unsigned int c) const;
 
-  /**
-   * \returns The variable description object for variable \p c.
-   */
-  const Variable & variable (const unsigned int c) const;
+  const Variable & variable (const unsigned int c) const override;
 
   /**
    * \returns The approximation order for variable \p c.
@@ -627,12 +625,7 @@ class DofMap : public ReferenceCountedObject<DofMap>,
   unsigned int n_variable_groups() const
   { return cast_int<unsigned int>(_variable_groups.size()); }
 
-  /**
-   * \returns The number of variables in the global solution vector. Defaults
-   * to 1, should be 1 for a scalar equation, 3 for 2D incompressible Navier
-   * Stokes (u,v,p), etc...
-   */
-  unsigned int n_variables() const
+  unsigned int n_variables() const override
   { return cast_int<unsigned int>(_variables.size()); }
 
   /**
@@ -667,11 +660,7 @@ class DofMap : public ReferenceCountedObject<DofMap>,
     return (this->has_blocked_representation() ? this->n_variables() : 1);
   }
 
-  /**
-   * \returns The total number of degrees of freedom in the problem.
-   */
-  dof_id_type n_dofs() const { return _n_dfs; }
-
+  using DofMapBase::n_dofs;
   /**
    * \returns The total number of degrees of freedom for a particular
    * variable \p vn.
@@ -688,12 +677,7 @@ class DofMap : public ReferenceCountedObject<DofMap>,
    */
   dof_id_type n_SCALAR_dofs() const { return _n_SCALAR_dofs; }
 
-  /**
-   * \returns The number of degrees of freedom on this processor.
-   */
-  dof_id_type n_local_dofs () const
-  { return this->n_dofs_on_processor (this->processor_id()); }
-
+  using DofMapBase::n_local_dofs;
   /**
    * \returns The number of degrees of freedom on this processor for a
    * particular variable \p vn.  This is an O(N) operation on serialized or
@@ -706,15 +690,6 @@ class DofMap : public ReferenceCountedObject<DofMap>,
     return n;
   }
 
-  /**
-   * \returns The number of degrees of freedom on partition \p proc.
-   */
-  dof_id_type n_dofs_on_processor(const processor_id_type proc) const
-  {
-    libmesh_assert_less (proc, _first_df.size());
-    return cast_int<dof_id_type>(_end_df[proc] - _first_df[proc]);
-  }
-
   /**
    * \returns The number of degrees of freedom on each partition for a
    * particular variable \p vn.
@@ -726,39 +701,6 @@ class DofMap : public ReferenceCountedObject<DofMap>,
     return n_local_dofs;
   }
 
-  /**
-   * \returns The first dof index that is local to partition \p proc.
-   */
-  dof_id_type first_dof(const processor_id_type proc) const
-  { libmesh_assert_less (proc, _first_df.size()); return _first_df[proc]; }
-
-  dof_id_type first_dof() const
-  { return this->first_dof(this->processor_id()); }
-
-#ifdef LIBMESH_ENABLE_AMR
-  /**
-   * \returns The first old dof index that is local to partition \p proc.
-   */
-  dof_id_type first_old_dof(const processor_id_type proc) const
-  { libmesh_assert_less (proc, _first_old_df.size()); return _first_old_df[proc]; }
-
-  dof_id_type first_old_dof() const
-  { return this->first_old_dof(this->processor_id()); }
-
-#endif //LIBMESH_ENABLE_AMR
-
-  /**
-   * \returns The first dof index that is after all indices local to
-   * processor \p proc.
-   *
-   * Analogous to the end() member function of STL containers.
-   */
-  dof_id_type end_dof(const processor_id_type proc) const
-  { libmesh_assert_less (proc, _end_df.size()); return _end_df[proc]; }
-
-  dof_id_type end_dof() const
-  { return this->end_dof(this->processor_id()); }
-
   /**
    * \returns The processor id that owns the dof index \p dof
    */
@@ -769,25 +711,6 @@ class DofMap : public ReferenceCountedObject<DofMap>,
     return cast_int<processor_id_type>(ub - _end_df.begin());
   }
 
-#ifdef LIBMESH_ENABLE_AMR
-  /**
-   * \returns The first old dof index that is after all indices local
-   * to processor \p proc.
-   *
-   * Analogous to the end() member function of STL containers.
-   */
-  dof_id_type end_old_dof(const processor_id_type proc) const
-  { libmesh_assert_less (proc, _end_old_df.size()); return _end_old_df[proc]; }
-
-  dof_id_type end_old_dof() const
-  { return this->end_old_dof(this->processor_id()); }
-
-#endif //LIBMESH_ENABLE_AMR
-
-  /**
-   * Fills the vector \p di with the global degree of freedom indices
-   * for the element.
-   */
   void dof_indices (const Elem * const elem,
                     std::vector<dof_id_type> & di) const;
 
@@ -799,7 +722,7 @@ class DofMap : public ReferenceCountedObject<DofMap>,
   void dof_indices (const Elem * const elem,
                     std::vector<dof_id_type> & di,
                     const unsigned int vn,
-                    int p_level = -12345) const;
+                    int p_level = -12345) const override;
 
   /**
    * Retrieves degree of freedom indices for a given \p elem and then performs actions for these
@@ -851,7 +774,7 @@ class DofMap : public ReferenceCountedObject<DofMap>,
    */
   void dof_indices (const Node * const node,
                     std::vector<dof_id_type> & di,
-                    const unsigned int vn) const;
+                    const unsigned int vn) const override;
 
   /**
    * Appends to the vector \p di the global degree of freedom indices
@@ -972,6 +895,18 @@ class DofMap : public ReferenceCountedObject<DofMap>,
                               const MeshBase & mesh,
                               unsigned int var_num) const;
 
+  /**
+   * If T == dof_id_type, counts, if T == std::vector<dof_id_type>, fills an
+   * array of, those dof indices which belong to the given variable number and
+   * live on the current processor.
+   */
+  template <typename T,
+            std::enable_if_t<std::is_same_v<T, dof_id_type> ||
+                                 std::is_same_v<T, std::vector<dof_id_type>>,
+                             int> = 0>
+  void local_variable_indices(T & idx, unsigned int var_num) const
+  { this->local_variable_indices(idx, this->_mesh, var_num); }
+
 #ifdef LIBMESH_ENABLE_CONSTRAINTS
 
   //--------------------------------------------------------------------
@@ -1624,11 +1559,6 @@ class DofMap : public ReferenceCountedObject<DofMap>,
                         std::vector<dof_id_type> & di,
                         const unsigned int vn = libMesh::invalid_uint) const;
 
-  /**
-   * \returns The total number of degrees of freedom on old_dof_objects
-   */
-  dof_id_type n_old_dofs() const { return _n_old_dfs; }
-
 #endif // LIBMESH_ENABLE_AMR
 
   /**
@@ -1655,7 +1585,7 @@ class DofMap : public ReferenceCountedObject<DofMap>,
    * Free all new memory associated with the object, but restore its
    * original state, with the mesh pointer and any default ghosting.
    */
-  void clear ();
+  virtual void clear () override;
 
   /**
    * Prints summary info about the sparsity bandwidth and constraints.
@@ -1702,7 +1632,8 @@ class DofMap : public ReferenceCountedObject<DofMap>,
    * constraints or sufficiently many user constraints.
    */
   std::unique_ptr<SparsityPattern::Build> build_sparsity(const MeshBase & mesh,
-                                                         bool calculate_constrained = false) const;
+                                                         bool calculate_constrained = false,
+                                                         bool use_condensed_system = false) const;
 
   /**
    * Describe whether the given variable group should be p-refined. If this API is not called with
@@ -1729,7 +1660,7 @@ class DofMap : public ReferenceCountedObject<DofMap>,
   /**
    * Add a static condensation class
    */
-  void add_static_condensation(const StaticCondensation & sc) { _sc = &sc; }
+  void add_static_condensation(const StaticCondensationDofMap & sc) { _sc = &sc; }
 
   /**
    * Checks whether we have static condensation
@@ -1740,7 +1671,7 @@ class DofMap : public ReferenceCountedObject<DofMap>,
    * @returns the static condensation class. This should have been already added with a call to \p
    * add_static_condensation()
    */
-  const StaticCondensation & get_static_condensation() const;
+  const StaticCondensationDofMap & get_static_condensation() const;
 
 private:
 
@@ -2040,16 +1971,6 @@ class DofMap : public ReferenceCountedObject<DofMap>,
    */
   std::vector<SparseMatrix<Number> * > _matrices;
 
-  /**
-   * First DOF index on processor \p p.
-   */
-  std::vector<dof_id_type> _first_df;
-
-  /**
-   * Last DOF index (plus 1) on processor \p p.
-   */
-  std::vector<dof_id_type> _end_df;
-
   /**
    * First DOF index for SCALAR variable v, or garbage for non-SCALAR
    * variable v
@@ -2152,11 +2073,6 @@ class DofMap : public ReferenceCountedObject<DofMap>,
    */
   std::unique_ptr<SparsityPattern::Build> _sp;
 
-  /**
-   * Total number of degrees of freedom.
-   */
-  dof_id_type _n_dfs;
-
   /**
    * The total number of SCALAR dofs associated to
    * all SCALAR variables.
@@ -2165,21 +2081,6 @@ class DofMap : public ReferenceCountedObject<DofMap>,
 
 #ifdef LIBMESH_ENABLE_AMR
 
-  /**
-   * Total number of degrees of freedom on old dof objects
-   */
-  dof_id_type _n_old_dfs;
-
-  /**
-   * First old DOF index on processor \p p.
-   */
-  std::vector<dof_id_type> _first_old_df;
-
-  /**
-   * Last old DOF index (plus 1) on processor \p p.
-   */
-  std::vector<dof_id_type> _end_old_df;
-
   /**
    * First old DOF index for SCALAR variable v, or garbage for
    * non-SCALAR variable v
@@ -2256,7 +2157,7 @@ class DofMap : public ReferenceCountedObject<DofMap>,
   bool _verify_dirichlet_bc_consistency;
 
   /// Static condensation class
-  const StaticCondensation * _sc;
+  const StaticCondensationDofMap * _sc;
 };
 
 
@@ -2780,7 +2681,7 @@ void DofMap::dof_indices (const Elem * const elem,
 }
 
 inline
-const StaticCondensation & DofMap::get_static_condensation() const
+const StaticCondensationDofMap & DofMap::get_static_condensation() const
 {
   libmesh_assert(_sc);
   return *_sc;