Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "ExplicitDynamicsContactAction.h"
11 :
12 : #include "Factory.h"
13 : #include "FEProblem.h"
14 : #include "Conversion.h"
15 : #include "AddVariableAction.h"
16 : #include "NonlinearSystemBase.h"
17 : #include "Parser.h"
18 :
19 : #include "NanoflannMeshAdaptor.h"
20 : #include "PointListAdaptor.h"
21 :
22 : #include <set>
23 : #include <algorithm>
24 : #include <unordered_map>
25 : #include <limits>
26 :
27 : #include "libmesh/petsc_nonlinear_solver.h"
28 : #include "libmesh/string_to_enum.h"
29 :
30 : // Counter for distinct contact action objects
31 : static unsigned int ed_contact_action_counter = 0;
32 :
33 : registerMooseAction("ContactApp", ExplicitDynamicsContactAction, "add_aux_variable");
34 : registerMooseAction("ContactApp", ExplicitDynamicsContactAction, "add_contact_aux_variable");
35 : registerMooseAction("ContactApp", ExplicitDynamicsContactAction, "add_aux_kernel");
36 :
37 : registerMooseAction("ContactApp", ExplicitDynamicsContactAction, "add_constraint");
38 : registerMooseAction("ContactApp", ExplicitDynamicsContactAction, "add_user_object");
39 :
40 : InputParameters
41 91 : ExplicitDynamicsContactAction::validParams()
42 : {
43 91 : InputParameters params = Action::validParams();
44 91 : params += ExplicitDynamicsContactAction::commonParameters();
45 :
46 182 : params.addParam<std::vector<BoundaryName>>(
47 : "primary", "The list of boundary IDs referring to primary sidesets");
48 182 : params.addParam<std::vector<BoundaryName>>(
49 : "secondary", "The list of boundary IDs referring to secondary sidesets");
50 182 : params.addParam<std::vector<VariableName>>(
51 : "displacements",
52 : "The displacements appropriate for the simulation geometry and coordinate system");
53 182 : params.addParam<Real>("friction_coefficient", 0, "The friction coefficient");
54 182 : params.addParam<MooseEnum>(
55 182 : "model", ExplicitDynamicsContactAction::getModelEnum(), "The contact model to use");
56 182 : params.addParam<Real>("tangential_tolerance",
57 : "Tangential distance to extend edges of contact surfaces");
58 182 : params.addParam<Real>("penalty", 1e8, "Penalty factor for normal contact.");
59 91 : params.addClassDescription("Sets up all objects needed for mechanical contact enforcement in "
60 : "explicit dynamics simulations.");
61 182 : params.addParam<std::vector<TagName>>(
62 : "extra_vector_tags",
63 : "The tag names for extra vectors that residual data should be saved into");
64 182 : params.addParam<std::vector<TagName>>(
65 : "absolute_value_vector_tags",
66 : "The tags for the vectors this residual object should fill with the "
67 : "absolute value of the residual contribution");
68 :
69 182 : params.addParam<VariableName>("secondary_gap_offset",
70 : "Offset to gap distance from secondary side");
71 182 : params.addParam<VariableName>("mapped_primary_gap_offset",
72 : "Offset to gap distance mapped from primary side");
73 182 : params.addParam<bool>("verbose", false, "Verbose output. May increase runtime.");
74 :
75 91 : return params;
76 0 : }
77 :
78 91 : ExplicitDynamicsContactAction::ExplicitDynamicsContactAction(const InputParameters & params)
79 : : Action(params),
80 182 : _boundary_pairs(getParam<BoundaryName, BoundaryName>("primary", "secondary")),
81 182 : _model(getParam<MooseEnum>("model").getEnum<ExplicitDynamicsContactModel>()),
82 273 : _verbose(getParam<bool>("verbose"))
83 : {
84 : // The resulting velocity of the contact algorithm is applied, as the code stands, by modifying
85 : // the old position. This causes artifacts in the internal forces as old position, new position,
86 : // and velocities are not fully consistent. Results improve with shorter time steps, but it would
87 : // be best to modify all the solution arrays for consistency or replace the way the explicit
88 : // system is solved to obtain accelerations (a = M^{-1}F) such that only velocities on the
89 : // contacting boundary need to be updated with the contact algorithm output.
90 91 : mooseWarning("Verification of explicit dynamics capabilities is an ongoing effort.");
91 91 : }
92 :
93 : void
94 455 : ExplicitDynamicsContactAction::act()
95 : {
96 : // proform problem checks/corrections once during the first feasible task
97 455 : if (_current_task == "add_contact_aux_variable")
98 : {
99 182 : if (!_problem->getDisplacedProblem())
100 0 : mooseError(
101 : "Contact requires updated coordinates. Use the 'displacements = ...' parameter in the "
102 : "Mesh block.");
103 :
104 : // It is risky to apply this optimization to contact problems
105 : // since the problem configuration may be changed during Jacobian
106 : // evaluation. We therefore turn it off for all contact problems so that
107 : // PETSc-3.8.4 or higher will have the same behavior as PETSc-3.8.3.
108 91 : if (!_problem->isSNESMFReuseBaseSetbyUser())
109 : _problem->setSNESMFReuseBase(false, false);
110 : }
111 :
112 455 : addNodeFaceContact();
113 :
114 455 : if (_current_task == "add_aux_kernel" && _verbose)
115 : { // Add ContactPenetrationAuxAction.
116 146 : if (!_problem->getDisplacedProblem())
117 0 : mooseError("Contact requires updated coordinates. Use the 'displacements = ...' line in the "
118 : "Mesh block.");
119 : unsigned int pair_number(0);
120 : // Create auxiliary kernels for each contact pairs
121 146 : for (const auto & contact_pair : _boundary_pairs)
122 : {
123 : {
124 73 : InputParameters params = _factory.getValidParams("PenetrationAux");
125 73 : params.applyParameters(parameters(),
126 : {"secondary_gap_offset", "mapped_primary_gap_offset", "order"});
127 :
128 : std::vector<VariableName> displacements =
129 146 : getParam<std::vector<VariableName>>("displacements");
130 73 : const auto order = _problem->systemBaseNonlinear(/*nl_sys_num=*/0)
131 73 : .system()
132 73 : .variable_type(displacements[0])
133 : .order.get_order();
134 :
135 146 : params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
136 292 : params.set<ExecFlagEnum>("execute_on") = {EXEC_INITIAL, EXEC_LINEAR};
137 219 : params.set<std::vector<BoundaryName>>("boundary") = {contact_pair.second};
138 146 : params.set<BoundaryName>("paired_boundary") = contact_pair.first;
139 146 : params.set<AuxVariableName>("variable") = "penetration";
140 146 : if (isParamValid("secondary_gap_offset"))
141 0 : params.set<std::vector<VariableName>>("secondary_gap_offset") = {
142 0 : getParam<VariableName>("secondary_gap_offset")};
143 146 : if (isParamValid("mapped_primary_gap_offset"))
144 0 : params.set<std::vector<VariableName>>("mapped_primary_gap_offset") = {
145 0 : getParam<VariableName>("mapped_primary_gap_offset")};
146 73 : params.set<bool>("use_displaced_mesh") = true;
147 146 : std::string name = _name + "_contact_" + Moose::stringify(pair_number);
148 73 : pair_number++;
149 146 : _problem->addAuxKernel("PenetrationAux", name, params);
150 73 : }
151 : }
152 :
153 73 : addContactPressureAuxKernel();
154 : }
155 :
156 455 : if (_current_task == "add_contact_aux_variable")
157 : {
158 182 : std::vector<VariableName> displacements = getParam<std::vector<VariableName>>("displacements");
159 91 : const auto order = _problem->systemBaseNonlinear(/*nl_sys_num=*/0)
160 91 : .system()
161 91 : .variable_type(displacements[0])
162 : .order.get_order();
163 :
164 : // Add gap rate for output
165 : {
166 91 : auto var_params = _factory.getValidParams("MooseVariable");
167 182 : var_params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
168 182 : var_params.set<MooseEnum>("family") = "LAGRANGE";
169 182 : _problem->addAuxVariable("MooseVariable", "gap_rate", var_params);
170 91 : }
171 91 : if (_verbose)
172 : {
173 : // Add ContactPenetrationVarAction
174 : {
175 73 : auto var_params = _factory.getValidParams("MooseVariable");
176 146 : var_params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
177 146 : var_params.set<MooseEnum>("family") = "LAGRANGE";
178 :
179 146 : _problem->addAuxVariable("MooseVariable", "penetration", var_params);
180 73 : }
181 : {
182 73 : auto var_params = _factory.getValidParams("MooseVariable");
183 146 : var_params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
184 146 : var_params.set<MooseEnum>("family") = "LAGRANGE";
185 :
186 146 : _problem->addAuxVariable("MooseVariable", "contact_pressure", var_params);
187 73 : }
188 : // Add nodal area contact variable
189 : {
190 73 : auto var_params = _factory.getValidParams("MooseVariable");
191 146 : var_params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
192 146 : var_params.set<MooseEnum>("family") = "LAGRANGE";
193 :
194 146 : _problem->addAuxVariable("MooseVariable", "nodal_area", var_params);
195 73 : }
196 : // Add nodal density variable
197 : {
198 73 : auto var_params = _factory.getValidParams("MooseVariable");
199 146 : var_params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
200 146 : var_params.set<MooseEnum>("family") = "LAGRANGE";
201 :
202 146 : _problem->addAuxVariable("MooseVariable", "nodal_density", var_params);
203 73 : }
204 : // Add nodal wave speed
205 : {
206 73 : auto var_params = _factory.getValidParams("MooseVariable");
207 146 : var_params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
208 146 : var_params.set<MooseEnum>("family") = "LAGRANGE";
209 :
210 146 : _problem->addAuxVariable("MooseVariable", "nodal_wave_speed", var_params);
211 73 : }
212 : }
213 91 : }
214 :
215 455 : if (_current_task == "add_user_object" && _verbose)
216 : {
217 : {
218 73 : auto var_params = _factory.getValidParams("NodalArea");
219 :
220 : // Get secondary_boundary_vector from possibly updated set from the
221 : // ContactAction constructor cleanup
222 73 : const auto actions = _awh.getActions<ExplicitDynamicsContactAction>();
223 :
224 : std::vector<BoundaryName> secondary_boundary_vector;
225 146 : for (const auto * const action : actions)
226 146 : for (const auto j : index_range(action->_boundary_pairs))
227 73 : secondary_boundary_vector.push_back(action->_boundary_pairs[j].second);
228 :
229 146 : var_params.set<std::vector<BoundaryName>>("boundary") = secondary_boundary_vector;
230 219 : var_params.set<std::vector<VariableName>>("variable") = {"nodal_area"};
231 :
232 : mooseAssert(_problem, "Problem pointer is NULL");
233 292 : var_params.set<ExecFlagEnum>("execute_on", true) = {EXEC_INITIAL, EXEC_TIMESTEP_BEGIN};
234 73 : var_params.set<bool>("use_displaced_mesh") = true;
235 :
236 219 : _problem->addUserObject(
237 73 : "NodalArea", "nodal_area_object_" + Moose::stringify(name()), var_params);
238 73 : }
239 : // Add nodal density and nodal wave speed user
240 : {
241 : // Add nodal density
242 73 : auto var_params = _factory.getValidParams("NodalDensity");
243 :
244 : // Get secondary_boundary_vector from possibly updated set from the
245 : // ContactAction constructor cleanup
246 73 : const auto actions = _awh.getActions<ExplicitDynamicsContactAction>();
247 :
248 : std::vector<BoundaryName> secondary_boundary_vector;
249 146 : for (const auto * const action : actions)
250 146 : for (const auto j : index_range(action->_boundary_pairs))
251 73 : secondary_boundary_vector.push_back(action->_boundary_pairs[j].second);
252 :
253 146 : var_params.set<std::vector<BoundaryName>>("boundary") = secondary_boundary_vector;
254 219 : var_params.set<std::vector<VariableName>>("variable") = {"nodal_density"};
255 :
256 : mooseAssert(_problem, "Problem pointer is NULL");
257 292 : var_params.set<ExecFlagEnum>("execute_on", true) = {EXEC_INITIAL, EXEC_TIMESTEP_BEGIN};
258 73 : var_params.set<bool>("use_displaced_mesh") = true;
259 :
260 219 : _problem->addUserObject(
261 73 : "NodalDensity", "nodal_density_object_" + Moose::stringify(name()), var_params);
262 73 : }
263 : {
264 : // Add wave speed
265 73 : auto var_params = _factory.getValidParams("NodalWaveSpeed");
266 :
267 : // Get secondary_boundary_vector from possibly updated set from the
268 : // ContactAction constructor cleanup
269 73 : const auto actions = _awh.getActions<ExplicitDynamicsContactAction>();
270 :
271 : std::vector<BoundaryName> secondary_boundary_vector;
272 146 : for (const auto * const action : actions)
273 146 : for (const auto j : index_range(action->_boundary_pairs))
274 73 : secondary_boundary_vector.push_back(action->_boundary_pairs[j].second);
275 :
276 146 : var_params.set<std::vector<BoundaryName>>("boundary") = secondary_boundary_vector;
277 219 : var_params.set<std::vector<VariableName>>("variable") = {"nodal_wave_speed"};
278 :
279 : mooseAssert(_problem, "Problem pointer is NULL");
280 292 : var_params.set<ExecFlagEnum>("execute_on", true) = {EXEC_INITIAL, EXEC_TIMESTEP_BEGIN};
281 73 : var_params.set<bool>("use_displaced_mesh") = true;
282 :
283 219 : _problem->addUserObject(
284 73 : "NodalWaveSpeed", "nodal_wavespeed_object_" + Moose::stringify(name()), var_params);
285 73 : }
286 : }
287 747 : }
288 :
289 : void
290 73 : ExplicitDynamicsContactAction::addContactPressureAuxKernel()
291 : {
292 : // Add ContactPressureAux: Only one object for all contact pairs
293 73 : const auto actions = _awh.getActions<ExplicitDynamicsContactAction>();
294 :
295 : // Increment counter for contact action objects
296 73 : ed_contact_action_counter++;
297 :
298 : // Add auxiliary kernel if we are the last contact action object.
299 73 : if (ed_contact_action_counter == actions.size() && _verbose)
300 : {
301 : std::vector<BoundaryName> boundary_vector;
302 : std::vector<BoundaryName> pair_boundary_vector;
303 :
304 146 : for (const auto * const action : actions)
305 146 : for (const auto j : index_range(action->_boundary_pairs))
306 : {
307 73 : boundary_vector.push_back(action->_boundary_pairs[j].second);
308 73 : pair_boundary_vector.push_back(action->_boundary_pairs[j].first);
309 : }
310 :
311 73 : InputParameters params = _factory.getValidParams("ContactPressureAux");
312 73 : params.applyParameters(parameters(), {"order"});
313 :
314 146 : std::vector<VariableName> displacements = getParam<std::vector<VariableName>>("displacements");
315 73 : const auto order = _problem->systemBaseNonlinear(/*nl_sys_num=*/0)
316 73 : .system()
317 73 : .variable_type(displacements[0])
318 : .order.get_order();
319 :
320 146 : params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
321 73 : params.set<std::vector<BoundaryName>>("boundary") = boundary_vector;
322 146 : params.set<std::vector<BoundaryName>>("paired_boundary") = pair_boundary_vector;
323 146 : params.set<AuxVariableName>("variable") = "contact_pressure";
324 146 : params.addRequiredCoupledVar("nodal_area", "The nodal area");
325 219 : params.set<std::vector<VariableName>>("nodal_area") = {"nodal_area"};
326 73 : params.set<bool>("use_displaced_mesh") = true;
327 :
328 73 : std::string name = _name + "_contact_pressure";
329 219 : params.set<ExecFlagEnum>("execute_on", true) = {EXEC_TIMESTEP_END};
330 146 : _problem->addAuxKernel("ContactPressureAux", name, params);
331 73 : }
332 146 : }
333 :
334 : void
335 455 : ExplicitDynamicsContactAction::addNodeFaceContact()
336 : {
337 455 : if (_current_task != "add_constraint")
338 364 : return;
339 :
340 91 : std::string action_name = MooseUtils::shortName(name());
341 273 : std::vector<VariableName> displacements = getParam<std::vector<VariableName>>("displacements");
342 91 : const unsigned int ndisp = displacements.size();
343 :
344 : std::string constraint_type;
345 :
346 : constraint_type = "ExplicitDynamicsContactConstraint";
347 :
348 91 : InputParameters params = _factory.getValidParams(constraint_type);
349 :
350 91 : params.applyParameters(parameters(),
351 : {"displacements",
352 : "secondary_gap_offset",
353 : "mapped_primary_gap_offset",
354 : "primary",
355 : "secondary"});
356 :
357 91 : const auto order = _problem->systemBaseNonlinear(/*nl_sys_num=*/0)
358 91 : .system()
359 91 : .variable_type(displacements[0])
360 : .order.get_order();
361 :
362 91 : params.set<std::vector<VariableName>>("displacements") = displacements;
363 91 : params.set<bool>("use_displaced_mesh") = true;
364 182 : params.set<Real>("penalty") = getParam<Real>("penalty");
365 :
366 182 : params.set<MooseEnum>("order") = Utility::enum_to_string<Order>(OrderWrapper{order});
367 :
368 182 : for (const auto & contact_pair : _boundary_pairs)
369 : {
370 :
371 182 : if (isParamValid("vel_x"))
372 216 : params.set<std::vector<VariableName>>("vel_x") = getParam<std::vector<VariableName>>("vel_x");
373 182 : if (isParamValid("vel_y"))
374 216 : params.set<std::vector<VariableName>>("vel_y") = getParam<std::vector<VariableName>>("vel_y");
375 182 : if (isParamValid("vel_z"))
376 216 : params.set<std::vector<VariableName>>("vel_z") = getParam<std::vector<VariableName>>("vel_z");
377 :
378 273 : params.set<std::vector<VariableName>>("gap_rate") = {"gap_rate"};
379 :
380 91 : params.set<BoundaryName>("boundary") = contact_pair.first;
381 182 : if (isParamValid("secondary_gap_offset"))
382 0 : params.set<std::vector<VariableName>>("secondary_gap_offset") = {
383 0 : getParam<VariableName>("secondary_gap_offset")};
384 182 : if (isParamValid("mapped_primary_gap_offset"))
385 0 : params.set<std::vector<VariableName>>("mapped_primary_gap_offset") = {
386 0 : getParam<VariableName>("mapped_primary_gap_offset")};
387 :
388 364 : for (unsigned int i = 0; i < ndisp; ++i)
389 : {
390 819 : std::string name = action_name + "_constraint_" + Moose::stringify(contact_pair, "_") + "_" +
391 273 : Moose::stringify(i);
392 :
393 273 : params.set<unsigned int>("component") = i;
394 :
395 273 : params.set<BoundaryName>("primary") = contact_pair.first;
396 273 : params.set<BoundaryName>("secondary") = contact_pair.second;
397 546 : params.set<NonlinearVariableName>("variable") = displacements[i];
398 819 : params.set<std::vector<VariableName>>("primary_variable") = {displacements[i]};
399 273 : params.applySpecificParameters(parameters(),
400 : {"extra_vector_tags", "absolute_value_vector_tags"});
401 273 : _problem->addConstraint(constraint_type, name, params);
402 : }
403 : }
404 182 : }
405 :
406 : MooseEnum
407 546 : ExplicitDynamicsContactAction::getModelEnum()
408 : {
409 1092 : return MooseEnum("frictionless frictionless_balance", "frictionless");
410 : }
411 :
412 : InputParameters
413 455 : ExplicitDynamicsContactAction::commonParameters()
414 : {
415 455 : InputParameters params = emptyInputParameters();
416 :
417 910 : params.addParam<MooseEnum>(
418 910 : "model", ExplicitDynamicsContactAction::getModelEnum(), "The contact model to use");
419 :
420 : // Gap rate input
421 910 : params.addCoupledVar("vel_x", "x-component of velocity.");
422 910 : params.addCoupledVar("vel_y", "y-component of velocity.");
423 910 : params.addCoupledVar("vel_z", "z-component of velocity.");
424 : // Gap rate output
425 910 : params.addCoupledVar("gap_rate", "Gap rate for output.");
426 :
427 455 : return params;
428 0 : }
|
