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 224 : NodalStickConstraint::validParams()
31 : {
32 224 : InputParameters params = NodalConstraint::validParams();
33 224 : params.addClassDescription("Sticky nodal constraint for contact");
34 448 : params.addRequiredParam<BoundaryName>("boundary", "The primary boundary");
35 448 : params.addRequiredParam<BoundaryName>("secondary", "The secondary boundary");
36 448 : params.addRequiredParam<Real>("penalty", "Stiffness of the spring.");
37 224 : return params;
38 0 : }
39 :
40 112 : NodalStickConstraint::NodalStickConstraint(const InputParameters & parameters)
41 : : NodalConstraint(parameters),
42 224 : _primary_boundary_id(getParam<BoundaryName>("boundary")),
43 112 : _secondary_boundary_id(getParam<BoundaryName>("secondary")),
44 336 : _penalty(getParam<Real>("penalty"))
45 : {
46 112 : 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 112 : updateConstrainedNodes();
52 112 : }
53 :
54 : void
55 0 : NodalStickConstraint::meshChanged()
56 : {
57 0 : updateConstrainedNodes();
58 0 : }
59 :
60 : void
61 112 : NodalStickConstraint::updateConstrainedNodes()
62 : {
63 112 : _primary_node_vector.clear();
64 112 : _connected_nodes.clear();
65 112 : _primary_conn.clear();
66 :
67 : std::vector<dof_id_type> secondary_nodelist =
68 112 : _mesh.getNodeList(_mesh.getBoundaryID(_secondary_boundary_id));
69 : std::vector<dof_id_type> primary_nodelist =
70 112 : _mesh.getNodeList(_mesh.getBoundaryID(_primary_boundary_id));
71 :
72 : // Fill in _connected_nodes, which defines secondary nodes in the base class
73 224 : for (auto in : secondary_nodelist)
74 : {
75 112 : if (_mesh.nodeRef(in).processor_id() == _subproblem.processor_id())
76 80 : _connected_nodes.push_back(in);
77 : }
78 :
79 : // Fill in _primary_node_vector, which defines secondary nodes in the base class
80 336 : for (auto in : primary_nodelist)
81 224 : _primary_node_vector.push_back(in);
82 :
83 112 : const auto & node_to_elem_map = _mesh.nodeToElemMap();
84 112 : 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 112 : 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 336 : 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 224 : if (!found_elems)
156 0 : mooseError("Couldn't find any elements connected to primary node");
157 :
158 224 : elems[i] = node_to_elem_pair->second;
159 : }
160 : }
161 :
162 336 : for (unsigned int i = 0; i < _primary_node_vector.size(); ++i)
163 : {
164 224 : if (elems[i].size() == 0)
165 0 : mooseError("Couldn't find any elements connected to primary node");
166 :
167 448 : for (unsigned int j = 0; j < elems[i].size(); ++j)
168 224 : _subproblem.addGhostedElem(elems[i][j]);
169 : }
170 :
171 : // Cache map between secondary node and primary node
172 112 : _connected_nodes.clear();
173 112 : _primary_conn.clear();
174 224 : for (unsigned int j = 0; j < secondary_nodelist.size(); ++j)
175 : {
176 112 : if (_mesh.nodeRef(secondary_nodelist[j]).processor_id() == _subproblem.processor_id())
177 : {
178 80 : Node & secondary_node = _mesh.nodeRef(secondary_nodelist[j]);
179 160 : for (unsigned int i = 0; i < _primary_node_vector.size(); ++i)
180 : {
181 160 : Node & primary_node = _mesh.nodeRef(_primary_node_vector[i]);
182 160 : Real d = (secondary_node - primary_node).norm();
183 160 : if (MooseUtils::absoluteFuzzyEqual(d, 0.0))
184 : {
185 80 : _primary_conn.push_back(i);
186 80 : _connected_nodes.push_back(secondary_nodelist[j]);
187 : break;
188 : }
189 : }
190 : }
191 : }
192 112 : }
193 :
194 : void
195 1110 : NodalStickConstraint::computeJacobian(SparseMatrix<Number> & jacobian)
196 : {
197 : // Calculate Jacobian enteries and cache those entries along with the row and column indices
198 1110 : std::vector<dof_id_type> secondarydof = _var.dofIndicesNeighbor();
199 2220 : std::vector<dof_id_type> primarydof = _var.dofIndices();
200 :
201 1110 : DenseMatrix<Number> Kee(primarydof.size(), primarydof.size());
202 1110 : DenseMatrix<Number> Ken(primarydof.size(), secondarydof.size());
203 1110 : DenseMatrix<Number> Kne(secondarydof.size(), primarydof.size());
204 1110 : DenseMatrix<Number> Knn(secondarydof.size(), secondarydof.size());
205 :
206 : Kee.zero();
207 : Ken.zero();
208 : Kne.zero();
209 : Knn.zero();
210 :
211 2220 : for (_i = 0; _i < secondarydof.size(); ++_i)
212 : {
213 1110 : _j = _primary_conn[_i];
214 1110 : switch (_formulation)
215 : {
216 1110 : case Moose::Penalty:
217 1110 : Kee(_j, _j) += computeQpJacobian(Moose::PrimaryPrimary);
218 1110 : Ken(_j, _i) += computeQpJacobian(Moose::PrimarySecondary);
219 1110 : Kne(_i, _j) += computeQpJacobian(Moose::SecondaryPrimary);
220 1110 : Knn(_i, _i) += computeQpJacobian(Moose::SecondarySecondary);
221 1110 : 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 1110 : addJacobian(_assembly, Kee, primarydof, primarydof, _var.scalingFactor());
233 1110 : addJacobian(_assembly, Ken, primarydof, secondarydof, _var.scalingFactor());
234 1110 : addJacobian(_assembly, Kne, secondarydof, primarydof, _var.scalingFactor());
235 1110 : addJacobian(_assembly, Knn, secondarydof, secondarydof, _var.scalingFactor());
236 1110 : }
237 :
238 : void
239 5653 : NodalStickConstraint::computeResidual(NumericVector<Number> & residual)
240 : {
241 5653 : std::vector<dof_id_type> primarydof = _var.dofIndices();
242 5653 : std::vector<dof_id_type> secondarydof = _var.dofIndicesNeighbor();
243 5653 : DenseVector<Number> re(primarydof.size());
244 5653 : DenseVector<Number> neighbor_re(secondarydof.size());
245 :
246 : re.zero();
247 : neighbor_re.zero();
248 11306 : for (_i = 0; _i < secondarydof.size(); ++_i)
249 : {
250 5653 : _j = _primary_conn[_i];
251 5653 : switch (_formulation)
252 : {
253 5653 : case Moose::Penalty:
254 5653 : re(_j) += computeQpResidual(Moose::Primary) * _var.scalingFactor();
255 5653 : neighbor_re(_i) += computeQpResidual(Moose::Secondary) * _var.scalingFactor();
256 5653 : 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 5653 : addResiduals(_assembly, re, primarydof, /*scaling_factor=*/1);
267 5653 : addResiduals(_assembly, neighbor_re, secondarydof, /*scaling_fator=*/1);
268 5653 : }
269 :
270 : Real
271 11306 : NodalStickConstraint::computeQpResidual(Moose::ConstraintType type)
272 : {
273 11306 : switch (type)
274 : {
275 5653 : case Moose::Secondary:
276 5653 : return (_u_secondary[_i] - _u_primary[_j]) * _penalty;
277 5653 : case Moose::Primary:
278 5653 : return (_u_primary[_j] - _u_secondary[_i]) * _penalty;
279 : }
280 : return 0.;
281 : }
282 :
283 : Real
284 4440 : NodalStickConstraint::computeQpJacobian(Moose::ConstraintJacobianType type)
285 : {
286 4440 : switch (type)
287 : {
288 1110 : case Moose::SecondarySecondary:
289 1110 : return _penalty;
290 1110 : case Moose::SecondaryPrimary:
291 1110 : return -_penalty;
292 1110 : case Moose::PrimaryPrimary:
293 1110 : return _penalty;
294 1110 : case Moose::PrimarySecondary:
295 1110 : return -_penalty;
296 0 : default:
297 0 : mooseError("Invalid type");
298 : }
299 : return 0.;
300 : }
|
