Skip to content

Commit 98ea7ef

Browse files
dragos-anadragos-ana
authored
Merge pull request #1778 from dragos-ana/pass-evaluation-context-solid-surface-elements
Fix material evaluation in Nitsche trace estimation for solid-surface elements
2 parents 597ed12 + 756cf96 commit 98ea7ef

19 files changed

Lines changed: 595 additions & 66 deletions

‎src/contact/src/4C_contact_manager.cpp‎

Lines changed: 72 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,7 @@
3636
#include "4C_linear_solver_method.hpp"
3737
#include "4C_mortar_input.hpp"
3838
#include "4C_mortar_utils.hpp"
39+
#include "4C_utils_exceptions.hpp"
3940

4041
#include <Teuchos_StandardParameterEntryValidators.hpp>
4142

@@ -1083,20 +1084,79 @@ bool CONTACT::Manager::read_and_check_input(Teuchos::ParameterList& cparams) con
10831084
cparams.setParameters(contact);
10841085
cparams.setParameters(wearlist);
10851086
cparams.setParameters(tsic);
1086-
if (problemtype == Core::ProblemType::tsi)
1087-
cparams.set<double>(
1088-
"TIMESTEP", Global::Problem::instance()->tsi_dynamic_params().get<double>("TIMESTEP"));
1089-
else if (problemtype != Core::ProblemType::structure)
1087+
1088+
1089+
switch (problemtype)
10901090
{
1091-
// rauch 01/16
1092-
if (Core::Communication::my_mpi_rank(get_comm()) == 0)
1093-
std::cout << "\n \n Warning: CONTACT::Manager::read_and_check_input() reads TIMESTEP = "
1094-
<< stru.get<double>("TIMESTEP") << " from --STRUCTURAL DYNAMIC \n"
1095-
<< std::endl;
1096-
cparams.set<double>("TIMESTEP", stru.get<double>("TIMESTEP"));
1091+
case Core::ProblemType::tsi:
1092+
{
1093+
const double timestep =
1094+
Global::Problem::instance()->tsi_dynamic_params().get<double>("TIMESTEP");
1095+
1096+
cparams.set<double>("TIMESTEP", timestep);
1097+
1098+
break;
1099+
}
1100+
case Core::ProblemType::structure:
1101+
{
1102+
const double timestep = stru.get<double>("TIMESTEP");
1103+
1104+
cparams.set<double>("TIMESTEP", timestep);
1105+
break;
1106+
}
1107+
case Core::ProblemType::poroelast:
1108+
case Core::ProblemType::poroscatra:
1109+
{
1110+
const Teuchos::ParameterList& porodyn =
1111+
Global::Problem::instance()->poroelast_dynamic_params();
1112+
const double timestep = porodyn.get<double>("TIMESTEP");
1113+
1114+
cparams.set<double>("TIMESTEP", timestep);
1115+
break;
1116+
}
1117+
case Core::ProblemType::ssi:
1118+
{
1119+
const Teuchos::ParameterList& ssi = Global::Problem::instance()->ssi_control_params();
1120+
const double timestep = ssi.get<double>("TIMESTEP");
1121+
1122+
cparams.set<double>("TIMESTEP", timestep);
1123+
break;
1124+
}
1125+
case Core::ProblemType::fsi:
1126+
{
1127+
const Teuchos::ParameterList& fsi = Global::Problem::instance()->fsi_dynamic_params();
1128+
const double timestep = fsi.get<double>("TIMESTEP");
1129+
1130+
cparams.set<double>("TIMESTEP", timestep);
1131+
1132+
break;
1133+
}
1134+
case Core::ProblemType::fsi_xfem:
1135+
{
1136+
const Teuchos::ParameterList& fsi_xfem = Global::Problem::instance()->fluid_dynamic_params();
1137+
const double timestep = fsi_xfem.get<double>("TIMESTEP");
1138+
1139+
cparams.set<double>("TIMESTEP", timestep);
1140+
1141+
break;
1142+
}
1143+
1144+
case Core::ProblemType::fpsi:
1145+
case Core::ProblemType::fpsi_xfem:
1146+
{
1147+
const Teuchos::ParameterList& fpsi = Global::Problem::instance()->fpsi_dynamic_params();
1148+
const double timestep = fpsi.get<double>("TIMESTEP");
1149+
1150+
cparams.set<double>("TIMESTEP", timestep);
1151+
1152+
break;
1153+
}
1154+
1155+
default:
1156+
FOUR_C_THROW("Problem type {} not implemented for old contact formulation!",
1157+
EnumTools::enum_name(problemtype));
10971158
}
1098-
else
1099-
cparams.set<double>("TIMESTEP", stru.get<double>("TIMESTEP"));
1159+
11001160

