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 : // Mutex to use when accessing _penetration_info;
27 : Threads::spin_mutex pinfo_mutex;
28 :
29 170144 : PenetrationThread::PenetrationThread(
30 : SubProblem & subproblem,
31 : const MooseMesh & mesh,
32 : BoundaryID primary_boundary,
33 : BoundaryID secondary_boundary,
34 : std::map<dof_id_type, PenetrationInfo *> & penetration_info,
35 : bool check_whether_reasonable,
36 : bool update_location,
37 : Real tangential_tolerance,
38 : bool do_normal_smoothing,
39 : Real normal_smoothing_distance,
40 : PenetrationLocator::NORMAL_SMOOTHING_METHOD normal_smoothing_method,
41 : std::vector<std::vector<FEBase *>> & fes,
42 : FEType & fe_type,
43 : NearestNodeLocator & nearest_node,
44 : const std::map<dof_id_type, std::vector<dof_id_type>> & node_to_elem_map,
45 170144 : const std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>> & bc_tuples)
46 170144 : : _subproblem(subproblem),
47 170144 : _mesh(mesh),
48 170144 : _primary_boundary(primary_boundary),
49 170144 : _secondary_boundary(secondary_boundary),
50 170144 : _penetration_info(penetration_info),
51 170144 : _check_whether_reasonable(check_whether_reasonable),
52 170144 : _update_location(update_location),
53 170144 : _tangential_tolerance(tangential_tolerance),
54 170144 : _do_normal_smoothing(do_normal_smoothing),
55 170144 : _normal_smoothing_distance(normal_smoothing_distance),
56 170144 : _normal_smoothing_method(normal_smoothing_method),
57 170144 : _nodal_normal_x(NULL),
58 170144 : _nodal_normal_y(NULL),
59 170144 : _nodal_normal_z(NULL),
60 170144 : _fes(fes),
61 170144 : _fe_type(fe_type),
62 170144 : _nearest_node(nearest_node),
63 170144 : _node_to_elem_map(node_to_elem_map),
64 170144 : _bc_tuples(bc_tuples)
65 : {
66 170144 : }
67 :
68 : // Splitting Constructor
69 16090 : PenetrationThread::PenetrationThread(PenetrationThread & x, Threads::split /*split*/)
70 16090 : : _subproblem(x._subproblem),
71 16090 : _mesh(x._mesh),
72 16090 : _primary_boundary(x._primary_boundary),
73 16090 : _secondary_boundary(x._secondary_boundary),
74 16090 : _penetration_info(x._penetration_info),
75 16090 : _check_whether_reasonable(x._check_whether_reasonable),
76 16090 : _update_location(x._update_location),
77 16090 : _tangential_tolerance(x._tangential_tolerance),
78 16090 : _do_normal_smoothing(x._do_normal_smoothing),
79 16090 : _normal_smoothing_distance(x._normal_smoothing_distance),
80 16090 : _normal_smoothing_method(x._normal_smoothing_method),
81 16090 : _fes(x._fes),
82 16090 : _fe_type(x._fe_type),
83 16090 : _nearest_node(x._nearest_node),
84 16090 : _node_to_elem_map(x._node_to_elem_map),
85 16090 : _bc_tuples(x._bc_tuples)
86 : {
87 16090 : }
88 :
89 : void
90 186234 : PenetrationThread::operator()(const NodeIdRange & range)
91 : {
92 186234 : ParallelUniqueId puid;
93 186234 : _tid = puid.id;
94 :
95 : // Must get the variables every time this is run because _tid can change
96 186234 : if (_do_normal_smoothing &&
97 70064 : _normal_smoothing_method == PenetrationLocator::NSM_NODAL_NORMAL_BASED)
98 : {
99 14176 : _nodal_normal_x = &_subproblem.getStandardVariable(_tid, "nodal_normal_x");
100 14176 : _nodal_normal_y = &_subproblem.getStandardVariable(_tid, "nodal_normal_y");
101 14176 : _nodal_normal_z = &_subproblem.getStandardVariable(_tid, "nodal_normal_z");
102 : }
103 :
104 1334938 : for (const auto & node_id : range)
105 : {
106 1148704 : const Node & node = _mesh.nodeRef(node_id);
107 :
108 : // We're going to get a reference to the pointer for the pinfo for this node
109 : // This will allow us to manipulate this pointer without having to go through
110 : // the _penetration_info map... meaning this is the only mutex we'll have to do!
111 1148704 : pinfo_mutex.lock();
112 1148704 : PenetrationInfo *& info = _penetration_info[node.id()];
113 1148704 : pinfo_mutex.unlock();
114 :
115 1148704 : std::vector<PenetrationInfo *> p_info;
116 1148704 : bool info_set(false);
117 :
118 : // See if we already have info about this node
119 1148704 : if (info)
120 : {
121 651652 : FEBase * fe_elem = _fes[_tid][info->_elem->dim()];
122 651652 : FEBase * fe_side = _fes[_tid][info->_side->dim()];
123 :
124 651652 : if (!_update_location && (info->_distance >= 0 || info->isCaptured()))
125 : {
126 0 : const Point contact_ref = info->_closest_point_ref;
127 0 : bool contact_point_on_side(false);
128 :
129 : // Secondary position must be the previous contact point
130 : // Use the previous reference coordinates
131 0 : std::vector<Point> points(1);
132 0 : points[0] = contact_ref;
133 0 : const std::vector<Point> & secondary_pos = fe_side->get_xyz();
134 0 : bool search_succeeded = false;
135 :
136 0 : Moose::findContactPoint(*info,
137 : fe_elem,
138 : fe_side,
139 : _fe_type,
140 0 : secondary_pos[0],
141 : false,
142 : _tangential_tolerance,
143 : contact_point_on_side,
144 : search_succeeded);
145 :
146 : // Restore the original reference coordinates
147 0 : info->_closest_point_ref = contact_ref;
148 : // Just calculated as the distance of the contact point off the surface (0). Set to 0 to
149 : // avoid round-off.
150 0 : info->_distance = 0.0;
151 0 : info_set = true;
152 0 : }
153 : else
154 : {
155 651652 : Real old_tangential_distance(info->_tangential_distance);
156 651652 : bool contact_point_on_side(false);
157 651652 : bool search_succeeded = false;
158 :
159 651652 : Moose::findContactPoint(*info,
160 : fe_elem,
161 : fe_side,
162 : _fe_type,
163 : node,
164 : false,
165 : _tangential_tolerance,
166 : contact_point_on_side,
167 : search_succeeded);
168 :
169 651652 : if (contact_point_on_side)
170 : {
171 591663 : if (info->_tangential_distance <= 0.0) // on the face
172 : {
173 506635 : info_set = true;
174 : }
175 85028 : else if (info->_tangential_distance > 0.0 && old_tangential_distance > 0.0)
176 : { // off the face but within tolerance, was that way on the last step too
177 82098 : if (info->_side->dim() == 2 && info->_off_edge_nodes.size() < 2)
178 : { // Closest point on face is on a node rather than an edge. Another
179 : // face might be a better candidate.
180 : }
181 : else
182 : {
183 80540 : info_set = true;
184 : }
185 : }
186 : }
187 : }
188 : }
189 :
190 1148704 : if (!info_set)
191 : {
192 561529 : const Node * closest_node = _nearest_node.nearestNode(node.id());
193 561529 : auto node_to_elem_pair = _node_to_elem_map.find(closest_node->id());
194 : mooseAssert(node_to_elem_pair != _node_to_elem_map.end(),
195 : "Missing entry in node to elem map");
196 561529 : const std::vector<dof_id_type> & closest_elems = node_to_elem_pair->second;
197 :
198 1329522 : for (const auto & elem_id : closest_elems)
199 : {
200 767993 : const Elem * elem = _mesh.elemPtr(elem_id);
201 :
202 767993 : std::vector<PenetrationInfo *> thisElemInfo;
203 767993 : std::vector<const Node *> nodesThatMustBeOnSide;
204 767993 : nodesThatMustBeOnSide.push_back(closest_node);
205 767993 : createInfoForElem(
206 767993 : thisElemInfo, p_info, &node, elem, nodesThatMustBeOnSide, _check_whether_reasonable);
207 767993 : }
208 :
209 561529 : if (p_info.size() == 1)
210 : {
211 399343 : if (p_info[0]->_tangential_distance <= _tangential_tolerance)
212 : {
213 9419 : switchInfo(info, p_info[0]);
214 9419 : info_set = true;
215 : }
216 : }
217 162186 : else if (p_info.size() > 1)
218 : {
219 :
220 : // Loop through all pairs of faces, and check for contact on ridge between each face pair
221 157731 : std::vector<RidgeData> ridgeDataVec;
222 350218 : for (unsigned int i = 0; i + 1 < p_info.size(); ++i)
223 438567 : for (unsigned int j = i + 1; j < p_info.size(); ++j)
224 : {
225 246080 : Point closest_coor;
226 246080 : Real tangential_distance(0.0);
227 246080 : const Node * closest_node_on_ridge = NULL;
228 246080 : unsigned int index = 0;
229 246080 : Point closest_coor_ref;
230 246080 : bool found_ridge_contact_point = findRidgeContactPoint(closest_coor,
231 : tangential_distance,
232 : closest_node_on_ridge,
233 : index,
234 : closest_coor_ref,
235 : p_info,
236 : i,
237 : j);
238 246080 : if (found_ridge_contact_point)
239 : {
240 90595 : RidgeData rpd;
241 90595 : rpd._closest_coor = closest_coor;
242 90595 : rpd._tangential_distance = tangential_distance;
243 90595 : rpd._closest_node = closest_node_on_ridge;
244 90595 : rpd._index = index;
245 90595 : rpd._closest_coor_ref = closest_coor_ref;
246 90595 : ridgeDataVec.push_back(rpd);
247 : }
248 : }
249 :
250 157731 : if (ridgeDataVec.size() > 0) // Either find the ridge pair that is the best or find a peak
251 : {
252 : // Group together ridges for which we are off the edge of a common node.
253 : // Those are peaks.
254 72338 : std::vector<RidgeSetData> ridgeSetDataVec;
255 162933 : for (unsigned int i = 0; i < ridgeDataVec.size(); ++i)
256 : {
257 90595 : bool foundSetWithMatchingNode = false;
258 107761 : for (unsigned int j = 0; j < ridgeSetDataVec.size(); ++j)
259 : {
260 32241 : if (ridgeDataVec[i]._closest_node != NULL &&
261 9200 : ridgeDataVec[i]._closest_node == ridgeSetDataVec[j]._closest_node)
262 : {
263 5875 : foundSetWithMatchingNode = true;
264 5875 : ridgeSetDataVec[j]._ridge_data_vec.push_back(ridgeDataVec[i]);
265 5875 : break;
266 : }
267 : }
268 90595 : if (!foundSetWithMatchingNode)
269 : {
270 84720 : RidgeSetData rsd;
271 84720 : rsd._distance = std::numeric_limits<Real>::max();
272 84720 : rsd._ridge_data_vec.push_back(ridgeDataVec[i]);
273 84720 : rsd._closest_node = ridgeDataVec[i]._closest_node;
274 84720 : ridgeSetDataVec.push_back(rsd);
275 84720 : }
276 : }
277 : // Compute distance to each set of ridges
278 157058 : for (unsigned int i = 0; i < ridgeSetDataVec.size(); ++i)
279 : {
280 84720 : if (ridgeSetDataVec[i]._closest_node !=
281 : NULL) // Either a peak or off the edge of single ridge
282 : {
283 46390 : if (ridgeSetDataVec[i]._ridge_data_vec.size() == 1) // off edge of single ridge
284 : {
285 41299 : if (ridgeSetDataVec[i]._ridge_data_vec[0]._tangential_distance <=
286 41299 : _tangential_tolerance) // off within tolerance
287 : {
288 13252 : ridgeSetDataVec[i]._closest_coor =
289 13252 : ridgeSetDataVec[i]._ridge_data_vec[0]._closest_coor;
290 13252 : Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
291 13252 : ridgeSetDataVec[i]._distance = contact_point_vec.norm();
292 : }
293 : }
294 : else // several ridges join at common node to make peak. The common node is the
295 : // contact point
296 : {
297 5091 : ridgeSetDataVec[i]._closest_coor = *ridgeSetDataVec[i]._closest_node;
298 5091 : Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
299 5091 : ridgeSetDataVec[i]._distance = contact_point_vec.norm();
300 : }
301 : }
302 : else // on a single ridge
303 : {
304 38330 : ridgeSetDataVec[i]._closest_coor =
305 38330 : ridgeSetDataVec[i]._ridge_data_vec[0]._closest_coor;
306 38330 : Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
307 38330 : ridgeSetDataVec[i]._distance = contact_point_vec.norm();
308 : }
309 : }
310 : // Find the set of ridges closest to us.
311 72338 : unsigned int closest_ridge_set_index(0);
312 72338 : Real closest_distance(ridgeSetDataVec[0]._distance);
313 72338 : Point closest_point(ridgeSetDataVec[0]._closest_coor);
314 84720 : for (unsigned int i = 1; i < ridgeSetDataVec.size(); ++i)
315 : {
316 12382 : if (ridgeSetDataVec[i]._distance < closest_distance)
317 : {
318 1979 : closest_ridge_set_index = i;
319 1979 : closest_distance = ridgeSetDataVec[i]._distance;
320 1979 : closest_point = ridgeSetDataVec[i]._closest_coor;
321 : }
322 : }
323 :
324 72338 : if (closest_distance <
325 72338 : std::numeric_limits<Real>::max()) // contact point is on the closest ridge set
326 : {
327 : // find the face in the ridge set with the smallest index, assign that one to the
328 : // interaction
329 46369 : unsigned int face_index(std::numeric_limits<unsigned int>::max());
330 97671 : for (unsigned int i = 0;
331 97671 : i < ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size();
332 : ++i)
333 : {
334 51302 : if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[i]._index < face_index)
335 46369 : face_index = ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[i]._index;
336 : }
337 :
338 : mooseAssert(face_index < std::numeric_limits<unsigned int>::max(),
339 : "face_index invalid");
340 :
341 46369 : p_info[face_index]->_closest_point = closest_point;
342 92738 : p_info[face_index]->_distance =
343 46369 : (p_info[face_index]->_distance >= 0.0 ? 1.0 : -1.0) * closest_distance;
344 : // Calculate the normal as the vector from the ridge to the point only if we're not
345 : // doing normal
346 : // smoothing. Normal smoothing will average out the normals on its own.
347 46369 : if (!_do_normal_smoothing)
348 : {
349 12114 : Point normal(closest_point - node);
350 12114 : const Real len(normal.norm());
351 12114 : if (len > 0)
352 : {
353 12058 : normal /= len;
354 : }
355 12114 : const Real dot(normal * p_info[face_index]->_normal);
356 12114 : if (dot < 0)
357 : {
358 9781 : normal *= -1;
359 : }
360 12114 : p_info[face_index]->_normal = normal;
361 : }
362 46369 : p_info[face_index]->_tangential_distance = 0.0;
363 :
364 46369 : Point closest_point_ref;
365 46369 : if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size() ==
366 : 1) // contact with a single ridge rather than a peak
367 : {
368 84440 : p_info[face_index]->_tangential_distance =
369 42220 : ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[0]._tangential_distance;
370 42220 : p_info[face_index]->_closest_point_ref =
371 42220 : ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[0]._closest_coor_ref;
372 : }
373 : else
374 : { // peak
375 : const Node * closest_node_on_face;
376 4149 : bool restricted = restrictPointToFace(p_info[face_index]->_closest_point_ref,
377 : closest_node_on_face,
378 4149 : p_info[face_index]->_side);
379 4149 : if (restricted)
380 : {
381 4149 : if (closest_node_on_face != ridgeSetDataVec[closest_ridge_set_index]._closest_node)
382 : {
383 0 : mooseError("Closest node when restricting point to face != closest node from "
384 : "RidgeSetData");
385 : }
386 : }
387 : }
388 :
389 46369 : FEBase * fe = _fes[_tid][p_info[face_index]->_side->dim()];
390 46369 : std::vector<Point> points(1);
391 46369 : points[0] = p_info[face_index]->_closest_point_ref;
392 46369 : fe->reinit(p_info[face_index]->_side, &points);
393 46369 : p_info[face_index]->_side_phi = fe->get_phi();
394 46369 : p_info[face_index]->_side_grad_phi = fe->get_dphi();
395 46369 : p_info[face_index]->_dxyzdxi = fe->get_dxyzdxi();
396 46369 : p_info[face_index]->_dxyzdeta = fe->get_dxyzdeta();
397 46369 : p_info[face_index]->_d2xyzdxideta = fe->get_d2xyzdxideta();
398 :
399 46369 : switchInfo(info, p_info[face_index]);
400 46369 : info_set = true;
401 46369 : }
402 : else
403 : { // todo:remove invalid ridge cases so they don't mess up individual face competition????
404 : }
405 72338 : }
406 :
407 157731 : if (!info_set) // contact wasn't on a ridge -- compete individual interactions
408 : {
409 111362 : unsigned int best(0), i(1);
410 : do
411 : {
412 122761 : CompeteInteractionResult CIResult = competeInteractions(p_info[best], p_info[i]);
413 122761 : if (CIResult == FIRST_WINS)
414 : {
415 20161 : i++;
416 : }
417 102600 : else if (CIResult == SECOND_WINS)
418 : {
419 10438 : best = i;
420 10438 : i++;
421 : }
422 92162 : else if (CIResult == NEITHER_WINS)
423 : {
424 92162 : best = i + 1;
425 92162 : i += 2;
426 : }
427 122761 : } while (i < p_info.size() && best < p_info.size());
428 111362 : if (best < p_info.size())
429 : {
430 : // Ensure final info is within the tangential tolerance
431 23726 : if (p_info[best]->_tangential_distance <= _tangential_tolerance)
432 : {
433 23638 : switchInfo(info, p_info[best]);
434 23638 : info_set = true;
435 : }
436 : }
437 : }
438 157731 : }
439 : }
440 :
441 1148704 : if (!info_set)
442 : {
443 : // If penetration is not detected within the saved patch, it is possible
444 : // that the secondary node has moved outside the saved patch. So, the patch
445 : // for the secondary nodes saved in _recheck_secondary_nodes has to be updated
446 : // and the penetration detection has to be re-run on the updated patch.
447 :
448 482103 : _recheck_secondary_nodes.push_back(node_id);
449 :
450 482103 : delete info;
451 482103 : info = NULL;
452 : }
453 : else
454 : {
455 666601 : smoothNormal(info, p_info, node);
456 666601 : FEBase * fe = _fes[_tid][info->_side->dim()];
457 666601 : computeSlip(*fe, *info);
458 : }
459 :
460 2001476 : for (unsigned int j = 0; j < p_info.size(); ++j)
461 : {
462 852772 : if (p_info[j])
463 : {
464 773346 : delete p_info[j];
465 773346 : p_info[j] = NULL;
466 : }
467 : }
468 1148704 : }
469 186234 : }
470 :
471 : void
472 16090 : PenetrationThread::join(const PenetrationThread & other)
473 : {
474 16090 : _recheck_secondary_nodes.insert(_recheck_secondary_nodes.end(),
475 : other._recheck_secondary_nodes.begin(),
476 : other._recheck_secondary_nodes.end());
477 16090 : }
478 :
479 : void
480 79426 : PenetrationThread::switchInfo(PenetrationInfo *& info, PenetrationInfo *& infoNew)
481 : {
482 : mooseAssert(infoNew != NULL, "infoNew object is null");
483 79426 : if (info)
484 : {
485 57292 : infoNew->_starting_elem = info->_starting_elem;
486 57292 : infoNew->_starting_side_num = info->_starting_side_num;
487 57292 : infoNew->_starting_closest_point_ref = info->_starting_closest_point_ref;
488 57292 : infoNew->_incremental_slip = info->_incremental_slip;
489 57292 : infoNew->_accumulated_slip = info->_accumulated_slip;
490 57292 : infoNew->_accumulated_slip_old = info->_accumulated_slip_old;
491 57292 : infoNew->_frictional_energy = info->_frictional_energy;
492 57292 : infoNew->_frictional_energy_old = info->_frictional_energy_old;
493 57292 : infoNew->_contact_force = info->_contact_force;
494 57292 : infoNew->_contact_force_old = info->_contact_force_old;
495 57292 : infoNew->_lagrange_multiplier = info->_lagrange_multiplier;
496 57292 : infoNew->_lagrange_multiplier_slip = info->_lagrange_multiplier_slip;
497 57292 : infoNew->_locked_this_step = info->_locked_this_step;
498 57292 : infoNew->_stick_locked_this_step = info->_stick_locked_this_step;
499 57292 : infoNew->_mech_status = info->_mech_status;
500 57292 : infoNew->_mech_status_old = info->_mech_status_old;
501 : }
502 : else
503 : {
504 22134 : infoNew->_starting_elem = infoNew->_elem;
505 22134 : infoNew->_starting_side_num = infoNew->_side_num;
506 22134 : infoNew->_starting_closest_point_ref = infoNew->_closest_point_ref;
507 : }
508 79426 : delete info;
509 79426 : info = infoNew;
510 79426 : infoNew = NULL; // Set this to NULL so that we don't delete it (now owned by _penetration_info).
511 79426 : }
512 :
513 : PenetrationThread::CompeteInteractionResult
514 122761 : PenetrationThread::competeInteractions(PenetrationInfo * pi1, PenetrationInfo * pi2)
515 : {
516 :
517 122761 : CompeteInteractionResult result = NEITHER_WINS;
518 :
519 122761 : if (pi1->_tangential_distance > _tangential_tolerance &&
520 102160 : pi2->_tangential_distance > _tangential_tolerance) // out of tol on both faces
521 92162 : result = NEITHER_WINS;
522 :
523 30599 : else if (pi1->_tangential_distance == 0.0 &&
524 19459 : pi2->_tangential_distance > 0.0) // on face 1, off face 2
525 17359 : result = FIRST_WINS;
526 :
527 13240 : else if (pi2->_tangential_distance == 0.0 &&
528 11500 : pi1->_tangential_distance > 0.0) // on face 2, off face 1
529 9400 : result = SECOND_WINS;
530 :
531 3840 : else if (pi1->_tangential_distance <= _tangential_tolerance &&
532 3110 : pi2->_tangential_distance > _tangential_tolerance) // in face 1 tol, out of face 2 tol
533 300 : result = FIRST_WINS;
534 :
535 3540 : else if (pi2->_tangential_distance <= _tangential_tolerance &&
536 3540 : pi1->_tangential_distance > _tangential_tolerance) // in face 2 tol, out of face 1 tol
537 730 : result = SECOND_WINS;
538 :
539 2810 : else if (pi1->_tangential_distance == 0.0 && pi2->_tangential_distance == 0.0) // on both faces
540 2100 : result = competeInteractionsBothOnFace(pi1, pi2);
541 :
542 710 : else if (pi1->_tangential_distance <= _tangential_tolerance &&
543 710 : pi2->_tangential_distance <= _tangential_tolerance) // off but within tol of both faces
544 : {
545 710 : CommonEdgeResult cer = interactionsOffCommonEdge(pi1, pi2);
546 710 : if (cer == COMMON_EDGE || cer == COMMON_NODE) // ridge case.
547 : {
548 : // We already checked for ridges, and it got rejected, so neither face must be valid
549 0 : result = NEITHER_WINS;
550 : // mooseError("Erroneously encountered ridge case");
551 : }
552 710 : else if (cer == EDGE_AND_COMMON_NODE) // off side of face, off corner of another face. Favor
553 : // the off-side face
554 : {
555 454 : if (pi1->_off_edge_nodes.size() == pi2->_off_edge_nodes.size())
556 0 : mooseError("Invalid off_edge_nodes counts");
557 :
558 454 : else if (pi1->_off_edge_nodes.size() == 2)
559 146 : result = FIRST_WINS;
560 :
561 308 : else if (pi2->_off_edge_nodes.size() == 2)
562 308 : result = SECOND_WINS;
563 :
564 : else
565 0 : mooseError("Invalid off_edge_nodes counts");
566 : }
567 : else // The node projects to both faces within tangential tolerance.
568 256 : result = competeInteractionsBothOnFace(pi1, pi2);
569 : }
570 :
571 122761 : return result;
572 : }
573 :
574 : PenetrationThread::CompeteInteractionResult
575 2356 : PenetrationThread::competeInteractionsBothOnFace(PenetrationInfo * pi1, PenetrationInfo * pi2)
576 : {
577 2356 : CompeteInteractionResult result = NEITHER_WINS;
578 :
579 2356 : if (pi1->_distance >= 0.0 && pi2->_distance < 0.0)
580 0 : result = FIRST_WINS; // favor face with positive distance (penetrated) -- first in this case
581 :
582 2356 : else if (pi2->_distance >= 0.0 && pi1->_distance < 0.0)
583 0 : result = SECOND_WINS; // favor face with positive distance (penetrated) -- second in this case
584 :
585 : // TODO: This logic below could cause an abrupt jump from one face to the other with small mesh
586 : // movement. If there is some way to smooth the transition, we should do it.
587 2356 : else if (MooseUtils::relativeFuzzyLessThan(std::abs(pi1->_distance), std::abs(pi2->_distance)))
588 0 : result = FIRST_WINS; // otherwise, favor the closer face -- first in this case
589 :
590 2356 : else if (MooseUtils::relativeFuzzyLessThan(std::abs(pi2->_distance), std::abs(pi1->_distance)))
591 0 : result = SECOND_WINS; // otherwise, favor the closer face -- second in this case
592 :
593 : else // Equal within tolerance. Favor the one with a smaller element id (for repeatibility)
594 : {
595 2356 : if (pi1->_elem->id() < pi2->_elem->id())
596 2356 : result = FIRST_WINS;
597 :
598 : else
599 0 : result = SECOND_WINS;
600 : }
601 :
602 2356 : return result;
603 : }
604 :
605 : PenetrationThread::CommonEdgeResult
606 710 : PenetrationThread::interactionsOffCommonEdge(PenetrationInfo * pi1, PenetrationInfo * pi2)
607 : {
608 710 : CommonEdgeResult common_edge(NO_COMMON);
609 710 : const std::vector<const Node *> & off_edge_nodes1 = pi1->_off_edge_nodes;
610 710 : const std::vector<const Node *> & off_edge_nodes2 = pi2->_off_edge_nodes;
611 710 : const unsigned dim1 = pi1->_side->dim();
612 :
613 710 : if (dim1 == 1)
614 : {
615 : mooseAssert(pi2->_side->dim() == 1, "Incompatible dimensions.");
616 : mooseAssert(off_edge_nodes1.size() < 2 && off_edge_nodes2.size() < 2,
617 : "off_edge_nodes size should be <2 for 2D contact");
618 0 : if (off_edge_nodes1.size() == 1 && off_edge_nodes2.size() == 1 &&
619 0 : off_edge_nodes1[0] == off_edge_nodes2[0])
620 0 : common_edge = COMMON_EDGE;
621 : }
622 : else
623 : {
624 : mooseAssert(dim1 == 2 && pi2->_side->dim() == 2, "Incompatible dimensions.");
625 : mooseAssert(off_edge_nodes1.size() < 3 && off_edge_nodes2.size() < 3,
626 : "off_edge_nodes size should be <3 for 3D contact");
627 710 : if (off_edge_nodes1.size() == 1)
628 : {
629 308 : if (off_edge_nodes2.size() == 1)
630 : {
631 0 : if (off_edge_nodes1[0] == off_edge_nodes2[0])
632 0 : common_edge = COMMON_NODE;
633 : }
634 308 : else if (off_edge_nodes2.size() == 2)
635 : {
636 308 : if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[0] == off_edge_nodes2[1])
637 308 : common_edge = EDGE_AND_COMMON_NODE;
638 : }
639 : }
640 402 : else if (off_edge_nodes1.size() == 2)
641 : {
642 402 : if (off_edge_nodes2.size() == 1)
643 : {
644 146 : if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[1] == off_edge_nodes2[0])
645 146 : common_edge = EDGE_AND_COMMON_NODE;
646 : }
647 256 : else if (off_edge_nodes2.size() == 2)
648 : {
649 256 : if ((off_edge_nodes1[0] == off_edge_nodes2[0] &&
650 512 : off_edge_nodes1[1] == off_edge_nodes2[1]) ||
651 256 : (off_edge_nodes1[1] == off_edge_nodes2[0] && off_edge_nodes1[0] == off_edge_nodes2[1]))
652 0 : common_edge = COMMON_EDGE;
653 : }
654 : }
655 : }
656 710 : return common_edge;
657 : }
658 :
659 : bool
660 246080 : PenetrationThread::findRidgeContactPoint(Point & contact_point,
661 : Real & tangential_distance,
662 : const Node *& closest_node,
663 : unsigned int & index,
664 : Point & contact_point_ref,
665 : std::vector<PenetrationInfo *> & p_info,
666 : const unsigned int index1,
667 : const unsigned int index2)
668 : {
669 246080 : tangential_distance = 0.0;
670 246080 : closest_node = NULL;
671 246080 : PenetrationInfo * pi1 = p_info[index1];
672 246080 : PenetrationInfo * pi2 = p_info[index2];
673 246080 : const unsigned sidedim(pi1->_side->dim());
674 : mooseAssert(sidedim == pi2->_side->dim(), "Incompatible dimensionalities");
675 :
676 : // Nodes on faces for the two interactions
677 246080 : std::vector<const Node *> side1_nodes;
678 246080 : getSideCornerNodes(pi1->_side, side1_nodes);
679 246080 : std::vector<const Node *> side2_nodes;
680 246080 : getSideCornerNodes(pi2->_side, side2_nodes);
681 :
682 246080 : std::sort(side1_nodes.begin(), side1_nodes.end());
683 246080 : std::sort(side2_nodes.begin(), side2_nodes.end());
684 :
685 : // Find nodes shared by the two faces
686 246080 : std::vector<const Node *> common_nodes;
687 246080 : std::set_intersection(side1_nodes.begin(),
688 : side1_nodes.end(),
689 : side2_nodes.begin(),
690 : side2_nodes.end(),
691 : std::inserter(common_nodes, common_nodes.end()));
692 :
693 246080 : if (common_nodes.size() != sidedim)
694 36765 : return false;
695 :
696 : bool found_point1, found_point2;
697 209315 : Point closest_coor_ref1(pi1->_closest_point_ref);
698 : const Node * closest_node1;
699 209315 : found_point1 = restrictPointToSpecifiedEdgeOfFace(
700 : closest_coor_ref1, closest_node1, pi1->_side, common_nodes);
701 :
702 209315 : Point closest_coor_ref2(pi2->_closest_point_ref);
703 : const Node * closest_node2;
704 209315 : found_point2 = restrictPointToSpecifiedEdgeOfFace(
705 : closest_coor_ref2, closest_node2, pi2->_side, common_nodes);
706 :
707 209315 : if (!found_point1 || !found_point2)
708 118720 : return false;
709 :
710 : // if (sidedim == 2)
711 : // {
712 : // TODO:
713 : // We have the parametric coordinates of the closest intersection point for both faces.
714 : // We need to find a point somewhere in the middle of them so there's not an abrupt jump.
715 : // Find that point by taking dot products of vector from contact point to secondary node point
716 : // with face normal vectors to see which face we're closer to.
717 : // }
718 :
719 90595 : FEBase * fe = NULL;
720 90595 : std::vector<Point> points(1);
721 :
722 : // We have to pick one of the two faces to own the contact point. It doesn't really matter
723 : // which one we pick. For repeatibility, pick the face with the lowest index.
724 90595 : if (index1 < index2)
725 : {
726 90595 : fe = _fes[_tid][pi1->_side->dim()];
727 90595 : contact_point_ref = closest_coor_ref1;
728 90595 : points[0] = closest_coor_ref1;
729 90595 : fe->reinit(pi1->_side, &points);
730 90595 : index = index1;
731 : }
732 : else
733 : {
734 0 : fe = _fes[_tid][pi2->_side->dim()];
735 0 : contact_point_ref = closest_coor_ref2;
736 0 : points[0] = closest_coor_ref2;
737 0 : fe->reinit(pi2->_side, &points);
738 0 : index = index2;
739 : }
740 :
741 90595 : contact_point = fe->get_xyz()[0];
742 :
743 90595 : if (sidedim == 2)
744 : {
745 85939 : if (closest_node1) // point is off the ridge between the two elements
746 : {
747 : mooseAssert((closest_node1 == closest_node2 || closest_node2 == NULL),
748 : "If off edge of ridge, closest node must be the same on both elements");
749 52265 : closest_node = closest_node1;
750 :
751 52265 : RealGradient off_face = *closest_node1 - contact_point;
752 52265 : tangential_distance = off_face.norm();
753 : }
754 : }
755 :
756 90595 : return true;
757 246080 : }
758 :
759 : void
760 492160 : PenetrationThread::getSideCornerNodes(const Elem * side, std::vector<const Node *> & corner_nodes)
761 : {
762 492160 : const ElemType t(side->type());
763 492160 : corner_nodes.clear();
764 :
765 492160 : corner_nodes.push_back(side->node_ptr(0));
766 492160 : corner_nodes.push_back(side->node_ptr(1));
767 492160 : switch (t)
768 : {
769 33166 : case EDGE2:
770 : case EDGE3:
771 : case EDGE4:
772 : {
773 33166 : break;
774 : }
775 :
776 11022 : case TRI3:
777 : case TRI6:
778 : case TRI7:
779 : {
780 11022 : corner_nodes.push_back(side->node_ptr(2));
781 11022 : break;
782 : }
783 :
784 447972 : case QUAD4:
785 : case QUAD8:
786 : case QUAD9:
787 : {
788 447972 : corner_nodes.push_back(side->node_ptr(2));
789 447972 : corner_nodes.push_back(side->node_ptr(3));
790 447972 : break;
791 : }
792 :
793 0 : default:
794 : {
795 0 : mooseError("Unsupported face type: ", t);
796 : break;
797 : }
798 : }
799 492160 : }
800 :
801 : bool
802 418630 : PenetrationThread::restrictPointToSpecifiedEdgeOfFace(Point & p,
803 : const Node *& closest_node,
804 : const Elem * side,
805 : const std::vector<const Node *> & edge_nodes)
806 : {
807 418630 : const ElemType t = side->type();
808 418630 : Real & xi = p(0);
809 418630 : Real & eta = p(1);
810 418630 : closest_node = NULL;
811 :
812 418630 : std::vector<unsigned int> local_node_indices;
813 1222724 : for (const auto & edge_node : edge_nodes)
814 : {
815 804094 : unsigned int local_index = side->get_node_index(edge_node);
816 804094 : if (local_index == libMesh::invalid_uint)
817 0 : mooseError("Side does not contain node");
818 804094 : local_node_indices.push_back(local_index);
819 : }
820 : mooseAssert(local_node_indices.size() == side->dim(),
821 : "Number of edge nodes must match side dimensionality");
822 418630 : std::sort(local_node_indices.begin(), local_node_indices.end());
823 :
824 418630 : bool off_of_this_edge = false;
825 :
826 418630 : switch (t)
827 : {
828 33166 : case EDGE2:
829 : case EDGE3:
830 : case EDGE4:
831 : {
832 33166 : if (local_node_indices[0] == 0)
833 : {
834 16583 : if (xi <= -1.0)
835 : {
836 12488 : xi = -1.0;
837 12488 : off_of_this_edge = true;
838 12488 : closest_node = side->node_ptr(0);
839 : }
840 : }
841 16583 : else if (local_node_indices[0] == 1)
842 : {
843 16583 : if (xi >= 1.0)
844 : {
845 8740 : xi = 1.0;
846 8740 : off_of_this_edge = true;
847 8740 : closest_node = side->node_ptr(1);
848 : }
849 : }
850 : else
851 : {
852 0 : mooseError("Invalid local node indices");
853 : }
854 33166 : break;
855 : }
856 :
857 4580 : case TRI3:
858 : case TRI6:
859 : case TRI7:
860 : {
861 4580 : if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
862 : {
863 1194 : if (eta <= 0.0)
864 : {
865 600 : eta = 0.0;
866 600 : off_of_this_edge = true;
867 600 : if (xi < 0.0)
868 99 : closest_node = side->node_ptr(0);
869 501 : else if (xi > 1.0)
870 212 : closest_node = side->node_ptr(1);
871 : }
872 : }
873 3386 : else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
874 : {
875 1726 : if ((xi + eta) > 1.0)
876 : {
877 817 : Real delta = (xi + eta - 1.0) / 2.0;
878 817 : xi -= delta;
879 817 : eta -= delta;
880 817 : off_of_this_edge = true;
881 817 : if (xi > 1.0)
882 190 : closest_node = side->node_ptr(1);
883 627 : else if (xi < 0.0)
884 221 : closest_node = side->node_ptr(2);
885 : }
886 : }
887 1660 : else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 2))
888 : {
889 1660 : if (xi <= 0.0)
890 : {
891 835 : xi = 0.0;
892 835 : off_of_this_edge = true;
893 835 : if (eta > 1.0)
894 257 : closest_node = side->node_ptr(2);
895 578 : else if (eta < 0.0)
896 209 : closest_node = side->node_ptr(0);
897 : }
898 : }
899 : else
900 : {
901 0 : mooseError("Invalid local node indices");
902 : }
903 :
904 4580 : break;
905 : }
906 :
907 380884 : case QUAD4:
908 : case QUAD8:
909 : case QUAD9:
910 : {
911 380884 : if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
912 : {
913 165240 : if (eta <= -1.0)
914 : {
915 101256 : eta = -1.0;
916 101256 : off_of_this_edge = true;
917 101256 : if (xi < -1.0)
918 31737 : closest_node = side->node_ptr(0);
919 69519 : else if (xi > 1.0)
920 31825 : closest_node = side->node_ptr(1);
921 : }
922 : }
923 215644 : else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
924 : {
925 80070 : if (xi >= 1.0)
926 : {
927 66030 : xi = 1.0;
928 66030 : off_of_this_edge = true;
929 66030 : if (eta < -1.0)
930 15036 : closest_node = side->node_ptr(1);
931 50994 : else if (eta > 1.0)
932 12350 : closest_node = side->node_ptr(2);
933 : }
934 : }
935 135574 : else if ((local_node_indices[0] == 2) && (local_node_indices[1] == 3))
936 : {
937 111275 : if (eta >= 1.0)
938 : {
939 82596 : eta = 1.0;
940 82596 : off_of_this_edge = true;
941 82596 : if (xi < -1.0)
942 35910 : closest_node = side->node_ptr(3);
943 46686 : else if (xi > 1.0)
944 41909 : closest_node = side->node_ptr(2);
945 : }
946 : }
947 24299 : else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 3))
948 : {
949 24299 : if (xi <= -1.0)
950 : {
951 16550 : xi = -1.0;
952 16550 : off_of_this_edge = true;
953 16550 : if (eta < -1.0)
954 9892 : closest_node = side->node_ptr(0);
955 6658 : else if (eta > 1.0)
956 2637 : closest_node = side->node_ptr(3);
957 : }
958 : }
959 : else
960 : {
961 0 : mooseError("Invalid local node indices");
962 : }
963 380884 : break;
964 : }
965 :
966 0 : default:
967 : {
968 0 : mooseError("Unsupported face type: ", t);
969 : break;
970 : }
971 : }
972 418630 : return off_of_this_edge;
973 418630 : }
974 :
975 : bool
976 4149 : PenetrationThread::restrictPointToFace(Point & p, const Node *& closest_node, const Elem * side)
977 : {
978 4149 : const ElemType t(side->type());
979 4149 : Real & xi = p(0);
980 4149 : Real & eta = p(1);
981 4149 : closest_node = NULL;
982 :
983 4149 : bool off_of_this_face(false);
984 :
985 4149 : switch (t)
986 : {
987 0 : case EDGE2:
988 : case EDGE3:
989 : case EDGE4:
990 : {
991 0 : if (xi < -1.0)
992 : {
993 0 : xi = -1.0;
994 0 : off_of_this_face = true;
995 0 : closest_node = side->node_ptr(0);
996 : }
997 0 : else if (xi > 1.0)
998 : {
999 0 : xi = 1.0;
1000 0 : off_of_this_face = true;
1001 0 : closest_node = side->node_ptr(1);
1002 : }
1003 0 : break;
1004 : }
1005 :
1006 0 : case TRI3:
1007 : case TRI6:
1008 : case TRI7:
1009 : {
1010 0 : if (eta < 0.0)
1011 : {
1012 0 : eta = 0.0;
1013 0 : off_of_this_face = true;
1014 0 : if (xi < 0.5)
1015 : {
1016 0 : closest_node = side->node_ptr(0);
1017 0 : if (xi < 0.0)
1018 0 : xi = 0.0;
1019 : }
1020 : else
1021 : {
1022 0 : closest_node = side->node_ptr(1);
1023 0 : if (xi > 1.0)
1024 0 : xi = 1.0;
1025 : }
1026 : }
1027 0 : else if ((xi + eta) > 1.0)
1028 : {
1029 0 : Real delta = (xi + eta - 1.0) / 2.0;
1030 0 : xi -= delta;
1031 0 : eta -= delta;
1032 0 : off_of_this_face = true;
1033 0 : if (xi > 0.5)
1034 : {
1035 0 : closest_node = side->node_ptr(1);
1036 0 : if (xi > 1.0)
1037 : {
1038 0 : xi = 1.0;
1039 0 : eta = 0.0;
1040 : }
1041 : }
1042 : else
1043 : {
1044 0 : closest_node = side->node_ptr(2);
1045 0 : if (xi < 0.0)
1046 : {
1047 0 : xi = 0.0;
1048 0 : eta = 1.0;
1049 : }
1050 : }
1051 : }
1052 0 : else if (xi < 0.0)
1053 : {
1054 0 : xi = 0.0;
1055 0 : off_of_this_face = true;
1056 0 : if (eta > 0.5)
1057 : {
1058 0 : closest_node = side->node_ptr(2);
1059 0 : if (eta > 1.0)
1060 0 : eta = 1.0;
1061 : }
1062 : else
1063 : {
1064 0 : closest_node = side->node_ptr(0);
1065 0 : if (eta < 0.0)
1066 0 : eta = 0.0;
1067 : }
1068 : }
1069 0 : break;
1070 : }
1071 :
1072 4149 : case QUAD4:
1073 : case QUAD8:
1074 : case QUAD9:
1075 : {
1076 4149 : if (eta < -1.0)
1077 : {
1078 108 : eta = -1.0;
1079 108 : off_of_this_face = true;
1080 108 : if (xi < 0.0)
1081 : {
1082 14 : closest_node = side->node_ptr(0);
1083 14 : if (xi < -1.0)
1084 14 : xi = -1.0;
1085 : }
1086 : else
1087 : {
1088 94 : closest_node = side->node_ptr(1);
1089 94 : if (xi > 1.0)
1090 94 : xi = 1.0;
1091 : }
1092 : }
1093 4041 : else if (xi > 1.0)
1094 : {
1095 4041 : xi = 1.0;
1096 4041 : off_of_this_face = true;
1097 4041 : if (eta < 0.0)
1098 : {
1099 0 : closest_node = side->node_ptr(1);
1100 0 : if (eta < -1.0)
1101 0 : eta = -1.0;
1102 : }
1103 : else
1104 : {
1105 4041 : closest_node = side->node_ptr(2);
1106 4041 : if (eta > 1.0)
1107 1044 : eta = 1.0;
1108 : }
1109 : }
1110 0 : else if (eta > 1.0)
1111 : {
1112 0 : eta = 1.0;
1113 0 : off_of_this_face = true;
1114 0 : if (xi < 0.0)
1115 : {
1116 0 : closest_node = side->node_ptr(3);
1117 0 : if (xi < -1.0)
1118 0 : xi = -1.0;
1119 : }
1120 : else
1121 : {
1122 0 : closest_node = side->node_ptr(2);
1123 0 : if (xi > 1.0)
1124 0 : xi = 1.0;
1125 : }
1126 : }
1127 0 : else if (xi < -1.0)
1128 : {
1129 0 : xi = -1.0;
1130 0 : off_of_this_face = true;
1131 0 : if (eta < 0.0)
1132 : {
1133 0 : closest_node = side->node_ptr(0);
1134 0 : if (eta < -1.0)
1135 0 : eta = -1.0;
1136 : }
1137 : else
1138 : {
1139 0 : closest_node = side->node_ptr(3);
1140 0 : if (eta > 1.0)
1141 0 : eta = 1.0;
1142 : }
1143 : }
1144 4149 : break;
1145 : }
1146 :
1147 0 : default:
1148 : {
1149 0 : mooseError("Unsupported face type: ", t);
1150 : break;
1151 : }
1152 : }
1153 4149 : return off_of_this_face;
1154 : }
1155 :
1156 : bool
1157 757643 : PenetrationThread::isFaceReasonableCandidate(const Elem * primary_elem,
1158 : const Elem * side,
1159 : FEBase * fe,
1160 : const Point * secondary_point,
1161 : const Real tangential_tolerance)
1162 : {
1163 757643 : unsigned int dim = primary_elem->dim();
1164 :
1165 757643 : const std::vector<Point> & phys_point = fe->get_xyz();
1166 :
1167 757643 : const std::vector<RealGradient> & dxyz_dxi = fe->get_dxyzdxi();
1168 757643 : const std::vector<RealGradient> & dxyz_deta = fe->get_dxyzdeta();
1169 :
1170 757643 : Point ref_point;
1171 :
1172 757643 : std::vector<Point> points(1); // Default constructor gives us a point at 0,0,0
1173 :
1174 757643 : fe->reinit(side, &points);
1175 :
1176 757643 : RealGradient d = *secondary_point - phys_point[0];
1177 :
1178 757643 : const Real twosqrt2 = 2.8284; // way more precision than we actually need here
1179 757643 : Real max_face_length = side->hmax() + twosqrt2 * tangential_tolerance;
1180 :
1181 757643 : RealVectorValue normal;
1182 757643 : if (dim - 1 == 2)
1183 : {
1184 682055 : normal = dxyz_dxi[0].cross(dxyz_deta[0]);
1185 : }
1186 75588 : else if (dim - 1 == 1)
1187 : {
1188 75569 : const Node * const * elem_nodes = primary_elem->get_nodes();
1189 75569 : const Point in_plane_vector1 = *elem_nodes[1] - *elem_nodes[0];
1190 75569 : const Point in_plane_vector2 = *elem_nodes[2] - *elem_nodes[0];
1191 :
1192 75569 : Point out_of_plane_normal = in_plane_vector1.cross(in_plane_vector2);
1193 75569 : out_of_plane_normal /= out_of_plane_normal.norm();
1194 :
1195 75569 : normal = dxyz_dxi[0].cross(out_of_plane_normal);
1196 : }
1197 : else
1198 : {
1199 19 : return true;
1200 : }
1201 757624 : normal /= normal.norm();
1202 :
1203 757624 : const Real dot(d * normal);
1204 :
1205 757624 : const RealGradient normcomp = dot * normal;
1206 757624 : const RealGradient tangcomp = d - normcomp;
1207 :
1208 757624 : const Real tangdist = tangcomp.norm();
1209 :
1210 : // Increase the size of the zone that we consider if the vector from the face
1211 : // to the node has a larger normal component
1212 757624 : const Real faceExpansionFactor = 2.0 * (1.0 + normcomp.norm() / d.norm());
1213 :
1214 757624 : bool isReasonableCandidate = true;
1215 757624 : if (tangdist > faceExpansionFactor * max_face_length)
1216 : {
1217 8068 : isReasonableCandidate = false;
1218 : }
1219 757624 : return isReasonableCandidate;
1220 757643 : }
1221 :
1222 : void
1223 666601 : PenetrationThread::computeSlip(FEBase & fe, PenetrationInfo & info)
1224 : {
1225 : // Slip is current projected position of secondary node minus
1226 : // original projected position of secondary node
1227 666601 : std::vector<Point> points(1);
1228 666601 : points[0] = info._starting_closest_point_ref;
1229 666601 : const auto & side = _elem_side_builder(*info._starting_elem, info._starting_side_num);
1230 666601 : fe.reinit(&side, &points);
1231 666601 : const std::vector<Point> & starting_point = fe.get_xyz();
1232 666601 : info._incremental_slip = info._closest_point - starting_point[0];
1233 666601 : if (info.isCaptured())
1234 : {
1235 50848 : info._frictional_energy =
1236 50848 : info._frictional_energy_old + info._contact_force * info._incremental_slip;
1237 50848 : info._accumulated_slip = info._accumulated_slip_old + info._incremental_slip.norm();
1238 : }
1239 666601 : }
1240 :
1241 : void
1242 666601 : PenetrationThread::smoothNormal(PenetrationInfo * info,
1243 : std::vector<PenetrationInfo *> & p_info,
1244 : const Node & node)
1245 : {
1246 666601 : if (_do_normal_smoothing)
1247 : {
1248 324700 : if (_normal_smoothing_method == PenetrationLocator::NSM_EDGE_BASED)
1249 : {
1250 : // If we are within the smoothing distance of any edges or corners, find the
1251 : // corner nodes for those edges/corners, and weights from distance to edge/corner
1252 259380 : std::vector<Real> edge_face_weights;
1253 259380 : std::vector<PenetrationInfo *> edge_face_info;
1254 :
1255 259380 : getSmoothingFacesAndWeights(info, edge_face_info, edge_face_weights, p_info, node);
1256 :
1257 : mooseAssert(edge_face_info.size() == edge_face_weights.size(),
1258 : "edge_face_info.size() != edge_face_weights.size()");
1259 :
1260 259380 : if (edge_face_info.size() > 0)
1261 : {
1262 : // Smooth the normal using the weighting functions for all participating faces.
1263 116009 : RealVectorValue new_normal;
1264 116009 : Real this_face_weight = 1.0;
1265 :
1266 261062 : for (unsigned int efwi = 0; efwi < edge_face_weights.size(); ++efwi)
1267 : {
1268 145053 : PenetrationInfo * npi = edge_face_info[efwi];
1269 145053 : if (npi)
1270 145053 : new_normal += npi->_normal * edge_face_weights[efwi];
1271 :
1272 145053 : this_face_weight -= edge_face_weights[efwi];
1273 : }
1274 : mooseAssert(this_face_weight >= (0.25 - 1e-8),
1275 : "Sum of weights of other faces shouldn't exceed 0.75");
1276 116009 : new_normal += info->_normal * this_face_weight;
1277 :
1278 116009 : const Real len = new_normal.norm();
1279 116009 : if (len > 0)
1280 116009 : new_normal /= len;
1281 :
1282 116009 : info->_normal = new_normal;
1283 : }
1284 259380 : }
1285 65320 : else if (_normal_smoothing_method == PenetrationLocator::NSM_NODAL_NORMAL_BASED)
1286 : {
1287 : // params.addParam<VariableName>("var_name","description");
1288 : // getParam<VariableName>("var_name")
1289 65320 : info->_normal(0) = _nodal_normal_x->getValue(info->_side, info->_side_phi);
1290 65320 : info->_normal(1) = _nodal_normal_y->getValue(info->_side, info->_side_phi);
1291 65320 : info->_normal(2) = _nodal_normal_z->getValue(info->_side, info->_side_phi);
1292 65320 : const Real len(info->_normal.norm());
1293 65320 : if (len > 0)
1294 64572 : info->_normal /= len;
1295 : }
1296 : }
1297 666601 : }
1298 :
1299 : void
1300 259380 : PenetrationThread::getSmoothingFacesAndWeights(PenetrationInfo * info,
1301 : std::vector<PenetrationInfo *> & edge_face_info,
1302 : std::vector<Real> & edge_face_weights,
1303 : std::vector<PenetrationInfo *> & p_info,
1304 : const Node & secondary_node)
1305 : {
1306 259380 : const Elem * side = info->_side;
1307 259380 : const Point & p = info->_closest_point_ref;
1308 259380 : std::set<dof_id_type> elems_to_exclude;
1309 259380 : elems_to_exclude.insert(info->_elem->id());
1310 :
1311 259380 : std::vector<std::vector<const Node *>> edge_nodes;
1312 :
1313 : // Get the pairs of nodes along every edge that we are close enough to smooth with
1314 259380 : getSmoothingEdgeNodesAndWeights(p, side, edge_nodes, edge_face_weights);
1315 259380 : std::vector<Elem *> edge_neighbor_elems;
1316 259380 : edge_face_info.resize(edge_nodes.size(), NULL);
1317 :
1318 259380 : std::vector<unsigned int> edges_without_neighbors;
1319 :
1320 484890 : for (unsigned int i = 0; i < edge_nodes.size(); ++i)
1321 : {
1322 : // Sort all sets of edge nodes (needed for comparing edges)
1323 225510 : std::sort(edge_nodes[i].begin(), edge_nodes[i].end());
1324 :
1325 225510 : std::vector<PenetrationInfo *> face_info_comm_edge;
1326 225510 : getInfoForFacesWithCommonNodes(
1327 225510 : &secondary_node, elems_to_exclude, edge_nodes[i], face_info_comm_edge, p_info);
1328 :
1329 225510 : if (face_info_comm_edge.size() == 0)
1330 94979 : edges_without_neighbors.push_back(i);
1331 130531 : else if (face_info_comm_edge.size() > 1)
1332 0 : mooseError("Only one neighbor allowed per edge");
1333 : else
1334 130531 : edge_face_info[i] = face_info_comm_edge[0];
1335 225510 : }
1336 :
1337 : // Remove edges without neighbors from the vector, starting from end
1338 259380 : std::vector<unsigned int>::reverse_iterator rit;
1339 354359 : for (rit = edges_without_neighbors.rbegin(); rit != edges_without_neighbors.rend(); ++rit)
1340 : {
1341 94979 : unsigned int index = *rit;
1342 94979 : edge_nodes.erase(edge_nodes.begin() + index);
1343 94979 : edge_face_weights.erase(edge_face_weights.begin() + index);
1344 94979 : edge_face_info.erase(edge_face_info.begin() + index);
1345 : }
1346 :
1347 : // Handle corner case
1348 259380 : if (edge_nodes.size() > 1)
1349 : {
1350 14522 : if (edge_nodes.size() != 2)
1351 0 : mooseError("Invalid number of smoothing edges");
1352 :
1353 : // find common node
1354 14522 : std::vector<const Node *> common_nodes;
1355 58088 : std::set_intersection(edge_nodes[0].begin(),
1356 14522 : edge_nodes[0].end(),
1357 14522 : edge_nodes[1].begin(),
1358 14522 : edge_nodes[1].end(),
1359 : std::inserter(common_nodes, common_nodes.end()));
1360 :
1361 14522 : if (common_nodes.size() != 1)
1362 0 : mooseError("Invalid number of common nodes");
1363 :
1364 43566 : for (const auto & pinfo : edge_face_info)
1365 29044 : elems_to_exclude.insert(pinfo->_elem->id());
1366 :
1367 14522 : std::vector<PenetrationInfo *> face_info_comm_edge;
1368 14522 : getInfoForFacesWithCommonNodes(
1369 : &secondary_node, elems_to_exclude, common_nodes, face_info_comm_edge, p_info);
1370 :
1371 14522 : unsigned int num_corner_neighbors = face_info_comm_edge.size();
1372 :
1373 14522 : if (num_corner_neighbors > 0)
1374 : {
1375 14522 : Real fw0 = edge_face_weights[0];
1376 14522 : Real fw1 = edge_face_weights[1];
1377 :
1378 : // Corner weight is product of edge weights. Spread out over multiple neighbors.
1379 14522 : Real fw_corner = (fw0 * fw1) / static_cast<Real>(num_corner_neighbors);
1380 :
1381 : // Adjust original edge weights
1382 14522 : edge_face_weights[0] *= (1.0 - fw1);
1383 14522 : edge_face_weights[1] *= (1.0 - fw0);
1384 :
1385 29044 : for (unsigned int i = 0; i < num_corner_neighbors; ++i)
1386 : {
1387 14522 : edge_face_weights.push_back(fw_corner);
1388 14522 : edge_face_info.push_back(face_info_comm_edge[i]);
1389 : }
1390 : }
1391 14522 : }
1392 259380 : }
1393 :
1394 : void
1395 259380 : PenetrationThread::getSmoothingEdgeNodesAndWeights(
1396 : const Point & p,
1397 : const Elem * side,
1398 : std::vector<std::vector<const Node *>> & edge_nodes,
1399 : std::vector<Real> & edge_face_weights)
1400 : {
1401 259380 : const ElemType t(side->type());
1402 259380 : const Real & xi = p(0);
1403 259380 : const Real & eta = p(1);
1404 :
1405 259380 : Real smooth_limit = 1.0 - _normal_smoothing_distance;
1406 :
1407 259380 : switch (t)
1408 : {
1409 9936 : case EDGE2:
1410 : case EDGE3:
1411 : case EDGE4:
1412 : {
1413 9936 : if (xi < -smooth_limit)
1414 : {
1415 1988 : std::vector<const Node *> en;
1416 1988 : en.push_back(side->node_ptr(0));
1417 1988 : edge_nodes.push_back(en);
1418 1988 : Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
1419 1988 : if (fw > 0.5)
1420 132 : fw = 0.5;
1421 1988 : edge_face_weights.push_back(fw);
1422 1988 : }
1423 7948 : else if (xi > smooth_limit)
1424 : {
1425 1134 : std::vector<const Node *> en;
1426 1134 : en.push_back(side->node_ptr(1));
1427 1134 : edge_nodes.push_back(en);
1428 1134 : Real fw = 0.5 - (1.0 - xi) / (2.0 * _normal_smoothing_distance);
1429 1134 : if (fw > 0.5)
1430 180 : fw = 0.5;
1431 1134 : edge_face_weights.push_back(fw);
1432 1134 : }
1433 9936 : break;
1434 : }
1435 :
1436 0 : case TRI3:
1437 : case TRI6:
1438 : case TRI7:
1439 : {
1440 0 : if (eta < -smooth_limit)
1441 : {
1442 0 : std::vector<const Node *> en;
1443 0 : en.push_back(side->node_ptr(0));
1444 0 : en.push_back(side->node_ptr(1));
1445 0 : edge_nodes.push_back(en);
1446 0 : Real fw = 0.5 - (1.0 + eta) / (2.0 * _normal_smoothing_distance);
1447 0 : if (fw > 0.5)
1448 0 : fw = 0.5;
1449 0 : edge_face_weights.push_back(fw);
1450 0 : }
1451 0 : if ((xi + eta) > smooth_limit)
1452 : {
1453 0 : std::vector<const Node *> en;
1454 0 : en.push_back(side->node_ptr(1));
1455 0 : en.push_back(side->node_ptr(2));
1456 0 : edge_nodes.push_back(en);
1457 0 : Real fw = 0.5 - (1.0 - xi - eta) / (2.0 * _normal_smoothing_distance);
1458 0 : if (fw > 0.5)
1459 0 : fw = 0.5;
1460 0 : edge_face_weights.push_back(fw);
1461 0 : }
1462 0 : if (xi < -smooth_limit)
1463 : {
1464 0 : std::vector<const Node *> en;
1465 0 : en.push_back(side->node_ptr(2));
1466 0 : en.push_back(side->node_ptr(0));
1467 0 : edge_nodes.push_back(en);
1468 0 : Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
1469 0 : if (fw > 0.5)
1470 0 : fw = 0.5;
1471 0 : edge_face_weights.push_back(fw);
1472 0 : }
1473 0 : break;
1474 : }
1475 :
1476 249444 : case QUAD4:
1477 : case QUAD8:
1478 : case QUAD9:
1479 : {
1480 249444 : if (eta < -smooth_limit)
1481 : {
1482 41269 : std::vector<const Node *> en;
1483 41269 : en.push_back(side->node_ptr(0));
1484 41269 : en.push_back(side->node_ptr(1));
1485 41269 : edge_nodes.push_back(en);
1486 41269 : Real fw = 0.5 - (1.0 + eta) / (2.0 * _normal_smoothing_distance);
1487 41269 : if (fw > 0.5)
1488 14030 : fw = 0.5;
1489 41269 : edge_face_weights.push_back(fw);
1490 41269 : }
1491 249444 : if (xi > smooth_limit)
1492 : {
1493 77736 : std::vector<const Node *> en;
1494 77736 : en.push_back(side->node_ptr(1));
1495 77736 : en.push_back(side->node_ptr(2));
1496 77736 : edge_nodes.push_back(en);
1497 77736 : Real fw = 0.5 - (1.0 - xi) / (2.0 * _normal_smoothing_distance);
1498 77736 : if (fw > 0.5)
1499 17679 : fw = 0.5;
1500 77736 : edge_face_weights.push_back(fw);
1501 77736 : }
1502 249444 : if (eta > smooth_limit)
1503 : {
1504 79041 : std::vector<const Node *> en;
1505 79041 : en.push_back(side->node_ptr(2));
1506 79041 : en.push_back(side->node_ptr(3));
1507 79041 : edge_nodes.push_back(en);
1508 79041 : Real fw = 0.5 - (1.0 - eta) / (2.0 * _normal_smoothing_distance);
1509 79041 : if (fw > 0.5)
1510 22829 : fw = 0.5;
1511 79041 : edge_face_weights.push_back(fw);
1512 79041 : }
1513 249444 : if (xi < -smooth_limit)
1514 : {
1515 24342 : std::vector<const Node *> en;
1516 24342 : en.push_back(side->node_ptr(3));
1517 24342 : en.push_back(side->node_ptr(0));
1518 24342 : edge_nodes.push_back(en);
1519 24342 : Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
1520 24342 : if (fw > 0.5)
1521 13006 : fw = 0.5;
1522 24342 : edge_face_weights.push_back(fw);
1523 24342 : }
1524 249444 : break;
1525 : }
1526 :
1527 0 : default:
1528 : {
1529 0 : mooseError("Unsupported face type: ", t);
1530 : break;
1531 : }
1532 : }
1533 259380 : }
1534 :
1535 : void
1536 240032 : PenetrationThread::getInfoForFacesWithCommonNodes(
1537 : const Node * secondary_node,
1538 : const std::set<dof_id_type> & elems_to_exclude,
1539 : const std::vector<const Node *> edge_nodes,
1540 : std::vector<PenetrationInfo *> & face_info_comm_edge,
1541 : std::vector<PenetrationInfo *> & p_info)
1542 : {
1543 : // elems connected to a node on this edge, find one that has the same corners as this, and is not
1544 : // the current elem
1545 240032 : auto node_to_elem_pair = _node_to_elem_map.find(edge_nodes[0]->id()); // just need one of the
1546 : // nodes
1547 : mooseAssert(node_to_elem_pair != _node_to_elem_map.end(), "Missing entry in node to elem map");
1548 240032 : const std::vector<dof_id_type> & elems_connected_to_node = node_to_elem_pair->second;
1549 :
1550 240032 : std::vector<const Elem *> elems_connected_to_edge;
1551 :
1552 770689 : for (unsigned int ecni = 0; ecni < elems_connected_to_node.size(); ecni++)
1553 : {
1554 530657 : if (elems_to_exclude.find(elems_connected_to_node[ecni]) != elems_to_exclude.end())
1555 269076 : continue;
1556 261581 : const Elem * elem = _mesh.elemPtr(elems_connected_to_node[ecni]);
1557 :
1558 261581 : std::vector<const Node *> nodevec;
1559 4037689 : for (unsigned int ni = 0; ni < elem->n_nodes(); ++ni)
1560 3776108 : if (elem->is_vertex(ni))
1561 2084944 : nodevec.push_back(elem->node_ptr(ni));
1562 :
1563 261581 : std::vector<const Node *> common_nodes;
1564 261581 : std::sort(nodevec.begin(), nodevec.end());
1565 261581 : std::set_intersection(edge_nodes.begin(),
1566 : edge_nodes.end(),
1567 : nodevec.begin(),
1568 : nodevec.end(),
1569 : std::inserter(common_nodes, common_nodes.end()));
1570 :
1571 261581 : if (common_nodes.size() == edge_nodes.size())
1572 145053 : elems_connected_to_edge.push_back(elem);
1573 261581 : }
1574 :
1575 240032 : if (elems_connected_to_edge.size() > 0)
1576 : {
1577 :
1578 : // There are potentially multiple elements that share a common edge
1579 : // 2D:
1580 : // There can only be one element on the same surface
1581 : // 3D:
1582 : // If there are two edge nodes, there can only be one element on the same surface
1583 : // If there is only one edge node (a corner), there could be multiple elements on the same
1584 : // surface
1585 145053 : bool allowMultipleNeighbors = false;
1586 :
1587 145053 : if (elems_connected_to_edge[0]->dim() == 3)
1588 : {
1589 143127 : if (edge_nodes.size() == 1)
1590 : {
1591 14522 : allowMultipleNeighbors = true;
1592 : }
1593 : }
1594 :
1595 159575 : for (unsigned int i = 0; i < elems_connected_to_edge.size(); ++i)
1596 : {
1597 145053 : std::vector<PenetrationInfo *> thisElemInfo;
1598 145053 : getInfoForElem(thisElemInfo, p_info, elems_connected_to_edge[i]);
1599 145053 : if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
1600 : {
1601 34942 : if (thisElemInfo.size() > 1)
1602 0 : mooseError(
1603 : "Found multiple neighbors to current edge/face on surface when only one is allowed");
1604 34942 : face_info_comm_edge.push_back(thisElemInfo[0]);
1605 34942 : break;
1606 : }
1607 :
1608 110111 : createInfoForElem(
1609 110111 : thisElemInfo, p_info, secondary_node, elems_connected_to_edge[i], edge_nodes);
1610 110111 : if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
1611 : {
1612 95589 : if (thisElemInfo.size() > 1)
1613 0 : mooseError(
1614 : "Found multiple neighbors to current edge/face on surface when only one is allowed");
1615 95589 : face_info_comm_edge.push_back(thisElemInfo[0]);
1616 95589 : break;
1617 : }
1618 :
1619 29044 : for (unsigned int j = 0; j < thisElemInfo.size(); ++j)
1620 14522 : face_info_comm_edge.push_back(thisElemInfo[j]);
1621 145053 : }
1622 : }
1623 240032 : }
1624 :
1625 : void
1626 145053 : PenetrationThread::getInfoForElem(std::vector<PenetrationInfo *> & thisElemInfo,
1627 : std::vector<PenetrationInfo *> & p_info,
1628 : const Elem * elem)
1629 : {
1630 293091 : for (const auto & pi : p_info)
1631 : {
1632 148038 : if (!pi)
1633 41842 : continue;
1634 :
1635 106196 : if (pi->_elem == elem)
1636 41842 : thisElemInfo.push_back(pi);
1637 : }
1638 145053 : }
1639 :
1640 : void
1641 878104 : PenetrationThread::createInfoForElem(std::vector<PenetrationInfo *> & thisElemInfo,
1642 : std::vector<PenetrationInfo *> & p_info,
1643 : const Node * secondary_node,
1644 : const Elem * elem,
1645 : const std::vector<const Node *> & nodes_that_must_be_on_side,
1646 : const bool check_whether_reasonable)
1647 : {
1648 878104 : std::vector<unsigned int> sides;
1649 : // TODO: After libMesh update, add this line to MooseMesh.h, call sidesWithBoundaryID, delete
1650 : // getSidesOnPrimaryBoundary, and delete vectors used by it
1651 : // void sidesWithBoundaryID(std::vector<unsigned int>& sides, const Elem * const elem, const
1652 : // BoundaryID boundary_id) const
1653 : // {
1654 : // _mesh.get_boundary_info().sides_with_boundary_id(sides, elem, boundary_id);
1655 : // }
1656 878104 : getSidesOnPrimaryBoundary(sides, elem);
1657 : // _mesh.sidesWithBoundaryID(sides, elem, _primary_boundary);
1658 :
1659 1730890 : for (unsigned int i = 0; i < sides.size(); ++i)
1660 : {
1661 : // Don't create info for this side if one already exists
1662 867754 : bool already_have_info_this_side = false;
1663 867754 : for (const auto & pi : thisElemInfo)
1664 6900 : if (pi->_side_num == sides[i])
1665 : {
1666 6900 : already_have_info_this_side = true;
1667 6900 : break;
1668 : }
1669 :
1670 867754 : if (already_have_info_this_side)
1671 14968 : break;
1672 :
1673 860854 : const Elem * side = (elem->build_side_ptr(sides[i])).release();
1674 :
1675 : // Only continue with creating info for this side if the side contains
1676 : // all of the nodes in nodes_that_must_be_on_side
1677 860854 : std::vector<const Node *> nodevec;
1678 5514325 : for (unsigned int ni = 0; ni < side->n_nodes(); ++ni)
1679 4653471 : nodevec.push_back(side->node_ptr(ni));
1680 :
1681 860854 : std::sort(nodevec.begin(), nodevec.end());
1682 860854 : std::vector<const Node *> common_nodes;
1683 860854 : std::set_intersection(nodes_that_must_be_on_side.begin(),
1684 : nodes_that_must_be_on_side.end(),
1685 : nodevec.begin(),
1686 : nodevec.end(),
1687 : std::inserter(common_nodes, common_nodes.end()));
1688 860854 : if (common_nodes.size() != nodes_that_must_be_on_side.size())
1689 : {
1690 0 : delete side;
1691 0 : break;
1692 : }
1693 :
1694 860854 : FEBase * fe_elem = _fes[_tid][elem->dim()];
1695 860854 : FEBase * fe_side = _fes[_tid][side->dim()];
1696 :
1697 : // Optionally check to see whether face is reasonable candidate based on an
1698 : // estimate of how closely it is likely to project to the face
1699 860854 : if (check_whether_reasonable)
1700 757643 : if (!isFaceReasonableCandidate(elem, side, fe_side, secondary_node, _tangential_tolerance))
1701 : {
1702 8068 : delete side;
1703 8068 : break;
1704 : }
1705 :
1706 852786 : Point contact_phys;
1707 852786 : Point contact_ref;
1708 852786 : Point contact_on_face_ref;
1709 852786 : Real distance = 0.;
1710 852786 : Real tangential_distance = 0.;
1711 852786 : RealGradient normal;
1712 : bool contact_point_on_side;
1713 852786 : std::vector<const Node *> off_edge_nodes;
1714 852786 : std::vector<std::vector<Real>> side_phi;
1715 852786 : std::vector<std::vector<RealGradient>> side_grad_phi;
1716 852786 : std::vector<RealGradient> dxyzdxi;
1717 852786 : std::vector<RealGradient> dxyzdeta;
1718 852786 : std::vector<RealGradient> d2xyzdxideta;
1719 :
1720 : std::unique_ptr<PenetrationInfo> pen_info =
1721 : std::make_unique<PenetrationInfo>(elem,
1722 : side,
1723 852786 : sides[i],
1724 : normal,
1725 : distance,
1726 : tangential_distance,
1727 : contact_phys,
1728 : contact_ref,
1729 : contact_on_face_ref,
1730 : off_edge_nodes,
1731 : side_phi,
1732 : side_grad_phi,
1733 : dxyzdxi,
1734 : dxyzdeta,
1735 852786 : d2xyzdxideta);
1736 :
1737 852786 : bool search_succeeded = false;
1738 852786 : Moose::findContactPoint(*pen_info,
1739 : fe_elem,
1740 : fe_side,
1741 : _fe_type,
1742 : *secondary_node,
1743 : true,
1744 : _tangential_tolerance,
1745 : contact_point_on_side,
1746 : search_succeeded);
1747 :
1748 : // Do not add contact info from failed searches
1749 852786 : if (search_succeeded)
1750 : {
1751 852772 : thisElemInfo.push_back(pen_info.get());
1752 852772 : p_info.push_back(pen_info.release());
1753 : }
1754 868922 : }
1755 878104 : }
1756 :
1757 : // TODO: After libMesh update, replace this with a call to sidesWithBoundaryID, delete vectors used
1758 : // by this method
1759 : void
1760 878104 : PenetrationThread::getSidesOnPrimaryBoundary(std::vector<unsigned int> & sides,
1761 : const Elem * const elem)
1762 : {
1763 : // For each tuple, the fields are (0=elem_id, 1=side_id, 2=bc_id)
1764 878104 : sides.clear();
1765 : struct Comp
1766 : {
1767 2903074 : bool operator()(const libMesh::BoundaryInfo::BCTuple & tup, dof_id_type id) const
1768 : {
1769 2903074 : return std::get<0>(tup) < id;
1770 : }
1771 1869059 : bool operator()(dof_id_type id, const libMesh::BoundaryInfo::BCTuple & tup) const
1772 : {
1773 1869059 : return id < std::get<0>(tup);
1774 : }
1775 : };
1776 :
1777 878104 : auto range = std::equal_range(_bc_tuples.begin(), _bc_tuples.end(), elem->id(), Comp{});
1778 :
1779 1755094 : for (auto & t = range.first; t != range.second; ++t)
1780 876990 : if (std::get<2>(*t) == static_cast<boundary_id_type>(_primary_boundary))
1781 867754 : sides.push_back(std::get<1>(*t));
1782 878104 : }
|
