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