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
11 : #include "DenseMatrix.h"
12 : #include "FindContactPoint.h"
13 : #include "LineSegment.h"
14 : #include "PenetrationInfo.h"
15 :
16 : // libMesh
17 : #include "libmesh/fe_base.h"
18 : #include "libmesh/boundary_info.h"
19 : #include "libmesh/elem.h"
20 : #include "libmesh/plane.h"
21 : #include "libmesh/fe_interface.h"
22 : #include "libmesh/dense_vector.h"
23 : #include "libmesh/fe_base.h"
24 : #include "libmesh/vector_value.h"
25 :
26 : using namespace libMesh;
27 :
28 : namespace Moose
29 : {
30 :
31 : /**
32 : * Finds the closest point (called the contact point) on the primary_elem on side "side" to the
33 : * secondary_point.
34 : *
35 : * @param p_info The penetration info object, contains primary_elem, side, various other information
36 : * @param fe_elem FE object for the element
37 : * @param fe_side FE object for the side
38 : * @param fe_side_type Unused; was used for a now-deprecated
39 : * inverse_map overload.
40 : * @param start_with_centroid if true, start inverse mapping procedure from element centroid
41 : * @param tangential_tolerance 'tangential' tolerance for determining whether a contact point on a
42 : * side
43 : * @param secondary_point The physical space coordinates of the secondary node
44 : * @param contact_point_on_side whether or not the contact_point actually lies on _that_ side of the
45 : * element.
46 : * @param search_succeeded whether or not the search for the contact point succeeded. If not it
47 : * should likely be discarded
48 : */
49 : void
50 1504438 : findContactPoint(PenetrationInfo & p_info,
51 : FEBase * fe_elem,
52 : FEBase * fe_side,
53 : FEType & /* fe_side_type */,
54 : const libMesh::Point & secondary_point,
55 : bool start_with_centroid,
56 : const Real tangential_tolerance,
57 : bool & contact_point_on_side,
58 : bool & search_succeeded)
59 : {
60 : // Default to true and we'll switch on failures
61 1504438 : search_succeeded = true;
62 :
63 1504438 : const Elem * primary_elem = p_info._elem;
64 :
65 1504438 : unsigned int dim = primary_elem->dim();
66 :
67 1504438 : const Elem * side = p_info._side;
68 :
69 1504438 : const std::vector<libMesh::Point> & phys_point = fe_side->get_xyz();
70 :
71 1504438 : const std::vector<RealGradient> & dxyz_dxi = fe_side->get_dxyzdxi();
72 1504438 : const std::vector<RealGradient> & d2xyz_dxi2 = fe_side->get_d2xyzdxi2();
73 1504438 : const std::vector<RealGradient> & d2xyz_dxieta = fe_side->get_d2xyzdxideta();
74 :
75 1504438 : const std::vector<RealGradient> & dxyz_deta = fe_side->get_dxyzdeta();
76 1504438 : const std::vector<RealGradient> & d2xyz_deta2 = fe_side->get_d2xyzdeta2();
77 1504438 : const std::vector<RealGradient> & d2xyz_detaxi = fe_side->get_d2xyzdxideta();
78 :
79 1504438 : if (dim == 1)
80 : {
81 21 : const Node * nearest_node = side->node_ptr(0);
82 21 : p_info._closest_point = *nearest_node;
83 : p_info._closest_point_ref =
84 21 : primary_elem->master_point(primary_elem->get_node_index(nearest_node));
85 21 : std::vector<libMesh::Point> elem_points = {p_info._closest_point_ref};
86 :
87 21 : const std::vector<RealGradient> & elem_dxyz_dxi = fe_elem->get_dxyzdxi();
88 :
89 21 : fe_elem->reinit(primary_elem, &elem_points);
90 21 : fe_side->reinit(side, &elem_points);
91 :
92 21 : p_info._normal = elem_dxyz_dxi[0];
93 21 : if (nearest_node->id() == primary_elem->node_id(0))
94 21 : p_info._normal *= -1.0;
95 21 : p_info._normal /= p_info._normal.norm();
96 :
97 21 : libMesh::Point from_secondary_to_closest = p_info._closest_point - secondary_point;
98 21 : p_info._distance = from_secondary_to_closest * p_info._normal;
99 21 : libMesh::Point tangential = from_secondary_to_closest - p_info._distance * p_info._normal;
100 21 : p_info._tangential_distance = tangential.norm();
101 21 : p_info._dxyzdxi = dxyz_dxi;
102 21 : p_info._dxyzdeta = dxyz_deta;
103 21 : p_info._d2xyzdxideta = d2xyz_dxieta;
104 21 : p_info._side_phi = fe_side->get_phi();
105 21 : p_info._side_grad_phi = fe_side->get_dphi();
106 21 : contact_point_on_side = true;
107 21 : return;
108 21 : }
109 :
110 1504417 : libMesh::Point ref_point;
111 :
112 1504417 : if (start_with_centroid)
113 852767 : ref_point = FEMap::inverse_map(dim - 1, side, side->vertex_average(), TOLERANCE, false);
114 : else
115 651650 : ref_point = p_info._closest_point_ref;
116 :
117 1504417 : std::vector<libMesh::Point> points = {ref_point};
118 1504417 : fe_side->reinit(side, &points);
119 1504417 : RealGradient d = secondary_point - phys_point[0];
120 :
121 1504417 : Real update_size = std::numeric_limits<Real>::max();
122 :
123 : // Least squares
124 4051286 : for (unsigned int it = 0; it < 3 && update_size > TOLERANCE * 1e3; ++it)
125 : {
126 2546869 : DenseMatrix<Real> jac(dim - 1, dim - 1);
127 2546869 : jac(0, 0) = -(dxyz_dxi[0] * dxyz_dxi[0]);
128 :
129 2546869 : if (dim - 1 == 2)
130 : {
131 2179829 : jac(1, 0) = -(dxyz_dxi[0] * dxyz_deta[0]);
132 2179829 : jac(0, 1) = -(dxyz_deta[0] * dxyz_dxi[0]);
133 2179829 : jac(1, 1) = -(dxyz_deta[0] * dxyz_deta[0]);
134 : }
135 :
136 2546869 : DenseVector<Real> rhs(dim - 1);
137 2546869 : rhs(0) = dxyz_dxi[0] * d;
138 :
139 2546869 : if (dim - 1 == 2)
140 2179829 : rhs(1) = dxyz_deta[0] * d;
141 :
142 2546869 : DenseVector<Real> update(dim - 1);
143 2546869 : jac.lu_solve(rhs, update);
144 :
145 2546869 : ref_point(0) -= update(0);
146 :
147 2546869 : if (dim - 1 == 2)
148 2179829 : ref_point(1) -= update(1);
149 :
150 2546869 : points[0] = ref_point;
151 2546869 : fe_side->reinit(side, &points);
152 2546869 : d = secondary_point - phys_point[0];
153 :
154 2546869 : update_size = update.l2_norm();
155 2546869 : }
156 :
157 1504417 : update_size = std::numeric_limits<Real>::max();
158 :
159 1504417 : unsigned nit = 0;
160 :
161 : // Newton Loop
162 1504417 : const auto max_newton_its = 25;
163 1504417 : const auto tolerance_newton = 1e3 * TOLERANCE * TOLERANCE;
164 3015981 : for (; nit < max_newton_its && update_size > tolerance_newton; nit++)
165 : {
166 1511564 : d = secondary_point - phys_point[0];
167 :
168 1511564 : DenseMatrix<Real> jac(dim - 1, dim - 1);
169 1511564 : jac(0, 0) = (d2xyz_dxi2[0] * d) - (dxyz_dxi[0] * dxyz_dxi[0]);
170 :
171 1511564 : if (dim - 1 == 2)
172 : {
173 1273811 : jac(1, 0) = (d2xyz_dxieta[0] * d) - (dxyz_dxi[0] * dxyz_deta[0]);
174 :
175 1273811 : jac(0, 1) = (d2xyz_detaxi[0] * d) - (dxyz_deta[0] * dxyz_dxi[0]);
176 1273811 : jac(1, 1) = (d2xyz_deta2[0] * d) - (dxyz_deta[0] * dxyz_deta[0]);
177 : }
178 :
179 1511564 : DenseVector<Real> rhs(dim - 1);
180 1511564 : rhs(0) = -dxyz_dxi[0] * d;
181 :
182 1511564 : if (dim - 1 == 2)
183 1273811 : rhs(1) = -dxyz_deta[0] * d;
184 :
185 1511564 : DenseVector<Real> update(dim - 1);
186 1511564 : jac.lu_solve(rhs, update);
187 :
188 : // Improvised line search in case the update is too large and gets out of the element so bad
189 : // that we cannot reinit at the new point
190 1511564 : Real mult = 1;
191 : while (true)
192 : {
193 : try
194 : {
195 1512467 : ref_point(0) += mult * update(0);
196 :
197 1512467 : if (dim - 1 == 2)
198 1274714 : ref_point(1) += mult * update(1);
199 :
200 1512467 : points[0] = ref_point;
201 1512467 : fe_side->reinit(side, &points);
202 1511564 : d = secondary_point - phys_point[0];
203 :
204 : // we don't multiply by 'mult' because it is used for convergence
205 1511564 : update_size = update.l2_norm();
206 1511564 : break;
207 : }
208 903 : catch (libMesh::LogicError & e)
209 : {
210 903 : ref_point(0) -= mult * update(0);
211 903 : if (dim - 1 == 2)
212 903 : ref_point(1) -= mult * update(1);
213 :
214 903 : mult *= 0.5;
215 903 : if (mult < 1e-6)
216 : {
217 : #ifndef NDEBUG
218 : mooseWarning("We could not solve for the contact point.", e.what());
219 : #endif
220 0 : update_size = update.l2_norm();
221 0 : d = (secondary_point - phys_point[0]) * mult;
222 0 : break;
223 : }
224 903 : }
225 903 : }
226 : // We failed the line search, make sure to trigger the error
227 1511564 : if (mult < 1e-6)
228 : {
229 0 : nit = max_newton_its;
230 0 : update_size = 1;
231 0 : break;
232 : }
233 1511564 : }
234 :
235 1504417 : if (nit == max_newton_its && update_size > tolerance_newton)
236 : {
237 14 : search_succeeded = false;
238 : #ifndef NDEBUG
239 : const auto initial_point =
240 : start_with_centroid ? side->vertex_average() : ref_point = p_info._closest_point_ref;
241 : Moose::err << "Warning! Newton solve for contact point failed to converge!\nLast update "
242 : "distance was: "
243 : << update_size << "\nInitial point guess: " << initial_point
244 : << "\nLast considered point: " << phys_point[0]
245 : << "\nThis potential contact pair (face, point) will be discarded." << std::endl;
246 : #endif
247 14 : return;
248 : }
249 :
250 1504403 : p_info._closest_point_ref = ref_point;
251 1504403 : p_info._closest_point = phys_point[0];
252 1504403 : p_info._distance = d.norm();
253 :
254 1504403 : if (dim - 1 == 2)
255 : {
256 1266650 : p_info._normal = dxyz_dxi[0].cross(dxyz_deta[0]);
257 1266650 : if (!MooseUtils::absoluteFuzzyEqual(p_info._normal.norm(), 0))
258 1266447 : p_info._normal /= p_info._normal.norm();
259 : }
260 : else
261 : {
262 237753 : const Node * const * elem_nodes = primary_elem->get_nodes();
263 237753 : const libMesh::Point in_plane_vector1 = *elem_nodes[1] - *elem_nodes[0];
264 237753 : const libMesh::Point in_plane_vector2 = *elem_nodes[2] - *elem_nodes[0];
265 :
266 237753 : libMesh::Point out_of_plane_normal = in_plane_vector1.cross(in_plane_vector2);
267 237753 : out_of_plane_normal /= out_of_plane_normal.norm();
268 :
269 237753 : p_info._normal = dxyz_dxi[0].cross(out_of_plane_normal);
270 237753 : if (std::fabs(p_info._normal.norm()) > 1e-15)
271 237753 : p_info._normal /= p_info._normal.norm();
272 : }
273 :
274 : // If the point has not penetrated the face, make the distance negative
275 1504403 : const Real dot(d * p_info._normal);
276 1504403 : if (dot > 0.0)
277 1082880 : p_info._distance = -p_info._distance;
278 :
279 1504403 : contact_point_on_side = side->on_reference_element(ref_point);
280 :
281 1504403 : p_info._tangential_distance = 0.0;
282 :
283 1504403 : if (!contact_point_on_side)
284 : {
285 895713 : p_info._closest_point_on_face_ref = ref_point;
286 895713 : restrictPointToFace(p_info._closest_point_on_face_ref, side, p_info._off_edge_nodes);
287 :
288 895713 : points[0] = p_info._closest_point_on_face_ref;
289 895713 : fe_side->reinit(side, &points);
290 895713 : libMesh::Point closest_point_on_face(phys_point[0]);
291 :
292 895713 : RealGradient off_face = closest_point_on_face - p_info._closest_point;
293 895713 : Real tangential_distance = off_face.norm();
294 895713 : p_info._tangential_distance = tangential_distance;
295 895713 : if (tangential_distance <= tangential_tolerance)
296 : {
297 96782 : contact_point_on_side = true;
298 : }
299 : }
300 :
301 1504403 : const std::vector<std::vector<Real>> & phi = fe_side->get_phi();
302 1504403 : const std::vector<std::vector<RealGradient>> & grad_phi = fe_side->get_dphi();
303 :
304 1504403 : points[0] = p_info._closest_point_ref;
305 1504403 : fe_side->reinit(side, &points);
306 :
307 1504403 : p_info._side_phi = phi;
308 1504403 : p_info._side_grad_phi = grad_phi;
309 1504403 : p_info._dxyzdxi = dxyz_dxi;
310 1504403 : p_info._dxyzdeta = dxyz_deta;
311 1504403 : p_info._d2xyzdxideta = d2xyz_dxieta;
312 1504417 : }
313 :
314 : void
315 895713 : restrictPointToFace(libMesh::Point & p,
316 : const Elem * side,
317 : std::vector<const Node *> & off_edge_nodes)
318 : {
319 895713 : const ElemType t(side->type());
320 895713 : off_edge_nodes.clear();
321 895713 : Real & xi = p(0);
322 895713 : Real & eta = p(1);
323 :
324 895713 : switch (t)
325 : {
326 73640 : case EDGE2:
327 : case EDGE3:
328 : case EDGE4:
329 : {
330 : // The reference 1D element is [-1,1].
331 73640 : if (xi < -1.0)
332 : {
333 44743 : xi = -1.0;
334 44743 : off_edge_nodes.push_back(side->node_ptr(0));
335 : }
336 28897 : else if (xi > 1.0)
337 : {
338 28897 : xi = 1.0;
339 28897 : off_edge_nodes.push_back(side->node_ptr(1));
340 : }
341 73640 : break;
342 : }
343 :
344 2245 : case TRI3:
345 : case TRI6:
346 : case TRI7:
347 : {
348 : // The reference triangle is isosceles
349 : // and is bound by xi=0, eta=0, and xi+eta=1.
350 :
351 2245 : if (xi <= 0.0 && eta <= 0.0)
352 : {
353 235 : xi = 0.0;
354 235 : eta = 0.0;
355 235 : off_edge_nodes.push_back(side->node_ptr(0));
356 : }
357 2010 : else if (xi > 0.0 && xi < 1.0 && eta < 0.0)
358 : {
359 283 : eta = 0.0;
360 283 : off_edge_nodes.push_back(side->node_ptr(0));
361 283 : off_edge_nodes.push_back(side->node_ptr(1));
362 : }
363 1727 : else if (eta > 0.0 && eta < 1.0 && xi < 0.0)
364 : {
365 254 : xi = 0.0;
366 254 : off_edge_nodes.push_back(side->node_ptr(2));
367 254 : off_edge_nodes.push_back(side->node_ptr(0));
368 : }
369 1473 : else if (xi >= 1.0 && (eta - xi) <= -1.0)
370 : {
371 785 : xi = 1.0;
372 785 : eta = 0.0;
373 785 : off_edge_nodes.push_back(side->node_ptr(1));
374 : }
375 688 : else if (eta >= 1.0 && (eta - xi) >= 1.0)
376 : {
377 328 : xi = 0.0;
378 328 : eta = 1.0;
379 328 : off_edge_nodes.push_back(side->node_ptr(2));
380 : }
381 360 : else if ((xi + eta) > 1.0)
382 : {
383 360 : Real delta = (xi + eta - 1.0) / 2.0;
384 360 : xi -= delta;
385 360 : eta -= delta;
386 360 : off_edge_nodes.push_back(side->node_ptr(1));
387 360 : off_edge_nodes.push_back(side->node_ptr(2));
388 : }
389 2245 : break;
390 : }
391 :
392 819828 : case QUAD4:
393 : case QUAD8:
394 : case QUAD9:
395 : {
396 : // The reference quadrilateral element is [-1,1]^2.
397 819828 : if (xi < -1.0)
398 : {
399 283059 : xi = -1.0;
400 283059 : if (eta < -1.0)
401 : {
402 68099 : eta = -1.0;
403 68099 : off_edge_nodes.push_back(side->node_ptr(0));
404 : }
405 214960 : else if (eta > 1.0)
406 : {
407 48755 : eta = 1.0;
408 48755 : off_edge_nodes.push_back(side->node_ptr(3));
409 : }
410 : else
411 : {
412 166205 : off_edge_nodes.push_back(side->node_ptr(3));
413 166205 : off_edge_nodes.push_back(side->node_ptr(0));
414 : }
415 : }
416 536769 : else if (xi > 1.0)
417 : {
418 384334 : xi = 1.0;
419 384334 : if (eta < -1.0)
420 : {
421 77218 : eta = -1.0;
422 77218 : off_edge_nodes.push_back(side->node_ptr(1));
423 : }
424 307116 : else if (eta > 1.0)
425 : {
426 73392 : eta = 1.0;
427 73392 : off_edge_nodes.push_back(side->node_ptr(2));
428 : }
429 : else
430 : {
431 233724 : off_edge_nodes.push_back(side->node_ptr(1));
432 233724 : off_edge_nodes.push_back(side->node_ptr(2));
433 : }
434 : }
435 : else
436 : {
437 152435 : if (eta < -1.0)
438 : {
439 99122 : eta = -1.0;
440 99122 : off_edge_nodes.push_back(side->node_ptr(0));
441 99122 : off_edge_nodes.push_back(side->node_ptr(1));
442 : }
443 53313 : else if (eta > 1.0)
444 : {
445 53313 : eta = 1.0;
446 53313 : off_edge_nodes.push_back(side->node_ptr(2));
447 53313 : off_edge_nodes.push_back(side->node_ptr(3));
448 : }
449 : }
450 819828 : break;
451 : }
452 :
453 0 : default:
454 : {
455 0 : mooseError("Unsupported face type: ", t);
456 : break;
457 : }
458 : }
459 895713 : }
460 :
461 : } // namespace Moose
|