11011161
// *********************************************************************
11021162
// NURBS contact

‎src/contact/src/4C_contact_nitsche_strategy.cpp‎

Lines changed: 32 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@
1515
#include "4C_global_data.hpp"
1616
#include "4C_linalg_fevector.hpp"
1717
#include "4C_linalg_utils_sparse_algebra_manipulation.hpp"
18+
#include "4C_mat_so3_material.hpp"
1819

1920
#include <Teuchos_Time.hpp>
2021

@@ -98,7 +99,11 @@ void CONTACT::NitscheStrategy::do_read_restart(Core::IO::DiscretizationReader& r
9899
store_to_old(Mortar::StrategyBase::n_old);
99100
}
100101

101-
if (params().get<bool>("NITSCHE_PENALTY_ADAPTIVE")) update_trace_ineq_etimates();
102+
103+
save_time_step_size_and_total_time(cparams_ptr->get_delta_time(), cparams_ptr->get_total_time());
104+
105+
106+
if (params().get<bool>("NITSCHE_PENALTY_ADAPTIVE")) update_trace_ineq_estimates();
102107
}
103108

104109
void CONTACT::NitscheStrategy::set_state(
@@ -221,10 +226,13 @@ void CONTACT::NitscheStrategy::integrate(const CONTACT::ParamsInterface& cparams
221226
// time measurement (on each processor)
222227
const double t_start = Teuchos::Time::wallTime();
223228

229+
// set time step size and total time
230+
save_time_step_size_and_total_time(cparams.get_delta_time(), cparams.get_total_time());
231+
224232
// Evaluation for all interfaces
225233
for (const auto& interface : interface_)
226234
{
227-
interface->interface_params().set<double>("TIMESTEP", cparams.get_delta_time());
235+
interface->interface_params().set<double>("TIMESTEP", time_step_size_);
228236
interface->initialize();
229237
interface->evaluate(0, step_, iter_);
230238

@@ -384,9 +392,17 @@ void CONTACT::NitscheStrategy::setup(bool redistributed, bool init)
384392
curr_state_eval_ = false;
385393
}
386394

387-
void CONTACT::NitscheStrategy::update_trace_ineq_etimates()
395+
void CONTACT::NitscheStrategy::update_trace_ineq_estimates()
388396
{
389397
auto NitWgt = Teuchos::getIntegralValue<CONTACT::NitscheWeighting>(params(), "NITSCHE_WEIGHTING");
398+
399+
// create material evaluation context
400+
Mat::EvaluationContext mat_eval_context{.total_time = &total_time_,
401+
.time_step_size = &time_step_size_,
402+
.xi = nullptr,
403+
.ref_coords = nullptr};
404+
405+
390406
for (const auto& interface : interface_)
391407
{
392408
for (int e = 0; e < interface->discret().element_col_map()->num_my_elements(); ++e)
@@ -395,14 +411,14 @@ void CONTACT::NitscheStrategy::update_trace_ineq_etimates()
395411
interface->discret().g_element(interface->discret().element_col_map()->gid(e)));
396412
if (NitWgt == CONTACT::NitscheWeighting::slave && !mele->is_slave()) continue;
397413
if (NitWgt == CONTACT::NitscheWeighting::master && mele->is_slave()) continue;
398-
mele->estimate_nitsche_trace_max_eigenvalue();
414+
mele->estimate_nitsche_trace_max_eigenvalue(mat_eval_context);
399415
}
400416
}
401417
}
402418

403419
void CONTACT::NitscheStrategy::update(std::shared_ptr<const Core::LinAlg::Vector<double>> dis)
404420
{
405-
if (params().get<bool>("NITSCHE_PENALTY_ADAPTIVE")) update_trace_ineq_etimates();
421+
if (params().get<bool>("NITSCHE_PENALTY_ADAPTIVE")) update_trace_ineq_estimates();
406422
if (friction_)
407423
{
408424
store_to_old(Mortar::StrategyBase::n_old);
@@ -422,7 +438,8 @@ void CONTACT::NitscheStrategy::evaluate_reference_state()
422438
store_to_old(Mortar::StrategyBase::n_old);
423439
}
424440

425-
update_trace_ineq_etimates();
441+
442+
update_trace_ineq_estimates();
426443
}
427444

