Merge pull request #4368 from ChengHauYang/project_vector_solution

Allow specifying element range in project_vector/project_solution

Author: roystgnr

include/numerics/wrapped_functor.h

@@ -123,14 +123,28 @@ class WrappedFunctor : public FEMFunctionBase<Output>
 template <typename Output>
 void WrappedFunctor<Output>::init_context (const FEMContext & c)
 {
-  for (auto dim : c.elem_dimensions())
+  // If the context was constructed with an active var subset (e.g. projector
+  // projecting only selected variables), then only those FE objects are
+  // guaranteed to exist.
+  const std::vector<unsigned int> * active = c.active_vars();
+
+  auto check_and_init = [&c](unsigned int v, unsigned short dim)
+    {
+      FEAbstract * fe = nullptr;
+      c.get_element_fe(v, fe, dim);
+      libmesh_assert_msg(fe, "FEAbstract pointer is null for variable "
+                             << v << " in dimension " << dim);
+      fe->get_nothing();
+    };
+
+  for (const auto dim : c.elem_dimensions())
     {
-      for (auto v : make_range(c.n_vars()))
-        {
-          FEAbstract * fe;
-          c.get_element_fe(v, fe, dim);
-          fe->get_nothing();
-        }
+      if (active)
+        for (const auto v : *active)
+          check_and_init(v, dim);
+      else
+        for (const auto v : make_range(c.n_vars()))
+          check_and_init(v, dim);
     }
 }
 

include/systems/system.h

@@ -40,6 +40,7 @@
 #include <string_view>
 #include <vector>
 #include <memory>
+#include <optional>
 
 // This define may be useful in --disable-optional builds when it is
 // possible that libmesh will not have any solvers available.
@@ -514,19 +515,31 @@ class System : public ReferenceCountedObject<System>,
    * user-provided cloneable functors.
    * A gradient \p g is only required/used for projecting onto finite
    * element spaces with continuous derivatives.
