Merge pull request #3386 from roystgnr/code_dedup2

More code deduplication

Author: roystgnr

include/base/dof_map.h

@@ -1584,6 +1584,18 @@ class DofMap : public ReferenceCountedObject<DofMap>,
   void distribute_local_dofs_node_major (dof_id_type & next_free_dof,
                                          MeshBase & mesh);
 
+  /*
+   * Helper method for the above two to count + distriubte SCALAR dofs
+   */
+  void distribute_scalar_dofs (dof_id_type & next_free_dof);
+
+#ifdef DEBUG
+  /*
+   * Internal assertions for distribute_local_dofs_*
+   */
+  void assert_no_nodes_missed(MeshBase & mesh);
+#endif
+
   /*
    * A utility method for obtaining a set of elements to ghost along
    * with merged coupling matrices.

include/mesh/mesh_function.h

@@ -353,6 +353,14 @@ class MeshFunction : public FunctionBase<Number>,
   std::set<const Elem *> find_elements(const Point & p,
                                        const std::set<subdomain_id_type> * subdomain_ids = nullptr) const;
 
+  /**
+   * Helper function for finding a gradient as evaluated from a
+   * specific element
+   */
+  void _gradient_on_elem (const Point & p,
+                          const Elem * element,
+                          std::vector<Gradient> & output);
+
   /**
    * The equation systems handler, from which
    * the data are gathered.

include/numerics/dense_matrix.h

@@ -767,6 +767,20 @@ class DenseMatrix : public DenseMatrixBase<T>
                     DenseVector<T> & dest,
                     const DenseVector<T> & arg,
                     bool trans=false) const;
+
+  /**
+   * Left multiplies by the transpose of the matrix \p A which
+   * may contain a different numerical type.
+   */
+  template <typename T2>
+  void _left_multiply_transpose (const DenseMatrix<T2> & A);
+
+  /**
+   * Right multiplies by the transpose of the matrix \p A which
+   * may contain a different numerical type.
+   */
+  template <typename T2>
+  void _right_multiply_transpose (const DenseMatrix<T2> & A);
 };
 
 

include/numerics/dense_matrix_impl.h

@@ -95,68 +95,21 @@ void DenseMatrix<T>::left_multiply_transpose(const DenseMatrix<T> & A)
   if (this->use_blas_lapack)
     this->_multiply_blas(A, LEFT_MULTIPLY_TRANSPOSE);
   else
-    {
-      //Check to see if we are doing (A^T)*A
-      if (this == &A)
-        {
-          //libmesh_here();
-          DenseMatrix<T> B(*this);
-
-          // Simple but inefficient way
-          // return this->left_multiply_transpose(B);
-
-          // More efficient, but more code way
-          // If A is mxn, the result will be a square matrix of Size n x n.
-          const unsigned int n_rows = A.m();
-          const unsigned int n_cols = A.n();
-
-          // resize() *this and also zero out all entries.
-          this->resize(n_cols,n_cols);
-
-          // Compute the lower-triangular part
-          for (unsigned int i=0; i<n_cols; ++i)
-            for (unsigned int j=0; j<=i; ++j)
-              for (unsigned int k=0; k<n_rows; ++k) // inner products are over n_rows
-                (*this)(i,j) += B(k,i)*B(k,j);
-
-          // Copy lower-triangular part into upper-triangular part
-          for (unsigned int i=0; i<n_cols; ++i)
-            for (unsigned int j=i+1; j<n_cols; ++j)
-              (*this)(i,j) = (*this)(j,i);
-        }
+    this->_left_multiply_transpose(A);
+}
 
-      else
-        {
-          DenseMatrix<T> B(*this);
-
-          this->resize (A.n(), B.n());
-
-          libmesh_assert_equal_to (A.m(), B.m());
-          libmesh_assert_equal_to (this->m(), A.n());
-          libmesh_assert_equal_to (this->n(), B.n());
-
-          const unsigned int m_s = A.n();
-          const unsigned int p_s = A.m();
-          const unsigned int n_s = this->n();
-
-          // Do it this way because there is a
-          // decent chance (at least for constraint matrices)
-          // that A.transpose(i,k) = 0.
-          for (unsigned int i=0; i<m_s; i++)
-            for (unsigned int k=0; k<p_s; k++)
-              if (A.transpose(i,k) != 0.)
-                for (unsigned int j=0; j<n_s; j++)
-                  (*this)(i,j) += A.transpose(i,k)*B(k,j);
-        }
-    }
 
+template<typename T>
+template<typename T2>
+void DenseMatrix<T>::left_multiply_transpose(const DenseMatrix<T2> & A)
+{
+  this->_left_multiply_transpose(A);
 }
 
 
-
 template<typename T>
 template<typename T2>
