Skip to content

Commit 74949e9

Browse files
ebchinebchin
committed
make fields optional
1 parent 40b2699 commit 74949e9

3 files changed

Lines changed: 72 additions & 32 deletions

File tree

‎src/smith/physics/contact/contact_data.cpp‎

Lines changed: 38 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -74,42 +74,58 @@ void ContactData::reset()
7474
}
7575
}
7676

77-
void ContactData::updateGaps(int cycle, double time, double& dt, const mfem::Vector& u_shape, const mfem::Vector& u,
78-
bool eval_jacobian)
77+
void ContactData::updateGaps(int cycle, double time, double& dt,
78+
std::optional<std::reference_wrapper<const mfem::Vector>> u_shape,
79+
std::optional<std::reference_wrapper<const mfem::Vector>> u, bool eval_jacobian)
7980
{
8081
cycle_ = cycle;
8182
time_ = time;
8283
dt_ = dt;
8384

84-
setDisplacements(u_shape, u);
85+
if (u_shape && u) {
86+
setDisplacements(u_shape->get(), u->get());
87+
}
8588

8689
for (auto& interaction : interactions_) {
8790
interaction.evalJacobian(eval_jacobian);
8891
}
8992
// This updates the redecomposed surface mesh based on the current displacement, then transfers field quantities to
9093
// the updated mesh.
91-
tribol::updateMfemParallelDecomposition();
94+
if (u_shape && u) {
95+
tribol::updateMfemParallelDecomposition();
96+
}
9297
// This function computes gaps (and optionally geometric Jacobian blocks) based on the current mesh.
9398
tribol::update(cycle, time, dt);
9499
}
95100

96-
void ContactData::update(int cycle, double time, double& dt, const mfem::Vector& u_shape, const mfem::Vector& u,
97-
const mfem::Vector& p)
101+
void ContactData::update(int cycle, double time, double& dt,
102+
std::optional<std::reference_wrapper<const mfem::Vector>> u_shape,
103+
std::optional<std::reference_wrapper<const mfem::Vector>> u,
104+
std::optional<std::reference_wrapper<const mfem::Vector>> p)
98105
{
99-
// First pass: update gaps
100-
updateGaps(cycle, time, dt, u_shape, u);
106+
// First pass: update gaps if coordinates are provided
107+
if (u_shape && u) {
108+
updateGaps(cycle, time, dt, u_shape, u, false);
109+
} else {
110+
// Ensure internal timing is updated even if coordinates are not
111+
cycle_ = cycle;
112+
time_ = time;
113+
dt_ = dt;
114+
}
101115

102-
// with updated gaps, we can update pressure for contact interactions (active set detection and penalty)
103-
setPressures(p);
116+
// second pass: update pressures and compute forces/Jacobians if p is provided
117+
if (p) {
118+
// with updated gaps, we can update pressure for contact interactions (active set detection and penalty)
119+
setPressures(p->get());
104120

105-
// second pass: compute forces and Jacobians
106-
for (auto& interaction : interactions_) {
107-
interaction.evalJacobian(true);
121+
for (auto& interaction : interactions_) {
122+
interaction.evalJacobian(true);
123+
}
124+
// This second call is required to synchronize the updated pressures to Tribol's internal redecomposed surface mesh
125+
// and to ensure Tribol's internal state is correctly reset for the second pass.
126+
tribol::updateMfemParallelDecomposition();
127+
tribol::update(cycle, time, dt);
108128
}
109-
// This second call is required to synchronize the updated pressures to Tribol's internal redecomposed surface mesh
110-
// and to ensure Tribol's internal state is correctly reset for the second pass.
111-
tribol::updateMfemParallelDecomposition();
112-
tribol::update(cycle, time, dt);
113129
}
114130

115131
FiniteElementDual ContactData::forces() const
@@ -468,14 +484,16 @@ void ContactData::addContactInteraction([[maybe_unused]] int interaction_id,
468484
}
469485

470486
void ContactData::updateGaps([[maybe_unused]] int cycle, [[maybe_unused]] double time, [[maybe_unused]] double& dt,
471-
[[maybe_unused]] const mfem::Vector& u_shape, [[maybe_unused]] const mfem::Vector& u,
487+
[[maybe_unused]] std::optional<std::reference_wrapper<const mfem::Vector>> u_shape,
488+
[[maybe_unused]] std::optional<std::reference_wrapper<const mfem::Vector>> u,
472489
[[maybe_unused]] bool eval_jacobian)
473490
{
474491
}
475492

476493
void ContactData::update([[maybe_unused]] int cycle, [[maybe_unused]] double time, [[maybe_unused]] double& dt,
477-
[[maybe_unused]] const mfem::Vector& u_shape, [[maybe_unused]] const mfem::Vector& u,
478-
[[maybe_unused]] const mfem::Vector& p)
494+
[[maybe_unused]] std::optional<std::reference_wrapper<const mfem::Vector>> u_shape,
495+
[[maybe_unused]] std::optional<std::reference_wrapper<const mfem::Vector>> u,
496+
[[maybe_unused]] std::optional<std::reference_wrapper<const mfem::Vector>> p)
479497
{
480498
}
481499

‎src/smith/physics/contact/contact_data.hpp‎

