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 186010 : 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 186010 : const std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>> & bc_tuples)
46 186010 : : _subproblem(subproblem),
47 186010 : _mesh(mesh),
48 186010 : _primary_boundary(primary_boundary),
49 186010 : _secondary_boundary(secondary_boundary),
50 186010 : _penetration_info(penetration_info),
51 186010 : _check_whether_reasonable(check_whether_reasonable),
52 186010 : _update_location(update_location),
53 186010 : _tangential_tolerance(tangential_tolerance),
54 186010 : _do_normal_smoothing(do_normal_smoothing),
55 186010 : _normal_smoothing_distance(normal_smoothing_distance),
56 186010 : _normal_smoothing_method(normal_smoothing_method),
57 186010 : _nodal_normal_x(NULL),
58 186010 : _nodal_normal_y(NULL),
59 186010 : _nodal_normal_z(NULL),
60 186010 : _fes(fes),
61 186010 : _fe_type(fe_type),
62 186010 : _nearest_node(nearest_node),
63 186010 : _node_to_elem_map(node_to_elem_map),
64 186010 : _bc_tuples(bc_tuples)
65 : {
66 186010 : }
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 202100 : PenetrationThread::operator()(const NodeIdRange & range)
91 : {
92 202100 : ParallelUniqueId puid;
93 202100 : _tid = puid.id;
94 :
95 : // Must get the variables every time this is run because _tid can change
96 202100 : if (_do_normal_smoothing &&
97 76220 : _normal_smoothing_method == PenetrationLocator::NSM_NODAL_NORMAL_BASED)
98 : {
99 15388 : _nodal_normal_x = &_subproblem.getStandardVariable(_tid, "nodal_normal_x");
100 15388 : _nodal_normal_y = &_subproblem.getStandardVariable(_tid, "nodal_normal_y");
101 15388 : _nodal_normal_z = &_subproblem.getStandardVariable(_tid, "nodal_normal_z");
102 : }
103 :
104 1463847 : for (const auto & node_id : range)
105 : {
106 1261747 : 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 1261747 : pinfo_mutex.lock();
112 1261747 : PenetrationInfo *& info = _penetration_info[node.id()];
113 1261747 : pinfo_mutex.unlock();
114 :
115 1261747 : std::vector<PenetrationInfo *> p_info;
116 1261747 : bool info_set(false);
117 :
118 : // See if we already have info about this node
119 1261747 : if (info)
120 : {
121 715684 : FEBase * fe_elem = _fes[_tid][info->_elem->dim()];
122 715684 : FEBase * fe_side = _fes[_tid][info->_side->dim()];
123 :
124 715684 : 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 715684 : Real old_tangential_distance(info->_tangential_distance);
156 715684 : bool contact_point_on_side(false);
157 715684 : bool search_succeeded = false;
158 :
159 715684 : 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 715684 : if (contact_point_on_side)
170 : {
171 649815 : if (info->_tangential_distance <= 0.0) // on the face
172 : {
173 556331 : info_set = true;
174 : }
175 93484 : 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 90264 : 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 88548 : info_set = true;
184 : }
185 : }
186 : }
187 : }
188 : }
189 :
190 1261747 : if (!info_set)
191 : {
192 616868 : const Node * closest_node = _nearest_node.nearestNode(node.id());
193 616868 : 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 616868 : const std::vector<dof_id_type> & closest_elems = node_to_elem_pair->second;
197 :
198 1460902 : for (const auto & elem_id : closest_elems)
199 : {
200 844034 : const Elem * elem = _mesh.elemPtr(elem_id);
201 :
202 844034 : std::vector<PenetrationInfo *> thisElemInfo;
203 844034 : std::vector<const Node *> nodesThatMustBeOnSide;
204 844034 : nodesThatMustBeOnSide.push_back(closest_node);
205 844034 : createInfoForElem(
206 844034 : thisElemInfo, p_info, &node, elem, nodesThatMustBeOnSide, _check_whether_reasonable);
207 844034 : }
208 :
209 616868 : if (p_info.size() == 1)
210 : {
211 438717 : if (p_info[0]->_tangential_distance <= _tangential_tolerance)
212 : {
213 10334 : switchInfo(info, p_info[0]);
214 10334 : info_set = true;
215 : }
216 : }
217 178151 : 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 173258 : std::vector<RidgeData> ridgeDataVec;
222 384870 : for (unsigned int i = 0; i + 1 < p_info.size(); ++i)
223 482363 : for (unsigned int j = i + 1; j < p_info.size(); ++j)
224 : {
225 270751 : Point closest_coor;
226 270751 : Real tangential_distance(0.0);
227 270751 : const Node * closest_node_on_ridge = NULL;
228 270751 : unsigned int index = 0;
229 270751 : Point closest_coor_ref;
230 270751 : 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 270751 : if (found_ridge_contact_point)
239 : {
240 99340 : RidgeData rpd;
241 99340 : rpd._closest_coor = closest_coor;
242 99340 : rpd._tangential_distance = tangential_distance;
243 99340 : rpd._closest_node = closest_node_on_ridge;
244 99340 : rpd._index = index;
245 99340 : rpd._closest_coor_ref = closest_coor_ref;
246 99340 : ridgeDataVec.push_back(rpd);
247 : }
248 : }
249 :
250 173258 : 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 79318 : std::vector<RidgeSetData> ridgeSetDataVec;
255 178658 : for (unsigned int i = 0; i < ridgeDataVec.size(); ++i)
256 : {
257 99340 : bool foundSetWithMatchingNode = false;
258 118082 : for (unsigned int j = 0; j < ridgeSetDataVec.size(); ++j)
259 : {
260 35170 : if (ridgeDataVec[i]._closest_node != NULL &&
261 9970 : ridgeDataVec[i]._closest_node == ridgeSetDataVec[j]._closest_node)
262 : {
263 6458 : foundSetWithMatchingNode = true;
264 6458 : ridgeSetDataVec[j]._ridge_data_vec.push_back(ridgeDataVec[i]);
265 6458 : break;
266 : }
267 : }
268 99340 : if (!foundSetWithMatchingNode)
269 : {
270 92882 : RidgeSetData rsd;
271 92882 : rsd._distance = std::numeric_limits<Real>::max();
272 92882 : rsd._ridge_data_vec.push_back(ridgeDataVec[i]);
273 92882 : rsd._closest_node = ridgeDataVec[i]._closest_node;
274 92882 : ridgeSetDataVec.push_back(rsd);
275 92882 : }
276 : }
277 : // Compute distance to each set of ridges
278 172200 : for (unsigned int i = 0; i < ridgeSetDataVec.size(); ++i)
279 : {
280 92882 : if (ridgeSetDataVec[i]._closest_node !=
281 : NULL) // Either a peak or off the edge of single ridge
282 : {
283 50889 : if (ridgeSetDataVec[i]._ridge_data_vec.size() == 1) // off edge of single ridge
284 : {
285 45307 : if (ridgeSetDataVec[i]._ridge_data_vec[0]._tangential_distance <=
286 45307 : _tangential_tolerance) // off within tolerance
287 : {
288 14568 : ridgeSetDataVec[i]._closest_coor =
289 14568 : ridgeSetDataVec[i]._ridge_data_vec[0]._closest_coor;
290 14568 : Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
291 14568 : 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 5582 : ridgeSetDataVec[i]._closest_coor = *ridgeSetDataVec[i]._closest_node;
298 5582 : Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
299 5582 : ridgeSetDataVec[i]._distance = contact_point_vec.norm();
300 : }
301 : }
302 : else // on a single ridge
303 : {
304 41993 : ridgeSetDataVec[i]._closest_coor =
305 41993 : ridgeSetDataVec[i]._ridge_data_vec[0]._closest_coor;
306 41993 : Point contact_point_vec = node - ridgeSetDataVec[i]._closest_coor;
307 41993 : ridgeSetDataVec[i]._distance = contact_point_vec.norm();
308 : }
309 : }
310 : // Find the set of ridges closest to us.
311 79318 : unsigned int closest_ridge_set_index(0);
312 79318 : Real closest_distance(ridgeSetDataVec[0]._distance);
313 79318 : Point closest_point(ridgeSetDataVec[0]._closest_coor);
314 92882 : for (unsigned int i = 1; i < ridgeSetDataVec.size(); ++i)
315 : {
316 13564 : if (ridgeSetDataVec[i]._distance < closest_distance)
317 : {
318 2198 : closest_ridge_set_index = i;
319 2198 : closest_distance = ridgeSetDataVec[i]._distance;
320 2198 : closest_point = ridgeSetDataVec[i]._closest_coor;
321 : }
322 : }
323 :
324 79318 : if (closest_distance <
325 79318 : 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 50863 : unsigned int face_index(std::numeric_limits<unsigned int>::max());
330 107210 : for (unsigned int i = 0;
331 107210 : i < ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size();
332 : ++i)
333 : {
334 56347 : if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[i]._index < face_index)
335 50863 : 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 50863 : p_info[face_index]->_closest_point = closest_point;
342 101726 : p_info[face_index]->_distance =
343 50863 : (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 50863 : if (!_do_normal_smoothing)
348 : {
349 13249 : Point normal(closest_point - node);
350 13249 : const Real len(normal.norm());
351 13249 : if (len > 0)
352 : {
353 13186 : normal /= len;
354 : }
355 13249 : const Real dot(normal * p_info[face_index]->_normal);
356 13249 : if (dot < 0)
357 : {
358 10692 : normal *= -1;
359 : }
360 13249 : p_info[face_index]->_normal = normal;
361 : }
362 50863 : p_info[face_index]->_tangential_distance = 0.0;
363 :
364 50863 : Point closest_point_ref;
365 50863 : if (ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec.size() ==
366 : 1) // contact with a single ridge rather than a peak
367 : {
368 92510 : p_info[face_index]->_tangential_distance =
369 46255 : ridgeSetDataVec[closest_ridge_set_index]._ridge_data_vec[0]._tangential_distance;
370 46255 : p_info[face_index]->_closest_point_ref =
371 46255 : 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 4608 : bool restricted = restrictPointToFace(p_info[face_index]->_closest_point_ref,
377 : closest_node_on_face,
378 4608 : p_info[face_index]->_side);
379 4608 : if (restricted)
380 : {
381 4608 : 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 50863 : FEBase * fe = _fes[_tid][p_info[face_index]->_side->dim()];
390 50863 : std::vector<Point> points(1);
391 50863 : points[0] = p_info[face_index]->_closest_point_ref;
392 50863 : fe->reinit(p_info[face_index]->_side, &points);
393 50863 : p_info[face_index]->_side_phi = fe->get_phi();
394 50863 : p_info[face_index]->_side_grad_phi = fe->get_dphi();
395 50863 : p_info[face_index]->_dxyzdxi = fe->get_dxyzdxi();
396 50863 : p_info[face_index]->_dxyzdeta = fe->get_dxyzdeta();
397 50863 : p_info[face_index]->_d2xyzdxideta = fe->get_d2xyzdxideta();
398 :
399 50863 : switchInfo(info, p_info[face_index]);
400 50863 : info_set = true;
401 50863 : }
402 : else
403 : { // todo:remove invalid ridge cases so they don't mess up individual face competition????
404 : }
405 79318 : }
406 :
407 173258 : if (!info_set) // contact wasn't on a ridge -- compete individual interactions
408 : {
409 122395 : unsigned int best(0), i(1);
410 : do
411 : {
412 135048 : CompeteInteractionResult CIResult = competeInteractions(p_info[best], p_info[i]);
413 135048 : if (CIResult == FIRST_WINS)
414 : {
415 22344 : i++;
416 : }
417 112704 : else if (CIResult == SECOND_WINS)
418 : {
419 11522 : best = i;
420 11522 : i++;
421 : }
422 101182 : else if (CIResult == NEITHER_WINS)
423 : {
424 101182 : best = i + 1;
425 101182 : i += 2;
426 : }
427 135048 : } while (i < p_info.size() && best < p_info.size());
428 122395 : if (best < p_info.size())
429 : {
430 : // Ensure final info is within the tangential tolerance
431 26254 : if (p_info[best]->_tangential_distance <= _tangential_tolerance)
432 : {
433 26156 : switchInfo(info, p_info[best]);
434 26156 : info_set = true;
435 : }
436 : }
437 : }
438 173258 : }
439 : }
440 :
441 1261747 : 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 529515 : _recheck_secondary_nodes.push_back(node_id);
449 :
450 529515 : delete info;
451 529515 : info = NULL;
452 : }
453 : else
454 : {
455 732232 : smoothNormal(info, p_info, node);
456 732232 : FEBase * fe = _fes[_tid][info->_side->dim()];
457 732232 : computeSlip(*fe, *info);
458 : }
459 :
460 2198926 : for (unsigned int j = 0; j < p_info.size(); ++j)
461 : {
462 937179 : if (p_info[j])
463 : {
464 849826 : delete p_info[j];
465 849826 : p_info[j] = NULL;
466 : }
467 : }
468 1261747 : }
469 202100 : }
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 87353 : PenetrationThread::switchInfo(PenetrationInfo *& info, PenetrationInfo *& infoNew)
481 : {
482 : mooseAssert(infoNew != NULL, "infoNew object is null");
483 87353 : if (info)
484 : {
485 62920 : infoNew->_starting_elem = info->_starting_elem;
486 62920 : infoNew->_starting_side_num = info->_starting_side_num;
487 62920 : infoNew->_starting_closest_point_ref = info->_starting_closest_point_ref;
488 62920 : infoNew->_incremental_slip = info->_incremental_slip;
489 62920 : infoNew->_accumulated_slip = info->_accumulated_slip;
490 62920 : infoNew->_accumulated_slip_old = info->_accumulated_slip_old;
491 62920 : infoNew->_frictional_energy = info->_frictional_energy;
492 62920 : infoNew->_frictional_energy_old = info->_frictional_energy_old;
493 62920 : infoNew->_contact_force = info->_contact_force;
494 62920 : infoNew->_contact_force_old = info->_contact_force_old;
495 62920 : infoNew->_lagrange_multiplier = info->_lagrange_multiplier;
496 62920 : infoNew->_lagrange_multiplier_slip = info->_lagrange_multiplier_slip;
497 62920 : infoNew->_locked_this_step = info->_locked_this_step;
498 62920 : infoNew->_stick_locked_this_step = info->_stick_locked_this_step;
499 62920 : infoNew->_mech_status = info->_mech_status;
500 62920 : infoNew->_mech_status_old = info->_mech_status_old;
501 : }
502 : else
503 : {
504 24433 : infoNew->_starting_elem = infoNew->_elem;
505 24433 : infoNew->_starting_side_num = infoNew->_side_num;
506 24433 : infoNew->_starting_closest_point_ref = infoNew->_closest_point_ref;
507 : }
508 87353 : delete info;
509 87353 : info = infoNew;
510 87353 : infoNew = NULL; // Set this to NULL so that we don't delete it (now owned by _penetration_info).
511 87353 : }
512 :
513 : PenetrationThread::CompeteInteractionResult
514 135048 : PenetrationThread::competeInteractions(PenetrationInfo * pi1, PenetrationInfo * pi2)
515 : {
516 :
517 135048 : CompeteInteractionResult result = NEITHER_WINS;
518 :
519 135048 : if (pi1->_tangential_distance > _tangential_tolerance &&
520 112220 : pi2->_tangential_distance > _tangential_tolerance) // out of tol on both faces
521 101182 : result = NEITHER_WINS;
522 :
523 33866 : else if (pi1->_tangential_distance == 0.0 &&
524 21568 : pi2->_tangential_distance > 0.0) // on face 1, off face 2
525 19259 : result = FIRST_WINS;
526 :
527 14607 : else if (pi2->_tangential_distance == 0.0 &&
528 12691 : pi1->_tangential_distance > 0.0) // on face 2, off face 1
529 10382 : result = SECOND_WINS;
530 :
531 4225 : else if (pi1->_tangential_distance <= _tangential_tolerance &&
532 3425 : pi2->_tangential_distance > _tangential_tolerance) // in face 1 tol, out of face 2 tol
533 332 : result = FIRST_WINS;
534 :
535 3893 : else if (pi2->_tangential_distance <= _tangential_tolerance &&
536 3893 : pi1->_tangential_distance > _tangential_tolerance) // in face 2 tol, out of face 1 tol
537 800 : result = SECOND_WINS;
538 :
539 3093 : else if (pi1->_tangential_distance == 0.0 && pi2->_tangential_distance == 0.0) // on both faces
540 2309 : result = competeInteractionsBothOnFace(pi1, pi2);
541 :
542 784 : else if (pi1->_tangential_distance <= _tangential_tolerance &&
543 784 : pi2->_tangential_distance <= _tangential_tolerance) // off but within tol of both faces
544 : {
545 784 : CommonEdgeResult cer = interactionsOffCommonEdge(pi1, pi2);
546 784 : 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 784 : else if (cer == EDGE_AND_COMMON_NODE) // off side of face, off corner of another face. Favor
553 : // the off-side face
554 : {
555 500 : if (pi1->_off_edge_nodes.size() == pi2->_off_edge_nodes.size())
556 0 : mooseError("Invalid off_edge_nodes counts");
557 :
558 500 : else if (pi1->_off_edge_nodes.size() == 2)
559 160 : result = FIRST_WINS;
560 :
561 340 : else if (pi2->_off_edge_nodes.size() == 2)
562 340 : 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 284 : result = competeInteractionsBothOnFace(pi1, pi2);
569 : }
570 :
571 135048 : return result;
572 : }
573 :
574 : PenetrationThread::CompeteInteractionResult
575 2593 : PenetrationThread::competeInteractionsBothOnFace(PenetrationInfo * pi1, PenetrationInfo * pi2)
576 : {
577 2593 : CompeteInteractionResult result = NEITHER_WINS;
578 :
579 2593 : 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 2593 : 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 2593 : 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 2593 : 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 2593 : if (pi1->_elem->id() < pi2->_elem->id())
596 2593 : result = FIRST_WINS;
597 :
598 : else
599 0 : result = SECOND_WINS;
600 : }
601 :
602 2593 : return result;
603 : }
604 :
605 : PenetrationThread::CommonEdgeResult
606 784 : PenetrationThread::interactionsOffCommonEdge(PenetrationInfo * pi1, PenetrationInfo * pi2)
607 : {
608 784 : CommonEdgeResult common_edge(NO_COMMON);
609 784 : const std::vector<const Node *> & off_edge_nodes1 = pi1->_off_edge_nodes;
610 784 : const std::vector<const Node *> & off_edge_nodes2 = pi2->_off_edge_nodes;
611 784 : const unsigned dim1 = pi1->_side->dim();
612 :
613 784 : 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 784 : if (off_edge_nodes1.size() == 1)
628 : {
629 340 : 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 340 : else if (off_edge_nodes2.size() == 2)
635 : {
636 340 : if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[0] == off_edge_nodes2[1])
637 340 : common_edge = EDGE_AND_COMMON_NODE;
638 : }
639 : }
640 444 : else if (off_edge_nodes1.size() == 2)
641 : {
642 444 : if (off_edge_nodes2.size() == 1)
643 : {
644 160 : if (off_edge_nodes1[0] == off_edge_nodes2[0] || off_edge_nodes1[1] == off_edge_nodes2[0])
645 160 : common_edge = EDGE_AND_COMMON_NODE;
646 : }
647 284 : else if (off_edge_nodes2.size() == 2)
648 : {
649 284 : if ((off_edge_nodes1[0] == off_edge_nodes2[0] &&
650 568 : off_edge_nodes1[1] == off_edge_nodes2[1]) ||
651 284 : (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 784 : return common_edge;
657 : }
658 :
659 : bool
660 270751 : 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 270751 : tangential_distance = 0.0;
670 270751 : closest_node = NULL;
671 270751 : PenetrationInfo * pi1 = p_info[index1];
672 270751 : PenetrationInfo * pi2 = p_info[index2];
673 270751 : const unsigned sidedim(pi1->_side->dim());
674 : mooseAssert(sidedim == pi2->_side->dim(), "Incompatible dimensionalities");
675 :
676 : // Nodes on faces for the two interactions
677 270751 : std::vector<const Node *> side1_nodes;
678 270751 : getSideCornerNodes(pi1->_side, side1_nodes);
679 270751 : std::vector<const Node *> side2_nodes;
680 270751 : getSideCornerNodes(pi2->_side, side2_nodes);
681 :
682 270751 : std::sort(side1_nodes.begin(), side1_nodes.end());
683 270751 : std::sort(side2_nodes.begin(), side2_nodes.end());
684 :
685 : // Find nodes shared by the two faces
686 270751 : std::vector<const Node *> common_nodes;
687 270751 : 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 270751 : if (common_nodes.size() != sidedim)
694 40570 : return false;
695 :
696 : bool found_point1, found_point2;
697 230181 : Point closest_coor_ref1(pi1->_closest_point_ref);
698 : const Node * closest_node1;
699 230181 : found_point1 = restrictPointToSpecifiedEdgeOfFace(
700 : closest_coor_ref1, closest_node1, pi1->_side, common_nodes);
701 :
702 230181 : Point closest_coor_ref2(pi2->_closest_point_ref);
703 : const Node * closest_node2;
704 230181 : found_point2 = restrictPointToSpecifiedEdgeOfFace(
705 : closest_coor_ref2, closest_node2, pi2->_side, common_nodes);
706 :
707 230181 : if (!found_point1 || !found_point2)
708 130841 : 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 99340 : FEBase * fe = NULL;
720 99340 : 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 99340 : if (index1 < index2)
725 : {
726 99340 : fe = _fes[_tid][pi1->_side->dim()];
727 99340 : contact_point_ref = closest_coor_ref1;
728 99340 : points[0] = closest_coor_ref1;
729 99340 : fe->reinit(pi1->_side, &points);
730 99340 : 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 99340 : contact_point = fe->get_xyz()[0];
742 :
743 99340 : if (sidedim == 2)
744 : {
745 94243 : 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 57347 : closest_node = closest_node1;
750 :
751 57347 : RealGradient off_face = *closest_node1 - contact_point;
752 57347 : tangential_distance = off_face.norm();
753 : }
754 : }
755 :
756 99340 : return true;
757 270751 : }
758 :
759 : void
760 541502 : PenetrationThread::getSideCornerNodes(const Elem * side, std::vector<const Node *> & corner_nodes)
761 : {
762 541502 : const ElemType t(side->type());
763 541502 : corner_nodes.clear();
764 :
765 541502 : corner_nodes.push_back(side->node_ptr(0));
766 541502 : corner_nodes.push_back(side->node_ptr(1));
767 541502 : switch (t)
768 : {
769 36594 : case EDGE2:
770 : case EDGE3:
771 : case EDGE4:
772 : {
773 36594 : break;
774 : }
775 :
776 12112 : case TRI3:
777 : case TRI6:
778 : case TRI7:
779 : {
780 12112 : corner_nodes.push_back(side->node_ptr(2));
781 12112 : break;
782 : }
783 :
784 492796 : case QUAD4:
785 : case QUAD8:
786 : case QUAD9:
787 : {
788 492796 : corner_nodes.push_back(side->node_ptr(2));
789 492796 : corner_nodes.push_back(side->node_ptr(3));
790 492796 : break;
791 : }
792 :
793 0 : default:
794 : {
795 0 : mooseError("Unsupported face type: ", t);
796 : break;
797 : }
798 : }
799 541502 : }
800 :
801 : bool
802 460362 : PenetrationThread::restrictPointToSpecifiedEdgeOfFace(Point & p,
803 : const Node *& closest_node,
804 : const Elem * side,
805 : const std::vector<const Node *> & edge_nodes)
806 : {
807 460362 : const ElemType t = side->type();
808 460362 : Real & xi = p(0);
809 460362 : Real & eta = p(1);
810 460362 : closest_node = NULL;
811 :
812 460362 : std::vector<unsigned int> local_node_indices;
813 1344492 : for (const auto & edge_node : edge_nodes)
814 : {
815 884130 : unsigned int local_index = side->get_node_index(edge_node);
816 884130 : if (local_index == libMesh::invalid_uint)
817 0 : mooseError("Side does not contain node");
818 884130 : 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 460362 : std::sort(local_node_indices.begin(), local_node_indices.end());
823 :
824 460362 : bool off_of_this_edge = false;
825 :
826 460362 : switch (t)
827 : {
828 36594 : case EDGE2:
829 : case EDGE3:
830 : case EDGE4:
831 : {
832 36594 : if (local_node_indices[0] == 0)
833 : {
834 18297 : if (xi <= -1.0)
835 : {
836 13781 : xi = -1.0;
837 13781 : off_of_this_edge = true;
838 13781 : closest_node = side->node_ptr(0);
839 : }
840 : }
841 18297 : else if (local_node_indices[0] == 1)
842 : {
843 18297 : if (xi >= 1.0)
844 : {
845 9601 : xi = 1.0;
846 9601 : off_of_this_edge = true;
847 9601 : closest_node = side->node_ptr(1);
848 : }
849 : }
850 : else
851 : {
852 0 : mooseError("Invalid local node indices");
853 : }
854 36594 : break;
855 : }
856 :
857 5036 : case TRI3:
858 : case TRI6:
859 : case TRI7:
860 : {
861 5036 : if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
862 : {
863 1312 : if (eta <= 0.0)
864 : {
865 660 : eta = 0.0;
866 660 : off_of_this_edge = true;
867 660 : if (xi < 0.0)
868 109 : closest_node = side->node_ptr(0);
869 551 : else if (xi > 1.0)
870 233 : closest_node = side->node_ptr(1);
871 : }
872 : }
873 3724 : else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
874 : {
875 1898 : if ((xi + eta) > 1.0)
876 : {
877 896 : Real delta = (xi + eta - 1.0) / 2.0;
878 896 : xi -= delta;
879 896 : eta -= delta;
880 896 : off_of_this_edge = true;
881 896 : if (xi > 1.0)
882 209 : closest_node = side->node_ptr(1);
883 687 : else if (xi < 0.0)
884 242 : closest_node = side->node_ptr(2);
885 : }
886 : }
887 1826 : else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 2))
888 : {
889 1826 : if (xi <= 0.0)
890 : {
891 918 : xi = 0.0;
892 918 : off_of_this_edge = true;
893 918 : if (eta > 1.0)
894 281 : closest_node = side->node_ptr(2);
895 637 : else if (eta < 0.0)
896 231 : closest_node = side->node_ptr(0);
897 : }
898 : }
899 : else
900 : {
901 0 : mooseError("Invalid local node indices");
902 : }
903 :
904 5036 : break;
905 : }
906 :
907 418732 : case QUAD4:
908 : case QUAD8:
909 : case QUAD9:
910 : {
911 418732 : if ((local_node_indices[0] == 0) && (local_node_indices[1] == 1))
912 : {
913 181510 : if (eta <= -1.0)
914 : {
915 111216 : eta = -1.0;
916 111216 : off_of_this_edge = true;
917 111216 : if (xi < -1.0)
918 34897 : closest_node = side->node_ptr(0);
919 76319 : else if (xi > 1.0)
920 34965 : closest_node = side->node_ptr(1);
921 : }
922 : }
923 237222 : else if ((local_node_indices[0] == 1) && (local_node_indices[1] == 2))
924 : {
925 88092 : if (xi >= 1.0)
926 : {
927 72544 : xi = 1.0;
928 72544 : off_of_this_edge = true;
929 72544 : if (eta < -1.0)
930 16548 : closest_node = side->node_ptr(1);
931 55996 : else if (eta > 1.0)
932 13602 : closest_node = side->node_ptr(2);
933 : }
934 : }
935 149130 : else if ((local_node_indices[0] == 2) && (local_node_indices[1] == 3))
936 : {
937 122306 : if (eta >= 1.0)
938 : {
939 90633 : eta = 1.0;
940 90633 : off_of_this_edge = true;
941 90633 : if (xi < -1.0)
942 39377 : closest_node = side->node_ptr(3);
943 51256 : else if (xi > 1.0)
944 45953 : closest_node = side->node_ptr(2);
945 : }
946 : }
947 26824 : else if ((local_node_indices[0] == 0) && (local_node_indices[1] == 3))
948 : {
949 26824 : if (xi <= -1.0)
950 : {
951 18258 : xi = -1.0;
952 18258 : off_of_this_edge = true;
953 18258 : if (eta < -1.0)
954 10870 : closest_node = side->node_ptr(0);
955 7388 : else if (eta > 1.0)
956 2918 : closest_node = side->node_ptr(3);
957 : }
958 : }
959 : else
960 : {
961 0 : mooseError("Invalid local node indices");
962 : }
963 418732 : break;
964 : }
965 :
966 0 : default:
967 : {
968 0 : mooseError("Unsupported face type: ", t);
969 : break;
970 : }
971 : }
972 460362 : return off_of_this_edge;
973 460362 : }
974 :
975 : bool
976 4608 : PenetrationThread::restrictPointToFace(Point & p, const Node *& closest_node, const Elem * side)
977 : {
978 4608 : const ElemType t(side->type());
979 4608 : Real & xi = p(0);
980 4608 : Real & eta = p(1);
981 4608 : closest_node = NULL;
982 :
983 4608 : bool off_of_this_face(false);
984 :
985 4608 : 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 4608 : case QUAD4:
1073 : case QUAD8:
1074 : case QUAD9:
1075 : {
1076 4608 : 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 4500 : else if (xi > 1.0)
1094 : {
1095 4500 : xi = 1.0;
1096 4500 : off_of_this_face = true;
1097 4500 : 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 4500 : closest_node = side->node_ptr(2);
1106 4500 : if (eta > 1.0)
1107 1170 : 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 4608 : break;
1145 : }
1146 :
1147 0 : default:
1148 : {
1149 0 : mooseError("Unsupported face type: ", t);
1150 : break;
1151 : }
1152 : }
1153 4608 : return off_of_this_face;
1154 : }
1155 :
1156 : bool
1157 832484 : 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 832484 : unsigned int dim = primary_elem->dim();
1164 :
1165 832484 : const std::vector<Point> & phys_point = fe->get_xyz();
1166 :
1167 832484 : const std::vector<RealGradient> & dxyz_dxi = fe->get_dxyzdxi();
1168 832484 : const std::vector<RealGradient> & dxyz_deta = fe->get_dxyzdeta();
1169 :
1170 832484 : Point ref_point;
1171 :
1172 832484 : std::vector<Point> points(1); // Default constructor gives us a point at 0,0,0
1173 :
1174 832484 : fe->reinit(side, &points);
1175 :
1176 832484 : RealGradient d = *secondary_point - phys_point[0];
1177 :
1178 832484 : const Real twosqrt2 = 2.8284; // way more precision than we actually need here
1179 832484 : Real max_face_length = side->hmax() + twosqrt2 * tangential_tolerance;
1180 :
1181 832484 : RealVectorValue normal;
1182 832484 : if (dim - 1 == 2)
1183 : {
1184 749481 : normal = dxyz_dxi[0].cross(dxyz_deta[0]);
1185 : }
1186 83003 : else if (dim - 1 == 1)
1187 : {
1188 82982 : const Node * const * elem_nodes = primary_elem->get_nodes();
1189 82982 : const Point in_plane_vector1 = *elem_nodes[1] - *elem_nodes[0];
1190 82982 : const Point in_plane_vector2 = *elem_nodes[2] - *elem_nodes[0];
1191 :
1192 82982 : Point out_of_plane_normal = in_plane_vector1.cross(in_plane_vector2);
1193 82982 : out_of_plane_normal /= out_of_plane_normal.norm();
1194 :
1195 82982 : normal = dxyz_dxi[0].cross(out_of_plane_normal);
1196 : }
1197 : else
1198 : {
1199 21 : return true;
1200 : }
1201 832463 : normal /= normal.norm();
1202 :
1203 832463 : const Real dot(d * normal);
1204 :
1205 832463 : const RealGradient normcomp = dot * normal;
1206 832463 : const RealGradient tangcomp = d - normcomp;
1207 :
1208 832463 : 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 832463 : const Real faceExpansionFactor = 2.0 * (1.0 + normcomp.norm() / d.norm());
1213 :
1214 832463 : bool isReasonableCandidate = true;
1215 832463 : if (tangdist > faceExpansionFactor * max_face_length)
1216 : {
1217 8881 : isReasonableCandidate = false;
1218 : }
1219 832463 : return isReasonableCandidate;
1220 832484 : }
1221 :
1222 : void
1223 732232 : 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 732232 : std::vector<Point> points(1);
1228 732232 : points[0] = info._starting_closest_point_ref;
1229 732232 : const auto & side = _elem_side_builder(*info._starting_elem, info._starting_side_num);
1230 732232 : fe.reinit(&side, &points);
1231 732232 : const std::vector<Point> & starting_point = fe.get_xyz();
1232 732232 : info._incremental_slip = info._closest_point - starting_point[0];
1233 732232 : if (info.isCaptured())
1234 : {
1235 57248 : info._frictional_energy =
1236 57248 : info._frictional_energy_old + info._contact_force * info._incremental_slip;
1237 57248 : info._accumulated_slip = info._accumulated_slip_old + info._incremental_slip.norm();
1238 : }
1239 732232 : }
1240 :
1241 : void
1242 732232 : PenetrationThread::smoothNormal(PenetrationInfo * info,
1243 : std::vector<PenetrationInfo *> & p_info,
1244 : const Node & node)
1245 : {
1246 732232 : if (_do_normal_smoothing)
1247 : {
1248 356926 : 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 285510 : std::vector<Real> edge_face_weights;
1253 285510 : std::vector<PenetrationInfo *> edge_face_info;
1254 :
1255 285510 : 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 285510 : if (edge_face_info.size() > 0)
1261 : {
1262 : // Smooth the normal using the weighting functions for all participating faces.
1263 127626 : RealVectorValue new_normal;
1264 127626 : Real this_face_weight = 1.0;
1265 :
1266 287244 : for (unsigned int efwi = 0; efwi < edge_face_weights.size(); ++efwi)
1267 : {
1268 159618 : PenetrationInfo * npi = edge_face_info[efwi];
1269 159618 : if (npi)
1270 159618 : new_normal += npi->_normal * edge_face_weights[efwi];
1271 :
1272 159618 : 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 127626 : new_normal += info->_normal * this_face_weight;
1277 :
1278 127626 : const Real len = new_normal.norm();
1279 127626 : if (len > 0)
1280 127626 : new_normal /= len;
1281 :
1282 127626 : info->_normal = new_normal;
1283 : }
1284 285510 : }
1285 71416 : else if (_normal_smoothing_method == PenetrationLocator::NSM_NODAL_NORMAL_BASED)
1286 : {
1287 : // params.addParam<VariableName>("var_name","description");
1288 : // getParam<VariableName>("var_name")
1289 71416 : info->_normal(0) = _nodal_normal_x->getValue(info->_side, info->_side_phi);
1290 71416 : info->_normal(1) = _nodal_normal_y->getValue(info->_side, info->_side_phi);
1291 71416 : info->_normal(2) = _nodal_normal_z->getValue(info->_side, info->_side_phi);
1292 71416 : const Real len(info->_normal.norm());
1293 71416 : if (len > 0)
1294 70600 : info->_normal /= len;
1295 : }
1296 : }
1297 732232 : }
1298 :
1299 : void
1300 285510 : 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 285510 : const Elem * side = info->_side;
1307 285510 : const Point & p = info->_closest_point_ref;
1308 285510 : std::set<dof_id_type> elems_to_exclude;
1309 285510 : elems_to_exclude.insert(info->_elem->id());
1310 :
1311 285510 : 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 285510 : getSmoothingEdgeNodesAndWeights(p, side, edge_nodes, edge_face_weights);
1315 285510 : std::vector<Elem *> edge_neighbor_elems;
1316 285510 : edge_face_info.resize(edge_nodes.size(), NULL);
1317 :
1318 285510 : std::vector<unsigned int> edges_without_neighbors;
1319 :
1320 534093 : for (unsigned int i = 0; i < edge_nodes.size(); ++i)
1321 : {
1322 : // Sort all sets of edge nodes (needed for comparing edges)
1323 248583 : std::sort(edge_nodes[i].begin(), edge_nodes[i].end());
1324 :
1325 248583 : std::vector<PenetrationInfo *> face_info_comm_edge;
1326 248583 : getInfoForFacesWithCommonNodes(
1327 248583 : &secondary_node, elems_to_exclude, edge_nodes[i], face_info_comm_edge, p_info);
1328 :
1329 248583 : if (face_info_comm_edge.size() == 0)
1330 104961 : edges_without_neighbors.push_back(i);
1331 143622 : else if (face_info_comm_edge.size() > 1)
1332 0 : mooseError("Only one neighbor allowed per edge");
1333 : else
1334 143622 : edge_face_info[i] = face_info_comm_edge[0];
1335 248583 : }
1336 :
1337 : // Remove edges without neighbors from the vector, starting from end
1338 285510 : std::vector<unsigned int>::reverse_iterator rit;
1339 390471 : for (rit = edges_without_neighbors.rbegin(); rit != edges_without_neighbors.rend(); ++rit)
1340 : {
1341 104961 : unsigned int index = *rit;
1342 104961 : edge_nodes.erase(edge_nodes.begin() + index);
1343 104961 : edge_face_weights.erase(edge_face_weights.begin() + index);
1344 104961 : edge_face_info.erase(edge_face_info.begin() + index);
1345 : }
1346 :
1347 : // Handle corner case
1348 285510 : if (edge_nodes.size() > 1)
1349 : {
1350 15996 : if (edge_nodes.size() != 2)
1351 0 : mooseError("Invalid number of smoothing edges");
1352 :
1353 : // find common node
1354 15996 : std::vector<const Node *> common_nodes;
1355 63984 : std::set_intersection(edge_nodes[0].begin(),
1356 15996 : edge_nodes[0].end(),
1357 15996 : edge_nodes[1].begin(),
1358 15996 : edge_nodes[1].end(),
1359 : std::inserter(common_nodes, common_nodes.end()));
1360 :
1361 15996 : if (common_nodes.size() != 1)
1362 0 : mooseError("Invalid number of common nodes");
1363 :
1364 47988 : for (const auto & pinfo : edge_face_info)
1365 31992 : elems_to_exclude.insert(pinfo->_elem->id());
1366 :
1367 15996 : std::vector<PenetrationInfo *> face_info_comm_edge;
1368 15996 : getInfoForFacesWithCommonNodes(
1369 : &secondary_node, elems_to_exclude, common_nodes, face_info_comm_edge, p_info);
1370 :
1371 15996 : unsigned int num_corner_neighbors = face_info_comm_edge.size();
1372 :
1373 15996 : if (num_corner_neighbors > 0)
1374 : {
1375 15996 : Real fw0 = edge_face_weights[0];
1376 15996 : Real fw1 = edge_face_weights[1];
1377 :
1378 : // Corner weight is product of edge weights. Spread out over multiple neighbors.
1379 15996 : Real fw_corner = (fw0 * fw1) / static_cast<Real>(num_corner_neighbors);
1380 :
1381 : // Adjust original edge weights
1382 15996 : edge_face_weights[0] *= (1.0 - fw1);
1383 15996 : edge_face_weights[1] *= (1.0 - fw0);
1384 :
1385 31992 : for (unsigned int i = 0; i < num_corner_neighbors; ++i)
1386 : {
1387 15996 : edge_face_weights.push_back(fw_corner);
1388 15996 : edge_face_info.push_back(face_info_comm_edge[i]);
1389 : }
1390 : }
1391 15996 : }
1392 285510 : }
1393 :
1394 : void
1395 285510 : 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 285510 : const ElemType t(side->type());
1402 285510 : const Real & xi = p(0);
1403 285510 : const Real & eta = p(1);
1404 :
1405 285510 : Real smooth_limit = 1.0 - _normal_smoothing_distance;
1406 :
1407 285510 : switch (t)
1408 : {
1409 10870 : case EDGE2:
1410 : case EDGE3:
1411 : case EDGE4:
1412 : {
1413 10870 : if (xi < -smooth_limit)
1414 : {
1415 2178 : std::vector<const Node *> en;
1416 2178 : en.push_back(side->node_ptr(0));
1417 2178 : edge_nodes.push_back(en);
1418 2178 : Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
1419 2178 : if (fw > 0.5)
1420 144 : fw = 0.5;
1421 2178 : edge_face_weights.push_back(fw);
1422 2178 : }
1423 8692 : else if (xi > smooth_limit)
1424 : {
1425 1242 : std::vector<const Node *> en;
1426 1242 : en.push_back(side->node_ptr(1));
1427 1242 : edge_nodes.push_back(en);
1428 1242 : Real fw = 0.5 - (1.0 - xi) / (2.0 * _normal_smoothing_distance);
1429 1242 : if (fw > 0.5)
1430 198 : fw = 0.5;
1431 1242 : edge_face_weights.push_back(fw);
1432 1242 : }
1433 10870 : 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 274640 : case QUAD4:
1477 : case QUAD8:
1478 : case QUAD9:
1479 : {
1480 274640 : if (eta < -smooth_limit)
1481 : {
1482 45392 : std::vector<const Node *> en;
1483 45392 : en.push_back(side->node_ptr(0));
1484 45392 : en.push_back(side->node_ptr(1));
1485 45392 : edge_nodes.push_back(en);
1486 45392 : Real fw = 0.5 - (1.0 + eta) / (2.0 * _normal_smoothing_distance);
1487 45392 : if (fw > 0.5)
1488 15492 : fw = 0.5;
1489 45392 : edge_face_weights.push_back(fw);
1490 45392 : }
1491 274640 : if (xi > smooth_limit)
1492 : {
1493 85787 : std::vector<const Node *> en;
1494 85787 : en.push_back(side->node_ptr(1));
1495 85787 : en.push_back(side->node_ptr(2));
1496 85787 : edge_nodes.push_back(en);
1497 85787 : Real fw = 0.5 - (1.0 - xi) / (2.0 * _normal_smoothing_distance);
1498 85787 : if (fw > 0.5)
1499 19532 : fw = 0.5;
1500 85787 : edge_face_weights.push_back(fw);
1501 85787 : }
1502 274640 : if (eta > smooth_limit)
1503 : {
1504 87184 : std::vector<const Node *> en;
1505 87184 : en.push_back(side->node_ptr(2));
1506 87184 : en.push_back(side->node_ptr(3));
1507 87184 : edge_nodes.push_back(en);
1508 87184 : Real fw = 0.5 - (1.0 - eta) / (2.0 * _normal_smoothing_distance);
1509 87184 : if (fw > 0.5)
1510 25251 : fw = 0.5;
1511 87184 : edge_face_weights.push_back(fw);
1512 87184 : }
1513 274640 : if (xi < -smooth_limit)
1514 : {
1515 26800 : std::vector<const Node *> en;
1516 26800 : en.push_back(side->node_ptr(3));
1517 26800 : en.push_back(side->node_ptr(0));
1518 26800 : edge_nodes.push_back(en);
1519 26800 : Real fw = 0.5 - (1.0 + xi) / (2.0 * _normal_smoothing_distance);
1520 26800 : if (fw > 0.5)
1521 14352 : fw = 0.5;
1522 26800 : edge_face_weights.push_back(fw);
1523 26800 : }
1524 274640 : break;
1525 : }
1526 :
1527 0 : default:
1528 : {
1529 0 : mooseError("Unsupported face type: ", t);
1530 : break;
1531 : }
1532 : }
1533 285510 : }
1534 :
1535 : void
1536 264579 : 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 264579 : 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 264579 : const std::vector<dof_id_type> & elems_connected_to_node = node_to_elem_pair->second;
1549 :
1550 264579 : std::vector<const Elem *> elems_connected_to_edge;
1551 :
1552 846168 : for (unsigned int ecni = 0; ecni < elems_connected_to_node.size(); ecni++)
1553 : {
1554 581589 : if (elems_to_exclude.find(elems_connected_to_node[ecni]) != elems_to_exclude.end())
1555 296571 : continue;
1556 285018 : const Elem * elem = _mesh.elemPtr(elems_connected_to_node[ecni]);
1557 :
1558 285018 : std::vector<const Node *> nodevec;
1559 4510174 : for (unsigned int ni = 0; ni < elem->n_nodes(); ++ni)
1560 4225156 : if (elem->is_vertex(ni))
1561 2271696 : nodevec.push_back(elem->node_ptr(ni));
1562 :
1563 285018 : std::vector<const Node *> common_nodes;
1564 285018 : std::sort(nodevec.begin(), nodevec.end());
1565 285018 : 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 285018 : if (common_nodes.size() == edge_nodes.size())
1572 159618 : elems_connected_to_edge.push_back(elem);
1573 285018 : }
1574 :
1575 264579 : 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 159618 : bool allowMultipleNeighbors = false;
1586 :
1587 159618 : if (elems_connected_to_edge[0]->dim() == 3)
1588 : {
1589 157506 : if (edge_nodes.size() == 1)
1590 : {
1591 15996 : allowMultipleNeighbors = true;
1592 : }
1593 : }
1594 :
1595 175614 : for (unsigned int i = 0; i < elems_connected_to_edge.size(); ++i)
1596 : {
1597 159618 : std::vector<PenetrationInfo *> thisElemInfo;
1598 159618 : getInfoForElem(thisElemInfo, p_info, elems_connected_to_edge[i]);
1599 159618 : if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
1600 : {
1601 38438 : if (thisElemInfo.size() > 1)
1602 0 : mooseError(
1603 : "Found multiple neighbors to current edge/face on surface when only one is allowed");
1604 38438 : face_info_comm_edge.push_back(thisElemInfo[0]);
1605 38438 : break;
1606 : }
1607 :
1608 121180 : createInfoForElem(
1609 121180 : thisElemInfo, p_info, secondary_node, elems_connected_to_edge[i], edge_nodes);
1610 121180 : if (thisElemInfo.size() > 0 && !allowMultipleNeighbors)
1611 : {
1612 105184 : if (thisElemInfo.size() > 1)
1613 0 : mooseError(
1614 : "Found multiple neighbors to current edge/face on surface when only one is allowed");
1615 105184 : face_info_comm_edge.push_back(thisElemInfo[0]);
1616 105184 : break;
1617 : }
1618 :
1619 31992 : for (unsigned int j = 0; j < thisElemInfo.size(); ++j)
1620 15996 : face_info_comm_edge.push_back(thisElemInfo[j]);
1621 159618 : }
1622 : }
1623 264579 : }
1624 :
1625 : void
1626 159618 : PenetrationThread::getInfoForElem(std::vector<PenetrationInfo *> & thisElemInfo,
1627 : std::vector<PenetrationInfo *> & p_info,
1628 : const Elem * elem)
1629 : {
1630 322518 : for (const auto & pi : p_info)
1631 : {
1632 162900 : if (!pi)
1633 46026 : continue;
1634 :
1635 116874 : if (pi->_elem == elem)
1636 46026 : thisElemInfo.push_back(pi);
1637 : }
1638 159618 : }
1639 :
1640 : void
1641 965214 : 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 965214 : 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 965214 : getSidesOnPrimaryBoundary(sides, elem);
1657 : // _mesh.sidesWithBoundaryID(sides, elem, _primary_boundary);
1658 :
1659 1902409 : for (unsigned int i = 0; i < sides.size(); ++i)
1660 : {
1661 : // Don't create info for this side if one already exists
1662 953664 : bool already_have_info_this_side = false;
1663 953664 : for (const auto & pi : thisElemInfo)
1664 7588 : if (pi->_side_num == sides[i])
1665 : {
1666 7588 : already_have_info_this_side = true;
1667 7588 : break;
1668 : }
1669 :
1670 953664 : if (already_have_info_this_side)
1671 16469 : break;
1672 :
1673 946076 : 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 946076 : std::vector<const Node *> nodevec;
1678 6066307 : for (unsigned int ni = 0; ni < side->n_nodes(); ++ni)
1679 5120231 : nodevec.push_back(side->node_ptr(ni));
1680 :
1681 946076 : std::sort(nodevec.begin(), nodevec.end());
1682 946076 : std::vector<const Node *> common_nodes;
1683 946076 : 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 946076 : if (common_nodes.size() != nodes_that_must_be_on_side.size())
1689 : {
1690 0 : delete side;
1691 0 : break;
1692 : }
1693 :
1694 946076 : FEBase * fe_elem = _fes[_tid][elem->dim()];
1695 946076 : 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 946076 : if (check_whether_reasonable)
1700 832484 : if (!isFaceReasonableCandidate(elem, side, fe_side, secondary_node, _tangential_tolerance))
1701 : {
1702 8881 : delete side;
1703 8881 : break;
1704 : }
1705 :
1706 937195 : Point contact_phys;
1707 937195 : Point contact_ref;
1708 937195 : Point contact_on_face_ref;
1709 937195 : Real distance = 0.;
1710 937195 : Real tangential_distance = 0.;
1711 937195 : RealGradient normal;
1712 : bool contact_point_on_side;
1713 937195 : std::vector<const Node *> off_edge_nodes;
1714 937195 : std::vector<std::vector<Real>> side_phi;
1715 937195 : std::vector<std::vector<RealGradient>> side_grad_phi;
1716 937195 : std::vector<RealGradient> dxyzdxi;
1717 937195 : std::vector<RealGradient> dxyzdeta;
1718 937195 : std::vector<RealGradient> d2xyzdxideta;
1719 :
1720 : std::unique_ptr<PenetrationInfo> pen_info =
1721 : std::make_unique<PenetrationInfo>(elem,
1722 : side,
1723 937195 : 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 937195 : d2xyzdxideta);
1736 :
1737 937195 : bool search_succeeded = false;
1738 937195 : 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 937195 : if (search_succeeded)
1750 : {
1751 937179 : thisElemInfo.push_back(pen_info.get());
1752 937179 : p_info.push_back(pen_info.release());
1753 : }
1754 954957 : }
1755 965214 : }
1756 :
1757 : // TODO: After libMesh update, replace this with a call to sidesWithBoundaryID, delete vectors used
1758 : // by this method
1759 : void
1760 965214 : 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 965214 : sides.clear();
1765 : struct Comp
1766 : {
1767 3197751 : bool operator()(const libMesh::BoundaryInfo::BCTuple & tup, dof_id_type id) const
1768 : {
1769 3197751 : return std::get<0>(tup) < id;
1770 : }
1771 2051423 : bool operator()(dof_id_type id, const libMesh::BoundaryInfo::BCTuple & tup) const
1772 : {
1773 2051423 : return id < std::get<0>(tup);
1774 : }
1775 : };
1776 :
1777 965214 : auto range = std::equal_range(_bc_tuples.begin(), _bc_tuples.end(), elem->id(), Comp{});
1778 :
1779 1929192 : for (auto & t = range.first; t != range.second; ++t)
1780 963978 : if (std::get<2>(*t) == static_cast<boundary_id_type>(_primary_boundary))
1781 953664 : sides.push_back(std::get<1>(*t));
1782 965214 : }
|