-void DenseMatrix<T>::left_multiply_transpose(const DenseMatrix<T2> & A)
+void DenseMatrix<T>::_left_multiply_transpose(const DenseMatrix<T2> & A)
 {
   //Check to see if we are doing (A^T)*A
   if (this == &A)
@@ -267,67 +220,23 @@ void DenseMatrix<T>::right_multiply_transpose (const DenseMatrix<T> & B)
   if (this->use_blas_lapack)
     this->_multiply_blas(B, RIGHT_MULTIPLY_TRANSPOSE);
   else
-    {
-      //Check to see if we are doing B*(B^T)
-      if (this == &B)
-        {
-          //libmesh_here();
-          DenseMatrix<T> A(*this);
-
-          // Simple but inefficient way
-          // return this->right_multiply_transpose(A);
-
-          // More efficient, more code
-          // If B is mxn, the result will be a square matrix of Size m x m.
-          const unsigned int n_rows = B.m();
-          const unsigned int n_cols = B.n();
-
-          // resize() *this and also zero out all entries.
-          this->resize(n_rows,n_rows);
-
-          // Compute the lower-triangular part
-          for (unsigned int i=0; i<n_rows; ++i)
-            for (unsigned int j=0; j<=i; ++j)
-              for (unsigned int k=0; k<n_cols; ++k) // inner products are over n_cols
-                (*this)(i,j) += A(i,k)*A(j,k);
-
-          // Copy lower-triangular part into upper-triangular part
-          for (unsigned int i=0; i<n_rows; ++i)
-            for (unsigned int j=i+1; j<n_rows; ++j)
-              (*this)(i,j) = (*this)(j,i);
-        }
-
-      else
-        {
-          DenseMatrix<T> A(*this);
-
-          this->resize (A.m(), B.m());
-
-          libmesh_assert_equal_to (A.n(), B.n());
-          libmesh_assert_equal_to (this->m(), A.m());
-          libmesh_assert_equal_to (this->n(), B.m());
-
-          const unsigned int m_s = A.m();
-          const unsigned int p_s = A.n();
-          const unsigned int n_s = this->n();
-
-          // Do it this way because there is a
-          // decent chance (at least for constraint matrices)
-          // that B.transpose(k,j) = 0.
-          for (unsigned int j=0; j<n_s; j++)
-            for (unsigned int k=0; k<p_s; k++)
-              if (B.transpose(k,j) != 0.)
-                for (unsigned int i=0; i<m_s; i++)
-                  (*this)(i,j) += A(i,k)*B.transpose(k,j);
-        }
-    }
+    this->_right_multiply_transpose(B);
 }
 
 
 
 template<typename T>
 template<typename T2>
 void DenseMatrix<T>::right_multiply_transpose (const DenseMatrix<T2> & B)
+{
+  this->_right_multiply_transpose(B);
+}
+
+
+
+template<typename T>
+template<typename T2>
+void DenseMatrix<T>::_right_multiply_transpose (const DenseMatrix<T2> & B)
 {
   //Check to see if we are doing B*(B^T)
   if (this == &B)

include/solvers/petsc_linear_solver.h

@@ -279,12 +279,26 @@ class PetscLinearSolver : public LinearSolver<T>
    */
   virtual std::pair<unsigned int, Real>
   shell_solve_common (const ShellMatrix<T> & shell_matrix,
-                      const PetscMatrix<T> * precond_matrix,
+                      PetscMatrix<T> * precond_matrix,
                       NumericVector<T> & solution_in,
                       NumericVector<T> & rhs_in,
                       const double tol,
                       const unsigned int m_its);
 
+  /*
+   * Helper function for the helper functions
+   */
+  std::pair<unsigned int, Real>
+  solve_base (SparseMatrix<T> * matrix,
+              PetscMatrix<T> * precond,
+              Mat mat,
+              NumericVector<T> & solution_in,
+              NumericVector<T> & rhs_in,
+              const double tol,
+              const unsigned int m_its,
+              ksp_solve_func_type solve_func);
+
+
   /**
    * Tells PETSC to use the user-specified solver stored in
    * \p _solver_type

include/utils/xdr_cxx.h

@@ -36,12 +36,10 @@
 #include "libmesh/restore_warnings.h"
 #endif
 
+#include <complex>
 #include <iosfwd>
 #include <vector>
 #include <string>
-#ifdef LIBMESH_USE_COMPLEX_NUMBERS
-# include <complex>
-#endif
 
 const unsigned int xdr_MAX_STRING_LENGTH=256;
 
@@ -197,6 +195,27 @@ class Xdr
   template <typename T>
   void do_write(std::vector<std::complex<T>> & a);
 
+  /**
+   * Helper method for complex types
+   */
+  template <typename T>
+  void _complex_data_stream (std::complex<T> * val, const unsigned int len,
+                             const unsigned int line_break);
+
+  /**
+   * Helper method for extended FP types
+   */
+  template <typename XFP>
+  void _xfp_data_stream (XFP * val, const unsigned int len,
+#ifdef LIBMESH_HAVE_XDR
+                         xdrproc_t
+#else
+                         void *
+#endif
+                           xdr_proc,
+                         const unsigned int line_break,
+                         const int n_digits);
+
   /**
    * The mode used for accessing the file.
    */

src/base/dof_map.C

@@ -1274,46 +1274,10 @@ void DofMap::distribute_local_dofs_node_major(dof_id_type & next_free_dof,
             }
       }
 