428445

@@ -461,4 +478,13 @@ void CONTACT::NitscheStrategy::reconnect_parent_elements()
461478
}
462479
}
463480

481+
void CONTACT::NitscheStrategy::save_time_step_size_and_total_time(
482+
const double time_step_size, const double total_time)
483+
{
484+
// save time step size and total time
485+
time_step_size_ = time_step_size;
486+
total_time_ = total_time;
487+
}
488+
489+
464490
FOUR_C_NAMESPACE_CLOSE

‎src/contact/src/4C_contact_nitsche_strategy.hpp‎

Lines changed: 36 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,8 @@ namespace CONTACT
3939
NodeRowMap, params, dim, comm, alphaf, maxdof),
4040
interface_(std::move(interface)),
4141
curr_state_eval_(false)
42-
{ /* empty */
42+
{
43+
init_time_step_size_and_total_time(params);
4344
}
4445

4546
//! Shared data constructor
@@ -51,7 +52,8 @@ namespace CONTACT
5152
: AbstractStrategy(data_ptr, dof_row_map, NodeRowMap, params, dim, comm, alphaf, maxdof),
5253
interface_(std::move(interface)),
5354
curr_state_eval_(false)
54-
{ /* empty */
55+
{
56+
init_time_step_size_and_total_time(params);
5557
}
5658

5759
// don't want = operator and cctor
@@ -93,7 +95,7 @@ namespace CONTACT
9395
*/
9496
void setup(bool redistributed, bool init) override;
9597

96-
virtual void update_trace_ineq_etimates();
98+
virtual void update_trace_ineq_estimates();
9799

98100
/*! \brief Get dirichlet B.C. status and store into Nodes
99101
@@ -274,11 +276,42 @@ namespace CONTACT
274276
virtual std::shared_ptr<Core::LinAlg::SparseMatrix> create_matrix_block_ptr(
275277
const MatBlockType& bt);
276278

279+
280+
/*!
281+
* @brief Save time step size \f$ \Delta t \f$, and total time \f$ t_{n+1} \f$
282+
*
283+
* @param[in] time_step_size time step size \f$ \Delta t \f$
284+
* @param[in] total time total time \f$ t_{n+1} \f$
285+
*
286+
*/
287+
void save_time_step_size_and_total_time(const double time_step_size, const double total_time);
288+
289+
/*!
290+
* @brief Initialize time step size \f$ \Delta t \f$, and total time \f$ t_{n+1} \f$
291+
*
292+
* @param[in] params parameter list containing timestep for evaluating the reference state
293+
*/
294+
void init_time_step_size_and_total_time(const Teuchos::ParameterList& params)
295+
{
296+
const double time_step_size =
297+
params.isParameter("TIMESTEP") ? params.get<double>("TIMESTEP") : 0.0;
298+
const double dummy_total_time = 0.0;
299+
save_time_step_size_and_total_time(time_step_size, dummy_total_time);
300+
}
301+
302+
277303
std::vector<std::shared_ptr<CONTACT::Interface>> interface_;
278304

279305
std::shared_ptr<Core::LinAlg::Vector<double>> curr_state_;
280306
bool curr_state_eval_;
281307

308+
//! time step size
309+
double time_step_size_;
310+
311+
//! total time \f$ t_{n+1} \f$
312+
double total_time_;
313+
314+
282315
std::shared_ptr<Core::LinAlg::FEVector<double>> fc_;
283316
std::shared_ptr<Core::LinAlg::SparseMatrix> kc_;
284317
};

‎src/contact/src/4C_contact_nitsche_strategy_fpi.hpp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,8 @@ namespace CONTACT
4141
weighting_(
4242
Teuchos::getIntegralValue<CONTACT::NitscheWeighting>(params, "NITSCHE_WEIGHTING"))
4343
{
44+
init_time_step_size_and_total_time(params);
45+
4446
if (Teuchos::getIntegralValue<CONTACT::FrictionType>(params, "FRICTION") !=
4547
CONTACT::FrictionType::none)
4648
FOUR_C_THROW("NitscheStrategyFpi: No frictional contact implemented for Nitsche FPSCI!");
@@ -57,6 +59,8 @@ namespace CONTACT
5759
weighting_(
5860
Teuchos::getIntegralValue<CONTACT::NitscheWeighting>(params, "NITSCHE_WEIGHTING"))
5961
{
62+
init_time_step_size_and_total_time(params);
63+
6064
if (Teuchos::getIntegralValue<CONTACT::FrictionType>(params, "FRICTION") !=
6165
CONTACT::FrictionType::none)
6266
FOUR_C_THROW("NitscheStrategyFpi: No frictional contact implemented for Nitsche FPSCI!");

‎src/contact/src/4C_contact_nitsche_strategy_fsi.hpp‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,9 @@ namespace CONTACT
4141
weighting_(
4242
Teuchos::getIntegralValue<CONTACT::NitscheWeighting>(params, "NITSCHE_WEIGHTING"))
4343
{
44+
init_time_step_size_and_total_time(params);
45+
46+
4447
if (Teuchos::getIntegralValue<CONTACT::FrictionType>(params, "FRICTION") !=
4548
CONTACT::FrictionType::none)
4649
FOUR_C_THROW("NitscheStrategyFsi: No frictional contact implemented for Nitsche FSCI!");
@@ -57,6 +60,9 @@ namespace CONTACT
5760
weighting_(
5861
Teuchos::getIntegralValue<CONTACT::NitscheWeighting>(params, "NITSCHE_WEIGHTING"))
5962
{
63+
init_time_step_size_and_total_time(params);
64+
65+
6066
if (Teuchos::getIntegralValue<CONTACT::FrictionType>(params, "FRICTION") !=
6167
CONTACT::FrictionType::none)
6268
FOUR_C_THROW("NitscheStrategyFsi: No frictional contact implemented for Nitsche FSCI!");

‎src/contact/src/4C_contact_nitsche_strategy_poro.hpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@ namespace CONTACT
3737
dof_row_map, NodeRowMap, params, std::move(interface), dim, comm, alphaf, maxdof),
3838
no_penetration_(params.get<bool>("CONTACT_NO_PENETRATION"))
3939
{
40+
init_time_step_size_and_total_time(params);
4041
}
4142

4243
//! Shared data constructor
@@ -48,6 +49,7 @@ namespace CONTACT
4849
comm, alphaf, maxdof),
4950
no_penetration_(params.get<bool>("CONTACT_NO_PENETRATION"))
5051
{
52+
init_time_step_size_and_total_time(params);
5153
}
5254

5355
void apply_force_stiff_cmt(std::shared_ptr<Core::LinAlg::Vector<double>> dis,

‎src/contact/src/4C_contact_nitsche_strategy_ssi.hpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -35,6 +35,7 @@ namespace CONTACT
3535
: NitscheStrategy(
3636
dofRowMap, nodeRowMap, params, std::move(interface), dim, comm, alphaf, maxdof)
3737
{
38+
init_time_step_size_and_total_time(params);
3839
}
3940

4041
//! Shared data constructor
@@ -46,6 +47,7 @@ namespace CONTACT
4647
: NitscheStrategy(data_ptr, dofRowMap, nodeRowMap, params, std::move(interface), dim, comm,
4748
alphaf, maxdof)
4849
{
50+
init_time_step_size_and_total_time(params);
4951
}
5052

5153
void apply_force_stiff_cmt(std::shared_ptr<Core::LinAlg::Vector<double>> dis,

‎src/contact/src/4C_contact_nitsche_strategy_ssi_elch.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,7 @@ namespace CONTACT
3434
: NitscheStrategySsi(data_ptr, dof_row_map, NodeRowMap, params, std::move(interface), dim,
3535
comm, alphaf, maxdof)
3636
{
37+
init_time_step_size_and_total_time(params);
3738
}
3839

3940
void apply_force_stiff_cmt(std::shared_ptr<Core::LinAlg::Vector<double>> dis,

‎src/contact/src/4C_contact_paramsinterface.hpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -69,6 +69,9 @@ namespace CONTACT
6969
//! get the current time step [derived]
7070
virtual double get_delta_time() const = 0;
7171

72+
//! get the total time
73+
virtual double get_total_time() const = 0;
74+
7275
//! get a pointer to the contact model evaluator
7376
virtual const Solid::ModelEvaluator::Generic& get_model_evaluator() const = 0;
7477

0 commit comments

Comments
 (0)