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 "Conversion.h"
11 : #include "FEProblem.h"
12 : #include "Factory.h"
13 : #include "MooseMesh.h"
14 : #include "MooseObjectAction.h"
15 : #include "QuasiStaticSolidMechanicsPhysics.h"
16 : #include "Material.h"
17 : #include "CommonSolidMechanicsAction.h"
18 :
19 : #include "BlockRestrictable.h"
20 :
21 : #include "HomogenizationInterface.h"
22 : #include "AddVariableAction.h"
23 :
24 : #include "libmesh/string_to_enum.h"
25 : #include <algorithm>
26 :
27 : registerMooseAction("SolidMechanicsApp", QuasiStaticSolidMechanicsPhysics, "meta_action");
28 :
29 : registerMooseAction("SolidMechanicsApp", QuasiStaticSolidMechanicsPhysics, "setup_mesh_complete");
30 :
31 : registerMooseAction("SolidMechanicsApp",
32 : QuasiStaticSolidMechanicsPhysics,
33 : "validate_coordinate_systems");
34 :
35 : registerMooseAction("SolidMechanicsApp", QuasiStaticSolidMechanicsPhysics, "add_variable");
36 :
37 : registerMooseAction("SolidMechanicsApp", QuasiStaticSolidMechanicsPhysics, "add_aux_variable");
38 :
39 : registerMooseAction("SolidMechanicsApp", QuasiStaticSolidMechanicsPhysics, "add_kernel");
40 :
41 : registerMooseAction("SolidMechanicsApp", QuasiStaticSolidMechanicsPhysics, "add_aux_kernel");
42 :
43 : registerMooseAction("SolidMechanicsApp", QuasiStaticSolidMechanicsPhysics, "add_material");
44 :
45 : registerMooseAction("SolidMechanicsApp",
46 : QuasiStaticSolidMechanicsPhysics,
47 : "add_master_action_material");
48 :
49 : registerMooseAction("SolidMechanicsApp", QuasiStaticSolidMechanicsPhysics, "add_scalar_kernel");
50 :
51 : registerMooseAction("SolidMechanicsApp", QuasiStaticSolidMechanicsPhysics, "add_user_object");
52 :
53 : InputParameters
54 4511 : QuasiStaticSolidMechanicsPhysics::validParams()
55 : {
56 4511 : InputParameters params = QuasiStaticSolidMechanicsPhysicsBase::validParams();
57 4511 : params.addClassDescription("Set up stress divergence kernels with coordinate system aware logic");
58 :
59 : // parameters specified here only appear in the input file sub-blocks of the solid mechanics
60 : // action, not in the common parameters area
61 9022 : params.addParam<std::vector<SubdomainName>>("block",
62 : {},
63 : "The list of ids of the blocks (subdomain) "
64 : "that the stress divergence kernels will be "
65 : "applied to");
66 9022 : params.addParamNamesToGroup("block", "Advanced");
67 :
68 9022 : params.addParam<MultiMooseEnum>("additional_generate_output",
69 9022 : QuasiStaticSolidMechanicsPhysicsBase::outputPropertiesType(),
70 : "Add scalar quantity output for stress and/or strain (will be "
71 : "appended to the list in `generate_output`)");
72 9022 : params.addParam<MultiMooseEnum>(
73 : "additional_material_output_order",
74 9022 : QuasiStaticSolidMechanicsPhysicsBase::materialOutputOrders(),
75 : "Specifies the order of the FE shape function to use for this variable.");
76 :
77 9022 : params.addParam<MultiMooseEnum>(
78 : "additional_material_output_family",
79 9022 : QuasiStaticSolidMechanicsPhysicsBase::materialOutputFamilies(),
80 : "Specifies the family of FE shape functions to use for this variable.");
81 :
82 9022 : params.addParamNamesToGroup("additional_generate_output additional_material_output_order "
83 : "additional_material_output_family",
84 : "Output");
85 9022 : params.addParam<std::string>(
86 : "strain_base_name",
87 : "The base name used for the strain. If not provided, it will be set equal to base_name");
88 9022 : params.addParam<std::vector<TagName>>(
89 : "extra_vector_tags",
90 : "The tag names for extra vectors that residual data should be saved into");
91 9022 : params.addParam<std::vector<TagName>>("absolute_value_vector_tags",
92 : "The tag names for extra vectors that the absolute value "
93 : "of the residual should be accumulated into");
94 9022 : params.addParam<Real>("scaling", "The scaling to apply to the displacement variables");
95 9022 : params.addParam<Point>(
96 : "cylindrical_axis_point1",
97 : "Starting point for direction of axis of rotation for cylindrical stress/strain.");
98 9022 : params.addParam<Point>(
99 : "cylindrical_axis_point2",
100 : "Ending point for direction of axis of rotation for cylindrical stress/strain.");
101 9022 : params.addParam<Point>("spherical_center_point",
102 : "Center point of the spherical coordinate system.");
103 9022 : params.addParam<Point>("direction", "Direction stress/strain is calculated in");
104 9022 : params.addParam<bool>("automatic_eigenstrain_names",
105 9022 : false,
106 : "Collects all material eigenstrains and passes to required strain "
107 : "calculator within TMA internally.");
108 :
109 : // Homogenization system input
110 9022 : params.addParam<MultiMooseEnum>(
111 : "constraint_types",
112 : Homogenization::constraintType,
113 : "Type of each constraint: strain, stress, or none. The types are specified in the "
114 : "column-major order, and there must be 9 entries in total.");
115 9022 : params.addParam<std::vector<FunctionName>>(
116 : "targets", {}, "Functions giving the targets to hit for constraint types that are not none.");
117 9022 : params.addParam<bool>(
118 : "homogenized_off_diagonal_jacobian",
119 9022 : true,
120 : "Whether to include the off-diagonal scalar contributions to the homogenized jacobian");
121 :
122 9022 : params.addParamNamesToGroup("scaling", "Variables");
123 9022 : params.addParamNamesToGroup("strain_base_name automatic_eigenstrain_names", "Strain");
124 9022 : params.addParamNamesToGroup(
125 : "cylindrical_axis_point1 cylindrical_axis_point2 spherical_center_point direction",
126 : "Coordinate system");
127 9022 : params.addParamNamesToGroup("constraint_types targets homogenized_off_diagonal_jacobian",
128 : "Homogenization");
129 :
130 4511 : return params;
131 0 : }
132 :
133 4511 : QuasiStaticSolidMechanicsPhysics::QuasiStaticSolidMechanicsPhysics(const InputParameters & params)
134 : : QuasiStaticSolidMechanicsPhysicsBase(params),
135 9022 : _displacements(getParam<std::vector<VariableName>>("displacements")),
136 4511 : _ndisp(_displacements.size()),
137 4511 : _coupled_displacements(_ndisp),
138 9022 : _save_in(getParam<std::vector<AuxVariableName>>("save_in")),
139 9022 : _diag_save_in(getParam<std::vector<AuxVariableName>>("diag_save_in")),
140 13533 : _subdomain_names(getParam<std::vector<SubdomainName>>("block")),
141 : _subdomain_ids(),
142 9022 : _strain(getParam<MooseEnum>("strain").getEnum<Strain>()),
143 9022 : _planar_formulation(getParam<MooseEnum>("planar_formulation").getEnum<PlanarFormulation>()),
144 4511 : _out_of_plane_direction(
145 4511 : getParam<MooseEnum>("out_of_plane_direction").getEnum<OutOfPlaneDirection>()),
146 9106 : _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
147 9022 : _material_output_order(getParam<MultiMooseEnum>("material_output_order")),
148 9022 : _material_output_family(getParam<MultiMooseEnum>("material_output_family")),
149 4511 : _cylindrical_axis_point1_valid(params.isParamSetByUser("cylindrical_axis_point1")),
150 4511 : _cylindrical_axis_point2_valid(params.isParamSetByUser("cylindrical_axis_point2")),
151 4511 : _direction_valid(params.isParamSetByUser("direction")),
152 9022 : _verbose(getParam<bool>("verbose")),
153 4511 : _spherical_center_point_valid(params.isParamSetByUser("spherical_center_point")),
154 9022 : _auto_eigenstrain(getParam<bool>("automatic_eigenstrain_names")),
155 9022 : _compatibility_mode(getParam<bool>("compatibility_mode")),
156 9022 : _lagrangian_kernels(getParam<bool>("new_system") || _compatibility_mode),
157 4511 : _lk_large_kinematics(_strain == Strain::Finite),
158 : // compatibility_mode mandates the total Lagrangian kernel; otherwise honor the user choice.
159 9022 : _lk_formulation(_compatibility_mode
160 4511 : ? LKFormulation::Total
161 9010 : : getParam<MooseEnum>("formulation").getEnum<LKFormulation>()),
162 9022 : _lk_locking(getParam<bool>("volumetric_locking_correction")),
163 4511 : _lk_homogenization(false),
164 4511 : _constraint_types(getParam<MultiMooseEnum>("constraint_types")),
165 9022 : _targets(getParam<std::vector<FunctionName>>("targets")),
166 18044 : _lk_h_off_jac(getParam<bool>("homogenized_off_diagonal_jacobian"))
167 : {
168 : // determine if incremental strains are to be used
169 9022 : if (isParamValid("incremental"))
170 : {
171 2478 : const bool incremental = getParam<bool>("incremental");
172 1239 : if (!incremental && _strain == Strain::Small)
173 61 : _strain_and_increment = StrainAndIncrement::SmallTotal;
174 0 : else if (!incremental && _strain == Strain::Finite)
175 0 : _strain_and_increment = StrainAndIncrement::FiniteTotal;
176 1178 : else if (incremental && _strain == Strain::Small)
177 570 : _strain_and_increment = StrainAndIncrement::SmallIncremental;
178 608 : else if (incremental && _strain == Strain::Finite)
179 608 : _strain_and_increment = StrainAndIncrement::FiniteIncremental;
180 : else
181 0 : mooseError("Internal error");
182 : }
183 : else
184 : {
185 3272 : if (_strain == Strain::Small)
186 : {
187 1842 : _strain_and_increment = StrainAndIncrement::SmallTotal;
188 1842 : mooseInfo("SolidMechanics Action: selecting 'total small strain' formulation. Use "
189 : "`incremental = true` to select 'incremental small strain' instead.");
190 : }
191 1430 : else if (_strain == Strain::Finite)
192 : {
193 1430 : _strain_and_increment = StrainAndIncrement::FiniteIncremental;
194 1430 : mooseInfo("SolidMechanics Action: selecting 'incremental finite strain' formulation.");
195 : }
196 : else
197 0 : mooseError("Internal error");
198 : }
199 :
200 : // determine if displaced mesh is to be used
201 4511 : _use_displaced_mesh = (_strain == Strain::Finite);
202 9022 : if (params.isParamSetByUser("use_displaced_mesh") && !_lagrangian_kernels)
203 : {
204 516 : bool use_displaced_mesh_param = getParam<bool>("use_displaced_mesh");
205 481 : if (use_displaced_mesh_param != _use_displaced_mesh && params.isParamSetByUser("strain"))
206 0 : mooseError("Wrong combination of use displaced mesh and strain model");
207 258 : _use_displaced_mesh = use_displaced_mesh_param;
208 : }
209 :
210 : // convert vector of VariableName to vector of VariableName
211 16741 : for (unsigned int i = 0; i < _ndisp; ++i)
212 12230 : _coupled_displacements[i] = _displacements[i];
213 :
214 4511 : if (_save_in.size() != 0 && _save_in.size() != _ndisp)
215 0 : mooseError("Number of save_in variables should equal to the number of displacement variables ",
216 0 : _ndisp);
217 :
218 4511 : if (_diag_save_in.size() != 0 && _diag_save_in.size() != _ndisp)
219 0 : mooseError(
220 : "Number of diag_save_in variables should equal to the number of displacement variables ",
221 0 : _ndisp);
222 :
223 : // plane strain consistency check
224 4511 : if (_planar_formulation != PlanarFormulation::None)
225 : {
226 368 : if (_lagrangian_kernels)
227 0 : mooseDocumentedError(
228 : "moose",
229 : 27340,
230 : "Planar formulations are not yet available through the Physics syntax with new_system = "
231 : "true. They can be enabled by manually setting up the appropriate objects. Please refer "
232 : "to the documentation and regression tests for examples.");
233 736 : if (params.isParamSetByUser("out_of_plane_strain") &&
234 48 : _planar_formulation != PlanarFormulation::WeakPlaneStress)
235 0 : mooseError(
236 : "out_of_plane_strain should only be specified with planar_formulation=WEAK_PLANE_STRESS");
237 736 : else if (!params.isParamSetByUser("out_of_plane_strain") &&
238 320 : _planar_formulation == PlanarFormulation::WeakPlaneStress)
239 0 : mooseError("out_of_plane_strain must be specified with planar_formulation=WEAK_PLANE_STRESS");
240 : }
241 :
242 : // convert output variable names to lower case
243 19647 : for (const auto & out : getParam<MultiMooseEnum>("generate_output"))
244 : {
245 : std::string lower(out);
246 : std::transform(lower.begin(), lower.end(), lower.begin(), ::tolower);
247 10625 : _generate_output.push_back(lower);
248 : }
249 :
250 4511 : if (!_generate_output.empty())
251 2282 : verifyOrderAndFamilyOutputs();
252 :
253 : // Error if volumetric locking correction is true for 1D problems
254 4782 : if (_ndisp == 1 && getParam<bool>("volumetric_locking_correction"))
255 0 : mooseError("Volumetric locking correction should be set to false for 1D problems.");
256 :
257 10490 : if (!getParam<bool>("add_variables") && params.isParamSetByUser("scaling"))
258 0 : paramError("scaling",
259 : "The scaling parameter has no effect unless add_variables is set to true. Did you "
260 : "mean to set 'add_variables = true'?");
261 :
262 : // Get cylindrical axis points if set by user
263 4510 : if (_cylindrical_axis_point1_valid && _cylindrical_axis_point2_valid)
264 : {
265 24 : _cylindrical_axis_point1 = getParam<Point>("cylindrical_axis_point1");
266 24 : _cylindrical_axis_point2 = getParam<Point>("cylindrical_axis_point2");
267 : }
268 :
269 : // Get spherical center point if set by user
270 4510 : if (_spherical_center_point_valid)
271 66 : _spherical_center_point = getParam<Point>("spherical_center_point");
272 :
273 : // Get direction for tensor component if set by user
274 4510 : if (_direction_valid)
275 34 : _direction = getParam<Point>("direction");
276 :
277 : // Get eigenstrain names if passed by user
278 9020 : _eigenstrain_names = getParam<std::vector<MaterialPropertyName>>("eigenstrain_names");
279 :
280 : // Determine if we're going to use the homogenization system for the new
281 : // lagrangian kernels
282 4510 : bool ctype_set = params.isParamSetByUser("constraint_types");
283 4510 : bool targets_set = params.isParamSetByUser("targets");
284 4510 : if (ctype_set || targets_set)
285 : {
286 8 : if (!(ctype_set && targets_set))
287 0 : mooseError("To use the Lagrangian kernel homogenization system you "
288 : "most provide both the constraint_types and the targets "
289 : "parameters");
290 8 : _lk_homogenization = true;
291 : // Do consistency checking on the kernel options
292 8 : if (_lk_homogenization && (_lk_formulation == LKFormulation::Updated))
293 0 : mooseError("The Lagrangian kernel homogenization system requires the "
294 : "use of formulation = TOTAL");
295 : }
296 :
297 : // The action drives two mutually exclusive kernel systems: the legacy
298 : // StressDivergenceTensors + ComputeFiniteStrain system, and the new Lagrangian kernel system
299 : // (new_system / compatibility_mode). Many parameters belong exclusively to one system; setting
300 : // one that does not match the active system is almost always a copy-paste error, so reject it
301 : // explicitly with a clear message instead of silently ignoring it. The lists below are the
302 : // single place to register a parameter as system-specific.
303 : static const std::vector<std::string> new_system_only = {
304 4510 : "formulation", "constraint_types", "targets", "homogenized_off_diagonal_jacobian"};
305 : static const std::vector<std::string> legacy_only = {"use_finite_deform_jacobian",
306 4510 : "global_strain"};
307 :
308 4510 : if (_lagrangian_kernels)
309 : {
310 818 : for (const auto & p : legacy_only)
311 546 : if (params.isParamSetByUser(p))
312 1 : paramError(p,
313 : "'",
314 : p,
315 : "' belongs to the legacy StressDivergenceTensors kernel system and cannot be "
316 : "combined with new_system = true or compatibility_mode = true. The Lagrangian "
317 : "kernel system produces the exact Jacobian natively; homogenization replaces "
318 : "global_strain.");
319 :
320 : // use_automatic_differentiation is a value (defaults false); reject only when requested.
321 272 : if (_use_ad)
322 0 : paramError("use_automatic_differentiation",
323 : "The Lagrangian kernel system is not yet compatible with automatic "
324 : "differentiation.");
325 :
326 : // decomposition_method is legal only in compatibility_mode, where it maps to the strain
327 : // calculator's kinematic_approximation; plain new_system exposes those options directly on the
328 : // ComputeLagrangianStrain material via `kinematic_approximation`.
329 544 : if (params.isParamSetByUser("decomposition_method") && !_compatibility_mode)
330 0 : paramError("decomposition_method",
331 : "decomposition_method is not used by the Lagrangian kernel system. Set the "
332 : "equivalent option on the ComputeLagrangianStrain material via "
333 : "kinematic_approximation instead.");
334 :
335 : // compatibility_mode-specific consistency checks.
336 272 : if (_compatibility_mode)
337 : {
338 24 : if (params.isParamSetByUser("formulation") &&
339 12 : getParam<MooseEnum>("formulation").getEnum<LKFormulation>() != LKFormulation::Total)
340 0 : paramError("formulation",
341 : "compatibility_mode requires formulation = TOTAL (the OLD kernel system is "
342 : "always total Lagrangian on the displaced mesh).");
343 :
344 20 : if (_lk_large_kinematics && params.isParamSetByUser("kinematic_approximation"))
345 0 : paramError("kinematic_approximation",
346 : "compatibility_mode derives kinematic_approximation from decomposition_method; "
347 : "do not set it explicitly.");
348 12 : if (_lk_locking && _lk_large_kinematics &&
349 20 : params.isParamSetByUser("volumetric_locking_correction_mode"))
350 0 : paramError("volumetric_locking_correction_mode",
351 : "compatibility_mode derives the F-bar mode from volumetric_locking_correction; "
352 : "do not set volumetric_locking_correction_mode explicitly.");
353 24 : if (params.isParamSetByUser("generalized_midpoint_alpha"))
354 0 : paramError("generalized_midpoint_alpha",
355 : "compatibility_mode reproduces the legacy backward-Euler update; do not set "
356 : "generalized_midpoint_alpha explicitly.");
357 :
358 : // Validate decomposition_method <-> kinematic_approximation mapping. The NEW strain
359 : // calculator's `rashid_approximate` and `rashid_eigen` reproduce OLD's `TaylorExpansion`
360 : // and `EigenSolution` exactly; OLD's `HughesWinget` has no NEW counterpart.
361 12 : if (_lk_large_kinematics)
362 : {
363 16 : const std::string decomp = getParam<MooseEnum>("decomposition_method");
364 8 : if (decomp != "TaylorExpansion" && decomp != "EigenSolution")
365 0 : paramError("decomposition_method",
366 : "compatibility_mode + strain = FINITE supports only "
367 : "decomposition_method = TaylorExpansion or EigenSolution; '",
368 : decomp,
369 : "' has no Lagrangian kinematic_approximation counterpart.");
370 : }
371 : }
372 : }
373 : else
374 : {
375 21184 : for (const auto & p : new_system_only)
376 16948 : if (params.isParamSetByUser(p))
377 1 : paramError(p,
378 : "'",
379 : p,
380 : "' belongs to the new Lagrangian kernel system. Set new_system = true (or "
381 : "compatibility_mode = true) to use it.");
382 : }
383 4508 : }
384 :
385 : void
386 38175 : QuasiStaticSolidMechanicsPhysics::act()
387 : {
388 38175 : std::string ad_prepend = "";
389 38175 : if (_use_ad)
390 : ad_prepend = "AD";
391 :
392 : // Consistency checks across subdomains
393 38175 : actSubdomainChecks();
394 :
395 : // Gather info from all other QuasiStaticSolidMechanicsPhysics
396 38173 : actGatherActionParameters();
397 :
398 : // Deal with the optional AuxVariable based tensor quantity output
399 38169 : actOutputGeneration();
400 :
401 : // Meta action which optionally spawns other actions
402 38169 : if (_current_task == "meta_action")
403 : {
404 3283 : if (_planar_formulation == PlanarFormulation::GeneralizedPlaneStrain)
405 : {
406 89 : if (_use_ad)
407 0 : paramError("use_automatic_differentiation", "AD not setup for use with PlaneStrain");
408 : // Set the action parameters
409 89 : const std::string type = "GeneralizedPlaneStrainAction";
410 89 : auto action_params = _action_factory.getValidParams(type);
411 89 : action_params.set<bool>("_built_by_moose") = true;
412 89 : action_params.set<std::string>("registered_identifier") = "(AutoBuilt)";
413 :
414 : // Skipping selected parameters in applyParameters() and then manually setting them only if
415 : // they are set by the user is just to prevent both the current and deprecated variants of
416 : // these parameters from both getting passed to the UserObject. Once we get rid of the
417 : // deprecated versions, we can just set them all with applyParameters().
418 89 : action_params.applyParameters(parameters(),
419 : {"use_displaced_mesh",
420 : "out_of_plane_pressure",
421 : "out_of_plane_pressure_function",
422 : "factor",
423 : "pressure_factor"});
424 89 : action_params.set<bool>("use_displaced_mesh") = _use_displaced_mesh;
425 :
426 178 : if (parameters().isParamSetByUser("out_of_plane_pressure"))
427 0 : action_params.set<FunctionName>("out_of_plane_pressure") =
428 0 : getParam<FunctionName>("out_of_plane_pressure");
429 178 : if (parameters().isParamSetByUser("out_of_plane_pressure_function"))
430 0 : action_params.set<FunctionName>("out_of_plane_pressure_function") =
431 0 : getParam<FunctionName>("out_of_plane_pressure_function");
432 178 : if (parameters().isParamSetByUser("factor"))
433 0 : action_params.set<Real>("factor") = getParam<Real>("factor");
434 178 : if (parameters().isParamSetByUser("pressure_factor"))
435 0 : action_params.set<Real>("pressure_factor") = getParam<Real>("pressure_factor");
436 :
437 : // Create and add the action to the warehouse
438 : auto action = MooseSharedNamespace::static_pointer_cast<MooseObjectAction>(
439 89 : _action_factory.create(type, name() + "_gps", action_params));
440 267 : _awh.addActionBlock(action);
441 89 : }
442 : }
443 :
444 : // Add variables
445 34886 : else if (_current_task == "add_variable")
446 : {
447 : // Add variables here only if the CommonSolidMechanicsAction does not exist.
448 : // This happens notably if the QuasiStaticSolidMechanics was created by a meta_action
449 3271 : const auto common_actions = _awh.getActions<CommonSolidMechanicsAction>();
450 3271 : if (common_actions.empty() && getParam<bool>("add_variables"))
451 : {
452 0 : auto params = _factory.getValidParams("MooseVariable");
453 : // determine necessary order
454 0 : const bool second = _problem->mesh().hasSecondOrderElements();
455 :
456 0 : params.set<MooseEnum>("order") = second ? "SECOND" : "FIRST";
457 0 : params.set<MooseEnum>("family") = "LAGRANGE";
458 0 : if (isParamValid("scaling"))
459 0 : params.set<std::vector<Real>>("scaling") = {getParam<Real>("scaling")};
460 :
461 : // Note how we do not add the block restriction because BISON's meta-actions
462 : // currently rely on them not being added.
463 :
464 : // Loop through the displacement variables
465 0 : for (const auto & disp : _displacements)
466 : {
467 : // Create displacement variables
468 0 : _problem->addVariable("MooseVariable", disp, params);
469 : }
470 0 : }
471 :
472 : // Homogenization scalar
473 3271 : if (_lk_homogenization)
474 : {
475 8 : InputParameters params = _factory.getValidParams("MooseVariable");
476 : const std::map<bool, std::vector<unsigned int>> mg_order_max{{true, {1, 4, 9}},
477 24 : {false, {1, 3, 6}}};
478 : std::size_t mg_order = 0;
479 80 : for (auto i : index_range(_constraint_types))
480 : {
481 72 : const auto ctype = static_cast<Homogenization::ConstraintType>(_constraint_types.get(i));
482 72 : if (ctype != Homogenization::ConstraintType::None)
483 40 : mg_order++;
484 : }
485 8 : if (mg_order > mg_order_max.at(_lk_large_kinematics)[_ndisp - 1])
486 0 : paramError("constraint_types",
487 : "Number of non-none constraint types must not be greater than ",
488 0 : mg_order_max.at(_lk_large_kinematics)[_ndisp - 1],
489 : ", but ",
490 : mg_order,
491 : " are provided.");
492 16 : params.set<MooseEnum>("family") = "SCALAR";
493 8 : params.set<MooseEnum>("order") = mg_order;
494 8 : auto fe_type = AddVariableAction::feType(params);
495 8 : auto var_type = AddVariableAction::variableType(fe_type);
496 8 : _problem->addVariable(var_type, _hname, params);
497 8 : }
498 3271 : }
499 : // Add Materials
500 31615 : else if (_current_task == "add_master_action_material")
501 : {
502 : // Automatic eigenstrain names
503 3234 : if (_auto_eigenstrain)
504 31 : actEigenstrainNames();
505 :
506 : // Easiest just to branch on type here, as the strain systems are completely
507 : // different
508 3234 : if (_lagrangian_kernels)
509 272 : actLagrangianKernelStrain();
510 : else
511 2962 : actStressDivergenceTensorsStrain();
512 : }
513 :
514 : // Add Stress Divergence (and optionally WeakPlaneStress) Kernels
515 28381 : else if (_current_task == "add_kernel")
516 : {
517 16466 : for (unsigned int i = 0; i < _ndisp; ++i)
518 : {
519 12028 : auto tensor_kernel_type = getKernelType();
520 24056 : auto params = getKernelParameters(ad_prepend + tensor_kernel_type);
521 :
522 24056 : std::string kernel_name = "TM_" + name() + Moose::stringify(i);
523 :
524 : // Set appropriate components for kernels, including in the cases where a planar model is
525 : // running in planes other than the x-y plane (defined by _out_of_plane_strain_direction).
526 12028 : if (_out_of_plane_direction == OutOfPlaneDirection::x && i == 0)
527 30 : continue;
528 11998 : else if (_out_of_plane_direction == OutOfPlaneDirection::y && i == 1)
529 30 : continue;
530 :
531 11968 : params.set<unsigned int>("component") = i;
532 :
533 23936 : params.set<NonlinearVariableName>("variable") = _displacements[i];
534 :
535 11968 : if (_save_in.size() == _ndisp)
536 564 : params.set<std::vector<AuxVariableName>>("save_in") = {_save_in[i]};
537 11968 : if (_diag_save_in.size() == _ndisp)
538 0 : params.set<std::vector<AuxVariableName>>("diag_save_in") = {_diag_save_in[i]};
539 35904 : if (isParamValid("out_of_plane_strain") && !_lagrangian_kernels)
540 192 : params.set<std::vector<VariableName>>("out_of_plane_strain") = {
541 384 : getParam<VariableName>("out_of_plane_strain")};
542 :
543 11968 : if (_lk_homogenization)
544 : {
545 60 : params.set<std::vector<VariableName>>("scalar_variable") = {_hname};
546 20 : params.set<MultiMooseEnum>("constraint_types") = _constraint_types;
547 40 : params.set<std::vector<FunctionName>>("targets") = _targets;
548 : }
549 :
550 23936 : _problem->addKernel(ad_prepend + tensor_kernel_type, kernel_name, params);
551 12028 : }
552 :
553 4438 : if (_planar_formulation == PlanarFormulation::WeakPlaneStress)
554 : {
555 96 : auto params = getKernelParameters(ad_prepend + "WeakPlaneStress");
556 48 : std::string wps_kernel_name = "TM_WPS_" + name();
557 144 : params.set<NonlinearVariableName>("variable") = getParam<VariableName>("out_of_plane_strain");
558 :
559 96 : _problem->addKernel(ad_prepend + "WeakPlaneStress", wps_kernel_name, params);
560 48 : }
561 : }
562 38384 : }
563 :
564 : void
565 38175 : QuasiStaticSolidMechanicsPhysics::actSubdomainChecks()
566 : {
567 : // Do the coordinate system check only once the problem is created
568 38175 : if (_current_task == "setup_mesh_complete")
569 : {
570 : // get subdomain IDs
571 3605 : for (auto & name : _subdomain_names)
572 : {
573 318 : auto id = _mesh->getSubdomainID(name);
574 318 : if (id == Moose::INVALID_BLOCK_ID)
575 0 : paramError("block", "Subdomain \"" + name + "\" not found in mesh.");
576 : else
577 318 : _subdomain_ids.insert(id);
578 : }
579 : }
580 :
581 38175 : if (_current_task == "validate_coordinate_systems")
582 : {
583 : // use either block restriction list or list of all subdomains in the mesh
584 : const auto & check_subdomains =
585 4502 : _subdomain_ids.empty() ? _problem->mesh().meshSubdomains() : _subdomain_ids;
586 4502 : if (check_subdomains.empty())
587 0 : mooseError("No subdomains found");
588 :
589 : // make sure all subdomains are using the same coordinate system
590 4502 : _coord_system = _problem->getCoordSystem(*check_subdomains.begin());
591 9787 : for (auto subdomain : check_subdomains)
592 5287 : if (_problem->getCoordSystem(subdomain) != _coord_system)
593 2 : mooseError("The SolidMechanics action requires all subdomains to have the same coordinate "
594 : "system.");
595 :
596 4500 : if (_coord_system == Moose::COORD_RZ)
597 : {
598 235 : if (_out_of_plane_direction != OutOfPlaneDirection::z)
599 0 : mooseError("'out_of_plane_direction' must be 'z' for axisymmetric simulations");
600 : }
601 4265 : else if (_planar_formulation != PlanarFormulation::None)
602 : {
603 344 : if (_out_of_plane_direction == OutOfPlaneDirection::z && _ndisp != 2)
604 0 : mooseError(
605 : "Must specify two displacements for plane strain when the out of plane direction is z");
606 344 : else if (_out_of_plane_direction != OutOfPlaneDirection::z && _ndisp != 3)
607 0 : mooseError("Must specify three displacements for plane strain when the out of plane "
608 : "direction is x or y");
609 : }
610 : }
611 38173 : }
612 :
613 : void
614 38169 : QuasiStaticSolidMechanicsPhysics::actOutputGeneration()
615 : {
616 38169 : if (_current_task == "add_material")
617 3268 : actOutputMatProp();
618 :
619 : // Add variables (optional)
620 38169 : if (_current_task == "add_aux_variable")
621 : {
622 : unsigned int index = 0;
623 13881 : for (auto out : _generate_output)
624 : {
625 10610 : const auto & order = _material_output_order[index];
626 10610 : const auto & family = _material_output_family[index];
627 :
628 10588 : std::string type = (order == "CONSTANT" && family == "MONOMIAL")
629 : ? "MooseVariableConstMonomial"
630 21220 : : "MooseVariable";
631 :
632 : // Create output helper aux variables
633 10610 : auto params = _factory.getValidParams(type);
634 10610 : params.set<MooseEnum>("order") = order;
635 10610 : params.set<MooseEnum>("family") = family;
636 :
637 10610 : if (family == "MONOMIAL")
638 21194 : _problem->addAuxVariable(type, _base_name + out, params);
639 : else
640 26 : _problem->addVariable(type, _base_name + out, params);
641 :
642 10610 : index++;
643 10610 : }
644 : }
645 :
646 : // Add output AuxKernels
647 34898 : else if (_current_task == "add_aux_kernel")
648 : {
649 5775 : std::string ad_prepend = _use_ad ? "AD" : "";
650 : // Loop through output aux variables
651 : unsigned int index = 0;
652 13743 : for (auto out : _generate_output)
653 : {
654 10513 : if (_material_output_family[index] == "MONOMIAL")
655 : {
656 10500 : InputParameters params = emptyInputParameters();
657 :
658 21000 : params = _factory.getValidParams("MaterialRealAux");
659 10500 : params.applyParameters(parameters());
660 31500 : params.set<MaterialPropertyName>("property") = _base_name + out;
661 31500 : params.set<AuxVariableName>("variable") = _base_name + out;
662 10500 : params.set<ExecFlagEnum>("execute_on") = EXEC_TIMESTEP_END;
663 :
664 21000 : _problem->addAuxKernel(
665 21000 : ad_prepend + "MaterialRealAux", _base_name + out + '_' + name(), params);
666 10500 : }
667 10513 : index++;
668 : }
669 : }
670 31668 : else if (_current_task == "add_kernel")
671 : {
672 8170 : std::string ad_prepend = _use_ad ? "AD" : "";
673 : // Loop through output aux variables
674 : unsigned int index = 0;
675 14943 : for (auto out : _generate_output)
676 : {
677 10505 : if (_material_output_family[index] != "MONOMIAL")
678 : {
679 13 : InputParameters params = emptyInputParameters();
680 :
681 26 : params = _factory.getValidParams("MaterialPropertyValue");
682 13 : params.applyParameters(parameters());
683 39 : params.set<MaterialPropertyName>("prop_name") = _base_name + out;
684 39 : params.set<NonlinearVariableName>("variable") = _base_name + out;
685 :
686 26 : _problem->addKernel(
687 26 : ad_prepend + "MaterialPropertyValue", _base_name + out + '_' + name(), params);
688 13 : }
689 10505 : index++;
690 : }
691 : }
692 38169 : }
693 :
694 : void
695 31 : QuasiStaticSolidMechanicsPhysics::actEigenstrainNames()
696 : {
697 : // Create containers for collecting blockIDs and eigenstrain names from materials
698 : std::map<std::string, std::set<SubdomainID>> material_eigenstrain_map;
699 : std::set<std::string> eigenstrain_set;
700 :
701 : std::set<MaterialPropertyName> verified_eigenstrain_names;
702 :
703 : std::map<std::string, std::string> remove_add_map;
704 : std::set<std::string> remove_reduced_set;
705 :
706 : // Loop over all the materials(eigenstrains) already created
707 31 : auto materials = _problem->getMaterialWarehouse().getObjects();
708 300 : for (auto & mat : materials)
709 : {
710 269 : std::shared_ptr<BlockRestrictable> blk = std::dynamic_pointer_cast<BlockRestrictable>(mat);
711 : const InputParameters & mat_params = mat->parameters();
712 : auto & mat_name = mat->type();
713 :
714 : // Check for eigenstrain names, only deal with those materials
715 538 : if (mat_params.isParamValid("eigenstrain_name"))
716 : {
717 : std::shared_ptr<MaterialData> mat_dat;
718 59 : auto name = mat_params.get<std::string>("eigenstrain_name");
719 :
720 : // Check for base_name prefix
721 118 : if (mat_params.isParamValid("base_name"))
722 0 : name = mat_params.get<std::string>("base_name") + '_' + name;
723 :
724 : // Check block restrictions
725 59 : if (!blk)
726 0 : mooseError("Internal error, Material object that does not inherit form BlockRestricted");
727 : const std::set<SubdomainID> & blocks =
728 59 : blk->blockRestricted() ? blk->blockIDs() : blk->meshBlockIDs();
729 :
730 59 : if (std::includes(blocks.begin(), blocks.end(), _subdomain_ids.begin(), _subdomain_ids.end()))
731 : {
732 47 : material_eigenstrain_map[name].insert(blocks.begin(), blocks.end());
733 47 : eigenstrain_set.insert(name);
734 : }
735 : }
736 :
737 : // Account for reduced eigenstrains and CompositeEigenstrains
738 269 : if (mat_name == "ComputeReducedOrderEigenstrain")
739 : {
740 : auto input_eigenstrain_names =
741 3 : mat_params.get<std::vector<MaterialPropertyName>>("input_eigenstrain_names");
742 : remove_reduced_set.insert(input_eigenstrain_names.begin(), input_eigenstrain_names.end());
743 3 : }
744 : // Account for CompositeEigenstrains
745 269 : if (mat_name == "CompositeEigenstrain")
746 : {
747 3 : auto remove_list = mat_params.get<std::vector<MaterialPropertyName>>("tensors");
748 9 : for (auto i : remove_list)
749 6 : remove_reduced_set.insert(i);
750 3 : }
751 :
752 : // Account for MaterialADConverter , add or remove later
753 269 : if (mat_name == "RankTwoTensorMaterialADConverter")
754 : {
755 : std::vector<MaterialPropertyName> remove_list;
756 : std::vector<MaterialPropertyName> add_list;
757 :
758 18 : if (mat_params.isParamValid("ad_props_out") && mat_params.isParamValid("reg_props_in") &&
759 6 : _use_ad)
760 : {
761 6 : remove_list = mat_params.get<std::vector<MaterialPropertyName>>("reg_props_in");
762 6 : add_list = mat_params.get<std::vector<MaterialPropertyName>>("ad_props_out");
763 : }
764 18 : if (mat_params.isParamValid("ad_props_in") && mat_params.isParamValid("reg_props_out") &&
765 6 : !_use_ad)
766 : {
767 0 : remove_list = mat_params.get<std::vector<MaterialPropertyName>>("ad_props_in");
768 0 : add_list = mat_params.get<std::vector<MaterialPropertyName>>("reg_props_out");
769 : }
770 :
771 : // These vectors are the same size as checked in MaterialADConverter
772 12 : for (unsigned int index = 0; index < remove_list.size(); index++)
773 : remove_add_map.emplace(remove_list[index], add_list[index]);
774 6 : }
775 : }
776 : // All the materials have been accounted for, now remove or add parts
777 :
778 : // Remove names which aren't eigenstrains (converter properties)
779 37 : for (auto remove_add_index : remove_add_map)
780 : {
781 : const bool is_in = eigenstrain_set.find(remove_add_index.first) != eigenstrain_set.end();
782 6 : if (is_in)
783 : {
784 : eigenstrain_set.erase(remove_add_index.first);
785 3 : eigenstrain_set.insert(remove_add_index.second);
786 : }
787 : }
788 40 : for (auto index : remove_reduced_set)
789 : eigenstrain_set.erase(index);
790 :
791 : // Compare the blockIDs set of eigenstrain names with the vector of _eigenstrain_names for the
792 : // current subdomainID
793 31 : std::set_union(eigenstrain_set.begin(),
794 : eigenstrain_set.end(),
795 : _eigenstrain_names.begin(),
796 : _eigenstrain_names.end(),
797 : std::inserter(verified_eigenstrain_names, verified_eigenstrain_names.begin()));
798 :
799 : // Ensure the eigenstrain names previously passed include any missing names
800 31 : _eigenstrain_names.resize(verified_eigenstrain_names.size());
801 : std::copy(verified_eigenstrain_names.begin(),
802 : verified_eigenstrain_names.end(),
803 : _eigenstrain_names.begin());
804 :
805 31 : Moose::out << COLOR_CYAN << "*** Automatic Eigenstrain Names ***"
806 : << "\n"
807 93 : << _name << ": " << Moose::stringify(_eigenstrain_names) << "\n"
808 31 : << COLOR_DEFAULT << std::flush;
809 62 : }
810 :
811 : void
812 2282 : QuasiStaticSolidMechanicsPhysics::verifyOrderAndFamilyOutputs()
813 : {
814 : // Ensure material output order and family vectors are same size as generate output
815 :
816 : // check number of supplied orders and families
817 2282 : if (_material_output_order.size() > 1 && _material_output_order.size() < _generate_output.size())
818 0 : paramError("material_output_order",
819 : "The number of orders assigned to material outputs must be: 0 to be assigned "
820 : "CONSTANT; 1 to assign all outputs the same value, or the same size as the number "
821 : "of generate outputs listed.");
822 :
823 2289 : if (_material_output_family.size() > 1 &&
824 7 : _material_output_family.size() < _generate_output.size())
825 1 : paramError("material_output_family",
826 : "The number of families assigned to material outputs must be: 0 to be assigned "
827 : "MONOMIAL; 1 to assign all outputs the same value, or the same size as the number "
828 : "of generate outputs listed.");
829 :
830 : // if no value was provided, chose the default CONSTANT
831 2281 : if (_material_output_order.size() == 0)
832 4538 : _material_output_order.setAdditionalValue("CONSTANT");
833 :
834 : // For only one order, make all orders the same magnitude
835 2281 : if (_material_output_order.size() == 1)
836 : _material_output_order =
837 2275 : std::vector<std::string>(_generate_output.size(), _material_output_order[0]);
838 :
839 2281 : if (_verbose)
840 9 : Moose::out << COLOR_CYAN << "*** Automatic applied material output orders ***"
841 : << "\n"
842 27 : << _name << ": " << Moose::stringify(_material_output_order) << "\n"
843 9 : << COLOR_DEFAULT << std::flush;
844 :
845 : // if no value was provided, chose the default MONOMIAL
846 2281 : if (_material_output_family.size() == 0)
847 4544 : _material_output_family.setAdditionalValue("MONOMIAL");
848 :
849 : // For only one family, make all families that value
850 2281 : if (_material_output_family.size() == 1)
851 : _material_output_family =
852 2275 : std::vector<std::string>(_generate_output.size(), _material_output_family[0]);
853 :
854 2281 : if (_verbose)
855 9 : Moose::out << COLOR_CYAN << "*** Automatic applied material output families ***"
856 : << "\n"
857 27 : << _name << ": " << Moose::stringify(_material_output_family) << "\n"
858 9 : << COLOR_DEFAULT << std::flush;
859 2281 : }
860 :
861 : std::string
862 10598 : QuasiStaticSolidMechanicsPhysics::remapCompatOutputProp(const std::string & prop_name) const
863 : {
864 : // The OLD pipeline publishes `stress` and `mechanical_strain` directly. The NEW pipeline
865 : // publishes Cauchy stress under `cauchy_stress` and (when the strain calc is configured for
866 : // OLD compatibility) the OLD-equivalent mechanical strain under `rotated_mechanical_strain`.
867 : // Redirect the underlying property name so user output requests like `stress_xx` keep their
868 : // OLD-style aux variable name but read from the NEW-style source.
869 10598 : if (!_compatibility_mode)
870 10574 : return prop_name;
871 24 : if (prop_name == "stress")
872 12 : return "cauchy_stress";
873 12 : if (prop_name == "mechanical_strain")
874 12 : return "rotated_mechanical_strain";
875 0 : return prop_name;
876 : }
877 :
878 : void
879 3268 : QuasiStaticSolidMechanicsPhysics::actOutputMatProp()
880 : {
881 5849 : std::string ad_prepend = _use_ad ? "AD" : "";
882 :
883 3268 : if (_current_task == "add_material")
884 : {
885 : // Add output Materials
886 13866 : for (auto out : _generate_output)
887 : {
888 10598 : InputParameters params = emptyInputParameters();
889 :
890 : // RankTwoCartesianComponent
891 20155 : if (
892 21196 : [&]()
893 : {
894 117617 : for (const auto & r2q : _rank_two_cartesian_component_table)
895 445311 : for (unsigned int a = 0; a < 3; ++a)
896 1335353 : for (unsigned int b = 0; b < 3; ++b)
897 2013236 : if (r2q.first + '_' + _component_table[a] + _component_table[b] == out)
898 : {
899 9557 : auto type = ad_prepend + "RankTwoCartesianComponent";
900 9557 : params = _factory.getValidParams(type);
901 19114 : params.set<MaterialPropertyName>("rank_two_tensor") =
902 9557 : _base_name + remapCompatOutputProp(r2q.second);
903 9557 : params.set<unsigned int>("index_i") = a;
904 9557 : params.set<unsigned int>("index_j") = b;
905 :
906 9557 : params.applyParameters(parameters());
907 28671 : params.set<MaterialPropertyName>("property_name") = _base_name + out;
908 19114 : _problem->addMaterial(type, _base_name + out + '_' + name(), params);
909 : return true;
910 : }
911 : return false;
912 10598 : }())
913 9557 : continue;
914 :
915 : // RankTwoDirectionalComponent
916 2082 : if (setupOutput(out,
917 : _rank_two_directional_component_table,
918 0 : [&](std::string prop_name, std::string invariant)
919 : {
920 0 : auto type = ad_prepend + "RankTwoDirectionalComponent";
921 0 : params = _factory.getValidParams(type);
922 0 : params.set<MaterialPropertyName>("rank_two_tensor") =
923 0 : _base_name + remapCompatOutputProp(prop_name);
924 0 : params.set<MooseEnum>("invariant") = invariant;
925 0 : params.applyParameters(parameters());
926 0 : params.set<MaterialPropertyName>("property_name") = _base_name + out;
927 0 : _problem->addMaterial(type, _base_name + out + '_' + name(), params);
928 0 : }))
929 0 : continue;
930 :
931 : // RankTwoInvariant
932 2082 : if (setupOutput(out,
933 : _rank_two_invariant_table,
934 962 : [&](std::string prop_name, std::string invariant)
935 : {
936 962 : auto type = ad_prepend + "RankTwoInvariant";
937 962 : params = _factory.getValidParams(type);
938 1924 : params.set<MaterialPropertyName>("rank_two_tensor") =
939 962 : _base_name + remapCompatOutputProp(prop_name);
940 962 : params.set<MooseEnum>("invariant") = invariant;
941 962 : params.applyParameters(parameters());
942 2886 : params.set<MaterialPropertyName>("property_name") = _base_name + out;
943 1924 : _problem->addMaterial(type, _base_name + out + '_' + name(), params);
944 962 : }))
945 962 : continue;
946 :
947 : // RankTwoCylindricalComponent
948 158 : if (setupOutput(
949 : out,
950 : _rank_two_cylindrical_component_table,
951 24 : [&](std::string prop_name, std::string component)
952 : {
953 24 : if (_coord_system == Moose::COORD_RSPHERICAL)
954 0 : mooseError(
955 : "Cannot use cylindrical component output in a spherical coordinate system.");
956 24 : auto type = ad_prepend + "RankTwoCylindricalComponent";
957 24 : params = _factory.getValidParams(type);
958 48 : params.set<MaterialPropertyName>("rank_two_tensor") =
959 24 : _base_name + remapCompatOutputProp(prop_name);
960 24 : params.set<MooseEnum>("cylindrical_component") = component;
961 24 : params.applyParameters(parameters());
962 72 : params.set<MaterialPropertyName>("property_name") = _base_name + out;
963 48 : _problem->addMaterial(type, _base_name + out + '_' + name(), params);
964 24 : }))
965 24 : continue;
966 :
967 : // RankTwoSphericalComponent
968 110 : if (setupOutput(out,
969 : _rank_two_spherical_component_table,
970 55 : [&](std::string prop_name, std::string component)
971 : {
972 55 : auto type = ad_prepend + "RankTwoSphericalComponent";
973 55 : params = _factory.getValidParams(type);
974 110 : params.set<MaterialPropertyName>("rank_two_tensor") =
975 55 : _base_name + remapCompatOutputProp(prop_name);
976 55 : params.set<MooseEnum>("spherical_component") = component;
977 55 : params.applyParameters(parameters());
978 165 : params.set<MaterialPropertyName>("property_name") = _base_name + out;
979 110 : _problem->addMaterial(type, _base_name + out + '_' + name(), params);
980 55 : }))
981 55 : continue;
982 :
983 0 : paramError("generate_output", "Unable to add output Material for '", out, "'");
984 10598 : }
985 : }
986 3268 : }
987 :
988 : void
989 38173 : QuasiStaticSolidMechanicsPhysics::actGatherActionParameters()
990 : {
991 : // Gather info about all other solid mechanics physics when we add variables
992 47173 : if (_current_task == "validate_coordinate_systems" && getParam<bool>("add_variables"))
993 : {
994 3030 : auto actions = _awh.getActions<QuasiStaticSolidMechanicsPhysics>();
995 6174 : for (const auto & action : actions)
996 : {
997 : const auto size_before = _subdomain_id_union.size();
998 3148 : const auto added_size = action->_subdomain_ids.size();
999 3148 : _subdomain_id_union.insert(action->_subdomain_ids.begin(), action->_subdomain_ids.end());
1000 : const auto size_after = _subdomain_id_union.size();
1001 :
1002 3148 : if (size_after != size_before + added_size)
1003 2 : mooseError("The block restrictions in the SolidMechanics/QuasiStatic actions must be "
1004 : "non-overlapping.");
1005 :
1006 3146 : if (added_size == 0 && actions.size() > 1)
1007 2 : mooseError(
1008 : "No SolidMechanics/QuasiStatic action can be block unrestricted if more than one "
1009 : "SolidMechanics/QuasiStatic action is specified.");
1010 : }
1011 3026 : }
1012 38169 : }
1013 :
1014 : void
1015 272 : QuasiStaticSolidMechanicsPhysics::actLagrangianKernelStrain()
1016 : {
1017 : std::string type;
1018 272 : if (_coord_system == Moose::COORD_XYZ)
1019 : type = "ComputeLagrangianStrain";
1020 28 : else if (_coord_system == Moose::COORD_RZ)
1021 : type = "ComputeLagrangianStrainAxisymmetricCylindrical";
1022 8 : else if (_coord_system == Moose::COORD_RSPHERICAL)
1023 : type = "ComputeLagrangianStrainCentrosymmetricSpherical";
1024 : else
1025 0 : mooseError("Unsupported coordinate system");
1026 :
1027 272 : auto params = _factory.getValidParams(type);
1028 :
1029 544 : if (isParamValid("strain_base_name"))
1030 0 : params.set<std::string>("base_name") = getParam<std::string>("strain_base_name");
1031 :
1032 272 : params.set<std::vector<VariableName>>("displacements") = _coupled_displacements;
1033 272 : params.set<std::vector<MaterialPropertyName>>("eigenstrain_names") = _eigenstrain_names;
1034 272 : params.set<bool>("large_kinematics") = _lk_large_kinematics;
1035 272 : params.set<std::vector<SubdomainName>>("block") = _subdomain_names;
1036 :
1037 : // Error if volumetric locking correction is on for higher-order elements
1038 272 : if (_problem->mesh().hasSecondOrderElements() && _lk_locking)
1039 1 : mooseError("Volumetric locking correction should not be used for "
1040 : "higher-order elements.");
1041 :
1042 271 : params.set<bool>("stabilize_strain") = _lk_locking;
1043 :
1044 : // Forward the strain time-integration options. In compatibility mode kinematic_approximation
1045 : // and F_bar_mode are overridden below (and the action errors if the user also set them).
1046 813 : params.set<MooseEnum>("kinematic_approximation") = getParam<MooseEnum>("kinematic_approximation");
1047 542 : params.set<Real>("alpha") = getParam<Real>("generalized_midpoint_alpha");
1048 813 : params.set<MooseEnum>("F_bar_mode") = getParam<MooseEnum>("volumetric_locking_correction_mode");
1049 :
1050 271 : if (_lk_homogenization)
1051 : {
1052 16 : params.set<std::vector<MaterialPropertyName>>("homogenization_gradient_names") = {
1053 32 : _homogenization_strain_name};
1054 : }
1055 :
1056 : // OLD-compat configuration: map decomposition_method -> kinematic_approximation, switch the
1057 : // F-bar to incremental mode (matches OLD `ComputeFiniteStrain`'s `_Fhat` F-bar), and have
1058 : // the strain calc publish a real rotation_increment so the wrapped ComputeStressBase
1059 : // material's FSR rotates correctly across steps.
1060 271 : if (_compatibility_mode)
1061 : {
1062 12 : if (_lk_large_kinematics)
1063 : {
1064 8 : const std::string decomp = getParam<MooseEnum>("decomposition_method");
1065 8 : if (decomp == "TaylorExpansion")
1066 8 : params.set<MooseEnum>("kinematic_approximation") = "rashid_approximate";
1067 4 : else if (decomp == "EigenSolution")
1068 8 : params.set<MooseEnum>("kinematic_approximation") = "rashid_eigen";
1069 : // HughesWinget rejected in the ctor; SMALL strain ignores kinematic_approximation.
1070 : }
1071 :
1072 : // OLD's `volumetric_locking_correction = true` in incremental finite strain maps to
1073 : // `F_bar_mode = incremental` on the NEW strain calc. The strain calc rejects
1074 : // `incremental` for small kinematics -- for SMALL + locking we leave `F_bar_mode` at its
1075 : // default (`total`), which is the additive trace correction (matches OLD's small-strain
1076 : // B-bar in the locked-flag-equivalent form).
1077 12 : if (_lk_locking && _lk_large_kinematics)
1078 16 : params.set<MooseEnum>("F_bar_mode") = "incremental";
1079 :
1080 : // Only meaningful with FINITE -- in SMALL the wrap runs in passthrough-off mode and the
1081 : // rotation increment isn't consumed downstream.
1082 12 : params.set<bool>("publish_rotation_increment") = _lk_large_kinematics;
1083 : }
1084 :
1085 271 : _problem->addMaterial(type, name() + "_strain", params);
1086 :
1087 : // OLD-compat configuration: auto-add ComputeLagrangianWrappedStress around the user's
1088 : // ComputeStressBase-style stress material. `input_stress = "stress"` and
1089 : // `input_jacobian = "Jacobian_mult"` default to the property names every
1090 : // ComputeStressBase descendant publishes.
1091 271 : if (_compatibility_mode)
1092 : {
1093 12 : const std::string wrap_type = "ComputeLagrangianWrappedStress";
1094 12 : auto wrap_params = _factory.getValidParams(wrap_type);
1095 24 : if (isParamValid("base_name"))
1096 0 : wrap_params.set<std::string>("base_name") = getParam<std::string>("base_name");
1097 12 : wrap_params.set<std::vector<SubdomainName>>("block") = _subdomain_names;
1098 : // large_kinematics is derived from the strain calculator's LARGE_KINEMATICS guarantee.
1099 24 : wrap_params.set<MooseEnum>("objective_rate") = "rashid";
1100 : // FINITE: the wrap runs in passthrough mode (the wrapped material's FSR already produced
1101 : // the correctly-rotated cumulative Cauchy stress, fed by our published rotation_increment).
1102 : // SMALL: nothing to rotate, plain Cauchy passthrough via the rate's small-kinematics branch.
1103 12 : wrap_params.set<bool>("rotate_old_stress") = _lk_large_kinematics;
1104 12 : _problem->addMaterial(wrap_type, name() + "_compatibility_wrap", wrap_params);
1105 12 : }
1106 :
1107 : // Add the homogenization strain calculator
1108 271 : if (_lk_homogenization)
1109 : {
1110 8 : std::string type = "ComputeHomogenizedLagrangianStrain";
1111 8 : auto params = _factory.getValidParams(type);
1112 :
1113 16 : params.set<MaterialPropertyName>("homogenization_gradient_name") = _homogenization_strain_name;
1114 24 : params.set<std::vector<VariableName>>("macro_gradient") = {_hname};
1115 8 : params.set<MultiMooseEnum>("constraint_types") = _constraint_types;
1116 8 : params.set<std::vector<FunctionName>>("targets") = _targets;
1117 :
1118 8 : _problem->addMaterial(type, name() + "_compute_" + _homogenization_strain_name, params);
1119 8 : }
1120 542 : }
1121 :
1122 : void
1123 2962 : QuasiStaticSolidMechanicsPhysics::actStressDivergenceTensorsStrain()
1124 : {
1125 5239 : std::string ad_prepend = _use_ad ? "AD" : "";
1126 :
1127 : std::string type;
1128 :
1129 : // no plane strain
1130 2962 : if (_planar_formulation == PlanarFormulation::None)
1131 : {
1132 : std::map<std::pair<Moose::CoordinateSystemType, StrainAndIncrement>, std::string> type_map = {
1133 : {{Moose::COORD_XYZ, StrainAndIncrement::SmallTotal}, "ComputeSmallStrain"},
1134 : {{Moose::COORD_XYZ, StrainAndIncrement::SmallIncremental}, "ComputeIncrementalStrain"},
1135 : {{Moose::COORD_XYZ, StrainAndIncrement::FiniteIncremental}, "ComputeFiniteStrain"},
1136 : {{Moose::COORD_RZ, StrainAndIncrement::SmallTotal}, "ComputeAxisymmetricRZSmallStrain"},
1137 : {{Moose::COORD_RZ, StrainAndIncrement::SmallIncremental},
1138 : "ComputeAxisymmetricRZIncrementalStrain"},
1139 : {{Moose::COORD_RZ, StrainAndIncrement::FiniteIncremental},
1140 : "ComputeAxisymmetricRZFiniteStrain"},
1141 : {{Moose::COORD_RSPHERICAL, StrainAndIncrement::SmallTotal}, "ComputeRSphericalSmallStrain"},
1142 : {{Moose::COORD_RSPHERICAL, StrainAndIncrement::SmallIncremental},
1143 : "ComputeRSphericalIncrementalStrain"},
1144 : {{Moose::COORD_RSPHERICAL, StrainAndIncrement::FiniteIncremental},
1145 25950 : "ComputeRSphericalFiniteStrain"}};
1146 :
1147 2595 : auto type_it = type_map.find(std::make_pair(_coord_system, _strain_and_increment));
1148 2595 : if (type_it != type_map.end())
1149 2595 : type = type_it->second;
1150 : else
1151 0 : mooseError("Unsupported strain formulation");
1152 : }
1153 367 : else if (_planar_formulation == PlanarFormulation::WeakPlaneStress ||
1154 367 : _planar_formulation == PlanarFormulation::PlaneStrain ||
1155 : _planar_formulation == PlanarFormulation::GeneralizedPlaneStrain)
1156 : {
1157 367 : if (_use_ad && (_planar_formulation == PlanarFormulation::PlaneStrain ||
1158 : _planar_formulation == PlanarFormulation::GeneralizedPlaneStrain))
1159 0 : paramError("use_automatic_differentiation",
1160 : "AD not setup for use with PlaneStrain or GeneralizedPlaneStrain");
1161 :
1162 : std::map<std::pair<Moose::CoordinateSystemType, StrainAndIncrement>, std::string> type_map = {
1163 : {{Moose::COORD_XYZ, StrainAndIncrement::SmallTotal}, "ComputePlaneSmallStrain"},
1164 : {{Moose::COORD_XYZ, StrainAndIncrement::SmallIncremental}, "ComputePlaneIncrementalStrain"},
1165 : {{Moose::COORD_XYZ, StrainAndIncrement::FiniteIncremental}, "ComputePlaneFiniteStrain"},
1166 : {{Moose::COORD_RZ, StrainAndIncrement::SmallTotal}, "ComputeAxisymmetric1DSmallStrain"},
1167 : {{Moose::COORD_RZ, StrainAndIncrement::SmallIncremental},
1168 : "ComputeAxisymmetric1DIncrementalStrain"},
1169 : {{Moose::COORD_RZ, StrainAndIncrement::FiniteIncremental},
1170 2569 : "ComputeAxisymmetric1DFiniteStrain"}};
1171 :
1172 : // choose kernel type based on coordinate system
1173 367 : auto type_it = type_map.find(std::make_pair(_coord_system, _strain_and_increment));
1174 367 : if (type_it != type_map.end())
1175 367 : type = type_it->second;
1176 : else
1177 0 : mooseError("Unsupported coordinate system for plane strain.");
1178 : }
1179 : else
1180 0 : mooseError("Unsupported planar formulation");
1181 :
1182 : // set material parameters
1183 2962 : auto params = _factory.getValidParams(ad_prepend + type);
1184 2962 : params.applyParameters(
1185 : parameters(),
1186 : {"displacements", "use_displaced_mesh", "out_of_plane_strain", "scalar_out_of_plane_strain"});
1187 :
1188 5924 : if (isParamValid("strain_base_name"))
1189 0 : params.set<std::string>("base_name") = getParam<std::string>("strain_base_name");
1190 :
1191 2962 : params.set<std::vector<VariableName>>("displacements") = _coupled_displacements;
1192 2962 : params.set<bool>("use_displaced_mesh") = false;
1193 :
1194 5924 : if (isParamValid("scalar_out_of_plane_strain"))
1195 186 : params.set<std::vector<VariableName>>("scalar_out_of_plane_strain") = {
1196 372 : getParam<VariableName>("scalar_out_of_plane_strain")};
1197 :
1198 5924 : if (isParamValid("out_of_plane_strain"))
1199 96 : params.set<std::vector<VariableName>>("out_of_plane_strain") = {
1200 192 : getParam<VariableName>("out_of_plane_strain")};
1201 :
1202 2962 : params.set<std::vector<MaterialPropertyName>>("eigenstrain_names") = _eigenstrain_names;
1203 :
1204 5924 : _problem->addMaterial(ad_prepend + type, name() + "_strain", params);
1205 34725 : }
1206 :
1207 : std::string
1208 11549 : QuasiStaticSolidMechanicsPhysics::getKernelType()
1209 : {
1210 11549 : if (_lagrangian_kernels)
1211 : {
1212 : std::string type;
1213 742 : if (_coord_system == Moose::COORD_XYZ)
1214 : {
1215 694 : if (_lk_homogenization)
1216 : type = "HomogenizedTotalLagrangianStressDivergence";
1217 674 : else if (_lk_formulation == LKFormulation::Total)
1218 : type = "TotalLagrangianStressDivergence";
1219 52 : else if (_lk_formulation == LKFormulation::Updated)
1220 : type = "UpdatedLagrangianStressDivergence";
1221 : else
1222 0 : mooseError("Unknown formulation type");
1223 : }
1224 48 : else if (_coord_system == Moose::COORD_RZ)
1225 : {
1226 40 : if (_lk_homogenization)
1227 0 : mooseError("The Lagrangian mechanics kernels do not yet support homogenization in "
1228 : "coordinate systems other than Cartesian.");
1229 40 : else if (_lk_formulation == LKFormulation::Total)
1230 : type = "TotalLagrangianStressDivergenceAxisymmetricCylindrical";
1231 0 : else if (_lk_formulation == LKFormulation::Updated)
1232 0 : mooseError("The Lagrangian mechanics kernels do not yet support the updated Lagrangian "
1233 : "formulation in RZ coordinates.");
1234 : }
1235 8 : else if (_coord_system == Moose::COORD_RSPHERICAL)
1236 : {
1237 8 : if (_lk_homogenization)
1238 0 : mooseError("The Lagrangian mechanics kernels do not yet support homogenization in "
1239 : "coordinate systems other than Cartesian.");
1240 8 : else if (_lk_formulation == LKFormulation::Total)
1241 : type = "TotalLagrangianStressDivergenceCentrosymmetricSpherical";
1242 0 : else if (_lk_formulation == LKFormulation::Updated)
1243 0 : mooseError("The Lagrangian mechanics kernels do not yet support the updated Lagrangian "
1244 : "formulation in RZ coordinates.");
1245 : }
1246 : else
1247 0 : mooseError("Unsupported coordinate system");
1248 742 : return type;
1249 : }
1250 : else
1251 : {
1252 : std::map<Moose::CoordinateSystemType, std::string> type_map = {
1253 : {Moose::COORD_XYZ, "StressDivergenceTensors"},
1254 : {Moose::COORD_RZ, "StressDivergenceRZTensors"},
1255 43228 : {Moose::COORD_RSPHERICAL, "StressDivergenceRSphericalTensors"}};
1256 :
1257 : // choose kernel type based on coordinate system
1258 10807 : auto type_it = type_map.find(_coord_system);
1259 10807 : if (type_it != type_map.end())
1260 10807 : return type_it->second;
1261 : else
1262 0 : mooseError("Unsupported coordinate system");
1263 : }
1264 0 : }
1265 :
1266 : InputParameters
1267 12076 : QuasiStaticSolidMechanicsPhysics::getKernelParameters(std::string type)
1268 : {
1269 12076 : InputParameters params = _factory.getValidParams(type);
1270 12076 : params.applyParameters(
1271 : parameters(),
1272 : {"displacements", "use_displaced_mesh", "save_in", "diag_save_in", "out_of_plane_strain"});
1273 :
1274 12076 : params.set<std::vector<VariableName>>("displacements") = _coupled_displacements;
1275 :
1276 12076 : if (_lagrangian_kernels)
1277 : {
1278 742 : params.set<bool>("use_displaced_mesh") =
1279 742 : (_lk_large_kinematics && (_lk_formulation == LKFormulation::Updated));
1280 : // large_kinematics is derived from the strain calculator's LARGE_KINEMATICS guarantee.
1281 742 : params.set<bool>("stabilize_strain") = _lk_locking;
1282 2226 : params.set<MooseEnum>("F_bar_mode") = getParam<MooseEnum>("volumetric_locking_correction_mode");
1283 742 : if (_lk_homogenization)
1284 20 : params.set<bool>("off_diagonal_jacobian") = _lk_h_off_jac;
1285 : // Match the strain calc's F-bar mode on the kernel -- the TL kernel's OLD-compat B-bar
1286 : // residual is gated on `F_bar_mode = incremental && stabilize_strain`. Setting it on
1287 : // both keeps the kernel and the strain material on the same F-bar formulation. Only
1288 : // applicable to large kinematics (the strain calc rejects `incremental` for small).
1289 742 : if (_compatibility_mode && _lk_locking && _lk_large_kinematics)
1290 48 : params.set<MooseEnum>("F_bar_mode") = "incremental";
1291 : }
1292 : else
1293 11334 : params.set<bool>("use_displaced_mesh") = _use_displaced_mesh;
1294 :
1295 12076 : return params;
1296 0 : }
|
