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