Cleanup output system for writing truly constant monomial data

include/systems/equation_systems.h

@@ -36,6 +36,7 @@
 
 // C++ includes
 #include <cstddef>
+#include <functional>
 #include <map>
 #include <set>
 #include <string>
@@ -280,6 +281,24 @@ class EquationSystems : public ReferenceCountedObject<EquationSystems>,
                              const FEType * type=nullptr,
                              const std::set<std::string> * system_names=nullptr) const;
 
+  /**
+   * \returns Whether \p type can be represented as elemental data.
+   *
+   * Elemental data variables are CONSTANT MONOMIALs and CONSTANT
+   * MONOMIAL_VECs, regardless of their p_refinement flag.
+   */
+  static bool is_elemental_data_fe_type (const FEType & type);
+
+  /**
+   * Filter \p var_names to names of variables that can be represented as
+   * elemental data. If \p var_names is empty, all eligible variable names are
+   * returned. Vector-valued variables are decomposed into component names.
+   * If \p system_names!=nullptr, only include names from the specified systems.
+   */
+  void build_elemental_data_variable_names
+    (std::vector<std::string> & var_names,
+     const std::set<std::string> * system_names=nullptr) const;
+
   /**
    * Fill the input vector \p soln with the solution values for the
    * system named \p name.
@@ -294,7 +313,7 @@ class EquationSystems : public ReferenceCountedObject<EquationSystems>,
 
   /**
    * Fill the input vector \p soln with solution values.  The
-   * entries will be in variable-major format (corresponding to
+   * entries will be in node-major format (corresponding to
    * the names from \p build_variable_names()).
    *
    * If systems_names!=nullptr, only include data from the
@@ -365,6 +384,14 @@ class EquationSystems : public ReferenceCountedObject<EquationSystems>,
                          const FEType * type=nullptr,
                          const std::vector<FEType> * types=nullptr) const;
 
+  /**
+   * Finds system and variable numbers for variables that can be represented
+   * as elemental data. See find_variable_numbers() for name filtering,
+   * component decomposition, and sorting details.
+   */
+  std::vector<std::pair<unsigned int, unsigned int>>
+  find_elemental_data_variable_numbers (std::vector<std::string> & names) const;
+
   /**
    * Builds a parallel vector of CONSTANT MONOMIAL and/or components of
    * CONSTANT MONOMIAL_VEC solution values corresponding to the entries
@@ -388,8 +415,8 @@ class EquationSystems : public ReferenceCountedObject<EquationSystems>,
 
   /**
    * Fill the input vector \p soln with solution values.  The
-   * entries will be in variable-major format (corresponding to
-   * the names from \p build_variable_names()).
+   * entries will be in geometry-major format: element, then local
+   * node/vertex, then variable.
 
    * If systems_names!=nullptr, only include data from the
    * specified systems.
@@ -609,6 +636,14 @@ class EquationSystems : public ReferenceCountedObject<EquationSystems>,
   bool _enable_default_ghosting;
 
 private:
+  /**
+   * Implementation detail for find_variable_numbers() variants.
+   */
+  std::vector<std::pair<unsigned int, unsigned int>>
+  find_variable_numbers_by_predicate
+    (std::vector<std::string> & names,
+     const std::function<bool(const FEType &)> & type_filter) const;
+
   /**
    * This function is used in the implementation of add_system,
    * it loops over the nodes and elements of the Mesh, adding the

src/mesh/exodusII_io.C

@@ -72,40 +72,62 @@ namespace
   }
 
 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
+  enum class SolutionOrdering
+  {
+    variable_major,
+    geometry_major
+  };
+
+  std::size_t ordered_solution_index (const std::size_t var,
+                                      const std::size_t geometry_index,
+                                      const std::size_t num_vars,
+                                      const std::size_t num_geometric_entries,
+                                      const SolutionOrdering ordering)
+  {
+    if (ordering == SolutionOrdering::variable_major)
+      return var * num_geometric_entries + geometry_index;
+    else
+      return geometry_index * num_vars + var;
+  }
+
   std::vector<Real>
   complex_soln_components (const std::vector<Number> & soln,
-                           const unsigned int num_vars,
-                           const bool write_complex_abs)
+                           const std::size_t num_vars,
+                           const bool write_complex_abs,
+                           const SolutionOrdering source_ordering,
+                           const SolutionOrdering destination_ordering)
   {
-    unsigned int num_values = soln.size();
-    unsigned int num_elems = num_values / num_vars;
+    const std::size_t num_values = soln.size();
+    libmesh_assert_greater(num_vars, 0);
+    libmesh_assert_equal_to(num_values % num_vars, 0);
+    const std::size_t num_geometric_entries = num_values / num_vars;
 
     // This will contain the real and imaginary parts and the magnitude
     // of the values in soln
-    int nco = write_complex_abs ? 3 : 2;
+    const std::size_t nco = write_complex_abs ? 3 : 2;
+    const std::size_t num_complex_vars = nco * num_vars;
     std::vector<Real> complex_soln(nco * num_values);
 
-    for (unsigned i=0; i<num_vars; ++i)
-      {
-        for (unsigned int j=0; j<num_elems; ++j)
-          {
-            Number value = soln[i*num_vars + j];
-            complex_soln[nco*i*num_elems + j] = value.real();
-          }
-        for (unsigned int j=0; j<num_elems; ++j)
-          {
-            Number value = soln[i*num_vars + j];
-            complex_soln[nco*i*num_elems + num_elems + j] = value.imag();
-          }
-        if (write_complex_abs)
-          {
-            for (unsigned int j=0; j<num_elems; ++j)
-              {
-                Number value = soln[i*num_vars + j];
-                complex_soln[3*i*num_elems + 2*num_elems + j] = std::abs(value);
-              }
-          }
-      }
+    for (const auto var : make_range(num_vars))
+      for (const auto geometry_index : make_range(num_geometric_entries))
+        {
+          const Number value =
+            soln[ordered_solution_index
+                 (var, geometry_index, num_vars, num_geometric_entries, source_ordering)];
+
+          const std::size_t complex_var = nco * var;
+          complex_soln[ordered_solution_index
+                       (complex_var, geometry_index, num_complex_vars,
+                        num_geometric_entries, destination_ordering)] = value.real();
+          complex_soln[ordered_solution_index
+                       (complex_var + 1, geometry_index, num_complex_vars,
+                        num_geometric_entries, destination_ordering)] = value.imag();
+
+          if (write_complex_abs)
+            complex_soln[ordered_solution_index
+                         (complex_var + 2, geometry_index, num_complex_vars,
+                          num_geometric_entries, destination_ordering)] = std::abs(value);
+        }
 
     return complex_soln;
   }
@@ -1170,8 +1192,9 @@ void ExodusII_IO::copy_elemental_solution(System & system,
   //
   // NOTE: Currently, this reader is capable of reading only individual components of MONOMIAL_VEC
   //       types, and each must be written out to its own CONSTANT MONOMIAL variable
-  libmesh_error_msg_if((system.variable_type(var_num) != FEType(CONSTANT, MONOMIAL))
-                       && (system.variable_type(var_num) != FEType(CONSTANT, MONOMIAL_VEC)),
+  const auto & var_type = system.variable_type(var_num);
+  libmesh_error_msg_if(var_type.order != CONSTANT ||
+                       (var_type.family != MONOMIAL && var_type.family != MONOMIAL_VEC),
                        "Error! Trying to copy elemental solution into a variable that is not of CONSTANT MONOMIAL nor CONSTANT MONOMIAL_VEC type.");
 
   const MeshBase & mesh = MeshInput<MeshBase>::mesh();
@@ -1395,25 +1418,10 @@ void ExodusII_IO::write_element_data (const EquationSystems & es)
   // To be (possibly) filled with a filtered list of variable names to output.
   std::vector<std::string> names;
 
-  // If _output_variables is populated, only output the monomials which are
-  // also in the _output_variables vector.
+  // If _output_variables is populated, build_elemental_solution_vector() will filter this list to
+  // the monomial variables that can be written as elemental data.
   if (_output_variables.size() > 0)
-    {
-      // Create a list of CONSTANT MONOMIAL variable names
-      std::vector<std::string> monomials;
-      FEType type(CONSTANT, MONOMIAL);
-      es.build_variable_names(monomials, &type);
-
-      // Now concatenate a list of CONSTANT MONOMIAL_VEC variable names
-      type = FEType(CONSTANT, MONOMIAL_VEC);
-      es.build_variable_names(monomials, &type);
-
-      // Filter that list against the _output_variables list.  Note: if names is still empty after
-      // all this filtering, all the monomial variables will be gathered
-      for (const auto & var : monomials)
-        if (std::find(_output_variables.begin(), _output_variables.end(), var) != _output_variables.end())
-          names.push_back(var);
-    }
+    names.assign(_output_variables.begin(), _output_variables.end());
 
   // If we pass in a list of names to "build_elemental_solution_vector()"
   // it'll filter the variables coming back.
@@ -1446,7 +1454,9 @@ void ExodusII_IO::write_element_data (const EquationSystems & es)
   exio_helper->initialize_element_variables(complex_names, complex_vars_active_subdomains);
 
   const std::vector<Real> complex_soln =
-    complex_soln_components(soln, names.size(), _write_complex_abs);
+    complex_soln_components(soln, names.size(), _write_complex_abs,
+                            SolutionOrdering::variable_major,
+                            SolutionOrdering::variable_major);
 
   exio_helper->write_element_values(mesh, complex_soln, _timestep, complex_vars_active_subdomains);
 
@@ -1487,22 +1497,13 @@ ExodusII_IO::write_element_data_from_discontinuous_nodal_data
   std::vector<std::string> var_names;
   es.build_variable_names (var_names, /*fetype=*/nullptr, system_names);
 
