Merge pull request #1171 from LLNL/fix/atkins12/use-dual-and-state-in…

…terfaces Update contact and output to sum over shared components of dual vectors
LLNL · Jul 17, 2024 · 5b96ab5 · 5b96ab5
2 parents e9d8e68 + d32a49f
commit 5b96ab5
Show file tree

Hide file tree

Showing 2 changed files with 12 additions and 16 deletions.
diff --git a/src/serac/physics/base_physics.cpp b/src/serac/physics/base_physics.cpp
@@ -145,8 +145,6 @@ void BasePhysics::CreateParaviewDataCollection() const
   max_order_in_fields = std::max(max_order_in_fields, shape_displacement_.space().GetOrder(0));
 
   shape_sensitivity_grid_function_ = std::make_unique<mfem::ParGridFunction>(&shape_displacement_sensitivity_->space());
-  shape_displacement_sensitivity_->space().GetRestrictionMatrix()->MultTranspose(*shape_displacement_sensitivity_,
-                                                                                 *shape_sensitivity_grid_function_);
   max_order_in_fields = std::max(max_order_in_fields, shape_displacement_sensitivity_->space().GetOrder(0));
   paraview_dc_->RegisterField(shape_displacement_sensitivity_->name(), shape_sensitivity_grid_function_.get());
 
@@ -163,16 +161,16 @@ void BasePhysics::UpdateParaviewDataCollection(const std::string& paraview_outpu
     state->gridFunction();  // update grid function values
   }
   for (const FiniteElementDual* dual : duals_) {
-    // These are really const calls, but MFEM doesn't label them as such
     serac::FiniteElementDual* non_const_dual = const_cast<serac::FiniteElementDual*>(dual);
-    non_const_dual->space().GetRestrictionMatrix()->MultTranspose(*dual, *paraview_dual_grid_functions_[dual->name()]);
+    non_const_dual->linearForm().ParallelAssemble(paraview_dual_grid_functions_[dual->name()]->GetTrueVector());
+    paraview_dual_grid_functions_[dual->name()]->SetFromTrueVector();
   }
   for (auto& parameter : parameters_) {
     parameter.state->gridFunction();
   }
   shape_displacement_.gridFunction();
-  shape_displacement_sensitivity_->space().GetRestrictionMatrix()->MultTranspose(*shape_displacement_sensitivity_,
-                                                                                 *shape_sensitivity_grid_function_);
+  shape_displacement_sensitivity_->linearForm().ParallelAssemble(shape_sensitivity_grid_function_->GetTrueVector());
+  shape_sensitivity_grid_function_->SetFromTrueVector();
 
   // Set the current time, cycle, and requested paraview directory
   paraview_dc_->SetCycle(cycle_);

diff --git a/src/serac/physics/contact/contact_interaction.cpp b/src/serac/physics/contact/contact_interaction.cpp
@@ -67,28 +67,26 @@ ContactInteraction::ContactInteraction(int interaction_id, const mfem::ParMesh&
 FiniteElementDual ContactInteraction::forces() const
 {
   FiniteElementDual f(*current_coords_.ParFESpace());
-  mfem::Vector      f_tribol(current_coords_.ParFESpace()->GetVSize());
-  f_tribol = 0.0;
-  tribol::getMfemResponse(getInteractionId(), f_tribol);
-  current_coords_.ParFESpace()->GetRestrictionMatrix()->Mult(f_tribol, f);
+  auto&             f_loc = f.linearForm();
+  tribol::getMfemResponse(getInteractionId(), f_loc);
+  f.setFromLinearForm(f_loc);
   return f;
 }
 
 FiniteElementState ContactInteraction::pressure() const
 {
   auto&              p_tribol = tribol::getMfemPressure(getInteractionId());
   FiniteElementState p(*p_tribol.ParFESpace());
-  p_tribol.ParFESpace()->GetRestrictionMatrix()->Mult(p_tribol, p);
+  p.setFromGridFunction(p_tribol);
   return p;
 }
 
 FiniteElementDual ContactInteraction::gaps() const
 {
-  auto&             pressure_fes = *tribol::getMfemPressure(getInteractionId()).ParFESpace();
-  FiniteElementDual g(pressure_fes);
-  mfem::Vector      g_tribol;
-  tribol::getMfemGap(getInteractionId(), g_tribol);
-  pressure_fes.GetRestrictionMatrix()->Mult(g_tribol, g);
+  FiniteElementDual g(pressureSpace());
+  auto&             g_loc = g.linearForm();
+  tribol::getMfemGap(getInteractionId(), g_loc);
+  g.setFromLinearForm(g_loc);
   return g;
 }