www.mooseframework.org
MultiDContactConstraint.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
8 // MOOSE includes
10 #include "SystemBase.h"
11 #include "PenetrationLocator.h"
12 #include "MooseMesh.h"
13 
14 #include "libmesh/string_to_enum.h"
15 #include "libmesh/sparse_matrix.h"
16 
17 template <>
18 InputParameters
20 {
21  InputParameters params = validParams<NodeFaceConstraint>();
22  params.set<bool>("use_displaced_mesh") = true;
23  params.addParam<bool>("jacobian_update",
24  false,
25  "Whether or not to update the 'in contact' list "
26  "every jacobian evaluation (by default it will "
27  "happen once per timestep");
28 
29  params.addRequiredParam<unsigned int>("component",
30  "An integer corresponding to the direction "
31  "the variable this kernel acts in. (0 for x, "
32  "1 for y, 2 for z)");
33  params.addCoupledVar("disp_x", "The x displacement");
34  params.addCoupledVar("disp_y", "The y displacement");
35  params.addCoupledVar("disp_z", "The z displacement");
36 
37  params.addCoupledVar(
38  "displacements",
39  "The displacements appropriate for the simulation geometry and coordinate system");
40 
41  params.addParam<std::string>("model", "frictionless", "The contact model to use");
42  params.addParam<Real>(
43  "penalty",
44  1e8,
45  "The penalty to apply. This can vary depending on the stiffness of your materials");
46 
47  // TODO: Reenable this
48  // params.addParam<std::string>("order", "FIRST", "The finite element order");
49 
50  return params;
51 }
52 
53 MultiDContactConstraint::MultiDContactConstraint(const InputParameters & parameters)
54  : NodeFaceConstraint(parameters),
55  _residual_copy(_sys.residualGhosted()),
56  _jacobian_update(getParam<bool>("jacobian_update")),
57  _component(getParam<unsigned int>("component")),
58  _model(ContactMaster::contactModel(getParam<std::string>("model"))),
59  _penalty(getParam<Real>("penalty")),
60  _mesh_dimension(_mesh.dimension()),
61  _vars(3, libMesh::invalid_uint)
62 {
63  _overwrite_slave_residual = false;
64 
65  if (isParamValid("displacements"))
66  {
67  // modern parameter scheme for displacements
68  for (unsigned int i = 0; i < coupledComponents("displacements"); ++i)
69  _vars[i] = coupled("displacements", i);
70  }
71  else
72  {
73  // Legacy parameter scheme for displacements
74  if (isParamValid("disp_x"))
75  _vars[0] = coupled("disp_x");
76  if (isParamValid("disp_y"))
77  _vars[1] = coupled("disp_y");
78  if (isParamValid("disp_z"))
79  _vars[2] = coupled("disp_z");
80 
81  mooseDeprecated("use the `displacements` parameter rather than the `disp_*` parameters (those "
82  "will go away with the deprecation of the Solid Mechanics module).");
83  }
84 }
85 
86 void
88 {
89  if (_component == 0)
91 }
92 
93 void
95 {
96  if (_component == 0)
98 }
99 
100 void
102 {
103  auto & has_penetrated = _penetration_locator._has_penetrated;
104 
105  for (auto & pinfo_pair : _penetration_locator._penetration_info)
106  {
107  PenetrationInfo * pinfo = pinfo_pair.second;
108 
109  // Skip this pinfo if there are no DOFs on this node.
110  if (!pinfo || pinfo->_node->n_comp(_sys.number(), _vars[_component]) < 1)
111  continue;
112 
113  const Node * node = pinfo->_node;
114 
115  dof_id_type slave_node_num = pinfo_pair.first;
116  auto hpit = has_penetrated.find(slave_node_num);
117 
118  RealVectorValue res_vec;
119  // Build up residual vector
120  for (unsigned int i = 0; i < _mesh_dimension; ++i)
121  {
122  dof_id_type dof_number = node->dof_number(0, _vars[i], 0);
123  res_vec(i) = _residual_copy(dof_number);
124  }
125 
126  // Real resid = 0;
127  switch (_model)
128  {
129  case CM_FRICTIONLESS:
130  // resid = pinfo->_normal * res_vec;
131  break;
132 
133  case CM_GLUED:
134  // resid = pinfo->_normal * res_vec;
135  break;
136 
137  default:
138  mooseError("Invalid or unavailable contact model");
139  break;
140  }
141 
142  // if (hpit != has_penetrated.end() && resid < 0)
143  // Moose::err << resid << std::endl;
144  //
145  // if (hpit != has_penetrated.end() && resid < -.15)
146  // {
147  // Moose::err << std::endl
148  // << "Unlocking node " << node->id()
149  // << " because resid:
150  // "<<resid<<std::endl<<std::endl;
151  //
152  // has_penetrated.erase(hpit);
153  // unlocked_this_step[slave_node_num] = true;
154  // }
155  // else
156 
157  if (pinfo->_distance > 0 &&
158  hpit == has_penetrated.end()) // && !unlocked_this_step[slave_node_num])
159  {
160  // Moose::err << std::endl
161  // << "Locking node " << node->id() << " because distance:" << pinfo->_distance
162  // << std::endl
163  // << std::endl;
164 
165  has_penetrated.insert(slave_node_num);
166  }
167  }
168 }
169 
170 bool
172 {
173  std::set<dof_id_type>::iterator hpit =
174  _penetration_locator._has_penetrated.find(_current_node->id());
175  return (hpit != _penetration_locator._has_penetrated.end());
176 }
177 
178 Real
180 {
181  PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
182 
183  // Moose::err << std::endl
184  // << "Popping out node: " << _current_node->id() << std::endl
185  // << "Closest Point " << _component << ": " << pinfo->_closest_point(_component)
186  // << std::endl
187  // << "Current Node " << _component << ": " << (*_current_node)(_component) <<
188  // std::endl
189  // << "Current Value: " << _u_slave[_qp] << std::endl
190  // << "New Value: "
191  // << pinfo->_closest_point(_component) - ((*_current_node)(_component)-_u_slave[_qp])
192  // << std::endl
193  // << "Change: "
194  // << _u_slave[_qp] - (pinfo->_closest_point(_component) -
195  // ((*_current_node)(_component)-_u_slave[_qp]))
196  // << std::endl
197  // << std::endl;
198 
199  return pinfo->_closest_point(_component) - ((*_current_node)(_component)-_u_slave[_qp]);
200 }
201 
202 Real
204 {
205  PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
206  const Node * node = pinfo->_node;
207 
208  RealVectorValue res_vec;
209  // Build up residual vector
210  for (unsigned int i = 0; i < _mesh_dimension; ++i)
211  {
212  dof_id_type dof_number = node->dof_number(0, _vars[i], 0);
213  res_vec(i) = _residual_copy(dof_number);
214  }
215 
216  const RealVectorValue distance_vec(_mesh.nodeRef(node->id()) - pinfo->_closest_point);
217  const RealVectorValue pen_force(_penalty * distance_vec);
218  Real resid = 0.0;
219 
220  switch (type)
221  {
222  case Moose::Slave:
223  switch (_model)
224  {
225  case CM_FRICTIONLESS:
226  resid = pinfo->_normal(_component) * (pinfo->_normal * (pen_force - res_vec));
227  break;
228 
229  case CM_GLUED:
230  resid = pen_force(_component) - res_vec(_component);
231  break;
232 
233  default:
234  mooseError("Invalid or unavailable contact model");
235  break;
236  }
237  return _test_slave[_i][_qp] * resid;
238  case Moose::Master:
239  switch (_model)
240  {
241  case CM_FRICTIONLESS:
242  resid = pinfo->_normal(_component) * (pinfo->_normal * res_vec);
243  break;
244 
245  case CM_GLUED:
246  resid = res_vec(_component);
247  break;
248 
249  default:
250  mooseError("Invalid or unavailable contact model");
251  break;
252  }
253  return _test_master[_i][_qp] * resid;
254  }
255  return 0.0;
256 }
257 
258 Real
259 MultiDContactConstraint::computeQpJacobian(Moose::ConstraintJacobianType type)
260 {
261  PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
262 
263  Real slave_jac = 0.0;
264  switch (type)
265  {
266  case Moose::SlaveSlave:
267  switch (_model)
268  {
269  case CM_FRICTIONLESS:
270 
271  slave_jac = pinfo->_normal(_component) * pinfo->_normal(_component) *
272  (_penalty * _phi_slave[_j][_qp] -
273  (*_jacobian)(_current_node->dof_number(0, _var.number(), 0),
274  _connected_dof_indices[_j]));
275  break;
276 
277  case CM_GLUED:
278  // resid = pen_force(_component) - res_vec(_component);
279  break;
280 
281  default:
282  mooseError("Invalid or unavailable contact model");
283  break;
284  }
285  return _test_slave[_i][_qp] * slave_jac;
286 
287  case Moose::SlaveMaster:
288  switch (_model)
289  {
290  case CM_FRICTIONLESS:
291 
292  slave_jac = pinfo->_normal(_component) * pinfo->_normal(_component) *
293  (-_penalty * _phi_master[_j][_qp]);
294  break;
295 
296  case CM_GLUED:
297  /*
298  resid = pen_force(_component)
299  - res_vec(_component)
300  ;
301  */
302  break;
303 
304  default:
305  mooseError("Invalid or unavailable contact model");
306  break;
307  }
308  return _test_slave[_i][_qp] * slave_jac;
309 
310  case Moose::MasterSlave:
311  slave_jac =
312  (*_jacobian)(_current_node->dof_number(0, _var.number(), 0), _connected_dof_indices[_j]);
313  return slave_jac * _test_master[_i][_qp];
314 
315  case Moose::MasterMaster:
316  return 0.0;
317  }
318  return 0.0;
319 }
NumericVector< Number > & _residual_copy
virtual Real computeQpResidual(Moose::ConstraintType type)
MultiDContactConstraint(const InputParameters &parameters)
virtual Real computeQpJacobian(Moose::ConstraintJacobianType type)
InputParameters validParams< MultiDContactConstraint >()
std::vector< unsigned int > _vars
ContactModel contactModel(const std::string &the_name)
const unsigned int _mesh_dimension