Lines changed: 14 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,8 @@
1515
#include <memory>
1616
#include <set>
1717
#include <vector>
18+
#include <optional>
19+
#include <functional>
1820

1921
#include "mfem.hpp"
2022

@@ -75,11 +77,13 @@ class ContactData {
7577
* @param cycle The current simulation cycle
7678
* @param time The current time
7779
* @param dt The timestep size to attempt
78-
* @param u_shape Shape displacement vector
79-
* @param u Current displacement dof values
80+
* @param u_shape Optional shape displacement vector
81+
* @param u Optional current displacement dof values
8082
* @param eval_jacobian Whether to also evaluate the Jacobian contributions (default false)
8183
*/
82-
void updateGaps(int cycle, double time, double& dt, const mfem::Vector& u_shape, const mfem::Vector& u,
84+
void updateGaps(int cycle, double time, double& dt,
85+
std::optional<std::reference_wrapper<const mfem::Vector>> u_shape = std::nullopt,
86+
std::optional<std::reference_wrapper<const mfem::Vector>> u = std::nullopt,
8387
bool eval_jacobian = false);
8488

8589
/**
@@ -88,12 +92,14 @@ class ContactData {
8892
* @param cycle The current simulation cycle
8993
* @param time The current time
9094
* @param dt The timestep size to attempt
91-
* @param u_shape Shape displacement vector
92-
* @param u Current displacement dof values
93-
* @param p Current pressure true dof values
95+
* @param u_shape Optional shape displacement vector
96+
* @param u Optional current displacement dof values
97+
* @param p Optional current pressure true dof values
9498
*/
95-
void update(int cycle, double time, double& dt, const mfem::Vector& u_shape, const mfem::Vector& u,
96-
const mfem::Vector& p);
99+
void update(int cycle, double time, double& dt,
100+
std::optional<std::reference_wrapper<const mfem::Vector>> u_shape = std::nullopt,
101+
std::optional<std::reference_wrapper<const mfem::Vector>> u = std::nullopt,
102+
std::optional<std::reference_wrapper<const mfem::Vector>> p = std::nullopt);
97103

98104
/**
99105
* @brief Resets the contact pressures to zero

‎src/smith/physics/contact_constraint.hpp‎

Lines changed: 20 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -144,7 +144,11 @@ class ContactConstraint : public Constraint {
144144
// otherwise use previously cached Jacobian
145145
if (update_fields || fresh_derivative) {
146146
int cycle = 0;
147-
contact_.updateGaps(cycle, time, dt, *fields[ContactFields::SHAPE], *fields[ContactFields::DISP], true);
147+
if (update_fields) {
148+
contact_.updateGaps(cycle, time, dt, *fields[ContactFields::SHAPE], *fields[ContactFields::DISP], true);
149+
} else {
150+
contact_.updateGaps(cycle, time, dt, std::nullopt, std::nullopt, true);
151+
}
148152
J_contact_ = contact_.mergedJacobian();
149153
}
150154
// obtain (1, 0) block entry from the 2 x 2 block contact linear system
@@ -172,7 +176,11 @@ class ContactConstraint : public Constraint {
172176
SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative");
173177
int cycle = 0;
174178
if (update_fields || fresh_derivative) {
175-
contact_.update(cycle, time, dt, *fields[ContactFields::SHAPE], *fields[ContactFields::DISP], multipliers);
179+
if (update_fields) {
180+
contact_.update(cycle, time, dt, *fields[ContactFields::SHAPE], *fields[ContactFields::DISP], multipliers);
181+
} else {
182+
contact_.update(cycle, time, dt, std::nullopt, std::nullopt, multipliers);
183+
}
176184
}
177185
return contact_.forces();
178186
};
@@ -199,7 +207,11 @@ class ContactConstraint : public Constraint {
199207

200208
int cycle = 0;
201209
if (update_fields || fresh_derivative) {
202-
contact_.update(cycle, time, dt, *fields[ContactFields::SHAPE], *fields[ContactFields::DISP], multipliers);
210+
if (update_fields) {
211+
contact_.update(cycle, time, dt, *fields[ContactFields::SHAPE], *fields[ContactFields::DISP], multipliers);
212+
} else {
213+
contact_.update(cycle, time, dt, std::nullopt, std::nullopt, multipliers);
214+
}
203215
J_contact_ = contact_.mergedJacobian();
204216
}
205217
// obtain (0, 0) block entry from the 2 x 2 block contact linear system
@@ -226,7 +238,11 @@ class ContactConstraint : public Constraint {
226238
int cycle = 0;
227239
if (update_fields || fresh_derivative) {
228240
mfem::Vector p = contact_.mergedPressures();
229-
contact_.update(cycle, time, dt, *fields[ContactFields::SHAPE], *fields[ContactFields::DISP], p);
241+
if (update_fields) {
242+
contact_.update(cycle, time, dt, *fields[ContactFields::SHAPE], *fields[ContactFields::DISP], p);
243+
} else {
244+
contact_.update(cycle, time, dt, std::nullopt, std::nullopt, p);
245+
}
230246
J_contact_ = contact_.mergedJacobian();
231247
}
232248
// obtain (0, 1) block entry from the 2 x 2 block contact linear system

0 commit comments

Comments
 (0)