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 1652879 : 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 1652879 : search_succeeded = true;
62 :
63 1652879 : const Elem * primary_elem = p_info._elem;
64 :
65 1652879 : unsigned int dim = primary_elem->dim();
66 :
67 1652879 : const Elem * side = p_info._side;
68 :
69 1652879 : const std::vector<libMesh::Point> & phys_point = fe_side->get_xyz();
70 :
71 1652879 : const std::vector<RealGradient> & dxyz_dxi = fe_side->get_dxyzdxi();
72 1652879 : const std::vector<RealGradient> & d2xyz_dxi2 = fe_side->get_d2xyzdxi2();
73 1652879 : const std::vector<RealGradient> & d2xyz_dxieta = fe_side->get_d2xyzdxideta();
74 :
75 1652879 : const std::vector<RealGradient> & dxyz_deta = fe_side->get_dxyzdeta();
76 1652879 : const std::vector<RealGradient> & d2xyz_deta2 = fe_side->get_d2xyzdeta2();
77 1652879 : const std::vector<RealGradient> & d2xyz_detaxi = fe_side->get_d2xyzdxideta();
78 :
79 1652879 : if (dim == 1)
80 : {
81 23 : const Node * nearest_node = side->node_ptr(0);
82 23 : p_info._closest_point = *nearest_node;
83 : p_info._closest_point_ref =
84 23 : primary_elem->master_point(primary_elem->get_node_index(nearest_node));
85 23 : std::vector<libMesh::Point> elem_points = {p_info._closest_point_ref};
86 :
87 23 : const std::vector<RealGradient> & elem_dxyz_dxi = fe_elem->get_dxyzdxi();
88 :
89 23 : fe_elem->reinit(primary_elem, &elem_points);
90 23 : fe_side->reinit(side, &elem_points);
91 :
92 23 : p_info._normal = elem_dxyz_dxi[0];
93 23 : if (nearest_node->id() == primary_elem->node_id(0))
94 23 : p_info._normal *= -1.0;
95 23 : p_info._normal /= p_info._normal.norm();
96 :
97 23 : libMesh::Point from_secondary_to_closest = p_info._closest_point - secondary_point;
98 23 : p_info._distance = from_secondary_to_closest * p_info._normal;
99 23 : libMesh::Point tangential = from_secondary_to_closest - p_info._distance * p_info._normal;
100 23 : p_info._tangential_distance = tangential.norm();
101 23 : p_info._dxyzdxi = dxyz_dxi;
102 23 : p_info._dxyzdeta = dxyz_deta;
103 23 : p_info._d2xyzdxideta = d2xyz_dxieta;
104 23 : p_info._side_phi = fe_side->get_phi();
105 23 : p_info._side_grad_phi = fe_side->get_dphi();
106 23 : contact_point_on_side = true;
107 23 : return;
108 23 : }
109 :
110 1652856 : libMesh::Point ref_point;
111 :
112 1652856 : if (start_with_centroid)
113 937174 : ref_point = FEMap::inverse_map(dim - 1, side, side->vertex_average(), TOLERANCE, false);
114 : else
115 715682 : ref_point = p_info._closest_point_ref;
116 :
117 1652856 : std::vector<libMesh::Point> points = {ref_point};
118 1652856 : fe_side->reinit(side, &points);
119 1652856 : RealGradient d = secondary_point - phys_point[0];
120 :
121 1652856 : Real update_size = std::numeric_limits<Real>::max();
122 :
123 : // Least squares
124 4451874 : for (unsigned int it = 0; it < 3 && update_size > TOLERANCE * 1e3; ++it)
125 : {
126 2799018 : DenseMatrix<Real> jac(dim - 1, dim - 1);
127 2799018 : jac(0, 0) = -(dxyz_dxi[0] * dxyz_dxi[0]);
128 :
129 2799018 : if (dim - 1 == 2)
130 : {
131 2395713 : jac(1, 0) = -(dxyz_dxi[0] * dxyz_deta[0]);
132 2395713 : jac(0, 1) = -(dxyz_deta[0] * dxyz_dxi[0]);
133 2395713 : jac(1, 1) = -(dxyz_deta[0] * dxyz_deta[0]);
134 : }
135 :
136 2799018 : DenseVector<Real> rhs(dim - 1);
137 2799018 : rhs(0) = dxyz_dxi[0] * d;
138 :
139 2799018 : if (dim - 1 == 2)
140 2395713 : rhs(1) = dxyz_deta[0] * d;
141 :
142 2799018 : DenseVector<Real> update(dim - 1);
143 2799018 : jac.lu_solve(rhs, update);
144 :
145 2799018 : ref_point(0) -= update(0);
146 :
147 2799018 : if (dim - 1 == 2)
148 2395713 : ref_point(1) -= update(1);
149 :
150 2799018 : points[0] = ref_point;
151 2799018 : fe_side->reinit(side, &points);
152 2799018 : d = secondary_point - phys_point[0];
153 :
154 2799018 : update_size = update.l2_norm();
155 2799018 : }
156 :
157 1652856 : update_size = std::numeric_limits<Real>::max();
158 :
159 1652856 : unsigned nit = 0;
160 :
161 : // Newton Loop
162 1652856 : const auto max_newton_its = 25;
163 1652856 : const auto tolerance_newton = 1e3 * TOLERANCE * TOLERANCE;
164 3313880 : for (; nit < max_newton_its && update_size > tolerance_newton; nit++)
165 : {
166 1661024 : d = secondary_point - phys_point[0];
167 :
168 1661024 : DenseMatrix<Real> jac(dim - 1, dim - 1);
169 1661024 : jac(0, 0) = (d2xyz_dxi2[0] * d) - (dxyz_dxi[0] * dxyz_dxi[0]);
170 :
171 1661024 : if (dim - 1 == 2)
172 : {
173 1399649 : jac(1, 0) = (d2xyz_dxieta[0] * d) - (dxyz_dxi[0] * dxyz_deta[0]);
174 :
175 1399649 : jac(0, 1) = (d2xyz_detaxi[0] * d) - (dxyz_deta[0] * dxyz_dxi[0]);
176 1399649 : jac(1, 1) = (d2xyz_deta2[0] * d) - (dxyz_deta[0] * dxyz_deta[0]);
177 : }
178 :
179 1661024 : DenseVector<Real> rhs(dim - 1);
180 1661024 : rhs(0) = -dxyz_dxi[0] * d;
181 :
182 1661024 : if (dim - 1 == 2)
183 1399649 : rhs(1) = -dxyz_deta[0] * d;
184 :
185 1661024 : DenseVector<Real> update(dim - 1);
186 1661024 : 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 1661024 : Real mult = 1;
191 : while (true)
192 : {
193 : try
194 : {
195 1662056 : ref_point(0) += mult * update(0);
196 :
197 1662056 : if (dim - 1 == 2)
198 1400681 : ref_point(1) += mult * update(1);
199 :
200 1662056 : points[0] = ref_point;
201 1662056 : fe_side->reinit(side, &points);
202 1661024 : d = secondary_point - phys_point[0];
203 :
204 : // we don't multiply by 'mult' because it is used for convergence
205 1661024 : update_size = update.l2_norm();
206 1661024 : break;
207 : }
208 1032 : catch (libMesh::LogicError & e)
209 : {
210 1032 : ref_point(0) -= mult * update(0);
211 1032 : if (dim - 1 == 2)
212 1032 : ref_point(1) -= mult * update(1);
213 :
214 1032 : mult *= 0.5;
215 1032 : 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 1032 : }
225 1032 : }
226 : // We failed the line search, make sure to trigger the error
227 1661024 : if (mult < 1e-6)
228 : {
229 0 : nit = max_newton_its;
230 0 : update_size = 1;
231 0 : break;
232 : }
233 1661024 : }
234 :
235 1652856 : if (nit == max_newton_its && update_size > tolerance_newton)
236 : {
237 16 : 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 16 : return;
248 : }
249 :
250 1652840 : p_info._closest_point_ref = ref_point;
251 1652840 : p_info._closest_point = phys_point[0];
252 1652840 : p_info._distance = d.norm();
253 :
254 1652840 : if (dim - 1 == 2)
255 : {
256 1391465 : p_info._normal = dxyz_dxi[0].cross(dxyz_deta[0]);
257 1391465 : if (!MooseUtils::absoluteFuzzyEqual(p_info._normal.norm(), 0))
258 1391233 : p_info._normal /= p_info._normal.norm();
259 : }
260 : else
261 : {
262 261375 : const Node * const * elem_nodes = primary_elem->get_nodes();
263 261375 : const libMesh::Point in_plane_vector1 = *elem_nodes[1] - *elem_nodes[0];
264 261375 : const libMesh::Point in_plane_vector2 = *elem_nodes[2] - *elem_nodes[0];
265 :
266 261375 : libMesh::Point out_of_plane_normal = in_plane_vector1.cross(in_plane_vector2);
267 261375 : out_of_plane_normal /= out_of_plane_normal.norm();
268 :
269 261375 : p_info._normal = dxyz_dxi[0].cross(out_of_plane_normal);
270 261375 : if (std::fabs(p_info._normal.norm()) > 1e-15)
271 261375 : p_info._normal /= p_info._normal.norm();
272 : }
273 :
274 : // If the point has not penetrated the face, make the distance negative
275 1652840 : const Real dot(d * p_info._normal);
276 1652840 : if (dot > 0.0)
277 1189175 : p_info._distance = -p_info._distance;
278 :
279 1652840 : contact_point_on_side = side->on_reference_element(ref_point);
280 :
281 1652840 : p_info._tangential_distance = 0.0;
282 :
283 1652840 : if (!contact_point_on_side)
284 : {
285 984211 : p_info._closest_point_on_face_ref = ref_point;
286 984211 : restrictPointToFace(p_info._closest_point_on_face_ref, side, p_info._off_edge_nodes);
287 :
288 984211 : points[0] = p_info._closest_point_on_face_ref;
289 984211 : fe_side->reinit(side, &points);
290 984211 : libMesh::Point closest_point_on_face(phys_point[0]);
291 :
292 984211 : RealGradient off_face = closest_point_on_face - p_info._closest_point;
293 984211 : Real tangential_distance = off_face.norm();
294 984211 : p_info._tangential_distance = tangential_distance;
295 984211 : if (tangential_distance <= tangential_tolerance)
296 : {
297 106400 : contact_point_on_side = true;
298 : }
299 : }
300 :
301 1652840 : const std::vector<std::vector<Real>> & phi = fe_side->get_phi();
302 1652840 : const std::vector<std::vector<RealGradient>> & grad_phi = fe_side->get_dphi();
303 :
304 1652840 : points[0] = p_info._closest_point_ref;
305 1652840 : fe_side->reinit(side, &points);
306 :
307 1652840 : p_info._side_phi = phi;
308 1652840 : p_info._side_grad_phi = grad_phi;
309 1652840 : p_info._dxyzdxi = dxyz_dxi;
310 1652840 : p_info._dxyzdeta = dxyz_deta;
311 1652840 : p_info._d2xyzdxideta = d2xyz_dxieta;
312 1652856 : }
313 :
314 : void
315 984211 : restrictPointToFace(libMesh::Point & p,
316 : const Elem * side,
317 : std::vector<const Node *> & off_edge_nodes)
318 : {
319 984211 : const ElemType t(side->type());
320 984211 : off_edge_nodes.clear();
321 984211 : Real & xi = p(0);
322 984211 : Real & eta = p(1);
323 :
324 984211 : switch (t)
325 : {
326 80825 : case EDGE2:
327 : case EDGE3:
328 : case EDGE4:
329 : {
330 : // The reference 1D element is [-1,1].
331 80825 : if (xi < -1.0)
332 : {
333 49150 : xi = -1.0;
334 49150 : off_edge_nodes.push_back(side->node_ptr(0));
335 : }
336 31675 : else if (xi > 1.0)
337 : {
338 31675 : xi = 1.0;
339 31675 : off_edge_nodes.push_back(side->node_ptr(1));
340 : }
341 80825 : break;
342 : }
343 :
344 2474 : 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 2474 : if (xi <= 0.0 && eta <= 0.0)
352 : {
353 260 : xi = 0.0;
354 260 : eta = 0.0;
355 260 : off_edge_nodes.push_back(side->node_ptr(0));
356 : }
357 2214 : else if (xi > 0.0 && xi < 1.0 && eta < 0.0)
358 : {
359 312 : eta = 0.0;
360 312 : off_edge_nodes.push_back(side->node_ptr(0));
361 312 : off_edge_nodes.push_back(side->node_ptr(1));
362 : }
363 1902 : else if (eta > 0.0 && eta < 1.0 && xi < 0.0)
364 : {
365 280 : xi = 0.0;
366 280 : off_edge_nodes.push_back(side->node_ptr(2));
367 280 : off_edge_nodes.push_back(side->node_ptr(0));
368 : }
369 1622 : else if (xi >= 1.0 && (eta - xi) <= -1.0)
370 : {
371 868 : xi = 1.0;
372 868 : eta = 0.0;
373 868 : off_edge_nodes.push_back(side->node_ptr(1));
374 : }
375 754 : else if (eta >= 1.0 && (eta - xi) >= 1.0)
376 : {
377 359 : xi = 0.0;
378 359 : eta = 1.0;
379 359 : off_edge_nodes.push_back(side->node_ptr(2));
380 : }
381 395 : else if ((xi + eta) > 1.0)
382 : {
383 395 : Real delta = (xi + eta - 1.0) / 2.0;
384 395 : xi -= delta;
385 395 : eta -= delta;
386 395 : off_edge_nodes.push_back(side->node_ptr(1));
387 395 : off_edge_nodes.push_back(side->node_ptr(2));
388 : }
389 2474 : break;
390 : }
391 :
392 900912 : case QUAD4:
393 : case QUAD8:
394 : case QUAD9:
395 : {
396 : // The reference quadrilateral element is [-1,1]^2.
397 900912 : if (xi < -1.0)
398 : {
399 311377 : xi = -1.0;
400 311377 : if (eta < -1.0)
401 : {
402 74657 : eta = -1.0;
403 74657 : off_edge_nodes.push_back(side->node_ptr(0));
404 : }
405 236720 : else if (eta > 1.0)
406 : {
407 53751 : eta = 1.0;
408 53751 : off_edge_nodes.push_back(side->node_ptr(3));
409 : }
410 : else
411 : {
412 182969 : off_edge_nodes.push_back(side->node_ptr(3));
413 182969 : off_edge_nodes.push_back(side->node_ptr(0));
414 : }
415 : }
416 589535 : else if (xi > 1.0)
417 : {
418 422259 : xi = 1.0;
419 422259 : if (eta < -1.0)
420 : {
421 84745 : eta = -1.0;
422 84745 : off_edge_nodes.push_back(side->node_ptr(1));
423 : }
424 337514 : else if (eta > 1.0)
425 : {
426 80643 : eta = 1.0;
427 80643 : off_edge_nodes.push_back(side->node_ptr(2));
428 : }
429 : else
430 : {
431 256871 : off_edge_nodes.push_back(side->node_ptr(1));
432 256871 : off_edge_nodes.push_back(side->node_ptr(2));
433 : }
434 : }
435 : else
436 : {
437 167276 : if (eta < -1.0)
438 : {
439 108629 : eta = -1.0;
440 108629 : off_edge_nodes.push_back(side->node_ptr(0));
441 108629 : off_edge_nodes.push_back(side->node_ptr(1));
442 : }
443 58647 : else if (eta > 1.0)
444 : {
445 58647 : eta = 1.0;
446 58647 : off_edge_nodes.push_back(side->node_ptr(2));
447 58647 : off_edge_nodes.push_back(side->node_ptr(3));
448 : }
449 : }
450 900912 : break;
451 : }
452 :
453 0 : default:
454 : {
455 0 : mooseError("Unsupported face type: ", t);
456 : break;
457 : }
458 : }
459 984211 : }
460 :
461 : } // namespace Moose
|
