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 9972 : closest_point_to_edge(const Point & src, const Point & p0, const Point & p1)
32 : {
33 9972 : const Point line01 = p1 - p0;
34 9972 : const Real line0c_xi = ((src - p0) * line01) / line01.norm_sq();
35 : // The projection would be behind p0; p0 is closest
36 9972 : if (line0c_xi <= 0)
37 2004 : return p0;
38 : // The projection would be past p1; p1 is closest
39 7968 : if (line0c_xi >= 1)
40 2985 : return p1;
41 : // The projection is on the segment between p0 to p1.
42 4983 : return p0 + line0c_xi * line01;
43 : }
44 :
45 : Point
46 10848 : closest_point_to_side(const Point & src, const Elem & side)
47 : {
48 10848 : switch (side.type())
49 : {
50 560 : case EDGE2:
51 : case EDGE3:
52 : case EDGE4:
53 : mooseAssert(side.has_affine_map(),
54 : "Penetration of elements with curved sides not implemented");
55 560 : return closest_point_to_edge(src, side.point(0), side.point(1));
56 10288 : case TRI3:
57 : case TRI6:
58 : {
59 : mooseAssert(side.has_affine_map(),
60 : "Penetration of elements with curved sides not implemented");
61 10288 : const Point p0 = side.point(0), p1 = side.point(1), p2 = side.point(2);
62 10288 : const Point l01 = p1 - p0, l02 = p2 - p0;
63 10288 : const Point tri_normal = (l01.cross(l02)).unit();
64 10288 : const Point linecs = ((src - p0) * tri_normal) / tri_normal.norm_sq() * tri_normal;
65 10288 : const Point in_plane = src - linecs;
66 10288 : 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 10288 : if (planar_offset.cross(l01) * tri_normal > 0)
70 3753 : 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 6535 : if (planar_offset.cross(l02) * tri_normal < 0)
74 3186 : 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 3349 : if ((in_plane - p1).cross(p2 - p1) * tri_normal > 0)
78 2473 : return closest_point_to_edge(src, p1, p2);
79 : // We must be inside the triangle!
80 876 : 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 170008 : 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 170008 : const std::unordered_map<dof_id_type, std::vector<dof_id_type>> & node_to_elem_map)
115 170008 : : _subproblem(subproblem),
116 170008 : _mesh(mesh),
117 170008 : _primary_boundary(primary_boundary),
118 170008 : _secondary_boundary(secondary_boundary),
119 170008 : _penetration_info(penetration_info),
120 170008 : _check_whether_reasonable(check_whether_reasonable),
121 170008 : _update_location(update_location),
122 170008 : _tangential_tolerance(tangential_tolerance),
123 170008 : _do_normal_smoothing(do_normal_smoothing),
124 170008 : _normal_smoothing_distance(normal_smoothing_distance),
125 170008 : _normal_smoothing_method(normal_smoothing_method),
126 170008 : _use_point_locator(use_point_locator),
127 170008 : _nodal_normal_x(NULL),
128 170008 : _nodal_normal_y(NULL),
129 170008 : _nodal_normal_z(NULL),
130 170008 : _fes(fes),
131 170008 : _fe_type(fe_type),
132 170008 : _nearest_node(nearest_node),
133 170008 : _node_to_elem_map(node_to_elem_map)
134 : {
135 170008 : }
136 :
137 : // Splitting Constructor
138 16097 : PenetrationThread::PenetrationThread(PenetrationThread & x, Threads::split /*split*/)
139 16097 : : _subproblem(x._subproblem),
140 16097 : _mesh(x._mesh),
141 16097 : _primary_boundary(x._primary_boundary),
142 16097 : _secondary_boundary(x._secondary_boundary),
143 16097 : _penetration_info(x._penetration_info),
144 16097 : _check_whether_reasonable(x._check_whether_reasonable),
145 16097 : _update_location(x._update_location),
146 16097 : _tangential_tolerance(x._tangential_tolerance),
147 16097 : _do_normal_smoothing(x._do_normal_smoothing),
148 16097 : _normal_smoothing_distance(x._normal_smoothing_distance),
149 16097 : _normal_smoothing_method(x._normal_smoothing_method),
150 16097 : _use_point_locator(x._use_point_locator),
151 16097 : _fes(x._fes),
152 16097 : _fe_type(x._fe_type),
153 16097 : _nearest_node(x._nearest_node),
154 16097 : _node_to_elem_map(x._node_to_elem_map)
155 : {
156 16097 : }
157 :
158 : void
159 186105 : PenetrationThread::operator()(const NodeIdRange & range)
160 : {
161 186105 : ParallelUniqueId puid;
162 186105 : _tid = puid.id;
163 :
164 : // Must get the variables every time this is run because _tid can change
165 186105 : if (_do_normal_smoothing &&
166 70224 : _normal_smoothing_method == PenetrationLocator::NSM_NODAL_NORMAL_BASED)
167 : {
168 28416 : _nodal_normal_x = &_subproblem.getStandardVariable(_tid, "nodal_normal_x");
169 28416 : _nodal_normal_y = &_subproblem.getStandardVariable(_tid, "nodal_normal_y");
170 42624 : _nodal_normal_z = &_subproblem.getStandardVariable(_tid, "nodal_normal_z");
171 : }
172 :
173 186105 : const BoundaryInfo & boundary_info = _mesh.getMesh().get_boundary_info();
174 186105 : std::unique_ptr<PointLocatorBase> point_locator;
175 186105 : if (_use_point_locator)
176 39 : point_locator = _mesh.getPointLocator();
177 :
178 1339508 : for (const auto & node_id : range)
179 : {
180 1153403 : 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 1153403 : pinfo_mutex.lock();
186 1153403 : PenetrationInfo *& info = _penetration_info[node.id()];
187 1153403 : pinfo_mutex.unlock();
188 :
189 1153403 : std::vector<PenetrationInfo *> p_info;
190 1153403 : bool info_set(false);
191 :
192 : // See if we already have info about this node
193 1153403 : if (info)
194 : {
195 651398 : FEBase * fe_elem = _fes[_tid][info->_elem->dim()];
196 651398 : FEBase * fe_side = _fes[_tid][info->_side->dim()];
197 :
198 651398 : 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 651398 : Real old_tangential_distance(info->_tangential_distance);
230 651398 : bool contact_point_on_side(false);
231 651398 : bool search_succeeded = false;
232 :
233 651398 : 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 651398 : if (contact_point_on_side)
244 : {
245 595041 : if (info->_tangential_distance <= 0.0) // on the face
246 : {
247 508037 : info_set = true;
248 : }
249 87004 : 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 83534 : 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 81412 : info_set = true;
258 : }
259 : }
260 : }
261 : }
262 : }
263 :
264 1153403 : if (!info_set)
265 : {
266 563954 : const Node * closest_node = _nearest_node.nearestNode(node.id());
267 :
268 563954 : std::vector<dof_id_type> located_elem_ids;
269 : const std::vector<dof_id_type> * closest_elems;
270 :
271 563954 : if (_use_point_locator)
272 : {
273 3178 : std::set<const Elem *> candidate_elements;
274 3178 : (*point_locator)(*closest_node, candidate_elements);
275 :
276 3178 : 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 26538 : for (const Elem * elem : candidate_elements)
286 : {
287 75088 : for (auto s : elem->side_index_range())
288 62576 : if (boundary_info.has_boundary_id(elem, s, _primary_boundary))
289 : {
290 10848 : located_elem_ids.push_back(elem->id());
291 10848 : break;
292 : }
293 : }
294 :
295 3178 : 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 3178 : closest_elems = &located_elem_ids;
306 3178 : }
307 : else
308 : {
309 560776 : 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 560776 : closest_elems = &(node_to_elem_pair->second);
313 : }
314 :
315 1341776 : for (const auto & elem_id : *closest_elems)
316 : {
317 777822 : const Elem * elem = _mesh.elemPtr(elem_id);
318 :
319 777822 : std::vector<PenetrationInfo *> thisElemInfo;
320 :
321 777822 : 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 777822 : if (!_use_point_locator)
332 766974 : nodesThatMustBeOnSide.push_back(closest_node);
333 777822 : createInfoForElem(
334 777822 : thisElemInfo, p_info, &node, elem, nodesThatMustBeOnSide, _check_whether_reasonable);
335 777822 : }
336 :
337 563954 : if (_use_point_locator)
338 : {
339 3178 : Real min_distance_sq = std::numeric_limits<Real>::max();
340 3178 : Point best_point;
341 3178 : unsigned int best_i = invalid_uint;
342 :
343 : // Find closest point in all p_info to the node of interest
344 14026 : for (unsigned int i = 0; i < p_info.size(); ++i)
345 : {
346 10848 : const Point closest_point = closest_point_to_side(node, *p_info[i]->_side);
347 10848 : const Real distance_sq = (closest_point - node).norm_sq();
348 10848 : if (distance_sq < min_distance_sq)
349 : {
350 4998 : min_distance_sq = distance_sq;
351 4998 : best_point = closest_point;
352 4998 : best_i = i;
353 : }
354 : }
355 :
356 3178 : p_info[best_i]->_closest_point = best_point;
357 6356 : p_info[best_i]->_distance =
358 3178 : (p_info[best_i]->_distance >= 0.0 ? 1.0 : -1.0) * std::sqrt(min_distance_sq);
359 3178 : if (_do_normal_smoothing)
360 0 : mooseError("Normal smoothing not implemented with point locator code");
361 3178 : Point normal = (best_point - node).unit();
362 3178 : const Real dot = normal * p_info[best_i]->_normal;
363 3178 : if (dot < 0)
364 3133 : normal *= -1;
365 3178 : p_info[best_i]->_normal = normal;
366 :
367 3178 : switchInfo(info, p_info[best_i]);
368 3178 : info_set = true;
369 : }
370 : else
371 : {
372 560776 : if (p_info.size() == 1)
373 : {
374 400724 : if (p_info[0]->_tangential_distance <= _tangential_tolerance)
375 : {
376 9513 : switchInfo(info, p_info[0]);
377 9513 : info_set = true;
378 : }
379 : }
380 160052 : 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 155597 : std::vector<RidgeData> ridgeDataVec;
384 345859 : for (unsigned int i = 0; i + 1 < p_info.size(); ++i)
385 434043 : for (unsigned int j = i + 1; j < p_info.size(); ++j)
386 : {
387 243781 : Point closest_coor;
388 243781 : Real tangential_distance(0.0);
389 243781 : const Node * closest_node_on_ridge = NULL;
390 243781 : unsigned int index = 0;
391 243781 : Point closest_coor_ref;
392 243781 : 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 243781 : if (found_ridge_contact_point)
401 : {
402 89190 : RidgeData rpd;
403 89190 : rpd._closest_coor = closest_coor;
404 89190 : rpd._tangential_distance = tangential_distance;
405 89190 : rpd._closest_node = closest_node_on_ridge;
406 89190 : rpd._index = index;
407 89190 : rpd._closest_coor_ref = closest_coor_ref;
408 89190 : ridgeDataVec.push_back(rpd);
409 : }
410 : }
411 :
412 155597 : 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 71313 : std::vector<RidgeSetData> ridgeSetDataVec;
417 160503 : for (unsigned int i = 0; i < ridgeDataVec.size(); ++i)
418 : {
419 89190 : bool foundSetWithMatchingNode = false;
420 105836 : for (unsigned int j = 0; j < ridgeSetDataVec.size(); ++j)
421 : {
422 31343 : if (ridgeDataVec[i]._closest_node != NULL &&
423 8894 : ridgeDataVec[i]._closest_node == ridgeSetDataVec[j]._closest_node)
424 : {
425 5803 : foundSetWithMatchingNode = true;
426 5803 : ridgeSetDataVec[j]._ridge_data_vec.push_back(ridgeDataVec[i]);
427 5803 : break;
428 : }
429 : }
430 89190 : if (!foundSetWithMatchingNode)
431 : {
432 83387 : RidgeSetData rsd;
433 83387 : rsd._distance = std::numeric_limits<Real>::max();
434 83387 : rsd._ridge_data_vec.push_back(ridgeDataVec[i]);
435 83387 : rsd._closest_node = ridgeDataVec[i]._closest_node;
436 83387 : ridgeSetDataVec.push_back(rsd);
437 83387 : }
438 : }
439 : // Compute distance to each set of ridges
440 154700 : for (unsigned int i = 0; i < ridgeSetDataVec.size(); ++i)
441 : {
442 83387 : if (ridgeSetDataVec[i]._closest_node !=
443 : NULL) // Either a peak or off the edge of single ridge
444 : {
445 45948 : if (ridgeSetDataVec[i]._ridge_data_vec.size() == 1) // off edge of single ridge
446 : {
447 40903 : if (ridgeSetDataVec[i]._ridge_data_vec[0]._tangential_distance <=
448 40903 : _tangential_tolerance) // off within tolerance
449 : {
450 13262 : ridgeSetDataVec[i]._closest_coor =
451 13262 : ridgeSetDataVec[i]._ridge_data_vec[0]._closest_coor;
452 13262 : Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
453 13262 : ridgeSetDataVec[i]._distance = contact_point_vec.norm();
454 : }
455 : }
456 : else // several ridges join at common node to make a peak. The common node is the
457 : // contact point
458 : {
459 5045 : ridgeSetDataVec[i]._closest_coor = *ridgeSetDataVec[i]._closest_node;
460 5045 : Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
461 5045 : ridgeSetDataVec[i]._distance = contact_point_vec.norm();
462 : }
463 : }
464 : else // on a single ridge
465 : {
466 37439 : ridgeSetDataVec[i]._closest_coor =
467 37439 : ridgeSetDataVec[i]._ridge_data_vec[0]._closest_coor;
468 37439 : Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
469 37439 : ridgeSetDataVec[i]._distance = contact_point_vec.norm();
470 : }
471 : }
472 : // Find the set of ridges closest to us.
473 71313 : unsigned int closest_ridge_set_index(0);
474 71313 : Real closest_distance(ridgeSetDataVec[0]._distance);
475 71313 : Point closest_point(ridgeSetDataVec[0]._closest_coor);
476 83387 : for (unsigned int i = 1; i < ridgeSetDataVec.size(); ++i)
477 : {
478 12074 : if (ridgeSetDataVec[i]._distance < closest_distance)
479 : {
480 1529 : closest_ridge_set_index = i;
481 1529 : closest_distance = ridgeSetDataVec[i]._distance;
482 1529 : closest_point = ridgeSetDataVec[i]._closest_coor;
483 : }
484 : }
485 :
486 71313 : if (closest_distance <
487 71313 : 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 : // TODO: We may need to select the face with the largest projected distance
492 : // rather than the smallest index, similar to what is done when picking the
493 : // face in findRidgeContactPoint() to better condition corner-case checks for
494 : // sign changes. That's less likely to be a problem here, though, and the
495 : // code to do that would be messier.
496 45602 : unsigned int face_index(std::numeric_limits<unsigned int>::max());
497 95991 : for (unsigned int i = 0;
498 95991 : i < ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size();
499 : ++i)
500 : {
501 50389 : if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[i]._index < face_index)
502 45746 : face_index = ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[i]._index;
503 : }
504 :
505 : mooseAssert(face_index < std::numeric_limits<unsigned int>::max(),
506 : "face_index invalid");
507 :
508 45602 : p_info[face_index]->_closest_point = closest_point;
509 91204 : p_info[face_index]->_distance =
510 45602 : (p_info[face_index]->_distance >= 0.0 ? 1.0 : -1.0) * closest_distance;
511 : // Calculate the normal as the vector from the ridge to the point only if we're not
512 : // doing normal
513 : // smoothing. Normal smoothing will average out the normals on its own.
514 45602 : if (!_do_normal_smoothing)
515 : {
516 12135 : Point normal(closest_point - node);
517 12135 : const Real len(normal.norm());
518 12135 : if (len > 0)
519 : {
520 11958 : normal /= len;
521 : }
522 12135 : const Real dot(normal * p_info[face_index]->_normal);
523 12135 : if (dot < 0)
524 9782 : normal *= -1;
525 12135 : p_info[face_index]->_normal = normal;
526 : }
527 45602 : p_info[face_index]->_tangential_distance = 0.0;
528 :
529 45602 : Point closest_point_ref;
530 45602 : if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size() ==
531 : 1) // contact with a single ridge rather than a peak
532 : {
533 41573 : p_info[face_index]->_tangential_distance = ridgeSetDataVec[closest_ridge_set_index]
534 41573 : ._ridge_data_vec[0]
535 41573 : ._tangential_distance;
536 41573 : p_info[face_index]->_closest_point_ref =
537 41573 : ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[0]._closest_coor_ref;
538 : }
539 : else
540 : { // peak
541 : const Node * closest_node_on_face;
542 4029 : bool restricted = restrictPointToFace(p_info[face_index]->_closest_point_ref,
543 : closest_node_on_face,
544 4029 : p_info[face_index]->_side);
545 4029 : if (restricted)
546 : {
547 4029 : if (closest_node_on_face !=
548 4029 : ridgeSetDataVec[closest_ridge_set_index]._closest_node)
549 : {
550 0 : mooseError("Closest node when restricting point to face != closest node from "
551 : "RidgeSetData");
552 : }
553 : }
554 : }
555 :
556 45602 : FEBase * fe = _fes[_tid][p_info[face_index]->_side->dim()];
557 45602 : std::vector<Point> points(1);
558 45602 : points[0] = p_info[face_index]->_closest_point_ref;
559 45602 : fe->reinit(p_info[face_index]->_side, &points);
560 45602 : p_info[face_index]->_side_phi = fe->get_phi();
561 45602 : p_info[face_index]->_side_grad_phi = fe->get_dphi();
562 45602 : p_info[face_index]->_dxyzdxi = fe->get_dxyzdxi();
563 45602 : p_info[face_index]->_dxyzdeta = fe->get_dxyzdeta();
564 45602 : p_info[face_index]->_d2xyzdxideta = fe->get_d2xyzdxideta();
565 :
566 45602 : switchInfo(info, p_info[face_index]);
567 45602 : info_set = true;
568 45602 : }
569 : else
570 : { // todo:remove invalid ridge cases so they don't mess up individual face
571 : // competition????
572 : }
573 71313 : }
574 :
575 155597 : if (!info_set) // contact wasn't on a ridge -- compete individual interactions
576 : {
577 109995 : unsigned int best(0), i(1);
578 : do
579 : {
580 121532 : CompeteInteractionResult CIResult = competeInteractions(p_info[best], p_info[i]);
581 121532 : if (CIResult == FIRST_WINS)
582 : {
583 20190 : i++;
584 : }
585 101342 : else if (CIResult == SECOND_WINS)
586 : {
587 8786 : best = i;
588 8786 : i++;
589 : }
590 92556 : else if (CIResult == NEITHER_WINS)
591 : {
592 92556 : best = i + 1;
593 92556 : i += 2;
594 : }
595 121532 : } while (i < p_info.size() && best < p_info.size());
596 109995 : if (best < p_info.size())
597 : {
598 : // Ensure final info is within the tangential tolerance
599 22106 : if (p_info[best]->_tangential_distance <= _tangential_tolerance)
600 : {
601 21920 : switchInfo(info, p_info[best]);
602 21920 : info_set = true;
603 : }
604 : }
605 : }
606 155597 : }
607 : }
608 563954 : }
609 :
610 1153403 : if (!info_set)
611 : {
612 : // If penetration is not detected within the saved patch, it is possible
613 : // that the secondary node has moved outside the saved patch. So, the patch
614 : // for the secondary nodes saved in _recheck_secondary_nodes has to be updated
615 : // and the penetration detection has to be re-run on the updated patch.
616 :
617 483741 : _recheck_secondary_nodes.push_back(node_id);
618 :
619 483741 : delete info;
620 483741 : info = NULL;
621 : }
622 : else
623 : {
624 669662 : smoothNormal(info, p_info, node);
625 669662 : FEBase * fe = _fes[_tid][info->_side->dim()];
626 669662 : computeSlip(*fe, *info);
627 : }
628 :
629 2016396 : for (unsigned int j = 0; j < p_info.size(); ++j)
630 : {
631 862993 : if (p_info[j])
632 : {
633 782780 : delete p_info[j];
634 782780 : p_info[j] = NULL;
635 : }
636 : }
637 1153403 : }
638 186105 : }
639 :
640 : void
641 16097 : PenetrationThread::join(const PenetrationThread & other)
642 : {
643 16097 : _recheck_secondary_nodes.insert(_recheck_secondary_nodes.end(),
644 : other._recheck_secondary_nodes.begin(),
645 : other._recheck_secondary_nodes.end());
646 16097 : }
647 :
648 : void
649 80213 : PenetrationThread::switchInfo(PenetrationInfo *& info, PenetrationInfo *& infoNew)
650 : {
651 : mooseAssert(infoNew != NULL, "infoNew object is null");
652 80213 : if (info)
653 : {
654 54742 : infoNew->_starting_elem = info->_starting_elem;
655 54742 : infoNew->_starting_side_num = info->_starting_side_num;
656 54742 : infoNew->_starting_closest_point_ref = info->_starting_closest_point_ref;
657 54742 : infoNew->_incremental_slip = info->_incremental_slip;
658 54742 : infoNew->_accumulated_slip = info->_accumulated_slip;
659 54742 : infoNew->_accumulated_slip_old = info->_accumulated_slip_old;
660 54742 : infoNew->_frictional_energy = info->_frictional_energy;
661 54742 : infoNew->_frictional_energy_old = info->_frictional_energy_old;
662 54742 : infoNew->_contact_force = info->_contact_force;
663 54742 : infoNew->_contact_force_old = info->_contact_force_old;
664 54742 : infoNew->_lagrange_multiplier = info->_lagrange_multiplier;
665 54742 : infoNew->_lagrange_multiplier_slip = info->_lagrange_multiplier_slip;
666 54742 : infoNew->_locked_this_step = info->_locked_this_step;
667 54742 : infoNew->_stick_locked_this_step = info->_stick_locked_this_step;
668 54742 : infoNew->_mech_status = info->_mech_status;
669 54742 : infoNew->_mech_status_old = info->_mech_status_old;
670 : }
671 : else
672 : {
673 25471 : infoNew->_starting_elem = infoNew->_elem;
674 25471 : infoNew->_starting_side_num = infoNew->_side_num;
675 25471 : infoNew->_starting_closest_point_ref = infoNew->_closest_point_ref;
676 : }
677 80213 : delete info;
678 80213 : info = infoNew;
679 80213 : infoNew = NULL; // Set this to NULL so that we don't delete it (now owned by _penetration_info).
680 80213 : }
681 :
682 : PenetrationThread::CompeteInteractionResult
683 121532 : PenetrationThread::competeInteractions(PenetrationInfo * pi1, PenetrationInfo * pi2)
684 : {
685 :
686 121532 : CompeteInteractionResult result = NEITHER_WINS;
687 :
688 121532 : if (pi1->_tangential_distance > _tangential_tolerance &&
689 100846 : pi2->_tangential_distance > _tangential_tolerance) // out of tol on both faces
690 92556 : result = NEITHER_WINS;
691 :
692 28976 : else if (pi1->_tangential_distance == 0.0 &&
693 19536 : pi2->_tangential_distance > 0.0) // on face 1, off face 2
694 17411 : result = FIRST_WINS;
695 :
696 11565 : else if (pi2->_tangential_distance == 0.0 &&
697 10313 : pi1->_tangential_distance > 0.0) // on face 2, off face 1
698 8188 : result = SECOND_WINS;
699 :
700 3377 : else if (pi1->_tangential_distance <= _tangential_tolerance &&
701 3143 : pi2->_tangential_distance > _tangential_tolerance) // in face 1 tol, out of face 2 tol
702 304 : result = FIRST_WINS;
703 :
704 3073 : else if (pi2->_tangential_distance <= _tangential_tolerance &&
705 3073 : pi1->_tangential_distance > _tangential_tolerance) // in face 2 tol, out of face 1 tol
706 234 : result = SECOND_WINS;
707 :
708 2839 : else if (pi1->_tangential_distance == 0.0 && pi2->_tangential_distance == 0.0) // on both faces
709 2125 : result = competeInteractionsBothOnFace(pi1, pi2);
710 :
711 714 : else if (pi1->_tangential_distance <= _tangential_tolerance &&
712 714 : pi2->_tangential_distance <= _tangential_tolerance) // off but within tol of both faces
713 : {
714 714 : CommonEdgeResult cer = interactionsOffCommonEdge(pi1, pi2);
715 714 : if (cer == COMMON_EDGE || cer == COMMON_NODE) // ridge case.
716 : {
717 : // We already checked for ridges, and it got rejected, so neither face must be valid
718 0 : result = NEITHER_WINS;
719 : // mooseError("Erroneously encountered ridge case");
720 : }
721 714 : else if (cer == EDGE_AND_COMMON_NODE) // off side of face, off corner of another face. Favor
722 : // the off-side face
723 : {
724 510 : if (pi1->_off_edge_nodes.size() == pi2->_off_edge_nodes.size())
725 0 : mooseError("Invalid off_edge_nodes counts");
726 :
727 510 : else if (pi1->_off_edge_nodes.size() == 2)
728 146 : result = FIRST_WINS;
729 :
730 364 : else if (pi2->_off_edge_nodes.size() == 2)
731 364 : result = SECOND_WINS;
732 :
733 : else
734 0 : mooseError("Invalid off_edge_nodes counts");
735 : }
736 : else // The node projects to both faces within tangential tolerance.
737 204 : result = competeInteractionsBothOnFace(pi1, pi2);
738 : }
739 :
740 121532 : return result;
741 : }
742 :
743 : PenetrationThread::CompeteInteractionResult
744 2329 : PenetrationThread::competeInteractionsBothOnFace(PenetrationInfo * pi1, PenetrationInfo * pi2)
745 : {
746 2329 : CompeteInteractionResult result = NEITHER_WINS;
747 :
748 2329 : if (pi1->_distance >= 0.0 && pi2->_distance < 0.0)
749 0 : result = FIRST_WINS; // favor face with positive distance (penetrated) -- first in this case
750 :
751 2329 : else if (pi2->_distance >= 0.0 && pi1->_distance < 0.0)
752 0 : result = SECOND_WINS; // favor face with positive distance (penetrated) -- second in this case
753 :
754 : // TODO: This logic below could cause an abrupt jump from one face to the other with small mesh
755 : // movement. If there is some way to smooth the transition, we should do it.
756 2329 : else if (MooseUtils::relativeFuzzyLessThan(std::abs(pi1->_distance), std::abs(pi2->_distance)))
757 0 : result = FIRST_WINS; // otherwise, favor the closer face -- first in this case
758 :
759 2329 : else if (MooseUtils::relativeFuzzyLessThan(std::abs(pi2->_distance), std::abs(pi1->_distance)))
760 0 : result = SECOND_WINS; // otherwise, favor the closer face -- second in this case
761 :
762 : else // Equal within tolerance. Favor the one with a smaller element id (for repeatibility)
763 : {
764 2329 : if (pi1->_elem->id() < pi2->_elem->id())
765 2329 : result = FIRST_WINS;
766 :
767 : else
768 0 : result = SECOND_WINS;
769 : }
770 :
771 2329 : return result;
772 : }
773 :
774 : PenetrationThread::CommonEdgeResult
775 714 : PenetrationThread::interactionsOffCommonEdge(PenetrationInfo * pi1, PenetrationInfo * pi2)
776 : {
777 714 : CommonEdgeResult common_edge(NO_COMMON);
778 714 : const std::vector<const Node *> & off_edge_nodes1 = pi1->_off_edge_nodes;
779 714 : const std::vector<const Node *> & off_edge_nodes2 = pi2->_off_edge_nodes;
780 714 : const unsigned dim1 = pi1->_side->dim();
781 :
782 714 : if (dim1 == 1)
783 : {
784 : mooseAssert(pi2->_side->dim() == 1, "Incompatible dimensions.");
785 : mooseAssert(off_edge_nodes1.size() < 2 && off_edge_nodes2.size() < 2,
786 : "off_edge_nodes size should be <2 for 2D contact");
787 0 : if (off_edge_nodes1.size() == 1 && off_edge_nodes2.size() == 1 &&
788 0 : off_edge_nodes1[0] == off_edge_nodes2[0])
789 0 : common_edge = COMMON_EDGE;
790 : }
791 : else
792 : {
793 : mooseAssert(dim1 == 2 && pi2->_side->dim() == 2, "Incompatible dimensions.");
794 : mooseAssert(off_edge_nodes1.size() < 3 && off_edge_nodes2.size() < 3,
795 : "off_edge_nodes size should be <3 for 3D contact");
796 714 : if (off_edge_nodes1.size() == 1)
797 : {
798 364 : if (off_edge_nodes2.size() == 1)
799 : {
800 0 : if (off_edge_nodes1[0] == off_edge_nodes2[0])
801 0 : common_edge = COMMON_NODE;
802 : }
803 364 : else if (off_edge_nodes2.size() == 2)
804 : {
805 364 : if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[0] == off_edge_nodes2[1])
806 364 : common_edge = EDGE_AND_COMMON_NODE;
807 : }
808 : }
809 350 : else if (off_edge_nodes1.size() == 2)
810 : {
811 350 : if (off_edge_nodes2.size() == 1)
812 : {
813 146 : if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[1] == off_edge_nodes2[0])
814 146 : common_edge = EDGE_AND_COMMON_NODE;
815 : }
816 204 : else if (off_edge_nodes2.size() == 2)
817 : {
818 204 : if ((off_edge_nodes1[0] == off_edge_nodes2[0] &&
819 408 : off_edge_nodes1[1] == off_edge_nodes2[1]) ||
820 204 : (off_edge_nodes1[1] == off_edge_nodes2[0] && off_edge_nodes1[0] == off_edge_nodes2[1]))
821 0 : common_edge = COMMON_EDGE;
822 : }
823 : }
824 : }
825 714 : return common_edge;
826 : }
827 :
828 : bool
829 243781 : PenetrationThread::findRidgeContactPoint(Point & contact_point,
830 : Real & tangential_distance,
831 : const Node *& closest_node,
832 : unsigned int & index,
833 : Point & contact_point_ref,
834 : std::vector<PenetrationInfo *> & p_info,
835 : const unsigned int index1,
836 : const unsigned int index2)
837 : {
838 243781 : tangential_distance = 0.0;
839 243781 : closest_node = NULL;
840 243781 : PenetrationInfo * pi1 = p_info[index1];
841 243781 : PenetrationInfo * pi2 = p_info[index2];
842 243781 : const unsigned sidedim(pi1->_side->dim());
843 : mooseAssert(sidedim == pi2->_side->dim(), "Incompatible dimensionalities");
844 :
845 : // Nodes on faces for the two interactions
846 243781 : std::vector<const Node *> side1_nodes;
847 243781 : getSideCornerNodes(pi1->_side, side1_nodes);
848 243781 : std::vector<const Node *> side2_nodes;
849 243781 : getSideCornerNodes(pi2->_side, side2_nodes);
850 :
851 243781 : std::sort(side1_nodes.begin(), side1_nodes.end());
852 243781 : std::sort(side2_nodes.begin(), side2_nodes.end());
853 :
854 : // Find nodes shared by the two faces
855 243781 : std::vector<const Node *> common_nodes;
856 243781 : std::set_intersection(side1_nodes.begin(),
857 : side1_nodes.end(),
858 : side2_nodes.begin(),
859 : side2_nodes.end(),
860 : std::inserter(common_nodes, common_nodes.end()));
861 :
862 243781 : if (common_nodes.size() != sidedim)
863 36941 : return false;
864 :
865 : bool found_point1, found_point2;
866 206840 : Point closest_coor_ref1(pi1->_closest_point_ref);
867 : const Node * closest_node1;
868 206840 : found_point1 = restrictPointToSpecifiedEdgeOfFace(
869 : closest_coor_ref1, closest_node1, pi1->_side, common_nodes);
870 :
871 206840 : Point closest_coor_ref2(pi2->_closest_point_ref);
872 : const Node * closest_node2;
873 206840 : found_point2 = restrictPointToSpecifiedEdgeOfFace(
874 : closest_coor_ref2, closest_node2, pi2->_side, common_nodes);
875 :
876 206840 : if (!found_point1 || !found_point2)
877 117650 : return false;
878 :
879 : // if (sidedim == 2)
880 : // {
881 : // TODO:
882 : // We have the parametric coordinates of the closest intersection point for both faces.
883 : // We need to find a point somewhere in the middle of them so there's not an abrupt jump.
884 : // Find that point by taking dot products of vector from contact point to secondary node point
885 : // with face normal vectors to see which face we're closer to.
886 : // }
887 :
888 89190 : FEBase * fe = NULL;
889 178380 : std::vector<Point> points(1);
890 :
891 : // We have to pick one of the two faces to own the contact point. Either one would
892 : // generally work, but to avoid some corner-case numerical roundoff issues when
893 : // determining signs of the normal and displacement, pick the one that the point is
894 : // closer to projecting onto, which is the one with the larger projected distance.
895 : // If that distance is the same for both faces (within numerical precision), pick the
896 : // face with the lowest index for repeatability.
897 153156 : if (MooseUtils::absoluteFuzzyGreaterThan(std::abs(pi1->_distance), std::abs(pi2->_distance)) ||
898 153156 : (MooseUtils::absoluteFuzzyEqual(std::abs(pi1->_distance), std::abs(pi2->_distance)) &&
899 : index1 < index2))
900 : {
901 63936 : fe = _fes[_tid][pi1->_side->dim()];
902 63936 : contact_point_ref = closest_coor_ref1;
903 63936 : points[0] = closest_coor_ref1;
904 63936 : fe->reinit(pi1->_side, &points);
905 63936 : index = index1;
906 : }
907 : else
908 : {
909 25254 : fe = _fes[_tid][pi2->_side->dim()];
910 25254 : contact_point_ref = closest_coor_ref2;
911 25254 : points[0] = closest_coor_ref2;
912 25254 : fe->reinit(pi2->_side, &points);
913 25254 : index = index2;
914 : }
915 :
916 89190 : contact_point = fe->get_xyz()[0];
917 :
918 89190 : if (sidedim == 2)
919 : {
920 84570 : if (closest_node1) // point is off the ridge between the two elements
921 : {
922 : mooseAssert((closest_node1 == closest_node2 || closest_node2 == NULL),
923 : "If off edge of ridge, closest node must be the same on both elements");
924 51751 : closest_node = closest_node1;
925 :
926 51751 : RealGradient off_face = *closest_node1 - contact_point;
927 51751 : tangential_distance = off_face.norm();
928 : }
929 : }
930 :
931 89190 : return true;
932 243781 : }
933 :
934 : void
935 487562 : PenetrationThread::getSideCornerNodes(const Elem * side, std::vector<const Node *> & corner_nodes)
936 : {
937 487562 : const ElemType t(side->type());
938 487562 : corner_nodes.clear();
939 :
940 487562 : corner_nodes.push_back(side->node_ptr(0));
941 487562 : corner_nodes.push_back(side->node_ptr(1));
942 487562 : switch (t)
943 : {
944 31968 : case EDGE2:
945 : case EDGE3:
946 : case EDGE4:
947 : {
948 31968 : break;
949 : }
950 :
951 14180 : case TRI3:
952 : case TRI6:
953 : case TRI7:
954 : {
955 14180 : corner_nodes.push_back(side->node_ptr(2));
956 14180 : break;
957 : }
958 :
959 441414 : case QUAD4:
960 : case QUAD8:
961 : case QUAD9:
962 : {
963 441414 : corner_nodes.push_back(side->node_ptr(2));
964 441414 : corner_nodes.push_back(side->node_ptr(3));
965 441414 : break;
966 : }
967 :
968 0 : default:
969 : {
970 0 : mooseError("Unsupported face type: ", t);
971 : break;
972 : }
973 : }
974 487562 : }
975 :
976 : bool
977 413680 : PenetrationThread::restrictPointToSpecifiedEdgeOfFace(Point & p,
978 : const Node *& closest_node,
979 : const Elem * side,
980 : const std::vector<const Node *> & edge_nodes)
981 : {
982 413680 : const ElemType t = side->type();
983 413680 : Real & xi = p(0);
984 413680 : Real & eta = p(1);
985 413680 : closest_node = NULL;
986 :
987 413680 : std::vector<unsigned int> local_node_indices;
988 1209072 : for (const auto & edge_node : edge_nodes)
989 : {
990 795392 : unsigned int local_index = side->get_node_index(edge_node);
991 795392 : if (local_index == libMesh::invalid_uint)
992 0 : mooseError("Side does not contain node");
993 795392 : local_node_indices.push_back(local_index);
994 : }
995 : mooseAssert(local_node_indices.size() == side->dim(),
996 : "Number of edge nodes must match side dimensionality");
997 413680 : std::sort(local_node_indices.begin(), local_node_indices.end());
998 :
999 413680 : bool off_of_this_edge = false;
1000 :
1001 413680 : switch (t)
1002 : {
1003 31968 : case EDGE2:
1004 : case EDGE3:
1005 : case EDGE4:
1006 : {
1007 31968 : if (local_node_indices[0] == 0)
1008 : {
1009 15984 : if (xi <= -1.0)
1010 : {
1011 11931 : xi = -1.0;
1012 11931 : off_of_this_edge = true;
1013 11931 : closest_node = side->node_ptr(0);
1014 : }
1015 : }
1016 15984 : else if (local_node_indices[0] == 1)
1017 : {
1018 15984 : if (xi >= 1.0)
1019 : {
1020 8662 : xi = 1.0;
1021 8662 : off_of_this_edge = true;
1022 8662 : closest_node = side->node_ptr(1);
1023 : }
1024 : }
1025 : else
1026 : {
1027 0 : mooseError("Invalid local node indices");
1028 : }
1029 31968 : break;
1030 : }
1031 :
1032 6210 : case TRI3:
1033 : case TRI6:
1034 : case TRI7:
1035 : {
1036 6210 : if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
1037 : {
1038 1823 : if (eta <= 0.0)
1039 : {
1040 983 : eta = 0.0;
1041 983 : off_of_this_edge = true;
1042 983 : if (xi < 0.0)
1043 223 : closest_node = side->node_ptr(0);
1044 760 : else if (xi > 1.0)
1045 310 : closest_node = side->node_ptr(1);
1046 : }
1047 : }
1048 4387 : else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
1049 : {
1050 2169 : if ((xi + eta) > 1.0)
1051 : {
1052 1015 : Real delta = (xi + eta - 1.0) / 2.0;
1053 1015 : xi -= delta;
1054 1015 : eta -= delta;
1055 1015 : off_of_this_edge = true;
1056 1015 : if (xi > 1.0)
1057 256 : closest_node = side->node_ptr(1);
1058 759 : else if (xi < 0.0)
1059 254 : closest_node = side->node_ptr(2);
1060 : }
1061 : }
1062 2218 : else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 2))
1063 : {
1064 2218 : if (xi <= 0.0)
1065 : {
1066 1086 : xi = 0.0;
1067 1086 : off_of_this_edge = true;
1068 1086 : if (eta > 1.0)
1069 345 : closest_node = side->node_ptr(2);
1070 741 : else if (eta < 0.0)
1071 286 : closest_node = side->node_ptr(0);
1072 : }
1073 : }
1074 : else
1075 : {
1076 0 : mooseError("Invalid local node indices");
1077 : }
1078 :
1079 6210 : break;
1080 : }
1081 :
1082 375502 : case QUAD4:
1083 : case QUAD8:
1084 : case QUAD9:
1085 : {
1086 375502 : if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
1087 : {
1088 162843 : if (eta <= -1.0)
1089 : {
1090 99880 : eta = -1.0;
1091 99880 : off_of_this_edge = true;
1092 99880 : if (xi < -1.0)
1093 31578 : closest_node = side->node_ptr(0);
1094 68302 : else if (xi > 1.0)
1095 31293 : closest_node = side->node_ptr(1);
1096 : }
1097 : }
1098 212659 : else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
1099 : {
1100 77527 : if (xi >= 1.0)
1101 : {
1102 63715 : xi = 1.0;
1103 63715 : off_of_this_edge = true;
1104 63715 : if (eta < -1.0)
1105 14586 : closest_node = side->node_ptr(1);
1106 49129 : else if (eta > 1.0)
1107 12173 : closest_node = side->node_ptr(2);
1108 : }
1109 : }
1110 135132 : else if ((local_node_indices[0] == 2) && (local_node_indices[1] == 3))
1111 : {
1112 111127 : if (eta >= 1.0)
1113 : {
1114 82329 : eta = 1.0;
1115 82329 : off_of_this_edge = true;
1116 82329 : if (xi < -1.0)
1117 35800 : closest_node = side->node_ptr(3);
1118 46529 : else if (xi > 1.0)
1119 41771 : closest_node = side->node_ptr(2);
1120 : }
1121 : }
1122 24005 : else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 3))
1123 : {
1124 24005 : if (xi <= -1.0)
1125 : {
1126 16177 : xi = -1.0;
1127 16177 : off_of_this_edge = true;
1128 16177 : if (eta < -1.0)
1129 9600 : closest_node = side->node_ptr(0);
1130 6577 : else if (eta > 1.0)
1131 2637 : closest_node = side->node_ptr(3);
1132 : }
1133 : }
1134 : else
1135 : {
1136 0 : mooseError("Invalid local node indices");
1137 : }
1138 375502 : break;
1139 : }
1140 :
1141 0 : default:
1142 : {
1143 0 : mooseError("Unsupported face type: ", t);
1144 : break;
1145 : }
1146 : }
1147 413680 : return off_of_this_edge;
1148 413680 : }
1149 :
1150 : bool
1151 4029 : PenetrationThread::restrictPointToFace(Point & p, const Node *& closest_node, const Elem * side)
1152 : {
1153 4029 : const ElemType t(side->type());
1154 4029 : Real & xi = p(0);
1155 4029 : Real & eta = p(1);
1156 4029 : closest_node = NULL;
1157 :
1158 4029 : bool off_of_this_face(false);
1159 :
1160 4029 : switch (t)
1161 : {
1162 0 : case EDGE2:
1163 : case EDGE3:
1164 : case EDGE4:
1165 : {
1166 0 : if (xi < -1.0)
1167 : {
1168 0 : xi = -1.0;
1169 0 : off_of_this_face = true;
1170 0 : closest_node = side->node_ptr(0);
1171 : }
1172 0 : else if (xi > 1.0)
1173 : {
1174 0 : xi = 1.0;
1175 0 : off_of_this_face = true;
1176 0 : closest_node = side->node_ptr(1);
1177 : }
1178 0 : break;
1179 : }
1180 :
1181 0 : case TRI3:
1182 : case TRI6:
1183 : case TRI7:
1184 : {
1185 0 : if (eta < 0.0)
1186 : {
1187 0 : eta = 0.0;
1188 0 : off_of_this_face = true;
1189 0 : if (xi < 0.5)
1190 : {
1191 0 : closest_node = side->node_ptr(0);
1192 0 : if (xi < 0.0)
1193 0 : xi = 0.0;
1194 : }
1195 : else
1196 : {
1197 0 : closest_node = side->node_ptr(1);
1198 0 : if (xi > 1.0)
1199 0 : xi = 1.0;
1200 : }
1201 : }
1202 0 : else if ((xi + eta) > 1.0)
1203 : {
1204 0 : Real delta = (xi + eta - 1.0) / 2.0;
1205 0 : xi -= delta;
1206 0 : eta -= delta;
1207 0 : off_of_this_face = true;
1208 0 : if (xi > 0.5)
1209 : {
1210 0 : closest_node = side->node_ptr(1);
1211 0 : if (xi > 1.0)
1212 : {
1213 0 : xi = 1.0;
1214 0 : eta = 0.0;
1215 : }
1216 : }
1217 : else
1218 : {
1219 0 : closest_node = side->node_ptr(2);
1220 0 : if (xi < 0.0)
1221 : {
1222 0 : xi = 0.0;
1223 0 : eta = 1.0;
1224 : }
1225 : }
1226 : }
1227 0 : else if (xi < 0.0)
1228 : {
1229 0 : xi = 0.0;
1230 0 : off_of_this_face = true;
1231 0 : if (eta > 0.5)
1232 : {
1233 0 : closest_node = side->node_ptr(2);
1234 0 : if (eta > 1.0)
1235 0 : eta = 1.0;
1236 : }
1237 : else
1238 : {
1239 0 : closest_node = side->node_ptr(0);
1240 0 : if (eta < 0.0)
1241 0 : eta = 0.0;
1242 : }
1243 : }
1244 0 : break;
1245 : }
1246 :
1247 4029 : case QUAD4:
1248 : case QUAD8:
1249 : case QUAD9:
1250 : {
1251 4029 : if (eta < -1.0)
1252 : {
1253 198 : eta = -1.0;
1254 198 : off_of_this_face = true;
1255 198 : if (xi < 0.0)
1256 : {
1257 8 : closest_node = side->node_ptr(0);
1258 8 : if (xi < -1.0)
1259 8 : xi = -1.0;
1260 : }
1261 : else
1262 : {
1263 190 : closest_node = side->node_ptr(1);
1264 190 : if (xi > 1.0)
1265 190 : xi = 1.0;
1266 : }
1267 : }
1268 3831 : else if (xi > 1.0)
1269 : {
1270 3831 : xi = 1.0;
1271 3831 : off_of_this_face = true;
1272 3831 : if (eta < 0.0)
1273 : {
1274 0 : closest_node = side->node_ptr(1);
1275 0 : if (eta < -1.0)
1276 0 : eta = -1.0;
1277 : }
1278 : else
1279 : {
1280 3831 : closest_node = side->node_ptr(2);
1281 3831 : if (eta > 1.0)
1282 938 : eta = 1.0;
1283 : }
1284 : }
1285 0 : else if (eta > 1.0)
1286 : {
1287 0 : eta = 1.0;
1288 0 : off_of_this_face = true;
1289 0 : if (xi < 0.0)
1290 : {
1291 0 : closest_node = side->node_ptr(3);
1292 0 : if (xi < -1.0)
1293 0 : xi = -1.0;
1294 : }
1295 : else
1296 : {
1297 0 : closest_node = side->node_ptr(2);
1298 0 : if (xi > 1.0)
1299 0 : xi = 1.0;
1300 : }
1301 : }
1302 0 : else if (xi < -1.0)
1303 : {
1304 0 : xi = -1.0;
1305 0 : off_of_this_face = true;
1306 0 : if (eta < 0.0)
1307 : {
1308 0 : closest_node = side->node_ptr(0);
1309 0 : if (eta < -1.0)
1310 0 : eta = -1.0;
1311 : }
1312 : else
1313 : {
1314 0 : closest_node = side->node_ptr(3);
1315 0 : if (eta > 1.0)
1316 0 : eta = 1.0;
1317 : }
1318 : }
1319 4029 : break;
1320 : }
1321 :
1322 0 : default:
1323 : {
1324 0 : mooseError("Unsupported face type: ", t);
1325 : break;
1326 : }
1327 : }
1328 4029 : return off_of_this_face;
1329 : }
1330 :
1331 : bool
1332 765513 : PenetrationThread::isFaceReasonableCandidate(const Elem * primary_elem,
1333 : const Elem * side,
1334 : FEBase * fe,
1335 : const Point * secondary_point,
1336 : const Real tangential_tolerance)
1337 : {
1338 765513 : unsigned int dim = primary_elem->dim();
1339 :
1340 765513 : const std::vector<Point> & phys_point = fe->get_xyz();
1341 :
1342 765513 : const std::vector<RealGradient> & dxyz_dxi = fe->get_dxyzdxi();
1343 765513 : const std::vector<RealGradient> & dxyz_deta = fe->get_dxyzdeta();
1344 :
1345 765513 : Point ref_point;
1346 :
1347 765513 : std::vector<Point> points(1); // Default constructor gives us a point at 0,0,0
1348 :
1349 765513 : fe->reinit(side, &points);
1350 :
1351 765513 : RealGradient d = *secondary_point - phys_point[0];
1352 :
1353 765513 : const Real twosqrt2 = 2.8284; // way more precision than we actually need here
1354 765513 : Real max_face_length = side->hmax() + twosqrt2 * tangential_tolerance;
1355 :
1356 765513 : RealVectorValue normal;
1357 765513 : if (dim - 1 == 2)
1358 : {
1359 690363 : normal = dxyz_dxi[0].cross(dxyz_deta[0]);
1360 : }
1361 75150 : else if (dim - 1 == 1)
1362 : {
1363 75131 : const Node * const * elem_nodes = primary_elem->get_nodes();
1364 75131 : const Point in_plane_vector1 = *elem_nodes[1] - *elem_nodes[0];
1365 75131 : const Point in_plane_vector2 = *elem_nodes[2] - *elem_nodes[0];
1366 :
1367 75131 : Point out_of_plane_normal = in_plane_vector1.cross(in_plane_vector2);
1368 75131 : out_of_plane_normal /= out_of_plane_normal.norm();
1369 :
1370 75131 : normal = dxyz_dxi[0].cross(out_of_plane_normal);
1371 : }
1372 : else
1373 : {
1374 19 : return true;
1375 : }
1376 765494 : normal /= normal.norm();
1377 :
1378 765494 : const Real dot(d * normal);
1379 :
1380 765494 : const RealGradient normcomp = dot * normal;
1381 765494 : const RealGradient tangcomp = d - normcomp;
1382 :
1383 765494 : const Real tangdist = tangcomp.norm();
1384 :
1385 : // Increase the size of the zone that we consider if the vector from the face
1386 : // to the node has a larger normal component
1387 765494 : const Real faceExpansionFactor = 2.0 * (1.0 + normcomp.norm() / d.norm());
1388 :
1389 765494 : bool isReasonableCandidate = true;
1390 765494 : if (tangdist > faceExpansionFactor * max_face_length)
1391 : {
1392 8068 : isReasonableCandidate = false;
1393 : }
1394 765494 : return isReasonableCandidate;
1395 765513 : }
1396 :
1397 : void
1398 669662 : PenetrationThread::computeSlip(FEBase & fe, PenetrationInfo & info)
1399 : {
1400 : // Slip is current projected position of secondary node minus
1401 : // original projected position of secondary node
1402 669662 : std::vector<Point> points(1);
1403 669662 : points[0] = info._starting_closest_point_ref;
1404 669662 : const auto & side = _elem_side_builder(*info._starting_elem, info._starting_side_num);
1405 669662 : fe.reinit(&side, &points);
1406 669662 : const std::vector<Point> & starting_point = fe.get_xyz();
1407 669662 : info._incremental_slip = info._closest_point - starting_point[0];
1408 669662 : if (info.isCaptured())
1409 : {
1410 50896 : info._frictional_energy =
1411 50896 : info._frictional_energy_old + info._contact_force * info._incremental_slip;
1412 50896 : info._accumulated_slip = info._accumulated_slip_old + info._incremental_slip.norm();
1413 : }
1414 669662 : }
1415 :
1416 : void
1417 669662 : PenetrationThread::smoothNormal(PenetrationInfo * info,
1418 : std::vector<PenetrationInfo *> & p_info,
1419 : const Node & node)
1420 : {
1421 669662 : if (_do_normal_smoothing)
1422 : {
1423 325614 : if (_normal_smoothing_method == PenetrationLocator::NSM_EDGE_BASED)
1424 : {
1425 : // If we are within the smoothing distance of any edges or corners, find the
1426 : // corner nodes for those edges/corners, and weights from distance to edge/corner
1427 260130 : std::vector<Real> edge_face_weights;
1428 260130 : std::vector<PenetrationInfo *> edge_face_info;
1429 :
1430 260130 : getSmoothingFacesAndWeights(info, edge_face_info, edge_face_weights, p_info, node);
1431 :
1432 : mooseAssert(edge_face_info.size() == edge_face_weights.size(),
1433 : "edge_face_info.size() != edge_face_weights.size()");
1434 :
1435 260130 : if (edge_face_info.size() > 0)
1436 : {
1437 : // Smooth the normal using the weighting functions for all participating faces.
1438 116445 : RealVectorValue new_normal;
1439 116445 : Real this_face_weight = 1.0;
1440 :
1441 262116 : for (unsigned int efwi = 0; efwi < edge_face_weights.size(); ++efwi)
1442 : {
1443 145671 : PenetrationInfo * npi = edge_face_info[efwi];
1444 145671 : if (npi)
1445 145671 : new_normal += npi->_normal * edge_face_weights[efwi];
1446 :
1447 145671 : this_face_weight -= edge_face_weights[efwi];
1448 : }
1449 : mooseAssert(this_face_weight >= (0.25 - 1e-8),
1450 : "Sum of weights of other faces shouldn't exceed 0.75");
1451 116445 : new_normal += info->_normal * this_face_weight;
1452 :
1453 116445 : const Real len = new_normal.norm();
1454 116445 : if (len > 0)
1455 116445 : new_normal /= len;
1456 :
1457 116445 : info->_normal = new_normal;
1458 : }
1459 260130 : }
1460 65484 : else if (_normal_smoothing_method == PenetrationLocator::NSM_NODAL_NORMAL_BASED)
1461 : {
1462 : // params.addParam<VariableName>("var_name","description");
1463 : // getParam<VariableName>("var_name")
1464 65484 : info->_normal(0) = _nodal_normal_x->getValue(info->_side, info->_side_phi);
1465 65484 : info->_normal(1) = _nodal_normal_y->getValue(info->_side, info->_side_phi);
1466 65484 : info->_normal(2) = _nodal_normal_z->getValue(info->_side, info->_side_phi);
1467 65484 : const Real len(info->_normal.norm());
1468 65484 : if (len > 0)
1469 64736 : info->_normal /= len;
1470 : }
1471 : }
1472 669662 : }
1473 :
1474 : void
1475 260130 : PenetrationThread::getSmoothingFacesAndWeights(PenetrationInfo * info,
1476 : std::vector<PenetrationInfo *> & edge_face_info,
1477 : std::vector<Real> & edge_face_weights,
1478 : std::vector<PenetrationInfo *> & p_info,
1479 : const Node & secondary_node)
1480 : {
1481 260130 : const Elem * side = info->_side;
1482 260130 : const Point & p = info->_closest_point_ref;
1483 260130 : std::set<dof_id_type> elems_to_exclude;
1484 260130 : elems_to_exclude.insert(info->_elem->id());
1485 :
1486 260130 : std::vector<std::vector<const Node *>> edge_nodes;
1487 :
1488 : // Get the pairs of nodes along every edge that we are close enough to smooth with
1489 260130 : getSmoothingEdgeNodesAndWeights(p, side, edge_nodes, edge_face_weights);
1490 260130 : std::vector<Elem *> edge_neighbor_elems;
1491 260130 : edge_face_info.resize(edge_nodes.size(), NULL);
1492 :
1493 260130 : std::vector<unsigned int> edges_without_neighbors;
1494 :
1495 486408 : for (unsigned int i = 0; i < edge_nodes.size(); ++i)
1496 : {
1497 : // Sort all sets of edge nodes (needed for comparing edges)
1498 226278 : std::sort(edge_nodes[i].begin(), edge_nodes[i].end());
1499 :
1500 226278 : std::vector<PenetrationInfo *> face_info_comm_edge;
1501 226278 : getInfoForFacesWithCommonNodes(
1502 226278 : &secondary_node, elems_to_exclude, edge_nodes[i], face_info_comm_edge, p_info);
1503 :
1504 226278 : if (face_info_comm_edge.size() == 0)
1505 95220 : edges_without_neighbors.push_back(i);
1506 131058 : else if (face_info_comm_edge.size() > 1)
1507 0 : mooseError("Only one neighbor allowed per edge");
1508 : else
1509 131058 : edge_face_info[i] = face_info_comm_edge[0];
1510 226278 : }
1511 :
1512 : // Remove edges without neighbors from the vector, starting from end
1513 260130 : std::vector<unsigned int>::reverse_iterator rit;
1514 355350 : for (rit = edges_without_neighbors.rbegin(); rit != edges_without_neighbors.rend(); ++rit)
1515 : {
1516 95220 : unsigned int index = *rit;
1517 95220 : edge_nodes.erase(edge_nodes.begin() + index);
1518 95220 : edge_face_weights.erase(edge_face_weights.begin() + index);
1519 95220 : edge_face_info.erase(edge_face_info.begin() + index);
1520 : }
1521 :
1522 : // Handle corner case
1523 260130 : if (edge_nodes.size() > 1)
1524 : {
1525 14613 : if (edge_nodes.size() != 2)
1526 0 : mooseError("Invalid number of smoothing edges");
1527 :
1528 : // find common node
1529 14613 : std::vector<const Node *> common_nodes;
1530 58452 : std::set_intersection(edge_nodes[0].begin(),
1531 14613 : edge_nodes[0].end(),
1532 14613 : edge_nodes[1].begin(),
1533 14613 : edge_nodes[1].end(),
1534 : std::inserter(common_nodes, common_nodes.end()));
1535 :
1536 14613 : if (common_nodes.size() != 1)
1537 0 : mooseError("Invalid number of common nodes");
1538 :
1539 43839 : for (const auto & pinfo : edge_face_info)
1540 29226 : elems_to_exclude.insert(pinfo->_elem->id());
1541 :
1542 14613 : std::vector<PenetrationInfo *> face_info_comm_edge;
1543 14613 : getInfoForFacesWithCommonNodes(
1544 : &secondary_node, elems_to_exclude, common_nodes, face_info_comm_edge, p_info);
1545 :
1546 14613 : unsigned int num_corner_neighbors = face_info_comm_edge.size();
1547 :
1548 14613 : if (num_corner_neighbors > 0)
1549 : {
1550 14613 : Real fw0 = edge_face_weights[0];
1551 14613 : Real fw1 = edge_face_weights[1];
1552 :
1553 : // Corner weight is product of edge weights. Spread out over multiple neighbors.
1554 14613 : Real fw_corner = (fw0 * fw1) / static_cast<Real>(num_corner_neighbors);
1555 :
1556 : // Adjust original edge weights
1557 14613 : edge_face_weights[0] *= (1.0 - fw1);
1558 14613 : edge_face_weights[1] *= (1.0 - fw0);
1559 :
1560 29226 : for (unsigned int i = 0; i < num_corner_neighbors; ++i)
1561 : {
1562 14613 : edge_face_weights.push_back(fw_corner);
1563 14613 : edge_face_info.push_back(face_info_comm_edge[i]);
1564 : }
1565 : }
1566 14613 : }
1567 260130 : }
1568 :
1569 : void
1570 260130 : PenetrationThread::getSmoothingEdgeNodesAndWeights(
1571 : const Point & p,
1572 : const Elem * side,
1573 : std::vector<std::vector<const Node *>> & edge_nodes,
1574 : std::vector<Real> & edge_face_weights)
1575 : {
1576 260130 : const ElemType t(side->type());
1577 260130 : const Real & xi = p(0);
1578 260130 : const Real & eta = p(1);
1579 :
1580 260130 : Real smooth_limit = 1.0 - _normal_smoothing_distance;
1581 :
1582 260130 : switch (t)
1583 : {
1584 9970 : case EDGE2:
1585 : case EDGE3:
1586 : case EDGE4:
1587 : {
1588 9970 : if (xi < -smooth_limit)
1589 : {
1590 1600 : std::vector<const Node *> en;
1591 1600 : en.push_back(side->node_ptr(0));
1592 1600 : edge_nodes.push_back(en);
1593 1600 : Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
1594 1600 : if (fw > 0.5)
1595 132 : fw = 0.5;
1596 1600 : edge_face_weights.push_back(fw);
1597 1600 : }
1598 8370 : else if (xi > smooth_limit)
1599 : {
1600 1538 : std::vector<const Node *> en;
1601 1538 : en.push_back(side->node_ptr(1));
1602 1538 : edge_nodes.push_back(en);
1603 1538 : Real fw = 0.5 - (1.0 - xi) / (2.0 * _normal_smoothing_distance);
1604 1538 : if (fw > 0.5)
1605 248 : fw = 0.5;
1606 1538 : edge_face_weights.push_back(fw);
1607 1538 : }
1608 9970 : break;
1609 : }
1610 :
1611 0 : case TRI3:
1612 : case TRI6:
1613 : case TRI7:
1614 : {
1615 0 : if (eta < -smooth_limit)
1616 : {
1617 0 : std::vector<const Node *> en;
1618 0 : en.push_back(side->node_ptr(0));
1619 0 : en.push_back(side->node_ptr(1));
1620 0 : edge_nodes.push_back(en);
1621 0 : Real fw = 0.5 - (1.0 + eta) / (2.0 * _normal_smoothing_distance);
1622 0 : if (fw > 0.5)
1623 0 : fw = 0.5;
1624 0 : edge_face_weights.push_back(fw);
1625 0 : }
1626 0 : if ((xi + eta) > smooth_limit)
1627 : {
1628 0 : std::vector<const Node *> en;
1629 0 : en.push_back(side->node_ptr(1));
1630 0 : en.push_back(side->node_ptr(2));
1631 0 : edge_nodes.push_back(en);
1632 0 : Real fw = 0.5 - (1.0 - xi - eta) / (2.0 * _normal_smoothing_distance);
1633 0 : if (fw > 0.5)
1634 0 : fw = 0.5;
1635 0 : edge_face_weights.push_back(fw);
1636 0 : }
1637 0 : if (xi < -smooth_limit)
1638 : {
1639 0 : std::vector<const Node *> en;
1640 0 : en.push_back(side->node_ptr(2));
1641 0 : en.push_back(side->node_ptr(0));
1642 0 : edge_nodes.push_back(en);
1643 0 : Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
1644 0 : if (fw > 0.5)
1645 0 : fw = 0.5;
1646 0 : edge_face_weights.push_back(fw);
1647 0 : }
1648 0 : break;
1649 : }
1650 :
1651 250160 : case QUAD4:
1652 : case QUAD8:
1653 : case QUAD9:
1654 : {
1655 250160 : if (eta < -smooth_limit)
1656 : {
1657 51736 : std::vector<const Node *> en;
1658 51736 : en.push_back(side->node_ptr(0));
1659 51736 : en.push_back(side->node_ptr(1));
1660 51736 : edge_nodes.push_back(en);
1661 51736 : Real fw = 0.5 - (1.0 + eta) / (2.0 * _normal_smoothing_distance);
1662 51736 : if (fw > 0.5)
1663 13138 : fw = 0.5;
1664 51736 : edge_face_weights.push_back(fw);
1665 51736 : }
1666 250160 : if (xi > smooth_limit)
1667 : {
1668 66116 : std::vector<const Node *> en;
1669 66116 : en.push_back(side->node_ptr(1));
1670 66116 : en.push_back(side->node_ptr(2));
1671 66116 : edge_nodes.push_back(en);
1672 66116 : Real fw = 0.5 - (1.0 - xi) / (2.0 * _normal_smoothing_distance);
1673 66116 : if (fw > 0.5)
1674 17602 : fw = 0.5;
1675 66116 : edge_face_weights.push_back(fw);
1676 66116 : }
1677 250160 : if (eta > smooth_limit)
1678 : {
1679 73969 : std::vector<const Node *> en;
1680 73969 : en.push_back(side->node_ptr(2));
1681 73969 : en.push_back(side->node_ptr(3));
1682 73969 : edge_nodes.push_back(en);
1683 73969 : Real fw = 0.5 - (1.0 - eta) / (2.0 * _normal_smoothing_distance);
1684 73969 : if (fw > 0.5)
1685 21483 : fw = 0.5;
1686 73969 : edge_face_weights.push_back(fw);
1687 73969 : }
1688 250160 : if (xi < -smooth_limit)
1689 : {
1690 31319 : std::vector<const Node *> en;
1691 31319 : en.push_back(side->node_ptr(3));
1692 31319 : en.push_back(side->node_ptr(0));
1693 31319 : edge_nodes.push_back(en);
1694 31319 : Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
1695 31319 : if (fw > 0.5)
1696 16695 : fw = 0.5;
1697 31319 : edge_face_weights.push_back(fw);
1698 31319 : }
1699 250160 : break;
1700 : }
1701 :
1702 0 : default:
1703 : {
1704 0 : mooseError("Unsupported face type: ", t);
1705 : break;
1706 : }
1707 : }
1708 260130 : }
1709 :
1710 : void
1711 240891 : PenetrationThread::getInfoForFacesWithCommonNodes(
1712 : const Node * secondary_node,
1713 : const std::set<dof_id_type> & elems_to_exclude,
1714 : const std::vector<const Node *> edge_nodes,
1715 : std::vector<PenetrationInfo *> & face_info_comm_edge,
1716 : std::vector<PenetrationInfo *> & p_info)
1717 : {
1718 : // elems connected to a node on this edge, find one that has the same corners as this, and is not
1719 : // the current elem
1720 240891 : auto node_to_elem_pair = _node_to_elem_map.find(edge_nodes[0]->id()); // just need one of the
1721 : // nodes
1722 : mooseAssert(node_to_elem_pair != _node_to_elem_map.end(), "Missing entry in node to elem map");
1723 240891 : const std::vector<dof_id_type> & elems_connected_to_node = node_to_elem_pair->second;
1724 :
1725 240891 : std::vector<const Elem *> elems_connected_to_edge;
1726 :
1727 786140 : for (unsigned int ecni = 0; ecni < elems_connected_to_node.size(); ecni++)
1728 : {
1729 545249 : if (elems_to_exclude.find(elems_connected_to_node[ecni]) != elems_to_exclude.end())
1730 270117 : continue;
1731 275132 : const Elem * elem = _mesh.elemPtr(elems_connected_to_node[ecni]);
1732 :
1733 275132 : std::vector<const Node *> nodevec;
1734 4281816 : for (unsigned int ni = 0; ni < elem->n_nodes(); ++ni)
1735 4006684 : if (elem->is_vertex(ni))
1736 2193304 : nodevec.push_back(elem->node_ptr(ni));
1737 :
1738 275132 : std::vector<const Node *> common_nodes;
1739 275132 : std::sort(nodevec.begin(), nodevec.end());
1740 275132 : std::set_intersection(edge_nodes.begin(),
1741 : edge_nodes.end(),
1742 : nodevec.begin(),
1743 : nodevec.end(),
1744 : std::inserter(common_nodes, common_nodes.end()));
1745 :
1746 275132 : if (common_nodes.size() == edge_nodes.size())
1747 145671 : elems_connected_to_edge.push_back(elem);
1748 275132 : }
1749 :
1750 240891 : if (elems_connected_to_edge.size() > 0)
1751 : {
1752 :
1753 : // There are potentially multiple elements that share a common edge
1754 : // 2D:
1755 : // There can only be one element on the same surface
1756 : // 3D:
1757 : // If there are two edge nodes, there can only be one element on the same surface
1758 : // If there is only one edge node (a corner), there could be multiple elements on the same
1759 : // surface
1760 145671 : bool allowMultipleNeighbors = false;
1761 :
1762 145671 : if (elems_connected_to_edge[0]->dim() == 3)
1763 : {
1764 143733 : if (edge_nodes.size() == 1)
1765 : {
1766 14613 : allowMultipleNeighbors = true;
1767 : }
1768 : }
1769 :
1770 160284 : for (unsigned int i = 0; i < elems_connected_to_edge.size(); ++i)
1771 : {
1772 145671 : std::vector<PenetrationInfo *> thisElemInfo;
1773 145671 : getInfoForElem(thisElemInfo, p_info, elems_connected_to_edge[i]);
1774 145671 : if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
1775 : {
1776 33406 : if (thisElemInfo.size() > 1)
1777 0 : mooseError(
1778 : "Found multiple neighbors to current edge/face on surface when only one is allowed");
1779 33406 : face_info_comm_edge.push_back(thisElemInfo[0]);
1780 33406 : break;
1781 : }
1782 :
1783 112265 : createInfoForElem(
1784 112265 : thisElemInfo, p_info, secondary_node, elems_connected_to_edge[i], edge_nodes);
1785 112265 : if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
1786 : {
1787 97652 : if (thisElemInfo.size() > 1)
1788 0 : mooseError(
1789 : "Found multiple neighbors to current edge/face on surface when only one is allowed");
1790 97652 : face_info_comm_edge.push_back(thisElemInfo[0]);
1791 97652 : break;
1792 : }
1793 :
1794 29226 : for (unsigned int j = 0; j < thisElemInfo.size(); ++j)
1795 14613 : face_info_comm_edge.push_back(thisElemInfo[j]);
1796 145671 : }
1797 : }
1798 240891 : }
1799 :
1800 : void
1801 145671 : PenetrationThread::getInfoForElem(std::vector<PenetrationInfo *> & thisElemInfo,
1802 : std::vector<PenetrationInfo *> & p_info,
1803 : const Elem * elem)
1804 : {
1805 289925 : for (const auto & pi : p_info)
1806 : {
1807 144254 : if (!pi)
1808 40109 : continue;
1809 :
1810 104145 : if (pi->_elem == elem)
1811 40109 : thisElemInfo.push_back(pi);
1812 : }
1813 145671 : }
1814 :
1815 : void
1816 890087 : PenetrationThread::createInfoForElem(std::vector<PenetrationInfo *> & thisElemInfo,
1817 : std::vector<PenetrationInfo *> & p_info,
1818 : const Node * secondary_node,
1819 : const Elem * elem,
1820 : const std::vector<const Node *> & nodes_that_must_be_on_side,
1821 : const bool check_whether_reasonable)
1822 : {
1823 890087 : const BoundaryInfo & boundary_info = _mesh.getMesh().get_boundary_info();
1824 :
1825 5958607 : for (auto s : elem->side_index_range())
1826 : {
1827 5083291 : if (!boundary_info.has_boundary_id(elem, s, _primary_boundary))
1828 4205513 : continue;
1829 :
1830 : // Don't create info for this side if one already exists
1831 877778 : bool already_have_info_this_side = false;
1832 877778 : for (const auto & pi : thisElemInfo)
1833 6703 : if (pi->_side_num == s)
1834 : {
1835 6703 : already_have_info_this_side = true;
1836 6703 : break;
1837 : }
1838 :
1839 877778 : if (already_have_info_this_side)
1840 14771 : break;
1841 :
1842 871075 : const Elem * side = elem->build_side_ptr(s).release();
1843 :
1844 : // Only continue with creating info for this side if the side contains
1845 : // all of the nodes in nodes_that_must_be_on_side
1846 871075 : std::vector<const Node *> nodevec;
1847 5552998 : for (unsigned int ni = 0; ni < side->n_nodes(); ++ni)
1848 4681923 : nodevec.push_back(side->node_ptr(ni));
1849 :
1850 871075 : std::sort(nodevec.begin(), nodevec.end());
1851 871075 : std::vector<const Node *> common_nodes;
1852 871075 : std::set_intersection(nodes_that_must_be_on_side.begin(),
1853 : nodes_that_must_be_on_side.end(),
1854 : nodevec.begin(),
1855 : nodevec.end(),
1856 : std::inserter(common_nodes, common_nodes.end()));
1857 871075 : if (common_nodes.size() != nodes_that_must_be_on_side.size())
1858 : {
1859 0 : delete side;
1860 0 : break;
1861 : }
1862 :
1863 871075 : FEBase * fe_elem = _fes[_tid][elem->dim()];
1864 871075 : FEBase * fe_side = _fes[_tid][side->dim()];
1865 :
1866 : // Optionally check to see whether face is reasonable candidate based on an
1867 : // estimate of how closely it is likely to project to the face
1868 871075 : if (check_whether_reasonable)
1869 765513 : if (!isFaceReasonableCandidate(elem, side, fe_side, secondary_node, _tangential_tolerance))
1870 : {
1871 8068 : delete side;
1872 8068 : break;
1873 : }
1874 :
1875 863007 : Point contact_phys;
1876 863007 : Point contact_ref;
1877 863007 : Point contact_on_face_ref;
1878 863007 : Real distance = 0.;
1879 863007 : Real tangential_distance = 0.;
1880 863007 : RealGradient normal;
1881 : bool contact_point_on_side;
1882 863007 : std::vector<const Node *> off_edge_nodes;
1883 863007 : std::vector<std::vector<Real>> side_phi;
1884 863007 : std::vector<std::vector<RealGradient>> side_grad_phi;
1885 863007 : std::vector<RealGradient> dxyzdxi;
1886 863007 : std::vector<RealGradient> dxyzdeta;
1887 863007 : std::vector<RealGradient> d2xyzdxideta;
1888 :
1889 : std::unique_ptr<PenetrationInfo> pen_info =
1890 : std::make_unique<PenetrationInfo>(elem,
1891 : side,
1892 : s,
1893 : normal,
1894 : distance,
1895 : tangential_distance,
1896 : contact_phys,
1897 : contact_ref,
1898 : contact_on_face_ref,
1899 : off_edge_nodes,
1900 : side_phi,
1901 : side_grad_phi,
1902 : dxyzdxi,
1903 : dxyzdeta,
1904 863007 : d2xyzdxideta);
1905 :
1906 863007 : bool search_succeeded = false;
1907 863007 : Moose::findContactPoint(*pen_info,
1908 : fe_elem,
1909 : fe_side,
1910 : _fe_type,
1911 : *secondary_node,
1912 : true,
1913 : _tangential_tolerance,
1914 : contact_point_on_side,
1915 : search_succeeded);
1916 :
1917 : // Do not add contact info from failed searches
1918 863007 : if (search_succeeded)
1919 : {
1920 862993 : thisElemInfo.push_back(pen_info.get());
1921 862993 : p_info.push_back(pen_info.release());
1922 : }
1923 879143 : }
1924 890087 : }
|
