www.mooseframework.org
SlaveConstraint.C
Go to the documentation of this file.
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 #include "SlaveConstraint.h"
11 
12 // Moose includes
13 #include "SystemBase.h"
14 #include "MooseMesh.h"
15 #include "ContactAction.h"
16 
17 #include "libmesh/plane.h"
18 #include "libmesh/sparse_matrix.h"
19 #include "libmesh/string_to_enum.h"
20 
21 registerMooseObject("ContactApp", SlaveConstraint);
22 
23 template <>
24 InputParameters
25 validParams<SlaveConstraint>()
26 {
27  InputParameters params = validParams<DiracKernel>();
28  params += ContactAction::commonParameters();
29 
30  params.addRequiredParam<BoundaryName>("boundary", "The slave boundary");
31  params.addRequiredParam<BoundaryName>("master", "The master boundary");
32  params.addRequiredParam<unsigned int>("component",
33  "An integer corresponding to the direction "
34  "the variable this kernel acts in. (0 for x, "
35  "1 for y, 2 for z)");
36 
37  params.addCoupledVar("disp_x", "The x displacement");
38  params.addCoupledVar("disp_y", "The y displacement");
39  params.addCoupledVar("disp_z", "The z displacement");
40 
41  params.addCoupledVar(
42  "displacements",
43  "The displacements appropriate for the simulation geometry and coordinate system");
44 
45  params.addRequiredCoupledVar("nodal_area", "The nodal area");
46 
47  params.set<bool>("use_displaced_mesh") = true;
48  params.addParam<Real>(
49  "penalty",
50  1e8,
51  "The penalty to apply. This can vary depending on the stiffness of your materials");
52  params.addParam<Real>("friction_coefficient", 0, "The friction coefficient");
53  params.addParam<Real>("tangential_tolerance",
54  "Tangential distance to extend edges of contact surfaces");
55  params.addParam<bool>(
56  "normalize_penalty",
57  false,
58  "Whether to normalize the penalty parameter with the nodal area for penalty contact.");
59  return params;
60 }
61 
62 SlaveConstraint::SlaveConstraint(const InputParameters & parameters)
63  : DiracKernel(parameters),
64  _component(getParam<unsigned int>("component")),
65  _model(getParam<MooseEnum>("model").getEnum<ContactModel>()),
66  _formulation(getParam<MooseEnum>("formulation").getEnum<ContactFormulation>()),
67  _normalize_penalty(getParam<bool>("normalize_penalty")),
68  _penetration_locator(
69  getPenetrationLocator(getParam<BoundaryName>("master"),
70  getParam<BoundaryName>("boundary"),
71  Utility::string_to_enum<Order>(getParam<MooseEnum>("order")))),
72  _penalty(getParam<Real>("penalty")),
73  _friction_coefficient(getParam<Real>("friction_coefficient")),
74  _residual_copy(_sys.residualGhosted()),
75  _vars(3, libMesh::invalid_uint),
76  _mesh_dimension(_mesh.dimension()),
77  _nodal_area_var(getVar("nodal_area", 0)),
78  _aux_system(_nodal_area_var->sys()),
79  _aux_solution(_aux_system.currentSolution())
80 {
81  if (isParamValid("displacements"))
82  {
83  // modern parameter scheme for displacements
84  for (unsigned int i = 0; i < coupledComponents("displacements"); ++i)
85  _vars[i] = coupled("displacements", i);
86  }
87  else
88  {
89  // Legacy parameter scheme for displacements
90  if (isParamValid("disp_x"))
91  _vars[0] = coupled("disp_x");
92  if (isParamValid("disp_y"))
93  _vars[1] = coupled("disp_y");
94  if (isParamValid("disp_z"))
95  _vars[2] = coupled("disp_z");
96 
97  mooseDeprecated("use the `displacements` parameter rather than the `disp_*` parameters (those "
98  "will go away with the deprecation of the Solid Mechanics module).");
99  }
100 
101  if (parameters.isParamValid("tangential_tolerance"))
102  _penetration_locator.setTangentialTolerance(getParam<Real>("tangential_tolerance"));
103 
104  if (parameters.isParamValid("normal_smoothing_distance"))
105  _penetration_locator.setNormalSmoothingDistance(getParam<Real>("normal_smoothing_distance"));
106 
107  if (parameters.isParamValid("normal_smoothing_method"))
108  _penetration_locator.setNormalSmoothingMethod(
109  parameters.get<std::string>("normal_smoothing_method"));
110 }
111 
112 void
113 SlaveConstraint::addPoints()
114 {
115  _point_to_info.clear();
116 
117  std::map<dof_id_type, PenetrationInfo *>::iterator
118  it = _penetration_locator._penetration_info.begin(),
119  end = _penetration_locator._penetration_info.end();
120 
121  const auto & node_to_elem_map = _mesh.nodeToElemMap();
122  for (; it != end; ++it)
123  {
124  PenetrationInfo * pinfo = it->second;
125 
126  // Skip this pinfo if there are no DOFs on this node.
127  if (!pinfo || pinfo->_node->n_comp(_sys.number(), _vars[_component]) < 1)
128  continue;
129 
130  dof_id_type slave_node_num = it->first;
131  const Node * node = pinfo->_node;
132 
133  if (pinfo->isCaptured() && node->processor_id() == processor_id())
134  {
135  // Find an element that is connected to this node that and that is also on this processor
136  auto node_to_elem_pair = node_to_elem_map.find(slave_node_num);
137  mooseAssert(node_to_elem_pair != node_to_elem_map.end(), "Missing node in node to elem map");
138  const std::vector<dof_id_type> & connected_elems = node_to_elem_pair->second;
139 
140  Elem * elem = NULL;
141 
142  for (unsigned int i = 0; i < connected_elems.size() && !elem; ++i)
143  {
144  Elem * cur_elem = _mesh.elemPtr(connected_elems[i]);
145  if (cur_elem->processor_id() == processor_id())
146  elem = cur_elem;
147  }
148 
149  mooseAssert(elem,
150  "Couldn't find an element on this processor that is attached to the slave node!");
151 
152  addPoint(elem, *node);
153  _point_to_info[*node] = pinfo;
154  }
155  }
156 }
157 
158 Real
159 SlaveConstraint::computeQpResidual()
160 {
161  PenetrationInfo * pinfo = _point_to_info[_current_point];
162  const Node * node = pinfo->_node;
163 
164  Real resid = pinfo->_contact_force(_component);
165 
166  const Real area = nodalArea(*pinfo);
167 
168  if (_formulation == ContactFormulation::KINEMATIC)
169  {
170  RealVectorValue distance_vec(_mesh.nodeRef(node->id()) - pinfo->_closest_point);
171  RealVectorValue pen_force(_penalty * distance_vec);
172  if (_normalize_penalty)
173  pen_force *= area;
174 
175  if (_model == ContactModel::FRICTIONLESS)
176  resid += pinfo->_normal(_component) * pinfo->_normal * pen_force;
177 
178  else if (_model == ContactModel::GLUED || _model == ContactModel::COULOMB)
179  resid += pen_force(_component);
180  }
181 
182  return _test[_i][_qp] * resid;
183 }
184 
185 Real
186 SlaveConstraint::computeQpJacobian()
187 {
188 
189  // TODO: for the default formulation,
190  // we should subtract off the existing Jacobian weighted by the effect of the normal
191 
192  PenetrationInfo * pinfo = _point_to_info[_current_point];
193  const Node * node = pinfo->_node;
194  // long int dof_number = node->dof_number(0, _var.number(), 0);
195 
196  // RealVectorValue jac_vec;
197 
198  // Build up jac vector
199  // for (unsigned int i=0; i<_dim; i++)
200  // {
201  // unsigned int dof_row = _dof_data._var_dof_indices[_var_num][_i];
202  // unsigned int dof_col = _dof_data._var_dof_indices[_var_num][_j];
203 
204  // Real jac_value = _jacobian_copy(dof_row, dof_col);
205  // }
206 
207  // Real jac_mag = pinfo->_normal(_component) * jac_value;
208  /*
209  return _test[_i][_qp] * (
210  (1e8*-_phi[_j][_qp])
211  -_jacobian_copy(dof_number, dof_number)
212  );
213  */
214 
215  RealVectorValue normal(pinfo->_normal);
216 
217  Real penalty = _penalty;
218  if (_normalize_penalty)
219  penalty *= nodalArea(*pinfo);
220 
221  Real term(0);
222 
223  if (ContactModel::FRICTIONLESS == _model)
224  {
225 
226  const Real nnTDiag = normal(_component) * normal(_component);
227  term = penalty * nnTDiag;
228 
229  const RealGradient & A1(pinfo->_dxyzdxi[0]);
230  RealGradient A2;
231  RealGradient d2;
232  if (_mesh_dimension == 3)
233  {
234  A2 = pinfo->_dxyzdeta[0];
235  d2 = pinfo->_d2xyzdxideta[0];
236  }
237  else
238  {
239  A2.zero();
240  d2.zero();
241  }
242 
243  const RealVectorValue distance_vec(_mesh.nodeRef(node->id()) - pinfo->_closest_point);
244  const Real ATA11(A1 * A1);
245  const Real ATA12(A1 * A2);
246  const Real ATA22(A2 * A2);
247  const Real D11(-ATA11);
248  const Real D12(-ATA12 + d2 * distance_vec);
249  const Real D22(-ATA22);
250 
251  Real invD11(0);
252  Real invD12(0);
253  Real invD22(0);
254  if (_mesh_dimension == 3)
255  {
256  const Real detD(D11 * D22 - D12 * D12);
257  invD11 = D22 / detD;
258  invD12 = -D12 / detD;
259  invD22 = D11 / detD;
260  }
261  else if (_mesh_dimension == 2)
262  {
263  invD11 = 1 / D11;
264  }
265 
266  const Real AinvD11(A1(0) * invD11 + A2(0) * invD12);
267  const Real AinvD12(A1(0) * invD12 + A2(0) * invD22);
268  const Real AinvD21(A1(1) * invD11 + A2(1) * invD12);
269  const Real AinvD22(A1(1) * invD12 + A2(1) * invD22);
270  const Real AinvD31(A1(2) * invD11 + A2(2) * invD12);
271  const Real AinvD32(A1(2) * invD12 + A2(2) * invD22);
272 
273  const Real AinvDAT11(AinvD11 * A1(0) + AinvD12 * A2(0));
274  // const Real AinvDAT12( AinvD11*A1(1) + AinvD12*A2(1) );
275  // const Real AinvDAT13( AinvD11*A1(2) + AinvD12*A2(2) );
276  // const Real AinvDAT21( AinvD21*A1(0) + AinvD22*A2(0) );
277  const Real AinvDAT22(AinvD21 * A1(1) + AinvD22 * A2(1));
278  // const Real AinvDAT23( AinvD21*A1(2) + AinvD22*A2(2) );
279  // const Real AinvDAT31( AinvD31*A1(0) + AinvD32*A2(0) );
280  // const Real AinvDAT32( AinvD31*A1(1) + AinvD32*A2(1) );
281  const Real AinvDAT33(AinvD31 * A1(2) + AinvD32 * A2(2));
282 
283  if (_component == 0)
284  term += penalty * (1 - nnTDiag + AinvDAT11);
285 
286  else if (_component == 1)
287  term += penalty * (1 - nnTDiag + AinvDAT22);
288 
289  else
290  term += penalty * (1 - nnTDiag + AinvDAT33);
291  }
292  else if (ContactModel::GLUED == _model || ContactModel::COULOMB == _model)
293  {
294  normal.zero();
295  normal(_component) = 1;
296  term = penalty;
297  }
298  else
299  {
300  mooseError("Invalid or unavailable contact model");
301  }
302 
303  return _test[_i][_qp] * term * _phi[_j][_qp];
304 }
305 
306 Real
307 SlaveConstraint::nodalArea(PenetrationInfo & pinfo)
308 {
309  const Node * node = pinfo._node;
310 
311  dof_id_type dof = node->dof_number(_aux_system.number(), _nodal_area_var->number(), 0);
312 
313  Real area = (*_aux_solution)(dof);
314  if (area == 0)
315  {
316  if (_t_step > 1)
317  mooseError("Zero nodal area found");
318 
319  else
320  area = 1; // Avoid divide by zero during initialization
321  }
322  return area;
323 }
ContactModel
ContactModel
Definition: ContactAction.h:16
SlaveConstraint::nodalArea
Real nodalArea(PenetrationInfo &pinfo)
Definition: SlaveConstraint.C:307
SlaveConstraint::_aux_system
SystemBase & _aux_system
Definition: SlaveConstraint.h:54
SlaveConstraint::_component
const unsigned int _component
Definition: SlaveConstraint.h:36
SlaveConstraint::_point_to_info
std::map< Point, PenetrationInfo * > _point_to_info
Definition: SlaveConstraint.h:47
SlaveConstraint.h
libMesh::RealGradient
VectorValue< Real > RealGradient
Definition: GrainForceAndTorqueInterface.h:17
libMesh
Definition: RANFSNormalMechanicalContact.h:24
SlaveConstraint::_mesh_dimension
const unsigned int _mesh_dimension
Definition: SlaveConstraint.h:51
SlaveConstraint::_nodal_area_var
MooseVariable * _nodal_area_var
Definition: SlaveConstraint.h:53
ContactModel::GLUED
SlaveConstraint::_normalize_penalty
const bool _normalize_penalty
Definition: SlaveConstraint.h:39
validParams< SlaveConstraint >
InputParameters validParams< SlaveConstraint >()
Definition: SlaveConstraint.C:25
ContactFormulation::KINEMATIC
SlaveConstraint::_vars
std::vector< unsigned int > _vars
Definition: SlaveConstraint.h:49
ContactModel::COULOMB
SlaveConstraint::_model
const ContactModel _model
Definition: SlaveConstraint.h:37
ContactModel::FRICTIONLESS
SlaveConstraint
Definition: SlaveConstraint.h:24
ContactFormulation
ContactFormulation
Definition: ContactAction.h:23
SlaveConstraint::SlaveConstraint
SlaveConstraint(const InputParameters &parameters)
Definition: SlaveConstraint.C:62
SlaveConstraint::computeQpJacobian
virtual Real computeQpJacobian()
Definition: SlaveConstraint.C:186
SlaveConstraint::_penetration_locator
PenetrationLocator & _penetration_locator
Definition: SlaveConstraint.h:40
SlaveConstraint::_penalty
const Real _penalty
Definition: SlaveConstraint.h:42
SlaveConstraint::computeQpResidual
virtual Real computeQpResidual()
Definition: SlaveConstraint.C:159
SlaveConstraint::_formulation
const ContactFormulation _formulation
Definition: SlaveConstraint.h:38
ContactAction.h
ContactAction::commonParameters
static InputParameters commonParameters()
Definition: ContactAction.C:441
registerMooseObject
registerMooseObject("ContactApp", SlaveConstraint)
SlaveConstraint::addPoints
virtual void addPoints()
Definition: SlaveConstraint.C:113
Generated on Sat Jan 25 2020 12:04:36 for www.mooseframework.org by   doxygen 1.8.16