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