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