+   * elem_range \p active_local_range, if provided, indicates the range of elements
+   * over which to perform the projection.
+   * variable_numbers \p variable_numbers, if provided, indicates the variable numbers
+   * onto which to project.
    */
   void project_solution (FunctionBase<Number> * f,
-                         FunctionBase<Gradient> * g = nullptr) const;
+                         FunctionBase<Gradient> * g = nullptr,
+                         std::optional<ConstElemRange> active_local_range = std::nullopt,
+                         std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;
 
   /**
    * Projects arbitrary functions onto the current solution.
    * The function value \p f and its gradient \p g are
    * user-provided cloneable functors.
    * A gradient \p g is only required/used for projecting onto finite
    * element spaces with continuous derivatives.
+   * elem_range \p active_local_range, if provided, indicates the range of elements
+   * over which to perform the projection.
+   * variable_numbers \p variable_numbers, if provided, indicates the variable numbers
+   * onto which to project.
    */
   void project_solution (FEMFunctionBase<Number> * f,
-                         FEMFunctionBase<Gradient> * g = nullptr) const;
+                         FEMFunctionBase<Gradient> * g = nullptr,
+                         std::optional<ConstElemRange> active_local_range = std::nullopt,
+                         std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;
 
   /**
    * Projects arbitrary functions onto the current solution.
@@ -545,7 +558,9 @@ class System : public ReferenceCountedObject<System>,
                                               const std::string & unknown_name);
   void project_solution (ValueFunctionPointer fptr,
                          GradientFunctionPointer gptr,
-                         const Parameters & parameters) const;
+                         const Parameters & parameters,
+                         std::optional<ConstElemRange> active_local_range = std::nullopt,
+                         std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;
 
   /**
    * Projects arbitrary functions onto a vector of degree of freedom
@@ -554,14 +569,20 @@ class System : public ReferenceCountedObject<System>,
    * user-provided cloneable functors.
    * A gradient \p g is only required/used for projecting onto finite
    * element spaces with continuous derivatives.
+   * elem_range \p active_local_range, if provided, indicates the range of elements
+   * over which to perform the projection.
+   * variable_numbers \p variable_numbers, if provided, indicates the variable numbers
+   * onto which to project.
    *
    * Constrain the new vector using the requested adjoint rather than
    * primal constraints if is_adjoint is non-negative.
    */
   void project_vector (NumericVector<Number> & new_vector,
                        FunctionBase<Number> * f,
                        FunctionBase<Gradient> * g = nullptr,
-                       int is_adjoint = -1) const;
+                       int is_adjoint = -1,
+                       std::optional<ConstElemRange> active_local_range = std::nullopt,
+                       std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;
 
   /**
    * Projects arbitrary functions onto a vector of degree of freedom
@@ -570,14 +591,20 @@ class System : public ReferenceCountedObject<System>,
    * user-provided cloneable functors.
    * A gradient \p g is only required/used for projecting onto finite
    * element spaces with continuous derivatives.
+   * elem_range \p active_local_range, if provided, indicates the range of elements
+   * over which to perform the projection.
+   * variable_numbers \p variable_numbers, if provided, indicates the variable numbers
+   * onto which to project.
    *
    * Constrain the new vector using the requested adjoint rather than
    * primal constraints if is_adjoint is non-negative.
    */
   void project_vector (NumericVector<Number> & new_vector,
                        FEMFunctionBase<Number> * f,
                        FEMFunctionBase<Gradient> * g = nullptr,
-                       int is_adjoint = -1) const;
+                       int is_adjoint = -1,
+                       std::optional<ConstElemRange> active_local_range = std::nullopt,
+                       std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;
 
   /**
    * Projects arbitrary functions onto a vector of degree of freedom
@@ -586,6 +613,10 @@ class System : public ReferenceCountedObject<System>,
    * represented by function pointers.
    * A gradient \p gptr is only required/used for projecting onto
    * finite element spaces with continuous derivatives.
+   * elem_range \p active_local_range, if provided, indicates the range of elements
+   * over which to perform the projection.
+   * variable_numbers \p variable_numbers, if provided, indicates the variable numbers
+   * onto which to project.
    *
    * Constrain the new vector using the requested adjoint rather than
    * primal constraints if is_adjoint is non-negative.
@@ -594,7 +625,9 @@ class System : public ReferenceCountedObject<System>,
                        GradientFunctionPointer gptr,
                        const Parameters & parameters,
                        NumericVector<Number> & new_vector,
-                       int is_adjoint = -1) const;
+                       int is_adjoint = -1,
+                       std::optional<ConstElemRange> active_local_range = std::nullopt,
+                       std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;
 
   /**
    * Projects arbitrary boundary functions onto a vector of degree of
@@ -607,11 +640,14 @@ class System : public ReferenceCountedObject<System>,
    * user-provided cloneable functors.
    * A gradient \p g is only required/used for projecting onto finite
    * element spaces with continuous derivatives.
+   * elem_range \p active_local_range, if provided, indicates the range of elements
+   * over which to perform the projection.
    */
   void boundary_project_solution (const std::set<boundary_id_type> & b,
                                   const std::vector<unsigned int> & variables,
                                   FunctionBase<Number> * f,
-                                  FunctionBase<Gradient> * g = nullptr);
+                                  FunctionBase<Gradient> * g = nullptr,
+                                  std::optional<ConstElemRange> active_local_range = std::nullopt);
 
   /**
    * Projects arbitrary boundary functions onto a vector of degree of
@@ -624,12 +660,15 @@ class System : public ReferenceCountedObject<System>,
    * represented by function pointers.
    * A gradient \p gptr is only required/used for projecting onto
    * finite element spaces with continuous derivatives.
+   * elem_range \p active_local_range, if provided, indicates the range of elements
+   * over which to perform the projection.
    */
   void boundary_project_solution (const std::set<boundary_id_type> & b,
                                   const std::vector<unsigned int> & variables,
                                   ValueFunctionPointer fptr,
                                   GradientFunctionPointer gptr,
-                                  const Parameters & parameters);
+                                  const Parameters & parameters,
+                                  std::optional<ConstElemRange> active_local_range = std::nullopt);
 
   /**
    * Projects arbitrary boundary functions onto a vector of degree of
@@ -642,6 +681,8 @@ class System : public ReferenceCountedObject<System>,
    * user-provided cloneable functors.
    * A gradient \p g is only required/used for projecting onto finite
    * element spaces with continuous derivatives.
+   * elem_range \p active_local_range, if provided, indicates the range of elements
+   * over which to perform the projection.
    *
    * Constrain the new vector using the requested adjoint rather than
    * primal constraints if is_adjoint is non-negative.
@@ -651,7 +692,8 @@ class System : public ReferenceCountedObject<System>,
                                 NumericVector<Number> & new_vector,
                                 FunctionBase<Number> * f,
                                 FunctionBase<Gradient> * g = nullptr,
-                                int is_adjoint = -1) const;
+                                int is_adjoint = -1,
+                                std::optional<ConstElemRange> active_local_range = std::nullopt) const;
 
   /**
    * Projects arbitrary boundary functions onto a vector of degree of
@@ -664,6 +706,8 @@ class System : public ReferenceCountedObject<System>,
    * represented by function pointers.
    * A gradient \p gptr is only required/used for projecting onto
    * finite element spaces with continuous derivatives.
+   * elem_range \p active_local_range, if provided, indicates the range of elements
+   * over which to perform the projection.
    *
    * Constrain the new vector using the requested adjoint rather than
    * primal constraints if is_adjoint is non-negative.
@@ -674,7 +718,8 @@ class System : public ReferenceCountedObject<System>,
                                 GradientFunctionPointer gptr,
                                 const Parameters & parameters,
                                 NumericVector<Number> & new_vector,
-                                int is_adjoint = -1) const;
+                                int is_adjoint = -1,
+                                std::optional<ConstElemRange> active_local_range = std::nullopt) const;
 
   /**
    * \returns The system number.
@@ -1984,7 +2029,9 @@ class System : public ReferenceCountedObject<System>,
    * primal constraints if is_adjoint is non-negative.
    */
   void project_vector (NumericVector<Number> &,
-                       int is_adjoint = -1) const;
+                       int is_adjoint = -1,
+                       std::optional<ConstElemRange> active_local_range = std::nullopt,
+                       std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;
 
   /**
    * Projects the vector defined on the old mesh onto the
@@ -1996,7 +2043,9 @@ class System : public ReferenceCountedObject<System>,
    */
   void project_vector (const NumericVector<Number> &,
                        NumericVector<Number> &,
-                       int is_adjoint = -1) const;
+                       int is_adjoint = -1,
+                       std::optional<ConstElemRange> active_local_range = std::nullopt,
+                       std::optional<std::vector<unsigned int>> variable_numbers = std::nullopt) const;
 
   /**
    * Whether this object should condense out constrained degrees of freedom