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 : // MOOSE includes
11 : #include "NodalFrictionalConstraint.h"
12 : #include "NodalConstraintUtils.h"
13 : #include "MooseMesh.h"
14 : #include "Assembly.h"
15 : #include "SystemBase.h"
16 :
17 : #include "libmesh/null_output_iterator.h"
18 : #include "libmesh/parallel.h"
19 : #include "libmesh/parallel_elem.h"
20 : #include "libmesh/parallel_node.h"
21 :
22 : // C++ includes
23 : #include <limits.h>
24 :
25 : using namespace libMesh;
26 :
27 : registerMooseObject("SolidMechanicsApp", NodalFrictionalConstraint);
28 :
29 : InputParameters
30 28 : NodalFrictionalConstraint::validParams()
31 : {
32 28 : InputParameters params = NodalConstraint::validParams();
33 28 : params.addClassDescription("Frictional nodal constraint for contact");
34 56 : params.addRequiredParam<BoundaryName>("boundary", "The primary boundary");
35 56 : params.addRequiredParam<BoundaryName>("secondary", "The secondary boundary");
36 56 : params.addRequiredParam<Real>("friction_coefficient",
37 : "Friction coefficient for slippage in the normal direction");
38 56 : params.addRequiredParam<Real>("normal_force",
39 : "Normal force used together with friction_coefficient to compute "
40 : "the normal frictional capacity.");
41 56 : params.addRequiredParam<Real>("tangential_penalty",
42 : "Stiffness of the spring in the tangential direction.");
43 28 : return params;
44 0 : }
45 :
46 14 : NodalFrictionalConstraint::NodalFrictionalConstraint(const InputParameters & parameters)
47 : : NodalConstraint(parameters),
48 28 : _primary_boundary_id(getParam<BoundaryName>("boundary")),
49 14 : _secondary_boundary_id(getParam<BoundaryName>("secondary")),
50 28 : _normal_force(getParam<Real>("normal_force")),
51 28 : _tangential_penalty(getParam<Real>("tangential_penalty")),
52 28 : _friction_coefficient(getParam<Real>("friction_coefficient")),
53 14 : _u_secondary_old(_var.dofValuesOldNeighbor()),
54 28 : _u_primary_old(_var.dofValuesOld())
55 : {
56 14 : if (_var.number() != _var_secondary.number())
57 0 : paramError("variable_secondary",
58 : "Primary variable must be identical to secondary "
59 : "variable. Different variables are currently not supported.");
60 :
61 14 : updateConstrainedNodes();
62 :
63 28 : MooseEnum temp_formulation = getParam<MooseEnum>("formulation");
64 14 : if (temp_formulation == "penalty")
65 14 : _formulation = Moose::Penalty;
66 0 : else if (temp_formulation == "kinematic")
67 0 : mooseError("NodalFrictionalConstraint: Kinematic formulation is currently not supported for "
68 : "this constraint.");
69 : else
70 0 : mooseError("Formulation must be set to Penalty.");
71 14 : }
72 :
73 : void
74 0 : NodalFrictionalConstraint::meshChanged()
75 : {
76 0 : updateConstrainedNodes();
77 0 : }
78 :
79 : void
80 14 : NodalFrictionalConstraint::updateConstrainedNodes()
81 : {
82 14 : _primary_node_vector.clear();
83 14 : _connected_nodes.clear();
84 14 : _primary_conn.clear();
85 :
86 : std::vector<dof_id_type> secondary_nodelist =
87 14 : _mesh.getNodeList(_mesh.getBoundaryID(_secondary_boundary_id));
88 : std::vector<dof_id_type> primary_nodelist =
89 14 : _mesh.getNodeList(_mesh.getBoundaryID(_primary_boundary_id));
90 :
91 : // Fill in _connected_nodes, which defines secondary nodes in the base class
92 28 : for (auto in : secondary_nodelist)
93 : {
94 14 : if (_mesh.nodeRef(in).processor_id() == _subproblem.processor_id())
95 10 : _connected_nodes.push_back(in);
96 : }
97 :
98 : // Fill in _primary_node_vector, which defines secondary nodes in the base class
99 42 : for (auto in : primary_nodelist)
100 28 : _primary_node_vector.push_back(in);
101 :
102 14 : const auto & node_to_elem_map = _mesh.nodeToElemMap();
103 14 : std::vector<std::vector<dof_id_type>> elems(_primary_node_vector.size());
104 :
105 : // Add elements connected to primary node to Ghosted Elements.
106 :
107 : // On a distributed mesh, these elements might have already been
108 : // remoted, in which case we need to gather them back first.
109 14 : if (!_mesh.getMesh().is_serial())
110 : {
111 : std::set<Elem *, CompareElemsByLevel> primary_elems_to_ghost;
112 : std::set<Node *> nodes_to_ghost;
113 :
114 0 : for (unsigned int i = 0; i < primary_nodelist.size(); ++i)
115 : {
116 : auto node_to_elem_pair = node_to_elem_map.find(_primary_node_vector[i]);
117 :
118 : bool found_elems = (node_to_elem_pair != node_to_elem_map.end());
119 :
120 : #ifndef NDEBUG
121 : bool someone_found_elems = found_elems;
122 : _mesh.getMesh().comm().max(someone_found_elems);
123 : mooseAssert(someone_found_elems, "Missing entry in node to elem map");
124 : #endif
125 :
126 0 : if (found_elems)
127 : {
128 0 : for (auto id : node_to_elem_pair->second)
129 : {
130 0 : Elem * elem = _mesh.queryElemPtr(id);
131 0 : if (elem)
132 : {
133 0 : primary_elems_to_ghost.insert(elem);
134 :
135 0 : const unsigned int n_nodes = elem->n_nodes();
136 0 : for (unsigned int n = 0; n != n_nodes; ++n)
137 0 : nodes_to_ghost.insert(elem->node_ptr(n));
138 : }
139 : }
140 : }
141 : }
142 :
143 : // Send nodes first since elements need them
144 0 : _mesh.getMesh().comm().allgather_packed_range(&_mesh.getMesh(),
145 : nodes_to_ghost.begin(),
146 : nodes_to_ghost.end(),
147 : null_output_iterator<Node>());
148 :
149 0 : _mesh.getMesh().comm().allgather_packed_range(&_mesh.getMesh(),
150 : primary_elems_to_ghost.begin(),
151 : primary_elems_to_ghost.end(),
152 : null_output_iterator<Elem>());
153 :
154 0 : _mesh.update(); // Rebuild node_to_elem_map
155 :
156 : // Find elems again now that we know they're there
157 0 : const auto & new_node_to_elem_map = _mesh.nodeToElemMap();
158 : auto node_to_elem_pair = new_node_to_elem_map.find(_primary_node_vector[0]);
159 : bool found_elems = (node_to_elem_pair != new_node_to_elem_map.end());
160 :
161 0 : if (!found_elems)
162 0 : mooseError("Colundn't find any elements connected to primary node.");
163 :
164 0 : for (unsigned int i = 0; i < _primary_node_vector.size(); ++i)
165 0 : elems[i] = node_to_elem_pair->second;
166 : }
167 : else // serial mesh
168 : {
169 42 : for (unsigned int i = 0; i < _primary_node_vector.size(); ++i)
170 : {
171 : auto node_to_elem_pair = node_to_elem_map.find(_primary_node_vector[i]);
172 : bool found_elems = (node_to_elem_pair != node_to_elem_map.end());
173 :
174 28 : if (!found_elems)
175 0 : mooseError("Couldn't find any elements connected to primary node");
176 :
177 28 : elems[i] = node_to_elem_pair->second;
178 : }
179 : }
180 :
181 42 : for (unsigned int i = 0; i < _primary_node_vector.size(); ++i)
182 : {
183 28 : if (elems[i].size() == 0)
184 0 : mooseError("Couldn't find any elements connected to primary node");
185 :
186 56 : for (unsigned int j = 0; j < elems[i].size(); ++j)
187 28 : _subproblem.addGhostedElem(elems[i][j]);
188 : }
189 :
190 : // Cache map between secondary node and primary node
191 14 : _connected_nodes.clear();
192 14 : _primary_conn.clear();
193 28 : for (unsigned int j = 0; j < secondary_nodelist.size(); ++j)
194 : {
195 14 : if (_mesh.nodeRef(secondary_nodelist[j]).processor_id() == _subproblem.processor_id())
196 : {
197 10 : Node & secondary_node = _mesh.nodeRef(secondary_nodelist[j]);
198 20 : for (unsigned int i = 0; i < _primary_node_vector.size(); ++i)
199 : {
200 20 : Node & primary_node = _mesh.nodeRef(_primary_node_vector[i]);
201 20 : Real d = (secondary_node - primary_node).norm();
202 20 : if (MooseUtils::absoluteFuzzyEqual(d, 0.0))
203 : {
204 10 : _primary_conn.push_back(i);
205 10 : _connected_nodes.push_back(secondary_nodelist[j]);
206 : break;
207 : }
208 : }
209 : }
210 : }
211 :
212 14 : _console << "total secondary nodes, primary nodes: " << _primary_conn.size() << ", "
213 14 : << _primary_node_vector.size() << '\n';
214 14 : }
215 :
216 : void
217 713 : NodalFrictionalConstraint::computeResidual(NumericVector<Number> &
218 : /*residual*/)
219 : {
220 713 : const auto & primarydof = _var.dofIndices();
221 : const auto & secondarydof = _var.dofIndicesNeighbor();
222 713 : std::vector<Number> re(primarydof.size());
223 713 : std::vector<Number> neighbor_re(secondarydof.size());
224 :
225 713 : for (_i = 0; _i < secondarydof.size(); ++_i)
226 : {
227 713 : _j = _primary_conn[_i];
228 713 : re[_j] += computeQpResidual(Moose::Primary);
229 713 : neighbor_re[_i] += computeQpResidual(Moose::Secondary);
230 713 : break;
231 : }
232 713 : addResiduals(_assembly, re, primarydof, _var.scalingFactor());
233 713 : addResiduals(_assembly, neighbor_re, secondarydof, _var.scalingFactor());
234 713 : }
235 :
236 : Real
237 1426 : NodalFrictionalConstraint::computeQpResidual(Moose::ConstraintType type)
238 : {
239 : // check whether the tangential spring is already in the yielded state
240 1426 : Real old_force = (_u_secondary_old[_i] - _u_primary_old[_j]) * _tangential_penalty;
241 1426 : if (MooseUtils::absoluteFuzzyGreaterThan(std::abs(old_force),
242 1426 : _friction_coefficient * _normal_force))
243 620 : old_force = _friction_coefficient * _normal_force * old_force / std::abs(old_force);
244 :
245 : Real current_force =
246 1426 : ((_u_secondary[_i] - _u_secondary_old[_i]) - (_u_primary[_j] - _u_primary_old[_j])) *
247 : _tangential_penalty +
248 1426 : old_force;
249 1426 : if (MooseUtils::absoluteFuzzyGreaterThan(std::abs(current_force),
250 : _friction_coefficient * _normal_force))
251 1052 : current_force = _friction_coefficient * _normal_force * current_force / std::abs(current_force);
252 :
253 1426 : switch (type)
254 : {
255 : case Moose::Secondary:
256 : return current_force;
257 713 : case Moose::Primary:
258 713 : return -current_force;
259 : }
260 0 : return 0;
261 : }
262 :
263 : void
264 156 : NodalFrictionalConstraint::computeJacobian(SparseMatrix<Number> & /*jacobian*/)
265 : {
266 : // Calculate Jacobian entries and cache those entries along with the row and column indices
267 156 : std::vector<dof_id_type> secondarydof = _var.dofIndicesNeighbor();
268 156 : std::vector<dof_id_type> primarydof = _var.dofIndices();
269 :
270 156 : DenseMatrix<Number> Kee(primarydof.size(), primarydof.size());
271 156 : DenseMatrix<Number> Ken(primarydof.size(), secondarydof.size());
272 156 : DenseMatrix<Number> Kne(secondarydof.size(), primarydof.size());
273 156 : DenseMatrix<Number> Knn(secondarydof.size(), secondarydof.size());
274 :
275 : Kee.zero();
276 : Ken.zero();
277 : Kne.zero();
278 : Knn.zero();
279 :
280 312 : for (_i = 0; _i < secondarydof.size(); ++_i)
281 : {
282 156 : _j = _primary_conn[_i];
283 156 : Kee(_j, _j) += computeQpJacobian(Moose::PrimaryPrimary);
284 156 : Ken(_j, _i) += computeQpJacobian(Moose::PrimarySecondary);
285 156 : Kne(_i, _j) += computeQpJacobian(Moose::SecondaryPrimary);
286 156 : Knn(_i, _i) += computeQpJacobian(Moose::SecondarySecondary);
287 : }
288 156 : addJacobian(_assembly, Kee, primarydof, primarydof, _var.scalingFactor());
289 156 : addJacobian(_assembly, Ken, primarydof, secondarydof, _var.scalingFactor());
290 156 : addJacobian(_assembly, Kne, secondarydof, primarydof, _var.scalingFactor());
291 156 : addJacobian(_assembly, Knn, secondarydof, secondarydof, _var.scalingFactor());
292 156 : }
293 :
294 : Real
295 624 : NodalFrictionalConstraint::computeQpJacobian(Moose::ConstraintJacobianType type)
296 : {
297 624 : Real jac = _tangential_penalty;
298 :
299 : // set jacobian to zero if spring has yielded
300 624 : Real old_force = (_u_secondary_old[_i] - _u_primary_old[_j]) * _tangential_penalty;
301 624 : if (MooseUtils::absoluteFuzzyGreaterThan(std::abs(old_force),
302 624 : _friction_coefficient * _normal_force))
303 404 : old_force = _friction_coefficient * _normal_force * old_force / std::abs(old_force);
304 :
305 : Real current_force =
306 624 : ((_u_secondary[_i] - _u_secondary_old[_i]) - (_u_primary[_j] - _u_primary_old[_j])) *
307 : _tangential_penalty +
308 624 : old_force;
309 624 : if (MooseUtils::absoluteFuzzyGreaterThan(std::abs(current_force),
310 : _friction_coefficient * _normal_force))
311 : jac = 0.0;
312 :
313 624 : switch (type)
314 : {
315 : case Moose::SecondarySecondary:
316 : return jac;
317 156 : case Moose::SecondaryPrimary:
318 156 : return -jac;
319 : case Moose::PrimaryPrimary:
320 : return jac;
321 156 : case Moose::PrimarySecondary:
322 156 : return -jac;
323 0 : default:
324 0 : mooseError("Invalid type");
325 : }
326 : return 0.;
327 : }
|