-  // Get a subset of all variable names that are CONSTANT,
-  // MONOMIALs. We treat those slightly differently since they can
+  // Get a subset of all variable names that can be written directly
+  // as elemental data. We treat those slightly differently since they
   // truly only have a single value per Elem.
-  //
-  // Should the same apply here for CONSTANT MONOMIAL_VECs? [CW]
-  // That is, get rid of 'const' on 'fe_type' and rerun:
-  //    fe_type = FEType(CONSTANT, MONOMIAL_VEC);
-  //    es.build_variable_names(monomial_var_names, &fe_type);
-  // Then, es.find_variable_numbers() can be used without a type
-  // (since we know for sure they're monomials) like:
-  //    var_nums = es.find_variable_numbers(monomial_var_names)
-  // for which the DOF indices for 'var_nums' have to be resolved
-  // manually like in build_elemental_solution_vector()
   std::vector<std::string> monomial_var_names;
-  const FEType fe_type(CONSTANT, MONOMIAL);
-  es.build_variable_names(monomial_var_names, &fe_type);
+  if (!_output_variables.empty())
+    monomial_var_names.assign(_output_variables.begin(), _output_variables.end());
+  es.build_elemental_data_variable_names(monomial_var_names, system_names);
 
   // Remove all names from var_names that are not in _output_variables.
   // Note: This approach avoids errors when the user provides invalid