-  // Finally, count up the SCALAR dofs
-  this->_n_SCALAR_dofs = 0;
-  for (unsigned vg=0; vg<n_var_groups; vg++)
-    {
-      const VariableGroup & vg_description(this->variable_group(vg));
-
-      if (vg_description.type().family == SCALAR)
-        {
-          this->_n_SCALAR_dofs += (vg_description.n_variables()*
-                                   vg_description.type().order.get_order());
-          continue;
-        }
-    }
-
-  // Only increment next_free_dof if we're on the processor
-  // that holds this SCALAR variable
-  if (this->processor_id() == (this->n_processors()-1))
-    next_free_dof += _n_SCALAR_dofs;
+  this->distribute_scalar_dofs(next_free_dof);
 
 #ifdef DEBUG
-  {
-    // libMesh::out << "next_free_dof=" << next_free_dof << std::endl
-    //       << "_n_SCALAR_dofs=" << _n_SCALAR_dofs << std::endl;
-
-    // Make sure we didn't miss any nodes
-    MeshTools::libmesh_assert_valid_procids<Node>(mesh);
-
-    for (auto & node : mesh.local_node_ptr_range())
-      {
-        unsigned int n_var_g = node->n_var_groups(this->sys_number());
-        for (unsigned int vg=0; vg != n_var_g; ++vg)
-          {
-            unsigned int n_comp_g =
-              node->n_comp_group(this->sys_number(), vg);
-            dof_id_type my_first_dof = n_comp_g ?
-              node->vg_dof_base(this->sys_number(), vg) : 0;
-            libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
-          }
-      }
-  }
+  this->assert_no_nodes_missed(mesh);
 #endif // DEBUG
 }
 
@@ -1405,9 +1369,19 @@ void DofMap::distribute_local_dofs_var_major(dof_id_type & next_free_dof,
             }
     } // end loop on variable groups
 
-  // Finally, count up the SCALAR dofs
+  this->distribute_scalar_dofs(next_free_dof);
+
+#ifdef DEBUG
+  this->assert_no_nodes_missed(mesh);
+#endif
+}
+
+
+
+void DofMap::distribute_scalar_dofs(dof_id_type & next_free_dof)
+{
   this->_n_SCALAR_dofs = 0;
-  for (unsigned vg=0; vg<n_var_groups; vg++)
+  for (auto vg : make_range(this->n_variable_groups()))
     {
       const VariableGroup & vg_description(this->variable_group(vg));
 
@@ -1423,28 +1397,29 @@ void DofMap::distribute_local_dofs_var_major(dof_id_type & next_free_dof,
   // that holds this SCALAR variable
   if (this->processor_id() == (this->n_processors()-1))
     next_free_dof += _n_SCALAR_dofs;
+}
+
+
 
 #ifdef DEBUG
-  {
-    // Make sure we didn't miss any nodes
-    MeshTools::libmesh_assert_valid_procids<Node>(mesh);
+void DofMap::assert_no_nodes_missed(MeshBase & mesh)
+{
+  MeshTools::libmesh_assert_valid_procids<Node>(mesh);
 
-    for (auto & node : mesh.local_node_ptr_range())
-      {
-        unsigned int n_var_g = node->n_var_groups(this->sys_number());
-        for (unsigned int vg=0; vg != n_var_g; ++vg)
-          {
-            unsigned int n_comp_g =
-              node->n_comp_group(this->sys_number(), vg);
-            dof_id_type my_first_dof = n_comp_g ?
-              node->vg_dof_base(this->sys_number(), vg) : 0;
-            libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
-          }
-      }
-  }
-#endif // DEBUG
+  for (auto & node : mesh.local_node_ptr_range())
+    {
+      unsigned int n_var_g = node->n_var_groups(this->sys_number());
+      for (unsigned int vg=0; vg != n_var_g; ++vg)
+        {
+          unsigned int n_comp_g =
+            node->n_comp_group(this->sys_number(), vg);
+          dof_id_type my_first_dof = n_comp_g ?
+            node->vg_dof_base(this->sys_number(), vg) : 0;
+          libmesh_assert_not_equal_to (my_first_dof, DofObject::invalid_id);
+        }
+    }
 }
-
+#endif // DEBUG
 
 
 void