www.mooseframework.org
FindContactPoint.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
15 // Moose
16 #include "FindContactPoint.h"
17 #include "LineSegment.h"
18 #include "PenetrationInfo.h"
19 
20 // libMesh
21 #include "libmesh/boundary_info.h"
22 #include "libmesh/elem.h"
23 #include "libmesh/plane.h"
24 #include "libmesh/fe_interface.h"
25 #include "libmesh/dense_matrix.h"
26 #include "libmesh/dense_vector.h"
27 #include "libmesh/fe_base.h"
28 #include "libmesh/vector_value.h"
29 
30 namespace Moose
31 {
32 
48 void
50  FEBase * fe_elem,
51  FEBase * fe_side,
52  FEType & fe_side_type,
53  const Point & slave_point,
54  bool start_with_centroid,
55  const Real tangential_tolerance,
56  bool & contact_point_on_side)
57 {
58  const Elem * master_elem = p_info._elem;
59 
60  unsigned int dim = master_elem->dim();
61 
62  const Elem * side = p_info._side;
63 
64  const std::vector<Point> & phys_point = fe_side->get_xyz();
65 
66  const std::vector<RealGradient> & dxyz_dxi = fe_side->get_dxyzdxi();
67  const std::vector<RealGradient> & d2xyz_dxi2 = fe_side->get_d2xyzdxi2();
68  const std::vector<RealGradient> & d2xyz_dxieta = fe_side->get_d2xyzdxideta();
69 
70  const std::vector<RealGradient> & dxyz_deta = fe_side->get_dxyzdeta();
71  const std::vector<RealGradient> & d2xyz_deta2 = fe_side->get_d2xyzdeta2();
72  const std::vector<RealGradient> & d2xyz_detaxi = fe_side->get_d2xyzdxideta();
73 
74  if (dim == 1)
75  {
76  const Node * nearest_node = side->node_ptr(0);
77  p_info._closest_point = *nearest_node;
78  p_info._closest_point_ref =
79  master_elem->master_point(master_elem->get_node_index(nearest_node));
80  std::vector<Point> elem_points = {p_info._closest_point_ref};
81  fe_elem->reinit(master_elem, &elem_points);
82 
83  const std::vector<RealGradient> & elem_dxyz_dxi = fe_elem->get_dxyzdxi();
84  p_info._normal = elem_dxyz_dxi[0];
85  if (nearest_node->id() == master_elem->node_id(0))
86  p_info._normal *= -1.0;
87  p_info._normal /= p_info._normal.norm();
88 
89  Point from_slave_to_closest = p_info._closest_point - slave_point;
90  p_info._distance = from_slave_to_closest * p_info._normal;
91  Point tangential = from_slave_to_closest - p_info._distance * p_info._normal;
92  p_info._tangential_distance = tangential.norm();
93  p_info._dxyzdxi = dxyz_dxi;
94  p_info._dxyzdeta = dxyz_deta;
95  p_info._d2xyzdxideta = d2xyz_dxieta;
96  p_info._side_phi = fe_side->get_phi();
97  p_info._side_grad_phi = fe_side->get_dphi();
98  contact_point_on_side = true;
99  return;
100  }
101 
102  Point ref_point;
103 
104  if (start_with_centroid)
105  ref_point =
106  FEInterface::inverse_map(dim - 1, fe_side_type, side, side->centroid(), TOLERANCE, false);
107  else
108  ref_point = p_info._closest_point_ref;
109 
110  std::vector<Point> points = {ref_point};
111  fe_side->reinit(side, &points);
112  RealGradient d = slave_point - phys_point[0];
113 
114  Real update_size = std::numeric_limits<Real>::max();
115 
116  // Least squares
117  for (unsigned int it = 0; it < 3 && update_size > TOLERANCE * 1e3; ++it)
118  {
119  DenseMatrix<Real> jac(dim - 1, dim - 1);
120 
121  jac(0, 0) = -(dxyz_dxi[0] * dxyz_dxi[0]);
122 
123  if (dim - 1 == 2)
124  {
125  jac(1, 0) = -(dxyz_dxi[0] * dxyz_deta[0]);
126  jac(0, 1) = -(dxyz_deta[0] * dxyz_dxi[0]);
127  jac(1, 1) = -(dxyz_deta[0] * dxyz_deta[0]);
128  }
129 
130  DenseVector<Real> rhs(dim - 1);
131 
132  rhs(0) = dxyz_dxi[0] * d;
133 
134  if (dim - 1 == 2)
135  rhs(1) = dxyz_deta[0] * d;
136 
137  DenseVector<Real> update(dim - 1);
138 
139  jac.lu_solve(rhs, update);
140 
141  ref_point(0) -= update(0);
142 
143  if (dim - 1 == 2)
144  ref_point(1) -= update(1);
145 
146  points[0] = ref_point;
147  fe_side->reinit(side, &points);
148  d = slave_point - phys_point[0];
149 
150  update_size = update.l2_norm();
151  }
152 
153  update_size = std::numeric_limits<Real>::max();
154 
155  unsigned nit = 0;
156 
157  // Newton Loop
158  for (; nit < 12 && update_size > TOLERANCE * TOLERANCE; nit++)
159  {
160  d = slave_point - phys_point[0];
161 
162  DenseMatrix<Real> jac(dim - 1, dim - 1);
163 
164  jac(0, 0) = (d2xyz_dxi2[0] * d) - (dxyz_dxi[0] * dxyz_dxi[0]);
165 
166  if (dim - 1 == 2)
167  {
168  jac(1, 0) = (d2xyz_dxieta[0] * d) - (dxyz_dxi[0] * dxyz_deta[0]);
169 
170  jac(0, 1) = (d2xyz_detaxi[0] * d) - (dxyz_deta[0] * dxyz_dxi[0]);
171  jac(1, 1) = (d2xyz_deta2[0] * d) - (dxyz_deta[0] * dxyz_deta[0]);
172  }
173 
174  DenseVector<Real> rhs(dim - 1);
175 
176  rhs(0) = -dxyz_dxi[0] * d;
177 
178  if (dim - 1 == 2)
179  rhs(1) = -dxyz_deta[0] * d;
180 
181  DenseVector<Real> update(dim - 1);
182 
183  jac.lu_solve(rhs, update);
184 
185  ref_point(0) += update(0);
186 
187  if (dim - 1 == 2)
188  ref_point(1) += update(1);
189 
190  points[0] = ref_point;
191  fe_side->reinit(side, &points);
192  d = slave_point - phys_point[0];
193 
194  update_size = update.l2_norm();
195  }
196 
197  /*
198  if (nit == 12 && update_size > TOLERANCE*TOLERANCE)
199  Moose::err<<"Warning! Newton solve for contact point failed to converge!"<<std::endl;
200  */
201 
202  p_info._closest_point_ref = ref_point;
203  p_info._closest_point = phys_point[0];
204  p_info._distance = d.norm();
205 
206  if (dim - 1 == 2)
207  {
208  p_info._normal = dxyz_dxi[0].cross(dxyz_deta[0]);
209  p_info._normal /= p_info._normal.norm();
210  }
211  else
212  {
213  p_info._normal = RealGradient(dxyz_dxi[0](1), -dxyz_dxi[0](0));
214  if (std::fabs(p_info._normal.norm()) > 1e-15)
215  p_info._normal /= p_info._normal.norm();
216  }
217 
218  // If the point has not penetrated the face, make the distance negative
219  const Real dot(d * p_info._normal);
220  if (dot > 0.0)
221  p_info._distance = -p_info._distance;
222 
223  contact_point_on_side = FEInterface::on_reference_element(ref_point, side->type());
224 
225  p_info._tangential_distance = 0.0;
226 
227  if (!contact_point_on_side)
228  {
229  p_info._closest_point_on_face_ref = ref_point;
231 
232  points[0] = p_info._closest_point_on_face_ref;
233  fe_side->reinit(side, &points);
234  Point closest_point_on_face(phys_point[0]);
235 
236  RealGradient off_face = closest_point_on_face - p_info._closest_point;
237  Real tangential_distance = off_face.norm();
238  p_info._tangential_distance = tangential_distance;
239  if (tangential_distance <= tangential_tolerance)
240  {
241  contact_point_on_side = true;
242  }
243  }
244 
245  const std::vector<std::vector<Real>> & phi = fe_side->get_phi();
246  const std::vector<std::vector<RealGradient>> & grad_phi = fe_side->get_dphi();
247 
248  points[0] = p_info._closest_point_ref;
249  fe_side->reinit(side, &points);
250 
251  p_info._side_phi = phi;
252  p_info._side_grad_phi = grad_phi;
253  p_info._dxyzdxi = dxyz_dxi;
254  p_info._dxyzdeta = dxyz_deta;
255  p_info._d2xyzdxideta = d2xyz_dxieta;
256 }
257 
258 void
259 restrictPointToFace(Point & p, const Elem * side, std::vector<const Node *> & off_edge_nodes)
260 {
261  const ElemType t(side->type());
262  off_edge_nodes.clear();
263  Real & xi = p(0);
264  Real & eta = p(1);
265 
266  switch (t)
267  {
268  case EDGE2:
269  case EDGE3:
270  case EDGE4:
271  {
272  // The reference 1D element is [-1,1].
273  if (xi < -1.0)
274  {
275  xi = -1.0;
276  off_edge_nodes.push_back(side->node_ptr(0));
277  }
278  else if (xi > 1.0)
279  {
280  xi = 1.0;
281  off_edge_nodes.push_back(side->node_ptr(1));
282  }
283  break;
284  }
285 
286  case TRI3:
287  case TRI6:
288  {
289  // The reference triangle is isosceles
290  // and is bound by xi=0, eta=0, and xi+eta=1.
291 
292  if (xi <= 0.0 && eta <= 0.0)
293  {
294  xi = 0.0;
295  eta = 0.0;
296  off_edge_nodes.push_back(side->node_ptr(0));
297  }
298  else if (xi > 0.0 && xi < 1.0 && eta < 0.0)
299  {
300  eta = 0.0;
301  off_edge_nodes.push_back(side->node_ptr(0));
302  off_edge_nodes.push_back(side->node_ptr(1));
303  }
304  else if (eta > 0.0 && eta < 1.0 && xi < 0.0)
305  {
306  xi = 0.0;
307  off_edge_nodes.push_back(side->node_ptr(2));
308  off_edge_nodes.push_back(side->node_ptr(0));
309  }
310  else if (xi >= 1.0 && (eta - xi) <= -1.0)
311  {
312  xi = 1.0;
313  eta = 0.0;
314  off_edge_nodes.push_back(side->node_ptr(1));
315  }
316  else if (eta >= 1.0 && (eta - xi) >= 1.0)
317  {
318  xi = 0.0;
319  eta = 1.0;
320  off_edge_nodes.push_back(side->node_ptr(2));
321  }
322  else if ((xi + eta) > 1.0)
323  {
324  Real delta = (xi + eta - 1.0) / 2.0;
325  xi -= delta;
326  eta -= delta;
327  off_edge_nodes.push_back(side->node_ptr(1));
328  off_edge_nodes.push_back(side->node_ptr(2));
329  }
330  break;
331  }
332 
333  case QUAD4:
334  case QUAD8:
335  case QUAD9:
336  {
337  // The reference quadrilateral element is [-1,1]^2.
338  if (xi < -1.0)
339  {
340  xi = -1.0;
341  if (eta < -1.0)
342  {
343  eta = -1.0;
344  off_edge_nodes.push_back(side->node_ptr(0));
345  }
346  else if (eta > 1.0)
347  {
348  eta = 1.0;
349  off_edge_nodes.push_back(side->node_ptr(3));
350  }
351  else
352  {
353  off_edge_nodes.push_back(side->node_ptr(3));
354  off_edge_nodes.push_back(side->node_ptr(0));
355  }
356  }
357  else if (xi > 1.0)
358  {
359  xi = 1.0;
360  if (eta < -1.0)
361  {
362  eta = -1.0;
363  off_edge_nodes.push_back(side->node_ptr(1));
364  }
365  else if (eta > 1.0)
366  {
367  eta = 1.0;
368  off_edge_nodes.push_back(side->node_ptr(2));
369  }
370  else
371  {
372  off_edge_nodes.push_back(side->node_ptr(1));
373  off_edge_nodes.push_back(side->node_ptr(2));
374  }
375  }
376  else
377  {
378  if (eta < -1.0)
379  {
380  eta = -1.0;
381  off_edge_nodes.push_back(side->node_ptr(0));
382  off_edge_nodes.push_back(side->node_ptr(1));
383  }
384  else if (eta > 1.0)
385  {
386  eta = 1.0;
387  off_edge_nodes.push_back(side->node_ptr(2));
388  off_edge_nodes.push_back(side->node_ptr(3));
389  }
390  }
391  break;
392  }
393 
394  default:
395  {
396  mooseError("Unsupported face type: ", t);
397  break;
398  }
399  }
400 }
401 
402 } // namespace Moose
std::vector< RealGradient > _d2xyzdxideta
RealVectorValue RealGradient
Definition: Assembly.h:43
static PetscErrorCode Vec Mat jac
FEGenericBase< Real > FEBase
Definition: Assembly.h:34
RealVectorValue _normal
void restrictPointToFace(Point &p, const Elem *side, std::vector< const Node * > &off_edge_nodes)
Data structure used to hold penetration information.
void mooseError(Args &&...args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:182
Point _closest_point_on_face_ref
std::vector< std::vector< RealGradient > > _side_grad_phi
void findContactPoint(PenetrationInfo &p_info, FEBase *fe_elem, FEBase *fe_side, FEType &fe_side_type, const Point &slave_point, bool start_with_centroid, const Real tangential_tolerance, bool &contact_point_on_side)
Finds the closest point (called the contact point) on the master_elem on side "side" to the slave_poi...
std::vector< std::vector< Real > > _side_phi
const Elem * _elem
std::vector< RealGradient > _dxyzdxi
const Elem * _side
std::vector< const Node * > _off_edge_nodes
std::vector< RealGradient > _dxyzdeta
Definition: Moose.h:84