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 "PenetrationThread.h"
12 : #include "ParallelUniqueId.h"
13 : #include "FindContactPoint.h"
14 : #include "NearestNodeLocator.h"
15 : #include "SubProblem.h"
16 : #include "MooseVariableFE.h"
17 : #include "MooseMesh.h"
18 : #include "MooseUtils.h"
19 :
20 : #include "libmesh/threads.h"
21 :
22 : #include <algorithm>
23 :
24 : using namespace libMesh;
25 :
26 : // Anonymous namespace for helper functions that ought to be moved
27 : // into libMesh
28 : namespace
29 : {
30 : Point
31 10999 : closest_point_to_edge(const Point & src, const Point & p0, const Point & p1)
32 : {
33 10999 : const Point line01 = p1 - p0;
34 10999 : const Real line0c_xi = ((src - p0) * line01) / line01.norm_sq();
35 : // The projection would be behind p0; p0 is closest
36 10999 : if (line0c_xi <= 0)
37 2217 : return p0;
38 : // The projection would be past p1; p1 is closest
39 8782 : if (line0c_xi >= 1)
40 3279 : return p1;
41 : // The projection is on the segment between p0 to p1.
42 5503 : return p0 + line0c_xi * line01;
43 : }
44 :
45 : Point
46 11947 : closest_point_to_side(const Point & src, const Elem & side)
47 : {
48 11947 : switch (side.type())
49 : {
50 719 : case EDGE2:
51 : case EDGE3:
52 : case EDGE4:
53 : mooseAssert(side.has_affine_map(),
54 : "Penetration of elements with curved sides not implemented");
55 719 : return closest_point_to_edge(src, side.point(0), side.point(1));
56 11228 : case TRI3:
57 : case TRI6:
58 : {
59 : mooseAssert(side.has_affine_map(),
60 : "Penetration of elements with curved sides not implemented");
61 11228 : const Point p0 = side.point(0), p1 = side.point(1), p2 = side.point(2);
62 11228 : const Point l01 = p1 - p0, l02 = p2 - p0;
63 11228 : const Point tri_normal = (l01.cross(l02)).unit();
64 11228 : const Point linecs = ((src - p0) * tri_normal) / tri_normal.norm_sq() * tri_normal;
65 11228 : const Point in_plane = src - linecs;
66 11228 : const Point planar_offset = in_plane - p0;
67 : // If we're outside the triangle past line 01, our closest point
68 : // is on that line.
69 11228 : if (planar_offset.cross(l01) * tri_normal > 0)
70 4126 : return closest_point_to_edge(src, p0, p1);
71 : // If we're outside the triangle past line 02, our closest point
72 : // is on that line.
73 7102 : if (planar_offset.cross(l02) * tri_normal < 0)
74 3473 : return closest_point_to_edge(src, p0, p2);
75 : // If we're outside the triangle past line 12, our closest point
76 : // is on that line.
77 3629 : if ((in_plane - p1).cross(p2 - p1) * tri_normal > 0)
78 2681 : return closest_point_to_edge(src, p1, p2);
79 : // We must be inside the triangle!
80 948 : return in_plane;
81 : }
82 0 : case QUAD4:
83 : case QUAD8:
84 : case QUAD9:
85 : case C0POLYGON:
86 0 : mooseError("Not implemented");
87 0 : default:
88 0 : mooseError("Side type not recognized");
89 : break;
90 : }
91 : }
92 :
93 : } // anonymous namespace
94 :
95 : // Mutex to use when accessing _penetration_info;
96 : Threads::spin_mutex pinfo_mutex;
97 :
98 186064 : PenetrationThread::PenetrationThread(
99 : SubProblem & subproblem,
100 : const MooseMesh & mesh,
101 : BoundaryID primary_boundary,
102 : BoundaryID secondary_boundary,
103 : std::map<dof_id_type, PenetrationInfo *> & penetration_info,
104 : bool check_whether_reasonable,
105 : bool update_location,
106 : Real tangential_tolerance,
107 : bool do_normal_smoothing,
108 : Real normal_smoothing_distance,
109 : PenetrationLocator::NORMAL_SMOOTHING_METHOD normal_smoothing_method,
110 : bool use_point_locator,
111 : std::vector<std::vector<FEBase *>> & fes,
112 : FEType & fe_type,
113 : NearestNodeLocator & nearest_node,
114 186064 : const std::map<dof_id_type, std::vector<dof_id_type>> & node_to_elem_map)
115 186064 : : _subproblem(subproblem),
116 186064 : _mesh(mesh),
117 186064 : _primary_boundary(primary_boundary),
118 186064 : _secondary_boundary(secondary_boundary),
119 186064 : _penetration_info(penetration_info),
120 186064 : _check_whether_reasonable(check_whether_reasonable),
121 186064 : _update_location(update_location),
122 186064 : _tangential_tolerance(tangential_tolerance),
123 186064 : _do_normal_smoothing(do_normal_smoothing),
124 186064 : _normal_smoothing_distance(normal_smoothing_distance),
125 186064 : _normal_smoothing_method(normal_smoothing_method),
126 186064 : _use_point_locator(use_point_locator),
127 186064 : _nodal_normal_x(NULL),
128 186064 : _nodal_normal_y(NULL),
129 186064 : _nodal_normal_z(NULL),
130 186064 : _fes(fes),
131 186064 : _fe_type(fe_type),
132 186064 : _nearest_node(nearest_node),
133 186064 : _node_to_elem_map(node_to_elem_map)
134 : {
135 186064 : }
136 :
137 : // Splitting Constructor
138 16094 : PenetrationThread::PenetrationThread(PenetrationThread & x, Threads::split /*split*/)
139 16094 : : _subproblem(x._subproblem),
140 16094 : _mesh(x._mesh),
141 16094 : _primary_boundary(x._primary_boundary),
142 16094 : _secondary_boundary(x._secondary_boundary),
143 16094 : _penetration_info(x._penetration_info),
144 16094 : _check_whether_reasonable(x._check_whether_reasonable),
145 16094 : _update_location(x._update_location),
146 16094 : _tangential_tolerance(x._tangential_tolerance),
147 16094 : _do_normal_smoothing(x._do_normal_smoothing),
148 16094 : _normal_smoothing_distance(x._normal_smoothing_distance),
149 16094 : _normal_smoothing_method(x._normal_smoothing_method),
150 16094 : _use_point_locator(x._use_point_locator),
151 16094 : _fes(x._fes),
152 16094 : _fe_type(x._fe_type),
153 16094 : _nearest_node(x._nearest_node),
154 16094 : _node_to_elem_map(x._node_to_elem_map)
155 : {
156 16094 : }
157 :
158 : void
159 202158 : PenetrationThread::operator()(const NodeIdRange & range)
160 : {
161 202158 : ParallelUniqueId puid;
162 202158 : _tid = puid.id;
163 :
164 : // Must get the variables every time this is run because _tid can change
165 202158 : if (_do_normal_smoothing &&
166 76220 : _normal_smoothing_method == PenetrationLocator::NSM_NODAL_NORMAL_BASED)
167 : {
168 30776 : _nodal_normal_x = &_subproblem.getStandardVariable(_tid, "nodal_normal_x");
169 30776 : _nodal_normal_y = &_subproblem.getStandardVariable(_tid, "nodal_normal_y");
170 46164 : _nodal_normal_z = &_subproblem.getStandardVariable(_tid, "nodal_normal_z");
171 : }
172 :
173 202158 : const BoundaryInfo & boundary_info = _mesh.getMesh().get_boundary_info();
174 202158 : std::unique_ptr<PointLocatorBase> point_locator;
175 202158 : if (_use_point_locator)
176 44 : point_locator = _mesh.getPointLocator();
177 :
178 1468160 : for (const auto & node_id : range)
179 : {
180 1266002 : const Node & node = _mesh.nodeRef(node_id);
181 :
182 : // We're going to get a reference to the pointer for the pinfo for this node
183 : // This will allow us to manipulate this pointer without having to go through
184 : // the _penetration_info map... meaning this is the only mutex we'll have to do!
185 1266002 : pinfo_mutex.lock();
186 1266002 : PenetrationInfo *& info = _penetration_info[node.id()];
187 1266002 : pinfo_mutex.unlock();
188 :
189 1266002 : std::vector<PenetrationInfo *> p_info;
190 1266002 : bool info_set(false);
191 :
192 : // See if we already have info about this node
193 1266002 : if (info)
194 : {
195 715991 : FEBase * fe_elem = _fes[_tid][info->_elem->dim()];
196 715991 : FEBase * fe_side = _fes[_tid][info->_side->dim()];
197 :
198 715991 : if (!_update_location && (info->_distance >= 0 || info->isCaptured()))
199 : {
200 0 : const Point contact_ref = info->_closest_point_ref;
201 0 : bool contact_point_on_side(false);
202 :
203 : // Secondary position must be the previous contact point
204 : // Use the previous reference coordinates
205 0 : std::vector<Point> points(1);
206 0 : points[0] = contact_ref;
207 0 : const std::vector<Point> & secondary_pos = fe_side->get_xyz();
208 0 : bool search_succeeded = false;
209 :
210 0 : Moose::findContactPoint(*info,
211 : fe_elem,
212 : fe_side,
213 : _fe_type,
214 0 : secondary_pos[0],
215 : false,
216 : _tangential_tolerance,
217 : contact_point_on_side,
218 : search_succeeded);
219 :
220 : // Restore the original reference coordinates
221 0 : info->_closest_point_ref = contact_ref;
222 : // Just calculated as the distance of the contact point off the surface (0). Set to 0 to
223 : // avoid round-off.
224 0 : info->_distance = 0.0;
225 0 : info_set = true;
226 0 : }
227 : else
228 : {
229 715991 : Real old_tangential_distance(info->_tangential_distance);
230 715991 : bool contact_point_on_side(false);
231 715991 : bool search_succeeded = false;
232 :
233 715991 : Moose::findContactPoint(*info,
234 : fe_elem,
235 : fe_side,
236 : _fe_type,
237 : node,
238 : false,
239 : _tangential_tolerance,
240 : contact_point_on_side,
241 : search_succeeded);
242 :
243 715991 : if (contact_point_on_side)
244 : {
245 649958 : if (info->_tangential_distance <= 0.0) // on the face
246 : {
247 556474 : info_set = true;
248 : }
249 93484 : else if (info->_tangential_distance > 0.0 && old_tangential_distance > 0.0)
250 : { // off the face but within tolerance, was that way on the last step too
251 90264 : if (info->_side->dim() == 2 && info->_off_edge_nodes.size() < 2)
252 : { // Closest point on face is on a node rather than an edge. Another
253 : // face might be a better candidate.
254 : }
255 : else
256 : {
257 88548 : info_set = true;
258 : }
259 : }
260 : }
261 : }
262 : }
263 :
264 1266002 : if (!info_set)
265 : {
266 620980 : const Node * closest_node = _nearest_node.nearestNode(node.id());
267 :
268 620980 : std::vector<dof_id_type> located_elem_ids;
269 : const std::vector<dof_id_type> * closest_elems;
270 :
271 620980 : if (_use_point_locator)
272 : {
273 3532 : std::set<const Elem *> candidate_elements;
274 3532 : (*point_locator)(*closest_node, candidate_elements);
275 :
276 3532 : if (candidate_elements.empty())
277 0 : mooseError("No proximate elements found at node ",
278 0 : closest_node->id(),
279 : " at ",
280 0 : static_cast<const Point &>(*closest_node),
281 : " on boundary ",
282 0 : _nearest_node._boundary1,
283 : ". This should never happen.");
284 :
285 29241 : for (const Elem * elem : candidate_elements)
286 : {
287 82914 : for (auto s : elem->side_index_range())
288 69152 : if (boundary_info.has_boundary_id(elem, s, _primary_boundary))
289 : {
290 11947 : located_elem_ids.push_back(elem->id());
291 11947 : break;
292 : }
293 : }
294 :
295 3532 : if (located_elem_ids.empty())
296 0 : mooseError("No proximate elements found at node ",
297 0 : closest_node->id(),
298 : " at ",
299 0 : static_cast<const Point &>(*closest_node),
300 : " on boundary ",
301 0 : _nearest_node._boundary1,
302 : " share that boundary. This may happen if the mesh uses the same boundary id "
303 : "for a nodeset and an unrelated sideset.");
304 :
305 3532 : closest_elems = &located_elem_ids;
306 3532 : }
307 : else
308 : {
309 617448 : auto node_to_elem_pair = _node_to_elem_map.find(closest_node->id());
310 : mooseAssert(node_to_elem_pair != _node_to_elem_map.end(),
311 : "Missing entry in node to elem map");
312 617448 : closest_elems = &(node_to_elem_pair->second);
313 : }
314 :
315 1480534 : for (const auto & elem_id : *closest_elems)
316 : {
317 859554 : const Elem * elem = _mesh.elemPtr(elem_id);
318 :
319 859554 : std::vector<PenetrationInfo *> thisElemInfo;
320 :
321 859554 : std::vector<const Node *> nodesThatMustBeOnSide;
322 : // If we have a disconnected mesh, we might not have *any*
323 : // nodes that must be on a side we check; we'll rely on
324 : // boundary info to find valid sides, then rely on comparing
325 : // closest points from each to find the best.
326 : //
327 : // If we don't have a disconnected mesh, then for maximum
328 : // backwards compatibility we're still using the older ridge
329 : // and peak detection code, which depends on us ruling out
330 : // sides that don't touch closest_node.
331 859554 : if (!_use_point_locator)
332 847607 : nodesThatMustBeOnSide.push_back(closest_node);
333 859554 : createInfoForElem(
334 859554 : thisElemInfo, p_info, &node, elem, nodesThatMustBeOnSide, _check_whether_reasonable);
335 859554 : }
336 :
337 620980 : if (_use_point_locator)
338 : {
339 3532 : Real min_distance_sq = std::numeric_limits<Real>::max();
340 3532 : Point best_point;
341 3532 : unsigned int best_i = invalid_uint;
342 :
343 : // Find closest point in all p_info to the node of interest
344 15479 : for (unsigned int i = 0; i < p_info.size(); ++i)
345 : {
346 11947 : const Point closest_point = closest_point_to_side(node, *p_info[i]->_side);
347 11947 : const Real distance_sq = (closest_point - node).norm_sq();
348 11947 : if (distance_sq < min_distance_sq)
349 : {
350 5543 : min_distance_sq = distance_sq;
351 5543 : best_point = closest_point;
352 5543 : best_i = i;
353 : }
354 : }
355 :
356 3532 : p_info[best_i]->_closest_point = best_point;
357 7064 : p_info[best_i]->_distance =
358 3532 : (p_info[best_i]->_distance >= 0.0 ? 1.0 : -1.0) * std::sqrt(min_distance_sq);
359 3532 : if (_do_normal_smoothing)
360 0 : mooseError("Normal smoothing not implemented with point locator code");
361 3532 : Point normal = (best_point - node).unit();
362 3532 : const Real dot = normal * p_info[best_i]->_normal;
363 3532 : if (dot < 0)
364 3472 : normal *= -1;
365 3532 : p_info[best_i]->_normal = normal;
366 :
367 3532 : switchInfo(info, p_info[best_i]);
368 3532 : info_set = true;
369 : }
370 : else
371 : {
372 617448 : if (p_info.size() == 1)
373 : {
374 438998 : if (p_info[0]->_tangential_distance <= _tangential_tolerance)
375 : {
376 10334 : switchInfo(info, p_info[0]);
377 10334 : info_set = true;
378 : }
379 : }
380 178450 : else if (p_info.size() > 1)
381 : {
382 : // Loop through all pairs of faces, and check for contact on ridge between each face pair
383 173557 : std::vector<RidgeData> ridgeDataVec;
384 386010 : for (unsigned int i = 0; i + 1 < p_info.size(); ++i)
385 484926 : for (unsigned int j = i + 1; j < p_info.size(); ++j)
386 : {
387 272473 : Point closest_coor;
388 272473 : Real tangential_distance(0.0);
389 272473 : const Node * closest_node_on_ridge = NULL;
390 272473 : unsigned int index = 0;
391 272473 : Point closest_coor_ref;
392 272473 : bool found_ridge_contact_point = findRidgeContactPoint(closest_coor,
393 : tangential_distance,
394 : closest_node_on_ridge,
395 : index,
396 : closest_coor_ref,
397 : p_info,
398 : i,
399 : j);
400 272473 : if (found_ridge_contact_point)
401 : {
402 99396 : RidgeData rpd;
403 99396 : rpd._closest_coor = closest_coor;
404 99396 : rpd._tangential_distance = tangential_distance;
405 99396 : rpd._closest_node = closest_node_on_ridge;
406 99396 : rpd._index = index;
407 99396 : rpd._closest_coor_ref = closest_coor_ref;
408 99396 : ridgeDataVec.push_back(rpd);
409 : }
410 : }
411 :
412 173557 : if (ridgeDataVec.size() > 0) // Either find the ridge pair that is the best or find a peak
413 : {
414 : // Group together ridges for which we are off the edge of a common node.
415 : // Those are peaks.
416 79374 : std::vector<RidgeSetData> ridgeSetDataVec;
417 178770 : for (unsigned int i = 0; i < ridgeDataVec.size(); ++i)
418 : {
419 99396 : bool foundSetWithMatchingNode = false;
420 118138 : for (unsigned int j = 0; j < ridgeSetDataVec.size(); ++j)
421 : {
422 35170 : if (ridgeDataVec[i]._closest_node != NULL &&
423 9970 : ridgeDataVec[i]._closest_node == ridgeSetDataVec[j]._closest_node)
424 : {
425 6458 : foundSetWithMatchingNode = true;
426 6458 : ridgeSetDataVec[j]._ridge_data_vec.push_back(ridgeDataVec[i]);
427 6458 : break;
428 : }
429 : }
430 99396 : if (!foundSetWithMatchingNode)
431 : {
432 92938 : RidgeSetData rsd;
433 92938 : rsd._distance = std::numeric_limits<Real>::max();
434 92938 : rsd._ridge_data_vec.push_back(ridgeDataVec[i]);
435 92938 : rsd._closest_node = ridgeDataVec[i]._closest_node;
436 92938 : ridgeSetDataVec.push_back(rsd);
437 92938 : }
438 : }
439 : // Compute distance to each set of ridges
440 172312 : for (unsigned int i = 0; i < ridgeSetDataVec.size(); ++i)
441 : {
442 92938 : if (ridgeSetDataVec[i]._closest_node !=
443 : NULL) // Either a peak or off the edge of single ridge
444 : {
445 50912 : if (ridgeSetDataVec[i]._ridge_data_vec.size() == 1) // off edge of single ridge
446 : {
447 45330 : if (ridgeSetDataVec[i]._ridge_data_vec[0]._tangential_distance <=
448 45330 : _tangential_tolerance) // off within tolerance
449 : {
450 14568 : ridgeSetDataVec[i]._closest_coor =
451 14568 : ridgeSetDataVec[i]._ridge_data_vec[0]._closest_coor;
452 14568 : Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
453 14568 : ridgeSetDataVec[i]._distance = contact_point_vec.norm();
454 : }
455 : }
456 : else // several ridges join at common node to make peak. The common node is the
457 : // contact point
458 : {
459 5582 : ridgeSetDataVec[i]._closest_coor = *ridgeSetDataVec[i]._closest_node;
460 5582 : Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
461 5582 : ridgeSetDataVec[i]._distance = contact_point_vec.norm();
462 : }
463 : }
464 : else // on a single ridge
465 : {
466 42026 : ridgeSetDataVec[i]._closest_coor =
467 42026 : ridgeSetDataVec[i]._ridge_data_vec[0]._closest_coor;
468 42026 : Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
469 42026 : ridgeSetDataVec[i]._distance = contact_point_vec.norm();
470 : }
471 : }
472 : // Find the set of ridges closest to us.
473 79374 : unsigned int closest_ridge_set_index(0);
474 79374 : Real closest_distance(ridgeSetDataVec[0]._distance);
475 79374 : Point closest_point(ridgeSetDataVec[0]._closest_coor);
476 92938 : for (unsigned int i = 1; i < ridgeSetDataVec.size(); ++i)
477 : {
478 13564 : if (ridgeSetDataVec[i]._distance < closest_distance)
479 : {
480 2198 : closest_ridge_set_index = i;
481 2198 : closest_distance = ridgeSetDataVec[i]._distance;
482 2198 : closest_point = ridgeSetDataVec[i]._closest_coor;
483 : }
484 : }
485 :
486 79374 : if (closest_distance <
487 79374 : std::numeric_limits<Real>::max()) // contact point is on the closest ridge set
488 : {
489 : // find the face in the ridge set with the smallest index, assign that one to the
490 : // interaction
491 50896 : unsigned int face_index(std::numeric_limits<unsigned int>::max());
492 107276 : for (unsigned int i = 0;
493 107276 : i < ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size();
494 : ++i)
495 : {
496 56380 : if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[i]._index < face_index)
497 50896 : face_index = ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[i]._index;
498 : }
499 :
500 : mooseAssert(face_index < std::numeric_limits<unsigned int>::max(),
501 : "face_index invalid");
502 :
503 50896 : p_info[face_index]->_closest_point = closest_point;
504 101792 : p_info[face_index]->_distance =
505 50896 : (p_info[face_index]->_distance >= 0.0 ? 1.0 : -1.0) * closest_distance;
506 : // Calculate the normal as the vector from the ridge to the point only if we're not
507 : // doing normal
508 : // smoothing. Normal smoothing will average out the normals on its own.
509 50896 : if (!_do_normal_smoothing)
510 : {
511 13282 : Point normal(closest_point - node);
512 13282 : const Real len(normal.norm());
513 13282 : if (len > 0)
514 : {
515 13219 : normal /= len;
516 : }
517 13282 : const Real dot(normal * p_info[face_index]->_normal);
518 13282 : if (dot < 0)
519 : {
520 10725 : normal *= -1;
521 : }
522 13282 : p_info[face_index]->_normal = normal;
523 : }
524 50896 : p_info[face_index]->_tangential_distance = 0.0;
525 :
526 50896 : Point closest_point_ref;
527 50896 : if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size() ==
528 : 1) // contact with a single ridge rather than a peak
529 : {
530 46288 : p_info[face_index]->_tangential_distance = ridgeSetDataVec[closest_ridge_set_index]
531 46288 : ._ridge_data_vec[0]
532 46288 : ._tangential_distance;
533 46288 : p_info[face_index]->_closest_point_ref =
534 46288 : ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[0]._closest_coor_ref;
535 : }
536 : else
537 : { // peak
538 : const Node * closest_node_on_face;
539 4608 : bool restricted = restrictPointToFace(p_info[face_index]->_closest_point_ref,
540 : closest_node_on_face,
541 4608 : p_info[face_index]->_side);
542 4608 : if (restricted)
543 : {
544 4608 : if (closest_node_on_face !=
545 4608 : ridgeSetDataVec[closest_ridge_set_index]._closest_node)
546 : {
547 0 : mooseError("Closest node when restricting point to face != closest node from "
548 : "RidgeSetData");
549 : }
550 : }
551 : }
552 :
553 50896 : FEBase * fe = _fes[_tid][p_info[face_index]->_side->dim()];
554 50896 : std::vector<Point> points(1);
555 50896 : points[0] = p_info[face_index]->_closest_point_ref;
556 50896 : fe->reinit(p_info[face_index]->_side, &points);
557 50896 : p_info[face_index]->_side_phi = fe->get_phi();
558 50896 : p_info[face_index]->_side_grad_phi = fe->get_dphi();
559 50896 : p_info[face_index]->_dxyzdxi = fe->get_dxyzdxi();
560 50896 : p_info[face_index]->_dxyzdeta = fe->get_dxyzdeta();
561 50896 : p_info[face_index]->_d2xyzdxideta = fe->get_d2xyzdxideta();
562 :
563 50896 : switchInfo(info, p_info[face_index]);
564 50896 : info_set = true;
565 50896 : }
566 : else
567 : { // todo:remove invalid ridge cases so they don't mess up individual face
568 : // competition????
569 : }
570 79374 : }
571 :
572 173557 : if (!info_set) // contact wasn't on a ridge -- compete individual interactions
573 : {
574 122661 : unsigned int best(0), i(1);
575 : do
576 : {
577 135557 : CompeteInteractionResult CIResult = competeInteractions(p_info[best], p_info[i]);
578 135557 : if (CIResult == FIRST_WINS)
579 : {
580 22524 : i++;
581 : }
582 113033 : else if (CIResult == SECOND_WINS)
583 : {
584 11522 : best = i;
585 11522 : i++;
586 : }
587 101511 : else if (CIResult == NEITHER_WINS)
588 : {
589 101511 : best = i + 1;
590 101511 : i += 2;
591 : }
592 135557 : } while (i < p_info.size() && best < p_info.size());
593 122661 : if (best < p_info.size())
594 : {
595 : // Ensure final info is within the tangential tolerance
596 26433 : if (p_info[best]->_tangential_distance <= _tangential_tolerance)
597 : {
598 26228 : switchInfo(info, p_info[best]);
599 26228 : info_set = true;
600 : }
601 : }
602 : }
603 173557 : }
604 : }
605 620980 : }
606 :
607 1266002 : if (!info_set)
608 : {
609 : // If penetration is not detected within the saved patch, it is possible
610 : // that the secondary node has moved outside the saved patch. So, the patch
611 : // for the secondary nodes saved in _recheck_secondary_nodes has to be updated
612 : // and the penetration detection has to be re-run on the updated patch.
613 :
614 529990 : _recheck_secondary_nodes.push_back(node_id);
615 :
616 529990 : delete info;
617 529990 : info = NULL;
618 : }
619 : else
620 : {
621 736012 : smoothNormal(info, p_info, node);
622 736012 : FEBase * fe = _fes[_tid][info->_side->dim()];
623 736012 : computeSlip(*fe, *info);
624 : }
625 :
626 2216549 : for (unsigned int j = 0; j < p_info.size(); ++j)
627 : {
628 950547 : if (p_info[j])
629 : {
630 859557 : delete p_info[j];
631 859557 : p_info[j] = NULL;
632 : }
633 : }
634 1266002 : }
635 202158 : }
636 :
637 : void
638 16094 : PenetrationThread::join(const PenetrationThread & other)
639 : {
640 16094 : _recheck_secondary_nodes.insert(_recheck_secondary_nodes.end(),
641 : other._recheck_secondary_nodes.begin(),
642 : other._recheck_secondary_nodes.end());
643 16094 : }
644 :
645 : void
646 90990 : PenetrationThread::switchInfo(PenetrationInfo *& info, PenetrationInfo *& infoNew)
647 : {
648 : mooseAssert(infoNew != NULL, "infoNew object is null");
649 90990 : if (info)
650 : {
651 63084 : infoNew->_starting_elem = info->_starting_elem;
652 63084 : infoNew->_starting_side_num = info->_starting_side_num;
653 63084 : infoNew->_starting_closest_point_ref = info->_starting_closest_point_ref;
654 63084 : infoNew->_incremental_slip = info->_incremental_slip;
655 63084 : infoNew->_accumulated_slip = info->_accumulated_slip;
656 63084 : infoNew->_accumulated_slip_old = info->_accumulated_slip_old;
657 63084 : infoNew->_frictional_energy = info->_frictional_energy;
658 63084 : infoNew->_frictional_energy_old = info->_frictional_energy_old;
659 63084 : infoNew->_contact_force = info->_contact_force;
660 63084 : infoNew->_contact_force_old = info->_contact_force_old;
661 63084 : infoNew->_lagrange_multiplier = info->_lagrange_multiplier;
662 63084 : infoNew->_lagrange_multiplier_slip = info->_lagrange_multiplier_slip;
663 63084 : infoNew->_locked_this_step = info->_locked_this_step;
664 63084 : infoNew->_stick_locked_this_step = info->_stick_locked_this_step;
665 63084 : infoNew->_mech_status = info->_mech_status;
666 63084 : infoNew->_mech_status_old = info->_mech_status_old;
667 : }
668 : else
669 : {
670 27906 : infoNew->_starting_elem = infoNew->_elem;
671 27906 : infoNew->_starting_side_num = infoNew->_side_num;
672 27906 : infoNew->_starting_closest_point_ref = infoNew->_closest_point_ref;
673 : }
674 90990 : delete info;
675 90990 : info = infoNew;
676 90990 : infoNew = NULL; // Set this to NULL so that we don't delete it (now owned by _penetration_info).
677 90990 : }
678 :
679 : PenetrationThread::CompeteInteractionResult
680 135557 : PenetrationThread::competeInteractions(PenetrationInfo * pi1, PenetrationInfo * pi2)
681 : {
682 :
683 135557 : CompeteInteractionResult result = NEITHER_WINS;
684 :
685 135557 : if (pi1->_tangential_distance > _tangential_tolerance &&
686 112549 : pi2->_tangential_distance > _tangential_tolerance) // out of tol on both faces
687 101511 : result = NEITHER_WINS;
688 :
689 34046 : else if (pi1->_tangential_distance == 0.0 &&
690 21748 : pi2->_tangential_distance > 0.0) // on face 1, off face 2
691 19403 : result = FIRST_WINS;
692 :
693 14643 : else if (pi2->_tangential_distance == 0.0 &&
694 12727 : pi1->_tangential_distance > 0.0) // on face 2, off face 1
695 10382 : result = SECOND_WINS;
696 :
697 4261 : else if (pi1->_tangential_distance <= _tangential_tolerance &&
698 3461 : pi2->_tangential_distance > _tangential_tolerance) // in face 1 tol, out of face 2 tol
699 332 : result = FIRST_WINS;
700 :
701 3929 : else if (pi2->_tangential_distance <= _tangential_tolerance &&
702 3929 : pi1->_tangential_distance > _tangential_tolerance) // in face 2 tol, out of face 1 tol
703 800 : result = SECOND_WINS;
704 :
705 3129 : else if (pi1->_tangential_distance == 0.0 && pi2->_tangential_distance == 0.0) // on both faces
706 2345 : result = competeInteractionsBothOnFace(pi1, pi2);
707 :
708 784 : else if (pi1->_tangential_distance <= _tangential_tolerance &&
709 784 : pi2->_tangential_distance <= _tangential_tolerance) // off but within tol of both faces
710 : {
711 784 : CommonEdgeResult cer = interactionsOffCommonEdge(pi1, pi2);
712 784 : if (cer == COMMON_EDGE || cer == COMMON_NODE) // ridge case.
713 : {
714 : // We already checked for ridges, and it got rejected, so neither face must be valid
715 0 : result = NEITHER_WINS;
716 : // mooseError("Erroneously encountered ridge case");
717 : }
718 784 : else if (cer == EDGE_AND_COMMON_NODE) // off side of face, off corner of another face. Favor
719 : // the off-side face
720 : {
721 500 : if (pi1->_off_edge_nodes.size() == pi2->_off_edge_nodes.size())
722 0 : mooseError("Invalid off_edge_nodes counts");
723 :
724 500 : else if (pi1->_off_edge_nodes.size() == 2)
725 160 : result = FIRST_WINS;
726 :
727 340 : else if (pi2->_off_edge_nodes.size() == 2)
728 340 : result = SECOND_WINS;
729 :
730 : else
731 0 : mooseError("Invalid off_edge_nodes counts");
732 : }
733 : else // The node projects to both faces within tangential tolerance.
734 284 : result = competeInteractionsBothOnFace(pi1, pi2);
735 : }
736 :
737 135557 : return result;
738 : }
739 :
740 : PenetrationThread::CompeteInteractionResult
741 2629 : PenetrationThread::competeInteractionsBothOnFace(PenetrationInfo * pi1, PenetrationInfo * pi2)
742 : {
743 2629 : CompeteInteractionResult result = NEITHER_WINS;
744 :
745 2629 : if (pi1->_distance >= 0.0 && pi2->_distance < 0.0)
746 0 : result = FIRST_WINS; // favor face with positive distance (penetrated) -- first in this case
747 :
748 2629 : else if (pi2->_distance >= 0.0 && pi1->_distance < 0.0)
749 0 : result = SECOND_WINS; // favor face with positive distance (penetrated) -- second in this case
750 :
751 : // TODO: This logic below could cause an abrupt jump from one face to the other with small mesh
752 : // movement. If there is some way to smooth the transition, we should do it.
753 2629 : else if (MooseUtils::relativeFuzzyLessThan(std::abs(pi1->_distance), std::abs(pi2->_distance)))
754 0 : result = FIRST_WINS; // otherwise, favor the closer face -- first in this case
755 :
756 2629 : else if (MooseUtils::relativeFuzzyLessThan(std::abs(pi2->_distance), std::abs(pi1->_distance)))
757 0 : result = SECOND_WINS; // otherwise, favor the closer face -- second in this case
758 :
759 : else // Equal within tolerance. Favor the one with a smaller element id (for repeatibility)
760 : {
761 2629 : if (pi1->_elem->id() < pi2->_elem->id())
762 2629 : result = FIRST_WINS;
763 :
764 : else
765 0 : result = SECOND_WINS;
766 : }
767 :
768 2629 : return result;
769 : }
770 :
771 : PenetrationThread::CommonEdgeResult
772 784 : PenetrationThread::interactionsOffCommonEdge(PenetrationInfo * pi1, PenetrationInfo * pi2)
773 : {
774 784 : CommonEdgeResult common_edge(NO_COMMON);
775 784 : const std::vector<const Node *> & off_edge_nodes1 = pi1->_off_edge_nodes;
776 784 : const std::vector<const Node *> & off_edge_nodes2 = pi2->_off_edge_nodes;
777 784 : const unsigned dim1 = pi1->_side->dim();
778 :
779 784 : if (dim1 == 1)
780 : {
781 : mooseAssert(pi2->_side->dim() == 1, "Incompatible dimensions.");
782 : mooseAssert(off_edge_nodes1.size() < 2 && off_edge_nodes2.size() < 2,
783 : "off_edge_nodes size should be <2 for 2D contact");
784 0 : if (off_edge_nodes1.size() == 1 && off_edge_nodes2.size() == 1 &&
785 0 : off_edge_nodes1[0] == off_edge_nodes2[0])
786 0 : common_edge = COMMON_EDGE;
787 : }
788 : else
789 : {
790 : mooseAssert(dim1 == 2 && pi2->_side->dim() == 2, "Incompatible dimensions.");
791 : mooseAssert(off_edge_nodes1.size() < 3 && off_edge_nodes2.size() < 3,
792 : "off_edge_nodes size should be <3 for 3D contact");
793 784 : if (off_edge_nodes1.size() == 1)
794 : {
795 340 : if (off_edge_nodes2.size() == 1)
796 : {
797 0 : if (off_edge_nodes1[0] == off_edge_nodes2[0])
798 0 : common_edge = COMMON_NODE;
799 : }
800 340 : else if (off_edge_nodes2.size() == 2)
801 : {
802 340 : if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[0] == off_edge_nodes2[1])
803 340 : common_edge = EDGE_AND_COMMON_NODE;
804 : }
805 : }
806 444 : else if (off_edge_nodes1.size() == 2)
807 : {
808 444 : if (off_edge_nodes2.size() == 1)
809 : {
810 160 : if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[1] == off_edge_nodes2[0])
811 160 : common_edge = EDGE_AND_COMMON_NODE;
812 : }
813 284 : else if (off_edge_nodes2.size() == 2)
814 : {
815 284 : if ((off_edge_nodes1[0] == off_edge_nodes2[0] &&
816 568 : off_edge_nodes1[1] == off_edge_nodes2[1]) ||
817 284 : (off_edge_nodes1[1] == off_edge_nodes2[0] && off_edge_nodes1[0] == off_edge_nodes2[1]))
818 0 : common_edge = COMMON_EDGE;
819 : }
820 : }
821 : }
822 784 : return common_edge;
823 : }
824 :
825 : bool
826 272473 : PenetrationThread::findRidgeContactPoint(Point & contact_point,
827 : Real & tangential_distance,
828 : const Node *& closest_node,
829 : unsigned int & index,
830 : Point & contact_point_ref,
831 : std::vector<PenetrationInfo *> & p_info,
832 : const unsigned int index1,
833 : const unsigned int index2)
834 : {
835 272473 : tangential_distance = 0.0;
836 272473 : closest_node = NULL;
837 272473 : PenetrationInfo * pi1 = p_info[index1];
838 272473 : PenetrationInfo * pi2 = p_info[index2];
839 272473 : const unsigned sidedim(pi1->_side->dim());
840 : mooseAssert(sidedim == pi2->_side->dim(), "Incompatible dimensionalities");
841 :
842 : // Nodes on faces for the two interactions
843 272473 : std::vector<const Node *> side1_nodes;
844 272473 : getSideCornerNodes(pi1->_side, side1_nodes);
845 272473 : std::vector<const Node *> side2_nodes;
846 272473 : getSideCornerNodes(pi2->_side, side2_nodes);
847 :
848 272473 : std::sort(side1_nodes.begin(), side1_nodes.end());
849 272473 : std::sort(side2_nodes.begin(), side2_nodes.end());
850 :
851 : // Find nodes shared by the two faces
852 272473 : std::vector<const Node *> common_nodes;
853 272473 : std::set_intersection(side1_nodes.begin(),
854 : side1_nodes.end(),
855 : side2_nodes.begin(),
856 : side2_nodes.end(),
857 : std::inserter(common_nodes, common_nodes.end()));
858 :
859 272473 : if (common_nodes.size() != sidedim)
860 41403 : return false;
861 :
862 : bool found_point1, found_point2;
863 231070 : Point closest_coor_ref1(pi1->_closest_point_ref);
864 : const Node * closest_node1;
865 231070 : found_point1 = restrictPointToSpecifiedEdgeOfFace(
866 : closest_coor_ref1, closest_node1, pi1->_side, common_nodes);
867 :
868 231070 : Point closest_coor_ref2(pi2->_closest_point_ref);
869 : const Node * closest_node2;
870 231070 : found_point2 = restrictPointToSpecifiedEdgeOfFace(
871 : closest_coor_ref2, closest_node2, pi2->_side, common_nodes);
872 :
873 231070 : if (!found_point1 || !found_point2)
874 131674 : return false;
875 :
876 : // if (sidedim == 2)
877 : // {
878 : // TODO:
879 : // We have the parametric coordinates of the closest intersection point for both faces.
880 : // We need to find a point somewhere in the middle of them so there's not an abrupt jump.
881 : // Find that point by taking dot products of vector from contact point to secondary node point
882 : // with face normal vectors to see which face we're closer to.
883 : // }
884 :
885 99396 : FEBase * fe = NULL;
886 99396 : std::vector<Point> points(1);
887 :
888 : // We have to pick one of the two faces to own the contact point. It doesn't really matter
889 : // which one we pick. For repeatibility, pick the face with the lowest index.
890 99396 : if (index1 < index2)
891 : {
892 99396 : fe = _fes[_tid][pi1->_side->dim()];
893 99396 : contact_point_ref = closest_coor_ref1;
894 99396 : points[0] = closest_coor_ref1;
895 99396 : fe->reinit(pi1->_side, &points);
896 99396 : index = index1;
897 : }
898 : else
899 : {
900 0 : fe = _fes[_tid][pi2->_side->dim()];
901 0 : contact_point_ref = closest_coor_ref2;
902 0 : points[0] = closest_coor_ref2;
903 0 : fe->reinit(pi2->_side, &points);
904 0 : index = index2;
905 : }
906 :
907 99396 : contact_point = fe->get_xyz()[0];
908 :
909 99396 : if (sidedim == 2)
910 : {
911 94299 : if (closest_node1) // point is off the ridge between the two elements
912 : {
913 : mooseAssert((closest_node1 == closest_node2 || closest_node2 == NULL),
914 : "If off edge of ridge, closest node must be the same on both elements");
915 57370 : closest_node = closest_node1;
916 :
917 57370 : RealGradient off_face = *closest_node1 - contact_point;
918 57370 : tangential_distance = off_face.norm();
919 : }
920 : }
921 :
922 99396 : return true;
923 272473 : }
924 :
925 : void
926 544946 : PenetrationThread::getSideCornerNodes(const Elem * side, std::vector<const Node *> & corner_nodes)
927 : {
928 544946 : const ElemType t(side->type());
929 544946 : corner_nodes.clear();
930 :
931 544946 : corner_nodes.push_back(side->node_ptr(0));
932 544946 : corner_nodes.push_back(side->node_ptr(1));
933 544946 : switch (t)
934 : {
935 36594 : case EDGE2:
936 : case EDGE3:
937 : case EDGE4:
938 : {
939 36594 : break;
940 : }
941 :
942 15556 : case TRI3:
943 : case TRI6:
944 : case TRI7:
945 : {
946 15556 : corner_nodes.push_back(side->node_ptr(2));
947 15556 : break;
948 : }
949 :
950 492796 : case QUAD4:
951 : case QUAD8:
952 : case QUAD9:
953 : {
954 492796 : corner_nodes.push_back(side->node_ptr(2));
955 492796 : corner_nodes.push_back(side->node_ptr(3));
956 492796 : break;
957 : }
958 :
959 0 : default:
960 : {
961 0 : mooseError("Unsupported face type: ", t);
962 : break;
963 : }
964 : }
965 544946 : }
966 :
967 : bool
968 462140 : PenetrationThread::restrictPointToSpecifiedEdgeOfFace(Point & p,
969 : const Node *& closest_node,
970 : const Elem * side,
971 : const std::vector<const Node *> & edge_nodes)
972 : {
973 462140 : const ElemType t = side->type();
974 462140 : Real & xi = p(0);
975 462140 : Real & eta = p(1);
976 462140 : closest_node = NULL;
977 :
978 462140 : std::vector<unsigned int> local_node_indices;
979 1349826 : for (const auto & edge_node : edge_nodes)
980 : {
981 887686 : unsigned int local_index = side->get_node_index(edge_node);
982 887686 : if (local_index == libMesh::invalid_uint)
983 0 : mooseError("Side does not contain node");
984 887686 : local_node_indices.push_back(local_index);
985 : }
986 : mooseAssert(local_node_indices.size() == side->dim(),
987 : "Number of edge nodes must match side dimensionality");
988 462140 : std::sort(local_node_indices.begin(), local_node_indices.end());
989 :
990 462140 : bool off_of_this_edge = false;
991 :
992 462140 : switch (t)
993 : {
994 36594 : case EDGE2:
995 : case EDGE3:
996 : case EDGE4:
997 : {
998 36594 : if (local_node_indices[0] == 0)
999 : {
1000 18297 : if (xi <= -1.0)
1001 : {
1002 13781 : xi = -1.0;
1003 13781 : off_of_this_edge = true;
1004 13781 : closest_node = side->node_ptr(0);
1005 : }
1006 : }
1007 18297 : else if (local_node_indices[0] == 1)
1008 : {
1009 18297 : if (xi >= 1.0)
1010 : {
1011 9601 : xi = 1.0;
1012 9601 : off_of_this_edge = true;
1013 9601 : closest_node = side->node_ptr(1);
1014 : }
1015 : }
1016 : else
1017 : {
1018 0 : mooseError("Invalid local node indices");
1019 : }
1020 36594 : break;
1021 : }
1022 :
1023 6814 : case TRI3:
1024 : case TRI6:
1025 : case TRI7:
1026 : {
1027 6814 : if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
1028 : {
1029 1998 : if (eta <= 0.0)
1030 : {
1031 1078 : eta = 0.0;
1032 1078 : off_of_this_edge = true;
1033 1078 : if (xi < 0.0)
1034 244 : closest_node = side->node_ptr(0);
1035 834 : else if (xi > 1.0)
1036 340 : closest_node = side->node_ptr(1);
1037 : }
1038 : }
1039 4816 : else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
1040 : {
1041 2381 : if ((xi + eta) > 1.0)
1042 : {
1043 1112 : Real delta = (xi + eta - 1.0) / 2.0;
1044 1112 : xi -= delta;
1045 1112 : eta -= delta;
1046 1112 : off_of_this_edge = true;
1047 1112 : if (xi > 1.0)
1048 281 : closest_node = side->node_ptr(1);
1049 831 : else if (xi < 0.0)
1050 278 : closest_node = side->node_ptr(2);
1051 : }
1052 : }
1053 2435 : else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 2))
1054 : {
1055 2435 : if (xi <= 0.0)
1056 : {
1057 1192 : xi = 0.0;
1058 1192 : off_of_this_edge = true;
1059 1192 : if (eta > 1.0)
1060 377 : closest_node = side->node_ptr(2);
1061 815 : else if (eta < 0.0)
1062 315 : closest_node = side->node_ptr(0);
1063 : }
1064 : }
1065 : else
1066 : {
1067 0 : mooseError("Invalid local node indices");
1068 : }
1069 :
1070 6814 : break;
1071 : }
1072 :
1073 418732 : case QUAD4:
1074 : case QUAD8:
1075 : case QUAD9:
1076 : {
1077 418732 : if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
1078 : {
1079 181510 : if (eta <= -1.0)
1080 : {
1081 111216 : eta = -1.0;
1082 111216 : off_of_this_edge = true;
1083 111216 : if (xi < -1.0)
1084 34897 : closest_node = side->node_ptr(0);
1085 76319 : else if (xi > 1.0)
1086 34965 : closest_node = side->node_ptr(1);
1087 : }
1088 : }
1089 237222 : else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
1090 : {
1091 88092 : if (xi >= 1.0)
1092 : {
1093 72544 : xi = 1.0;
1094 72544 : off_of_this_edge = true;
1095 72544 : if (eta < -1.0)
1096 16548 : closest_node = side->node_ptr(1);
1097 55996 : else if (eta > 1.0)
1098 13602 : closest_node = side->node_ptr(2);
1099 : }
1100 : }
1101 149130 : else if ((local_node_indices[0] == 2) && (local_node_indices[1] == 3))
1102 : {
1103 122306 : if (eta >= 1.0)
1104 : {
1105 90633 : eta = 1.0;
1106 90633 : off_of_this_edge = true;
1107 90633 : if (xi < -1.0)
1108 39377 : closest_node = side->node_ptr(3);
1109 51256 : else if (xi > 1.0)
1110 45953 : closest_node = side->node_ptr(2);
1111 : }
1112 : }
1113 26824 : else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 3))
1114 : {
1115 26824 : if (xi <= -1.0)
1116 : {
1117 18258 : xi = -1.0;
1118 18258 : off_of_this_edge = true;
1119 18258 : if (eta < -1.0)
1120 10870 : closest_node = side->node_ptr(0);
1121 7388 : else if (eta > 1.0)
1122 2918 : closest_node = side->node_ptr(3);
1123 : }
1124 : }
1125 : else
1126 : {
1127 0 : mooseError("Invalid local node indices");
1128 : }
1129 418732 : break;
1130 : }
1131 :
1132 0 : default:
1133 : {
1134 0 : mooseError("Unsupported face type: ", t);
1135 : break;
1136 : }
1137 : }
1138 462140 : return off_of_this_edge;
1139 462140 : }
1140 :
1141 : bool
1142 4608 : PenetrationThread::restrictPointToFace(Point & p, const Node *& closest_node, const Elem * side)
1143 : {
1144 4608 : const ElemType t(side->type());
1145 4608 : Real & xi = p(0);
1146 4608 : Real & eta = p(1);
1147 4608 : closest_node = NULL;
1148 :
1149 4608 : bool off_of_this_face(false);
1150 :
1151 4608 : switch (t)
1152 : {
1153 0 : case EDGE2:
1154 : case EDGE3:
1155 : case EDGE4:
1156 : {
1157 0 : if (xi < -1.0)
1158 : {
1159 0 : xi = -1.0;
1160 0 : off_of_this_face = true;
1161 0 : closest_node = side->node_ptr(0);
1162 : }
1163 0 : else if (xi > 1.0)
1164 : {
1165 0 : xi = 1.0;
1166 0 : off_of_this_face = true;
1167 0 : closest_node = side->node_ptr(1);
1168 : }
1169 0 : break;
1170 : }
1171 :
1172 0 : case TRI3:
1173 : case TRI6:
1174 : case TRI7:
1175 : {
1176 0 : if (eta < 0.0)
1177 : {
1178 0 : eta = 0.0;
1179 0 : off_of_this_face = true;
1180 0 : if (xi < 0.5)
1181 : {
1182 0 : closest_node = side->node_ptr(0);
1183 0 : if (xi < 0.0)
1184 0 : xi = 0.0;
1185 : }
1186 : else
1187 : {
1188 0 : closest_node = side->node_ptr(1);
1189 0 : if (xi > 1.0)
1190 0 : xi = 1.0;
1191 : }
1192 : }
1193 0 : else if ((xi + eta) > 1.0)
1194 : {
1195 0 : Real delta = (xi + eta - 1.0) / 2.0;
1196 0 : xi -= delta;
1197 0 : eta -= delta;
1198 0 : off_of_this_face = true;
1199 0 : if (xi > 0.5)
1200 : {
1201 0 : closest_node = side->node_ptr(1);
1202 0 : if (xi > 1.0)
1203 : {
1204 0 : xi = 1.0;
1205 0 : eta = 0.0;
1206 : }
1207 : }
1208 : else
1209 : {
1210 0 : closest_node = side->node_ptr(2);
1211 0 : if (xi < 0.0)
1212 : {
1213 0 : xi = 0.0;
1214 0 : eta = 1.0;
1215 : }
1216 : }
1217 : }
1218 0 : else if (xi < 0.0)
1219 : {
1220 0 : xi = 0.0;
1221 0 : off_of_this_face = true;
1222 0 : if (eta > 0.5)
1223 : {
1224 0 : closest_node = side->node_ptr(2);
1225 0 : if (eta > 1.0)
1226 0 : eta = 1.0;
1227 : }
1228 : else
1229 : {
1230 0 : closest_node = side->node_ptr(0);
1231 0 : if (eta < 0.0)
1232 0 : eta = 0.0;
1233 : }
1234 : }
1235 0 : break;
1236 : }
1237 :
1238 4608 : case QUAD4:
1239 : case QUAD8:
1240 : case QUAD9:
1241 : {
1242 4608 : if (eta < -1.0)
1243 : {
1244 108 : eta = -1.0;
1245 108 : off_of_this_face = true;
1246 108 : if (xi < 0.0)
1247 : {
1248 14 : closest_node = side->node_ptr(0);
1249 14 : if (xi < -1.0)
1250 14 : xi = -1.0;
1251 : }
1252 : else
1253 : {
1254 94 : closest_node = side->node_ptr(1);
1255 94 : if (xi > 1.0)
1256 94 : xi = 1.0;
1257 : }
1258 : }
1259 4500 : else if (xi > 1.0)
1260 : {
1261 4500 : xi = 1.0;
1262 4500 : off_of_this_face = true;
1263 4500 : if (eta < 0.0)
1264 : {
1265 0 : closest_node = side->node_ptr(1);
1266 0 : if (eta < -1.0)
1267 0 : eta = -1.0;
1268 : }
1269 : else
1270 : {
1271 4500 : closest_node = side->node_ptr(2);
1272 4500 : if (eta > 1.0)
1273 1170 : eta = 1.0;
1274 : }
1275 : }
1276 0 : else if (eta > 1.0)
1277 : {
1278 0 : eta = 1.0;
1279 0 : off_of_this_face = true;
1280 0 : if (xi < 0.0)
1281 : {
1282 0 : closest_node = side->node_ptr(3);
1283 0 : if (xi < -1.0)
1284 0 : xi = -1.0;
1285 : }
1286 : else
1287 : {
1288 0 : closest_node = side->node_ptr(2);
1289 0 : if (xi > 1.0)
1290 0 : xi = 1.0;
1291 : }
1292 : }
1293 0 : else if (xi < -1.0)
1294 : {
1295 0 : xi = -1.0;
1296 0 : off_of_this_face = true;
1297 0 : if (eta < 0.0)
1298 : {
1299 0 : closest_node = side->node_ptr(0);
1300 0 : if (eta < -1.0)
1301 0 : eta = -1.0;
1302 : }
1303 : else
1304 : {
1305 0 : closest_node = side->node_ptr(3);
1306 0 : if (eta > 1.0)
1307 0 : eta = 1.0;
1308 : }
1309 : }
1310 4608 : break;
1311 : }
1312 :
1313 0 : default:
1314 : {
1315 0 : mooseError("Unsupported face type: ", t);
1316 : break;
1317 : }
1318 : }
1319 4608 : return off_of_this_face;
1320 : }
1321 :
1322 : bool
1323 845852 : PenetrationThread::isFaceReasonableCandidate(const Elem * primary_elem,
1324 : const Elem * side,
1325 : FEBase * fe,
1326 : const Point * secondary_point,
1327 : const Real tangential_tolerance)
1328 : {
1329 845852 : unsigned int dim = primary_elem->dim();
1330 :
1331 845852 : const std::vector<Point> & phys_point = fe->get_xyz();
1332 :
1333 845852 : const std::vector<RealGradient> & dxyz_dxi = fe->get_dxyzdxi();
1334 845852 : const std::vector<RealGradient> & dxyz_deta = fe->get_dxyzdeta();
1335 :
1336 845852 : Point ref_point;
1337 :
1338 845852 : std::vector<Point> points(1); // Default constructor gives us a point at 0,0,0
1339 :
1340 845852 : fe->reinit(side, &points);
1341 :
1342 845852 : RealGradient d = *secondary_point - phys_point[0];
1343 :
1344 845852 : const Real twosqrt2 = 2.8284; // way more precision than we actually need here
1345 845852 : Real max_face_length = side->hmax() + twosqrt2 * tangential_tolerance;
1346 :
1347 845852 : RealVectorValue normal;
1348 845852 : if (dim - 1 == 2)
1349 : {
1350 762130 : normal = dxyz_dxi[0].cross(dxyz_deta[0]);
1351 : }
1352 83722 : else if (dim - 1 == 1)
1353 : {
1354 83701 : const Node * const * elem_nodes = primary_elem->get_nodes();
1355 83701 : const Point in_plane_vector1 = *elem_nodes[1] - *elem_nodes[0];
1356 83701 : const Point in_plane_vector2 = *elem_nodes[2] - *elem_nodes[0];
1357 :
1358 83701 : Point out_of_plane_normal = in_plane_vector1.cross(in_plane_vector2);
1359 83701 : out_of_plane_normal /= out_of_plane_normal.norm();
1360 :
1361 83701 : normal = dxyz_dxi[0].cross(out_of_plane_normal);
1362 : }
1363 : else
1364 : {
1365 21 : return true;
1366 : }
1367 845831 : normal /= normal.norm();
1368 :
1369 845831 : const Real dot(d * normal);
1370 :
1371 845831 : const RealGradient normcomp = dot * normal;
1372 845831 : const RealGradient tangcomp = d - normcomp;
1373 :
1374 845831 : const Real tangdist = tangcomp.norm();
1375 :
1376 : // Increase the size of the zone that we consider if the vector from the face
1377 : // to the node has a larger normal component
1378 845831 : const Real faceExpansionFactor = 2.0 * (1.0 + normcomp.norm() / d.norm());
1379 :
1380 845831 : bool isReasonableCandidate = true;
1381 845831 : if (tangdist > faceExpansionFactor * max_face_length)
1382 : {
1383 8881 : isReasonableCandidate = false;
1384 : }
1385 845831 : return isReasonableCandidate;
1386 845852 : }
1387 :
1388 : void
1389 736012 : PenetrationThread::computeSlip(FEBase & fe, PenetrationInfo & info)
1390 : {
1391 : // Slip is current projected position of secondary node minus
1392 : // original projected position of secondary node
1393 736012 : std::vector<Point> points(1);
1394 736012 : points[0] = info._starting_closest_point_ref;
1395 736012 : const auto & side = _elem_side_builder(*info._starting_elem, info._starting_side_num);
1396 736012 : fe.reinit(&side, &points);
1397 736012 : const std::vector<Point> & starting_point = fe.get_xyz();
1398 736012 : info._incremental_slip = info._closest_point - starting_point[0];
1399 736012 : if (info.isCaptured())
1400 : {
1401 57248 : info._frictional_energy =
1402 57248 : info._frictional_energy_old + info._contact_force * info._incremental_slip;
1403 57248 : info._accumulated_slip = info._accumulated_slip_old + info._incremental_slip.norm();
1404 : }
1405 736012 : }
1406 :
1407 : void
1408 736012 : PenetrationThread::smoothNormal(PenetrationInfo * info,
1409 : std::vector<PenetrationInfo *> & p_info,
1410 : const Node & node)
1411 : {
1412 736012 : if (_do_normal_smoothing)
1413 : {
1414 356926 : if (_normal_smoothing_method == PenetrationLocator::NSM_EDGE_BASED)
1415 : {
1416 : // If we are within the smoothing distance of any edges or corners, find the
1417 : // corner nodes for those edges/corners, and weights from distance to edge/corner
1418 285510 : std::vector<Real> edge_face_weights;
1419 285510 : std::vector<PenetrationInfo *> edge_face_info;
1420 :
1421 285510 : getSmoothingFacesAndWeights(info, edge_face_info, edge_face_weights, p_info, node);
1422 :
1423 : mooseAssert(edge_face_info.size() == edge_face_weights.size(),
1424 : "edge_face_info.size() != edge_face_weights.size()");
1425 :
1426 285510 : if (edge_face_info.size() > 0)
1427 : {
1428 : // Smooth the normal using the weighting functions for all participating faces.
1429 127626 : RealVectorValue new_normal;
1430 127626 : Real this_face_weight = 1.0;
1431 :
1432 287244 : for (unsigned int efwi = 0; efwi < edge_face_weights.size(); ++efwi)
1433 : {
1434 159618 : PenetrationInfo * npi = edge_face_info[efwi];
1435 159618 : if (npi)
1436 159618 : new_normal += npi->_normal * edge_face_weights[efwi];
1437 :
1438 159618 : this_face_weight -= edge_face_weights[efwi];
1439 : }
1440 : mooseAssert(this_face_weight >= (0.25 - 1e-8),
1441 : "Sum of weights of other faces shouldn't exceed 0.75");
1442 127626 : new_normal += info->_normal * this_face_weight;
1443 :
1444 127626 : const Real len = new_normal.norm();
1445 127626 : if (len > 0)
1446 127626 : new_normal /= len;
1447 :
1448 127626 : info->_normal = new_normal;
1449 : }
1450 285510 : }
1451 71416 : else if (_normal_smoothing_method == PenetrationLocator::NSM_NODAL_NORMAL_BASED)
1452 : {
1453 : // params.addParam<VariableName>("var_name","description");
1454 : // getParam<VariableName>("var_name")
1455 71416 : info->_normal(0) = _nodal_normal_x->getValue(info->_side, info->_side_phi);
1456 71416 : info->_normal(1) = _nodal_normal_y->getValue(info->_side, info->_side_phi);
1457 71416 : info->_normal(2) = _nodal_normal_z->getValue(info->_side, info->_side_phi);
1458 71416 : const Real len(info->_normal.norm());
1459 71416 : if (len > 0)
1460 70600 : info->_normal /= len;
1461 : }
1462 : }
1463 736012 : }
1464 :
1465 : void
1466 285510 : PenetrationThread::getSmoothingFacesAndWeights(PenetrationInfo * info,
1467 : std::vector<PenetrationInfo *> & edge_face_info,
1468 : std::vector<Real> & edge_face_weights,
1469 : std::vector<PenetrationInfo *> & p_info,
1470 : const Node & secondary_node)
1471 : {
1472 285510 : const Elem * side = info->_side;
1473 285510 : const Point & p = info->_closest_point_ref;
1474 285510 : std::set<dof_id_type> elems_to_exclude;
1475 285510 : elems_to_exclude.insert(info->_elem->id());
1476 :
1477 285510 : std::vector<std::vector<const Node *>> edge_nodes;
1478 :
1479 : // Get the pairs of nodes along every edge that we are close enough to smooth with
1480 285510 : getSmoothingEdgeNodesAndWeights(p, side, edge_nodes, edge_face_weights);
1481 285510 : std::vector<Elem *> edge_neighbor_elems;
1482 285510 : edge_face_info.resize(edge_nodes.size(), NULL);
1483 :
1484 285510 : std::vector<unsigned int> edges_without_neighbors;
1485 :
1486 534093 : for (unsigned int i = 0; i < edge_nodes.size(); ++i)
1487 : {
1488 : // Sort all sets of edge nodes (needed for comparing edges)
1489 248583 : std::sort(edge_nodes[i].begin(), edge_nodes[i].end());
1490 :
1491 248583 : std::vector<PenetrationInfo *> face_info_comm_edge;
1492 248583 : getInfoForFacesWithCommonNodes(
1493 248583 : &secondary_node, elems_to_exclude, edge_nodes[i], face_info_comm_edge, p_info);
1494 :
1495 248583 : if (face_info_comm_edge.size() == 0)
1496 104961 : edges_without_neighbors.push_back(i);
1497 143622 : else if (face_info_comm_edge.size() > 1)
1498 0 : mooseError("Only one neighbor allowed per edge");
1499 : else
1500 143622 : edge_face_info[i] = face_info_comm_edge[0];
1501 248583 : }
1502 :
1503 : // Remove edges without neighbors from the vector, starting from end
1504 285510 : std::vector<unsigned int>::reverse_iterator rit;
1505 390471 : for (rit = edges_without_neighbors.rbegin(); rit != edges_without_neighbors.rend(); ++rit)
1506 : {
1507 104961 : unsigned int index = *rit;
1508 104961 : edge_nodes.erase(edge_nodes.begin() + index);
1509 104961 : edge_face_weights.erase(edge_face_weights.begin() + index);
1510 104961 : edge_face_info.erase(edge_face_info.begin() + index);
1511 : }
1512 :
1513 : // Handle corner case
1514 285510 : if (edge_nodes.size() > 1)
1515 : {
1516 15996 : if (edge_nodes.size() != 2)
1517 0 : mooseError("Invalid number of smoothing edges");
1518 :
1519 : // find common node
1520 15996 : std::vector<const Node *> common_nodes;
1521 63984 : std::set_intersection(edge_nodes[0].begin(),
1522 15996 : edge_nodes[0].end(),
1523 15996 : edge_nodes[1].begin(),
1524 15996 : edge_nodes[1].end(),
1525 : std::inserter(common_nodes, common_nodes.end()));
1526 :
1527 15996 : if (common_nodes.size() != 1)
1528 0 : mooseError("Invalid number of common nodes");
1529 :
1530 47988 : for (const auto & pinfo : edge_face_info)
1531 31992 : elems_to_exclude.insert(pinfo->_elem->id());
1532 :
1533 15996 : std::vector<PenetrationInfo *> face_info_comm_edge;
1534 15996 : getInfoForFacesWithCommonNodes(
1535 : &secondary_node, elems_to_exclude, common_nodes, face_info_comm_edge, p_info);
1536 :
1537 15996 : unsigned int num_corner_neighbors = face_info_comm_edge.size();
1538 :
1539 15996 : if (num_corner_neighbors > 0)
1540 : {
1541 15996 : Real fw0 = edge_face_weights[0];
1542 15996 : Real fw1 = edge_face_weights[1];
1543 :
1544 : // Corner weight is product of edge weights. Spread out over multiple neighbors.
1545 15996 : Real fw_corner = (fw0 * fw1) / static_cast<Real>(num_corner_neighbors);
1546 :
1547 : // Adjust original edge weights
1548 15996 : edge_face_weights[0] *= (1.0 - fw1);
1549 15996 : edge_face_weights[1] *= (1.0 - fw0);
1550 :
1551 31992 : for (unsigned int i = 0; i < num_corner_neighbors; ++i)
1552 : {
1553 15996 : edge_face_weights.push_back(fw_corner);
1554 15996 : edge_face_info.push_back(face_info_comm_edge[i]);
1555 : }
1556 : }
1557 15996 : }
1558 285510 : }
1559 :
1560 : void
1561 285510 : PenetrationThread::getSmoothingEdgeNodesAndWeights(
1562 : const Point & p,
1563 : const Elem * side,
1564 : std::vector<std::vector<const Node *>> & edge_nodes,
1565 : std::vector<Real> & edge_face_weights)
1566 : {
1567 285510 : const ElemType t(side->type());
1568 285510 : const Real & xi = p(0);
1569 285510 : const Real & eta = p(1);
1570 :
1571 285510 : Real smooth_limit = 1.0 - _normal_smoothing_distance;
1572 :
1573 285510 : switch (t)
1574 : {
1575 10870 : case EDGE2:
1576 : case EDGE3:
1577 : case EDGE4:
1578 : {
1579 10870 : if (xi < -smooth_limit)
1580 : {
1581 2178 : std::vector<const Node *> en;
1582 2178 : en.push_back(side->node_ptr(0));
1583 2178 : edge_nodes.push_back(en);
1584 2178 : Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
1585 2178 : if (fw > 0.5)
1586 144 : fw = 0.5;
1587 2178 : edge_face_weights.push_back(fw);
1588 2178 : }
1589 8692 : else if (xi > smooth_limit)
1590 : {
1591 1242 : std::vector<const Node *> en;
1592 1242 : en.push_back(side->node_ptr(1));
1593 1242 : edge_nodes.push_back(en);
1594 1242 : Real fw = 0.5 - (1.0 - xi) / (2.0 * _normal_smoothing_distance);
1595 1242 : if (fw > 0.5)
1596 198 : fw = 0.5;
1597 1242 : edge_face_weights.push_back(fw);
1598 1242 : }
1599 10870 : break;
1600 : }
1601 :
1602 0 : case TRI3:
1603 : case TRI6:
1604 : case TRI7:
1605 : {
1606 0 : if (eta < -smooth_limit)
1607 : {
1608 0 : std::vector<const Node *> en;
1609 0 : en.push_back(side->node_ptr(0));
1610 0 : en.push_back(side->node_ptr(1));
1611 0 : edge_nodes.push_back(en);
1612 0 : Real fw = 0.5 - (1.0 + eta) / (2.0 * _normal_smoothing_distance);
1613 0 : if (fw > 0.5)
1614 0 : fw = 0.5;
1615 0 : edge_face_weights.push_back(fw);
1616 0 : }
1617 0 : if ((xi + eta) > smooth_limit)
1618 : {
1619 0 : std::vector<const Node *> en;
1620 0 : en.push_back(side->node_ptr(1));
1621 0 : en.push_back(side->node_ptr(2));
1622 0 : edge_nodes.push_back(en);
1623 0 : Real fw = 0.5 - (1.0 - xi - eta) / (2.0 * _normal_smoothing_distance);
1624 0 : if (fw > 0.5)
1625 0 : fw = 0.5;
1626 0 : edge_face_weights.push_back(fw);
1627 0 : }
1628 0 : if (xi < -smooth_limit)
1629 : {
1630 0 : std::vector<const Node *> en;
1631 0 : en.push_back(side->node_ptr(2));
1632 0 : en.push_back(side->node_ptr(0));
1633 0 : edge_nodes.push_back(en);
1634 0 : Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
1635 0 : if (fw > 0.5)
1636 0 : fw = 0.5;
1637 0 : edge_face_weights.push_back(fw);
1638 0 : }
1639 0 : break;
1640 : }
1641 :
1642 274640 : case QUAD4:
1643 : case QUAD8:
1644 : case QUAD9:
1645 : {
1646 274640 : if (eta < -smooth_limit)
1647 : {
1648 45392 : std::vector<const Node *> en;
1649 45392 : en.push_back(side->node_ptr(0));
1650 45392 : en.push_back(side->node_ptr(1));
1651 45392 : edge_nodes.push_back(en);
1652 45392 : Real fw = 0.5 - (1.0 + eta) / (2.0 * _normal_smoothing_distance);
1653 45392 : if (fw > 0.5)
1654 15492 : fw = 0.5;
1655 45392 : edge_face_weights.push_back(fw);
1656 45392 : }
1657 274640 : if (xi > smooth_limit)
1658 : {
1659 85787 : std::vector<const Node *> en;
1660 85787 : en.push_back(side->node_ptr(1));
1661 85787 : en.push_back(side->node_ptr(2));
1662 85787 : edge_nodes.push_back(en);
1663 85787 : Real fw = 0.5 - (1.0 - xi) / (2.0 * _normal_smoothing_distance);
1664 85787 : if (fw > 0.5)
1665 19532 : fw = 0.5;
1666 85787 : edge_face_weights.push_back(fw);
1667 85787 : }
1668 274640 : if (eta > smooth_limit)
1669 : {
1670 87184 : std::vector<const Node *> en;
1671 87184 : en.push_back(side->node_ptr(2));
1672 87184 : en.push_back(side->node_ptr(3));
1673 87184 : edge_nodes.push_back(en);
1674 87184 : Real fw = 0.5 - (1.0 - eta) / (2.0 * _normal_smoothing_distance);
1675 87184 : if (fw > 0.5)
1676 25251 : fw = 0.5;
1677 87184 : edge_face_weights.push_back(fw);
1678 87184 : }
1679 274640 : if (xi < -smooth_limit)
1680 : {
1681 26800 : std::vector<const Node *> en;
1682 26800 : en.push_back(side->node_ptr(3));
1683 26800 : en.push_back(side->node_ptr(0));
1684 26800 : edge_nodes.push_back(en);
1685 26800 : Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
1686 26800 : if (fw > 0.5)
1687 14352 : fw = 0.5;
1688 26800 : edge_face_weights.push_back(fw);
1689 26800 : }
1690 274640 : break;
1691 : }
1692 :
1693 0 : default:
1694 : {
1695 0 : mooseError("Unsupported face type: ", t);
1696 : break;
1697 : }
1698 : }
1699 285510 : }
1700 :
1701 : void
1702 264579 : PenetrationThread::getInfoForFacesWithCommonNodes(
1703 : const Node * secondary_node,
1704 : const std::set<dof_id_type> & elems_to_exclude,
1705 : const std::vector<const Node *> edge_nodes,
1706 : std::vector<PenetrationInfo *> & face_info_comm_edge,
1707 : std::vector<PenetrationInfo *> & p_info)
1708 : {
1709 : // elems connected to a node on this edge, find one that has the same corners as this, and is not
1710 : // the current elem
1711 264579 : auto node_to_elem_pair = _node_to_elem_map.find(edge_nodes[0]->id()); // just need one of the
1712 : // nodes
1713 : mooseAssert(node_to_elem_pair != _node_to_elem_map.end(), "Missing entry in node to elem map");
1714 264579 : const std::vector<dof_id_type> & elems_connected_to_node = node_to_elem_pair->second;
1715 :
1716 264579 : std::vector<const Elem *> elems_connected_to_edge;
1717 :
1718 857370 : for (unsigned int ecni = 0; ecni < elems_connected_to_node.size(); ecni++)
1719 : {
1720 592791 : if (elems_to_exclude.find(elems_connected_to_node[ecni]) != elems_to_exclude.end())
1721 296571 : continue;
1722 296220 : const Elem * elem = _mesh.elemPtr(elems_connected_to_node[ecni]);
1723 :
1724 296220 : std::vector<const Node *> nodevec;
1725 4558564 : for (unsigned int ni = 0; ni < elem->n_nodes(); ++ni)
1726 4262344 : if (elem->is_vertex(ni))
1727 2361312 : nodevec.push_back(elem->node_ptr(ni));
1728 :
1729 296220 : std::vector<const Node *> common_nodes;
1730 296220 : std::sort(nodevec.begin(), nodevec.end());
1731 296220 : std::set_intersection(edge_nodes.begin(),
1732 : edge_nodes.end(),
1733 : nodevec.begin(),
1734 : nodevec.end(),
1735 : std::inserter(common_nodes, common_nodes.end()));
1736 :
1737 296220 : if (common_nodes.size() == edge_nodes.size())
1738 159618 : elems_connected_to_edge.push_back(elem);
1739 296220 : }
1740 :
1741 264579 : if (elems_connected_to_edge.size() > 0)
1742 : {
1743 :
1744 : // There are potentially multiple elements that share a common edge
1745 : // 2D:
1746 : // There can only be one element on the same surface
1747 : // 3D:
1748 : // If there are two edge nodes, there can only be one element on the same surface
1749 : // If there is only one edge node (a corner), there could be multiple elements on the same
1750 : // surface
1751 159618 : bool allowMultipleNeighbors = false;
1752 :
1753 159618 : if (elems_connected_to_edge[0]->dim() == 3)
1754 : {
1755 157506 : if (edge_nodes.size() == 1)
1756 : {
1757 15996 : allowMultipleNeighbors = true;
1758 : }
1759 : }
1760 :
1761 175614 : for (unsigned int i = 0; i < elems_connected_to_edge.size(); ++i)
1762 : {
1763 159618 : std::vector<PenetrationInfo *> thisElemInfo;
1764 159618 : getInfoForElem(thisElemInfo, p_info, elems_connected_to_edge[i]);
1765 159618 : if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
1766 : {
1767 38438 : if (thisElemInfo.size() > 1)
1768 0 : mooseError(
1769 : "Found multiple neighbors to current edge/face on surface when only one is allowed");
1770 38438 : face_info_comm_edge.push_back(thisElemInfo[0]);
1771 38438 : break;
1772 : }
1773 :
1774 121180 : createInfoForElem(
1775 121180 : thisElemInfo, p_info, secondary_node, elems_connected_to_edge[i], edge_nodes);
1776 121180 : if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
1777 : {
1778 105184 : if (thisElemInfo.size() > 1)
1779 0 : mooseError(
1780 : "Found multiple neighbors to current edge/face on surface when only one is allowed");
1781 105184 : face_info_comm_edge.push_back(thisElemInfo[0]);
1782 105184 : break;
1783 : }
1784 :
1785 31992 : for (unsigned int j = 0; j < thisElemInfo.size(); ++j)
1786 15996 : face_info_comm_edge.push_back(thisElemInfo[j]);
1787 159618 : }
1788 : }
1789 264579 : }
1790 :
1791 : void
1792 159618 : PenetrationThread::getInfoForElem(std::vector<PenetrationInfo *> & thisElemInfo,
1793 : std::vector<PenetrationInfo *> & p_info,
1794 : const Elem * elem)
1795 : {
1796 322518 : for (const auto & pi : p_info)
1797 : {
1798 162900 : if (!pi)
1799 46026 : continue;
1800 :
1801 116874 : if (pi->_elem == elem)
1802 46026 : thisElemInfo.push_back(pi);
1803 : }
1804 159618 : }
1805 :
1806 : void
1807 980734 : PenetrationThread::createInfoForElem(std::vector<PenetrationInfo *> & thisElemInfo,
1808 : std::vector<PenetrationInfo *> & p_info,
1809 : const Node * secondary_node,
1810 : const Elem * elem,
1811 : const std::vector<const Node *> & nodes_that_must_be_on_side,
1812 : const bool check_whether_reasonable)
1813 : {
1814 980734 : const BoundaryInfo & boundary_info = _mesh.getMesh().get_boundary_info();
1815 :
1816 6559431 : for (auto s : elem->side_index_range())
1817 : {
1818 5595166 : if (!boundary_info.has_boundary_id(elem, s, _primary_boundary))
1819 4628134 : continue;
1820 :
1821 : // Don't create info for this side if one already exists
1822 967032 : bool already_have_info_this_side = false;
1823 967032 : for (const auto & pi : thisElemInfo)
1824 7588 : if (pi->_side_num == s)
1825 : {
1826 7588 : already_have_info_this_side = true;
1827 7588 : break;
1828 : }
1829 :
1830 967032 : if (already_have_info_this_side)
1831 16469 : break;
1832 :
1833 959444 : const Elem * side = elem->build_side_ptr(s).release();
1834 :
1835 : // Only continue with creating info for this side if the side contains
1836 : // all of the nodes in nodes_that_must_be_on_side
1837 959444 : std::vector<const Node *> nodevec;
1838 6119060 : for (unsigned int ni = 0; ni < side->n_nodes(); ++ni)
1839 5159616 : nodevec.push_back(side->node_ptr(ni));
1840 :
1841 959444 : std::sort(nodevec.begin(), nodevec.end());
1842 959444 : std::vector<const Node *> common_nodes;
1843 959444 : std::set_intersection(nodes_that_must_be_on_side.begin(),
1844 : nodes_that_must_be_on_side.end(),
1845 : nodevec.begin(),
1846 : nodevec.end(),
1847 : std::inserter(common_nodes, common_nodes.end()));
1848 959444 : if (common_nodes.size() != nodes_that_must_be_on_side.size())
1849 : {
1850 0 : delete side;
1851 0 : break;
1852 : }
1853 :
1854 959444 : FEBase * fe_elem = _fes[_tid][elem->dim()];
1855 959444 : FEBase * fe_side = _fes[_tid][side->dim()];
1856 :
1857 : // Optionally check to see whether face is reasonable candidate based on an
1858 : // estimate of how closely it is likely to project to the face
1859 959444 : if (check_whether_reasonable)
1860 845852 : if (!isFaceReasonableCandidate(elem, side, fe_side, secondary_node, _tangential_tolerance))
1861 : {
1862 8881 : delete side;
1863 8881 : break;
1864 : }
1865 :
1866 950563 : Point contact_phys;
1867 950563 : Point contact_ref;
1868 950563 : Point contact_on_face_ref;
1869 950563 : Real distance = 0.;
1870 950563 : Real tangential_distance = 0.;
1871 950563 : RealGradient normal;
1872 : bool contact_point_on_side;
1873 950563 : std::vector<const Node *> off_edge_nodes;
1874 950563 : std::vector<std::vector<Real>> side_phi;
1875 950563 : std::vector<std::vector<RealGradient>> side_grad_phi;
1876 950563 : std::vector<RealGradient> dxyzdxi;
1877 950563 : std::vector<RealGradient> dxyzdeta;
1878 950563 : std::vector<RealGradient> d2xyzdxideta;
1879 :
1880 : std::unique_ptr<PenetrationInfo> pen_info =
1881 : std::make_unique<PenetrationInfo>(elem,
1882 : side,
1883 : s,
1884 : normal,
1885 : distance,
1886 : tangential_distance,
1887 : contact_phys,
1888 : contact_ref,
1889 : contact_on_face_ref,
1890 : off_edge_nodes,
1891 : side_phi,
1892 : side_grad_phi,
1893 : dxyzdxi,
1894 : dxyzdeta,
1895 950563 : d2xyzdxideta);
1896 :
1897 950563 : bool search_succeeded = false;
1898 950563 : Moose::findContactPoint(*pen_info,
1899 : fe_elem,
1900 : fe_side,
1901 : _fe_type,
1902 : *secondary_node,
1903 : true,
1904 : _tangential_tolerance,
1905 : contact_point_on_side,
1906 : search_succeeded);
1907 :
1908 : // Do not add contact info from failed searches
1909 950563 : if (search_succeeded)
1910 : {
1911 950547 : thisElemInfo.push_back(pen_info.get());
1912 950547 : p_info.push_back(pen_info.release());
1913 : }
1914 968325 : }
1915 980734 : }
|
