Recover boundary elements in zero algebraic range explicit NEML2

framework/include/neml2/userobjects/NEML2ModelExecutor.h

@@ -41,6 +41,7 @@ class NEML2ModelExecutor : public NEML2ModelInterface<GeneralUserObject>
   void finalize() override;
 
   void initialSetup() override;
+  void timestepSetup() override;
 
   /// Get the batch index for the given element ID
   std::size_t getBatchIndex(dof_id_type elem_id) const;
@@ -87,6 +88,9 @@ class NEML2ModelExecutor : public NEML2ModelInterface<GeneralUserObject>
   /// The NEML2BatchIndexGenerator used to generate the element-to-batch-index map
   const NEML2BatchIndexGenerator & _batch_index_generator;
 
+  /// Advance state on device (rather than via MOSOE material properties)
+  const bool _advance_step_on_device;
+
   /// flag that indicates if output data has been fully computed
   bool _output_ready;
 

framework/include/utils/MooseUtils.h

@@ -648,14 +648,10 @@ getUnion(const std::vector<T> & vector1, const std::vector<T> & vector2, std::ve
 {
   std::unordered_set<T> unique_elements;
   unique_elements.reserve(vector1.size() + vector2.size());
-
-  for (const T & entry : vector1)
-    unique_elements.insert(entry);
-  for (const T & entry : vector2)
-    unique_elements.insert(entry);
+  unique_elements.insert(vector1.begin(), vector1.end());
+  unique_elements.insert(vector2.begin(), vector2.end());
 
   // Now populate the common vector with the union
-  common.clear();
   common.assign(unique_elements.begin(), unique_elements.end());
 }
 

framework/src/neml2/userobjects/NEML2ModelExecutor.C

@@ -31,6 +31,12 @@ NEML2ModelExecutor::actionParams()
       "List of NEML2 variables to skip error checking when setting up the model input. If an "
       "input variable is skipped, its value will stay zero. If a required input variable is "
       "not skipped, an error will be raised.");
+  params.addParam<bool>(
+      "advance_step_on_device",
+      false,
+      "Keep state and forces on the device and advance it to old_state and old_forces without a "
+      "roundtrip through MOOSE materials. This is only recommended for explicit time integration "
+      "or when absolutely no restepping occurs (e.g. failed timesteps).");
   return params;
 }
 
@@ -68,6 +74,7 @@ NEML2ModelExecutor::NEML2ModelExecutor(const InputParameters & params)
 #ifdef NEML2_ENABLED
     ,
     _batch_index_generator(getUserObject<NEML2BatchIndexGenerator>("batch_index_generator")),
+    _advance_step_on_device(getParam<bool>("advance_step_on_device")),
     _output_ready(false),
     _error_message("")
 #endif
@@ -167,6 +174,29 @@ NEML2ModelExecutor::initialSetup()
           "those called `NEML2To*MOOSEMaterialProperty`.");
 }
 
+void
+NEML2ModelExecutor::timestepSetup()
+{
+  if (_advance_step_on_device && _t_step > 0)
+  {
+    // Set old state variables and old forces
+    if (model().input_axis().has_subaxis(neml2::OLD_STATE) &&
+        model().output_axis().has_subaxis(neml2::STATE))
+    {
+      const auto & input_old_state = model().input_axis().subaxis(neml2::OLD_STATE);
+      for (const auto & var : input_old_state.variable_names())
+        _in[var.prepend(neml2::OLD_STATE)] = _out[var.prepend(neml2::STATE)];
+    }
+    if (model().input_axis().has_subaxis(neml2::OLD_FORCES) &&
+        model().input_axis().has_subaxis(neml2::FORCES))
+    {
+      const auto & input_old_forces = model().input_axis().subaxis(neml2::OLD_FORCES);
+      for (const auto & var : input_old_forces.variable_names())
+        _in[var.prepend(neml2::OLD_FORCES)] = _in[var.prepend(neml2::FORCES)];
+    }
+  }
+}
+
 std::size_t
 NEML2ModelExecutor::getBatchIndex(dof_id_type elem_id) const
 {

modules/solid_mechanics/include/timeintegrators/NEML2CentralDifference.h

@@ -36,7 +36,7 @@ class NEML2CentralDifference : public ExplicitMixedOrder
 
 private:
   /// Empty element vector to help zero out the algebraic range
-  std::vector<const Elem *> _no_elem = {};
+  std::vector<const Elem *> _boundary_elems = {};
 
   /// Empty node vector to help zero out the algebraic range
   std::vector<const Node *> _no_node = {};

modules/solid_mechanics/src/timeintegrators/NEML2CentralDifference.C

@@ -7,6 +7,9 @@
 //* Licensed under LGPL 2.1, please see LICENSE for details
 //* https://www.gnu.org/licenses/lgpl-2.1.html
 
+#include "MooseTypes.h"
+#include "IntegratedBCBase.h"
+
 #ifdef NEML2_ENABLED
 
 // MOOSE includes
@@ -35,6 +38,27 @@ NEML2CentralDifference::NEML2CentralDifference(const InputParameters & parameter
 void
 NEML2CentralDifference::initialSetup()
 {
+  // build boundary element list by iterating over all integrated BCs
+  const auto & ibcs = _nl->getIntegratedBCWarehouse();
+  std::unordered_set<BoundaryID> bnds;
+  for (const auto & ibc : ibcs.getObjects())
+    bnds.insert(ibc->boundaryIDs().begin(), ibc->boundaryIDs().end());
+
+  if (!bnds.empty())
+  {
+    mooseInfo("Dectected BCs on ", bnds.size(), " boundaries.");
+
+    // deduplicate elements that have multiple boundaries
+    std::unordered_set<const Elem *> unique_elems;
+    const auto end = _fe_problem.mesh().bndElemsEnd();
+    for (auto it = _fe_problem.mesh().bndElemsBegin(); it != end; ++it)
+      if (bnds.find((*it)->_bnd_id) != bnds.end())
+        unique_elems.insert((*it)->_elem);
+
+    _boundary_elems.assign(unique_elems.begin(), unique_elems.end());
+    mooseInfo("Adding ", _boundary_elems.size(), " elements to the algebraic range.");
+  }
+
   ExplicitMixedOrder::initialSetup();
   _neml2_assembly = &_fe_problem.getUserObject<NEML2Assembly>("assembly", /*tid=*/0);
   _fe = &_fe_problem.getUserObject<NEML2FEInterpolation>("fe", /*tid=*/0);
@@ -52,8 +76,8 @@ NEML2CentralDifference::evaluateRHSResidual()
 {
   if (_fe->contextUpToDate() && _neml2_assembly->upToDate())
   {
-    libMesh::ConstElemRange null_elem_range(&_no_elem);
-    _fe_problem.setCurrentAlgebraicElementRange(&null_elem_range);
+    libMesh::ConstElemRange boundary_elem_range(&_boundary_elems);
+    _fe_problem.setCurrentAlgebraicElementRange(&boundary_elem_range);
 
     libMesh::ConstNodeRange null_node_range(&_no_node);
     _fe_problem.setCurrentAlgebraicNodeRange(&null_node_range);