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 "WeightedGapUserObject.h"
11 : #include "MooseVariableField.h"
12 : #include "SubProblem.h"
13 : #include "MortarUtils.h"
14 : #include "MooseUtils.h"
15 : #include "MortarContactUtils.h"
16 : #include "AutomaticMortarGeneration.h"
17 : #include "ADUtils.h"
18 :
19 : #include "libmesh/quadrature.h"
20 :
21 : #include <limits>
22 :
23 : InputParameters
24 2010 : WeightedGapUserObject::validParams()
25 : {
26 2010 : InputParameters params = MortarUserObject::validParams();
27 2010 : params += MortarConsumerInterface::validParams();
28 2010 : params += TwoMaterialPropertyInterface::validParams();
29 4020 : params.addRequiredCoupledVar("disp_x", "The x displacement variable");
30 4020 : params.addRequiredCoupledVar("disp_y", "The y displacement variable");
31 4020 : params.addCoupledVar("disp_z", "The z displacement variable");
32 2010 : params.set<bool>("use_displaced_mesh") = true;
33 2010 : params.set<bool>("interpolate_normals") = false;
34 2010 : params.addPrivateParam<bool>("allow_nodal_normal_derivatives", false);
35 2010 : params.addPrivateParam<bool>("use_nodal_normal_derivatives", false);
36 8040 : params.set<ExecFlagEnum>("execute_on") = {EXEC_LINEAR, EXEC_NONLINEAR};
37 2010 : params.suppressParameter<ExecFlagEnum>("execute_on");
38 2010 : return params;
39 2010 : }
40 :
41 886 : WeightedGapUserObject::WeightedGapUserObject(const InputParameters & parameters)
42 : : MortarUserObject(parameters),
43 1760 : _fe_problem(*getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
44 880 : _nodal(getVar("disp_x", 0)->feType().family == LAGRANGE),
45 880 : _disp_x_var(getVar("disp_x", 0)),
46 880 : _disp_y_var(getVar("disp_y", 0)),
47 880 : _has_disp_z(isCoupled("disp_z")),
48 880 : _disp_z_var(_has_disp_z ? getVar("disp_z", 0) : nullptr),
49 880 : _secondary_disp_x(_disp_x_var->adSln()),
50 880 : _primary_disp_x(_disp_x_var->adSlnNeighbor()),
51 880 : _secondary_disp_y(_disp_y_var->adSln()),
52 880 : _primary_disp_y(_disp_y_var->adSlnNeighbor()),
53 880 : _secondary_disp_z(_has_disp_z ? &_disp_z_var->adSln() : nullptr),
54 880 : _primary_disp_z(_has_disp_z ? &_disp_z_var->adSlnNeighbor() : nullptr),
55 880 : _coord(_assembly.mortarCoordTransformation()),
56 880 : _allow_nodal_normal_derivatives(getParam<bool>("allow_nodal_normal_derivatives")),
57 1766 : _use_nodal_normal_derivatives(false)
58 : {
59 1760 : if (!getParam<bool>("use_displaced_mesh"))
60 0 : paramError("use_displaced_mesh",
61 : "'use_displaced_mesh' must be true for the WeightedGapUserObject object");
62 :
63 : // A valid penetration_tolerance identifies augmented-Lagrange contact, which retains frozen
64 : // normal directions.
65 2139 : if (getParam<bool>("use_nodal_normal_derivatives") && !isParamValid("penetration_tolerance"))
66 170 : includeNodalNormalDerivatives();
67 880 : }
68 :
69 : void
70 869 : WeightedGapUserObject::includeNodalNormalDerivatives() const
71 : {
72 869 : if (_use_nodal_normal_derivatives)
73 161 : return;
74 :
75 708 : if (!_allow_nodal_normal_derivatives)
76 0 : mooseError("Nodal-normal derivatives are not supported by user object '", name(), "'.");
77 :
78 1416 : if (getParam<bool>("interpolate_normals"))
79 0 : paramError("interpolate_normals",
80 : "Nodal-normal derivatives require normalized secondary nodal normals and cannot be "
81 : "combined with quadrature-point normal interpolation.");
82 :
83 : const std::array<std::pair<const MooseVariable *, const char *>, 3> displacement_variables{
84 : {{_disp_x_var, "disp_x"}, {_disp_y_var, "disp_y"}, {_disp_z_var, "disp_z"}}};
85 2832 : for (const auto & [variable, parameter_name] : displacement_variables)
86 2124 : if (variable)
87 : {
88 1651 : if (!variable->isNodal())
89 0 : paramError(parameter_name,
90 : "Nodal-normal derivatives require a nodal displacement variable.");
91 1651 : if (&variable->sys() != &_sys)
92 0 : paramError(parameter_name,
93 : "Nodal-normal derivatives require displacement variables in the nonlinear "
94 : "system assembled by this contact object.");
95 : }
96 :
97 708 : _use_nodal_normal_derivatives = true;
98 : }
99 :
100 : bool
101 424734636 : WeightedGapUserObject::usesNodalNormalDerivatives() const
102 : {
103 424734636 : return _use_nodal_normal_derivatives && Moose::doDerivatives(_subproblem, _sys);
104 : }
105 :
106 : void
107 857 : WeightedGapUserObject::initialSetup()
108 : {
109 857 : MortarUserObject::initialSetup();
110 857 : _test = &test();
111 857 : }
112 :
113 : void
114 9047293 : WeightedGapUserObject::computeQpProperties()
115 : {
116 : // Trim interior node variable derivatives
117 : const auto & primary_ip_lowerd_map = amg().getPrimaryIpToLowerElementMap(
118 9047293 : *_lower_primary_elem, *_lower_primary_elem->interior_parent(), *_lower_secondary_elem);
119 : const auto & secondary_ip_lowerd_map =
120 9047293 : amg().getSecondaryIpToLowerElementMap(*_lower_secondary_elem);
121 :
122 9047293 : std::array<const MooseVariable *, 3> var_array{{_disp_x_var, _disp_y_var, _disp_z_var}};
123 : std::array<ADReal, 3> primary_disp{
124 9047293 : {_primary_disp_x[_qp], _primary_disp_y[_qp], _has_disp_z ? (*_primary_disp_z)[_qp] : 0}};
125 9047293 : std::array<ADReal, 3> secondary_disp{{_secondary_disp_x[_qp],
126 9047293 : _secondary_disp_y[_qp],
127 9047293 : _has_disp_z ? (*_secondary_disp_z)[_qp] : 0}};
128 :
129 9047293 : trimInteriorNodeDerivatives(primary_ip_lowerd_map, var_array, primary_disp, false);
130 9047293 : trimInteriorNodeDerivatives(secondary_ip_lowerd_map, var_array, secondary_disp, true);
131 :
132 : const ADReal & prim_x = primary_disp[0];
133 : const ADReal & prim_y = primary_disp[1];
134 : const ADReal * prim_z = nullptr;
135 9047293 : if (_has_disp_z)
136 : prim_z = &primary_disp[2];
137 :
138 : const ADReal & sec_x = secondary_disp[0];
139 : const ADReal & sec_y = secondary_disp[1];
140 : const ADReal * sec_z = nullptr;
141 9047293 : if (_has_disp_z)
142 : sec_z = &secondary_disp[2];
143 :
144 : // Compute gap vector
145 9047293 : ADRealVectorValue gap_vec = _phys_points_primary[_qp] - _phys_points_secondary[_qp];
146 :
147 : // Generic displacement for interface problems
148 9047293 : _qp_displacement_nodal(0) = prim_x - sec_x;
149 9047293 : _qp_displacement_nodal(1) = prim_y - sec_y;
150 9047293 : if (_has_disp_z)
151 7973790 : _qp_displacement_nodal(2) = *prim_z - *sec_z;
152 :
153 9047293 : _qp_displacement_nodal *= _JxW_msm[_qp] * _coord[_qp];
154 :
155 9047293 : gap_vec(0).derivatives() = prim_x.derivatives() - sec_x.derivatives();
156 9047293 : gap_vec(1).derivatives() = prim_y.derivatives() - sec_y.derivatives();
157 9047293 : if (_has_disp_z)
158 7973790 : gap_vec(2).derivatives() = prim_z->derivatives() - sec_z->derivatives();
159 :
160 : // Compute integration point quantities: Normals (geometry) is averaged at the node, but not
161 : // interpolated within the weak integration.
162 9047293 : _qp_gap_nodal = gap_vec * (_JxW_msm[_qp] * _coord[_qp]);
163 :
164 : // To do normalization of constraint coefficient (c_n)
165 9047293 : _qp_factor = _JxW_msm[_qp] * _coord[_qp];
166 9047293 : }
167 :
168 : void
169 34267874 : WeightedGapUserObject::computeQpIProperties()
170 : {
171 34267874 : if (usesNodalNormalDerivatives())
172 : mooseAssert(_lower_secondary_elem->n_nodes() > _i,
173 : "Making sure that the nodal normal index is valid");
174 : else
175 : mooseAssert(_normals.size() > _i, "Making sure that _normals is the expected size");
176 :
177 : // Get the _dof_to_weighted_gap map
178 34267874 : const auto * const dof = static_cast<const DofObject *>(_lower_secondary_elem->node_ptr(_i));
179 :
180 : auto & [weighted_gap, normalization] = _dof_to_weighted_gap[dof];
181 :
182 34267874 : if (usesNodalNormalDerivatives())
183 6185878 : weighted_gap += (*_test)[_i][_qp] * _qp_gap_nodal * contactNormal(*_lower_secondary_elem, _i);
184 : else
185 28081996 : weighted_gap += (*_test)[_i][_qp] * _qp_gap_nodal * _normals[_i];
186 :
187 34267874 : normalization += (*_test)[_i][_qp] * _qp_factor;
188 :
189 68535748 : _dof_to_weighted_displacements[dof] += (*_test)[_i][_qp] * _qp_displacement_nodal;
190 34267874 : }
191 :
192 : const ADRealVectorValue &
193 74307078 : WeightedGapUserObject::contactNormal(const Elem & lower_secondary_elem,
194 : const unsigned int nodal_index) const
195 : {
196 : mooseAssert(nodal_index < lower_secondary_elem.n_nodes(),
197 : "Nodal normal index must refer to a node on the secondary element.");
198 : mooseAssert(usesNodalNormalDerivatives(),
199 : "AD contact normals should only be requested while recording nodal-normal "
200 : "derivatives.");
201 74307078 : const Node * const node = lower_secondary_elem.node_ptr(nodal_index);
202 : const auto normal_it = _ad_nodal_normals.find(node);
203 74307078 : if (normal_it != _ad_nodal_normals.end())
204 74293743 : return normal_it->second;
205 :
206 13335 : amg().computeADNodalNormals(
207 : [this](const Node & coordinate_node, const Point & geometry_coordinate)
208 219377 : { return nodalCoordinate(coordinate_node, geometry_coordinate); },
209 13335 : _ad_nodal_normals);
210 13335 : return libmesh_map_find(_ad_nodal_normals, node);
211 : }
212 :
213 : ADPoint
214 219377 : WeightedGapUserObject::nodalCoordinate(const Node & node, const Point & geometry_coordinate) const
215 : {
216 : mooseAssert(usesNodalNormalDerivatives(),
217 : "AD nodal coordinates should only be requested while recording nodal-normal "
218 : "derivatives.");
219 219377 : ADPoint point = geometry_coordinate;
220 :
221 : const std::array<std::pair<const MooseVariable *, unsigned int>, 3> displacement_variables{
222 : {{_disp_x_var, 0}, {_disp_y_var, 1}, {_disp_z_var, 2}}};
223 877508 : for (const auto & [variable, component] : displacement_variables)
224 658131 : if (variable)
225 : {
226 505074 : const auto sys_num = variable->sys().number();
227 : const auto var_num = variable->number();
228 505074 : if (node.n_dofs(sys_num, var_num))
229 505074 : Moose::derivInsert(
230 505074 : point(component).derivatives(), node.dof_number(sys_num, var_num, 0), 1.0);
231 : }
232 :
233 219377 : return point;
234 : }
235 :
236 : void
237 114822 : WeightedGapUserObject::initialize()
238 : {
239 : _ad_nodal_normals.clear();
240 : _dof_to_weighted_gap.clear();
241 : _dof_to_weighted_displacements.clear();
242 114822 : }
243 :
244 : void
245 114822 : WeightedGapUserObject::finalize()
246 : {
247 : // If the constraint is performed by the owner, then we don't need any data sent back; the owner
248 : // will take care of it. But if the constraint is not performed by the owner and we might have to
249 : // do some of the constraining ourselves, then we need data sent back to us
250 114822 : const bool send_data_back = !constrainedByOwner();
251 114822 : Moose::Mortar::Contact::communicateGaps(_dof_to_weighted_gap,
252 114822 : _subproblem.mesh(),
253 114822 : _nodal,
254 : /*normalize_c*/ true,
255 : _communicator,
256 : send_data_back);
257 114822 : }
258 :
259 : void
260 2441208 : WeightedGapUserObject::execute()
261 : {
262 11488501 : for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
263 : {
264 9047293 : computeQpProperties();
265 43315167 : for (_i = 0; _i < _test->size(); ++_i)
266 34267874 : computeQpIProperties();
267 : }
268 2441208 : }
269 :
270 : Real
271 142883 : WeightedGapUserObject::getNormalGap(const Node * const node) const
272 : {
273 142883 : const auto it = _dof_to_weighted_gap.find(_subproblem.mesh().nodePtr(node->id()));
274 :
275 : // We are returning the physical weighted gap for analysis purposes
276 142883 : if (it != _dof_to_weighted_gap.end())
277 17872 : return physicalGap(it->second);
278 : else
279 : return 0.0;
280 : }
|
