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