@@ -1519,17 +1520,6 @@ ExodusII_IO::write_element_data_from_discontinuous_nodal_data
                               _output_variables.end(),
                               name);}),
          var_names.end());
-
-      // Also filter the monomial variable names.
-      monomial_var_names.erase
-        (std::remove_if
-         (monomial_var_names.begin(),
-          monomial_var_names.end(),
-          [this](const std::string & name)
-          {return !std::count(_output_variables.begin(),
-                              _output_variables.end(),
-                              name);}),
-         monomial_var_names.end());
     }
 
   // Build a solution vector, limiting the results to the variables in
@@ -2019,7 +2009,9 @@ void ExodusII_IO::write_global_data (const std::vector<Number> & soln,
   exio_helper->initialize_global_variables(complex_names);
 
   const std::vector<Real> complex_soln =
-    complex_soln_components(soln, names.size(), _write_complex_abs);
+    complex_soln_components(soln, names.size(), _write_complex_abs,
+                            SolutionOrdering::variable_major,
+                            SolutionOrdering::variable_major);
 
   exio_helper->write_global_values(complex_soln, _timestep);
 

src/mesh/nemesis_io.C

@@ -1471,30 +1471,16 @@ void Nemesis_IO::write_element_data (const EquationSystems & es)
   // To be (possibly) filled with a filtered list of variable names to output.
   std::vector<std::string> names;
 
-  // All of which should be low order monomials for now
-  const std::vector<FEType> type = {FEType(CONSTANT, MONOMIAL), FEType(CONSTANT, MONOMIAL_VEC)};
-
-  // If _output_variables is populated, only output the monomials which are
-  // also in the _output_variables vector.
+  // If _output_variables is populated, find_elemental_data_variable_numbers()
+  // will filter this list to the monomial variables that can be written as
+  // elemental data.
   if (_output_variables.size())
-    {
-      std::vector<std::string> monomials;
-
-      // Create a list of monomial variable names
-      es.build_variable_names(monomials, &type[0]); /*scalars*/
-      es.build_variable_names(monomials, &type[1]); /*vectors*/
-
-      // Filter that list against the _output_variables list.  Note: if names is still empty after
-      // all this filtering, all the monomial variables will be gathered
-      for (const auto & var : monomials)
-        if (std::find(_output_variables.begin(), _output_variables.end(), var) != _output_variables.end())
-          names.push_back(var);
-    }
+    names.assign(_output_variables.begin(), _output_variables.end());
 
   // The 'names' vector will here be updated with the variable's names
   // that are actually eligible to write
   std::vector<std::pair<unsigned int, unsigned int>> var_nums =
-    es.find_variable_numbers (names, /*type=*/nullptr, &type);
+    es.find_elemental_data_variable_numbers(names);
 
   // find_variable_numbers() can return a nullptr, in which case there are no constant monomial
   // variables to write, and we can just return.
@@ -1756,7 +1742,8 @@ void Nemesis_IO::copy_elemental_solution(System & system,
   parallel_object_only();
 
   const unsigned int var_num = system.variable_number(system_var_name);
-  libmesh_error_msg_if(system.variable_type(var_num) != FEType(CONSTANT, MONOMIAL),
+  const auto & var_type = system.variable_type(var_num);
+  libmesh_error_msg_if(var_type.order != CONSTANT || var_type.family != MONOMIAL,
                        "Error! Trying to copy elemental solution into a variable that is not of CONSTANT MONOMIAL type.");
 
   const MeshBase & mesh = MeshInput<MeshBase>::mesh();

src/mesh/nemesis_io_helper.C

@@ -2649,11 +2649,13 @@ Nemesis_IO_Helper::write_element_values(const MeshBase & mesh,
       const System & system = es.get_system(sys_num);
 
       // We need to check if the constant monomial is a scalar or a vector and set the number of
-      // components as the mesh spatial dimension for the latter as per es.find_variable_numbers().
+      // components for the latter as per es.find_variable_numbers().
       // Even for the case where a variable is not active on any subdomain belonging to the
       // processor, we still need to know this number to update 'var_ctr'.
+      const auto & var_type = system.variable_type(var);
       const unsigned int n_comps =
-        (system.variable_type(var) == FEType(CONSTANT, MONOMIAL_VEC)) ? mesh.spatial_dimension() : 1;
+        (FEInterface::field_type(var_type) == TYPE_VECTOR) ?
+        FEInterface::n_vec_dim(mesh, var_type) : 1;
 
       // Get list of active subdomains for variable v
       const auto & active_subdomains = vars_active_subdomains[v];

src/numerics/petsc_preconditioner.C

@@ -189,9 +189,15 @@ void PetscPreconditioner<T>::set_petsc_aux_data(PC & pc, System & sys, const uns
                            "variables");
       // If not a 1st order Nédélec or a 2d 1st order Raviart-Thomas system, we
       // error out as we do not support anything else at the moment
-      libmesh_error_msg_if(sys.variable(v).type() != FEType(1, NEDELEC_ONE) &&
-                          (sys.variable(v).type() != FEType(1, RAVIART_THOMAS) ||
-                           sys.get_mesh().mesh_dimension() != 2),
+      const FEType & var_type = sys.variable(v).type();
+      const bool first_order_nedelec =
+        var_type.order == FIRST && var_type.family == NEDELEC_ONE;
+      const bool first_order_raviart_thomas =
+        var_type.order == FIRST && var_type.family == RAVIART_THOMAS;
+
+      libmesh_error_msg_if(!first_order_nedelec &&
+                           (!first_order_raviart_thomas ||
+                            sys.get_mesh().mesh_dimension() != 2),
                            "Error applying hypre AMS to a system "
                            "whose variable is not 1st order Nedelec or 1st "
                            "order Raviart-Thomas on a 2d mesh");
@@ -205,7 +211,11 @@ void PetscPreconditioner<T>::set_petsc_aux_data(PC & pc, System & sys, const uns
                            "variables");
       // If not a 3d 1st order Raviart-Thomas system, we error out as we do not
       // support anything else at the moment
-      libmesh_error_msg_if(sys.variable(v).type() != FEType(1, RAVIART_THOMAS) ||
+      const FEType & var_type = sys.variable(v).type();
+      const bool first_order_raviart_thomas =
+        var_type.order == FIRST && var_type.family == RAVIART_THOMAS;
+
+      libmesh_error_msg_if(!first_order_raviart_thomas ||
                            sys.get_mesh().mesh_dimension() != 3,
                            "Error applying hypre ADS to a system "
                            "whose variable is not 1st "