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 : #include "MeshCut2DUserObjectBase.h"
11 : #include "MeshCut2DNucleationBase.h"
12 : #include "CrackFrontDefinition.h"
13 :
14 : #include "XFEMFuncs.h"
15 : #include "MooseError.h"
16 : #include "MooseMesh.h"
17 : #include "libmesh/edge_edge2.h"
18 : #include "libmesh/serial_mesh.h"
19 : #include "libmesh/mesh_tools.h"
20 :
21 : InputParameters
22 91 : MeshCut2DUserObjectBase::validParams()
23 : {
24 91 : InputParameters params = MeshCutUserObjectBase::validParams();
25 182 : params.addParam<UserObjectName>("nucleate_uo", "The MeshCutNucleation UO for nucleating cracks.");
26 182 : params.addParam<UserObjectName>("crack_front_definition",
27 : "crackFrontDefinition",
28 : "The CrackFrontDefinition user object name");
29 91 : params.addClassDescription("Creates a UserObject base class for a mesh cutter in 2D problems");
30 91 : return params;
31 0 : }
32 :
33 46 : MeshCut2DUserObjectBase::MeshCut2DUserObjectBase(const InputParameters & parameters)
34 : : MeshCutUserObjectBase(parameters),
35 92 : _mesh(_subproblem.mesh()),
36 46 : _nucleate_uo(isParamValid("nucleate_uo")
37 63 : ? &getUserObject<MeshCut2DNucleationBase>("nucleate_uo")
38 : : nullptr),
39 92 : _is_mesh_modified(false)
40 : {
41 138 : _depend_uo.insert(getParam<UserObjectName>("crack_front_definition"));
42 :
43 : // test element type; only line elements are allowed
44 326 : for (const auto & cut_elem : _cutter_mesh->element_ptr_range())
45 : {
46 117 : if (cut_elem->n_nodes() != 2)
47 0 : mooseError("The input cut mesh should include EDGE2 elements only!");
48 117 : if (cut_elem->dim() != 1)
49 0 : mooseError("The input cut mesh should have 1D elements (in a 2D space) only!");
50 46 : }
51 :
52 : // find node fronts of the original cutmesh. This is used to order EVERYTHING.
53 46 : findOriginalCrackFrontNodes();
54 46 : }
55 :
56 : void
57 45 : MeshCut2DUserObjectBase::initialSetup()
58 : {
59 45 : const auto uo_name = getParam<UserObjectName>("crack_front_definition");
60 45 : _crack_front_definition = &_fe_problem.getUserObject<CrackFrontDefinition>(uo_name);
61 45 : }
62 :
63 : bool
64 334904 : MeshCut2DUserObjectBase::cutElementByGeometry(const Elem * elem,
65 : std::vector<Xfem::CutEdge> & cut_edges,
66 : std::vector<Xfem::CutNode> & cut_nodes) const
67 : {
68 : // With the crack defined by a line, this method cuts a 2D elements by a line
69 : // Fixme lynn Copy and paste from InterfaceMeshCut2DUserObject::cutElementByGeometry
70 : mooseAssert(elem->dim() == 2, "Dimension of element to be cut must be 2");
71 :
72 : bool elem_cut = false;
73 :
74 4586838 : for (const auto & cut_elem : _cutter_mesh->element_ptr_range())
75 : {
76 1958515 : unsigned int n_sides = elem->n_sides();
77 :
78 9792575 : for (unsigned int i = 0; i < n_sides; ++i)
79 : {
80 7834060 : std::unique_ptr<const Elem> curr_side = elem->side_ptr(i);
81 :
82 : mooseAssert(curr_side->type() == EDGE2, "Element side type must be EDGE2.");
83 :
84 : const Node * node1 = curr_side->node_ptr(0);
85 : const Node * node2 = curr_side->node_ptr(1);
86 7834060 : Real seg_int_frac = 0.0;
87 :
88 7834060 : const std::pair<Point, Point> elem_endpoints(cut_elem->node_ref(0), cut_elem->node_ref(1));
89 :
90 7834060 : if (Xfem::intersectSegmentWithCutLine(*node1, *node2, elem_endpoints, 1.0, seg_int_frac))
91 : {
92 19929 : if (seg_int_frac > Xfem::tol && seg_int_frac < 1.0 - Xfem::tol)
93 : {
94 : elem_cut = true;
95 : Xfem::CutEdge mycut;
96 19929 : mycut._id1 = node1->id();
97 19929 : mycut._id2 = node2->id();
98 19929 : mycut._distance = seg_int_frac;
99 19929 : mycut._host_side_id = i;
100 19929 : cut_edges.push_back(mycut);
101 19929 : }
102 0 : else if (seg_int_frac < Xfem::tol)
103 : {
104 : elem_cut = true;
105 : Xfem::CutNode mycut;
106 0 : mycut._id = node1->id();
107 0 : mycut._host_id = i;
108 0 : cut_nodes.push_back(mycut);
109 : }
110 : }
111 7834060 : }
112 334904 : }
113 334904 : return elem_cut;
114 : }
115 :
116 : bool
117 0 : MeshCut2DUserObjectBase::cutElementByGeometry(const Elem * /*elem*/,
118 : std::vector<Xfem::CutFace> & /*cut_faces*/) const
119 : {
120 0 : mooseError("Invalid method for 2D mesh cutting.");
121 : return false;
122 : }
123 :
124 : bool
125 9527 : MeshCut2DUserObjectBase::cutFragmentByGeometry(std::vector<std::vector<Point>> & frag_edges,
126 : std::vector<Xfem::CutEdge> & cut_edges) const
127 : {
128 : bool cut_frag = false;
129 :
130 281284 : for (const auto & cut_elem : _cutter_mesh->element_ptr_range())
131 : {
132 131115 : const std::pair<Point, Point> elem_endpoints(cut_elem->node_ref(0), cut_elem->node_ref(1));
133 131115 : unsigned int n_sides = frag_edges.size();
134 665512 : for (unsigned int i = 0; i < n_sides; ++i)
135 : {
136 534397 : Real seg_int_frac = 0.0;
137 534397 : if (Xfem::intersectSegmentWithCutLine(
138 534397 : frag_edges[i][0], frag_edges[i][1], elem_endpoints, 1, seg_int_frac))
139 : {
140 : cut_frag = true;
141 : Xfem::CutEdge mycut;
142 17698 : mycut._id1 = i;
143 17698 : mycut._id2 = (i < (n_sides - 1) ? (i + 1) : 0);
144 17698 : mycut._distance = seg_int_frac;
145 17698 : mycut._host_side_id = i;
146 17698 : cut_edges.push_back(mycut);
147 : }
148 : }
149 9527 : }
150 9527 : return cut_frag;
151 : }
152 :
153 : bool
154 0 : MeshCut2DUserObjectBase::cutFragmentByGeometry(std::vector<std::vector<Point>> & /*frag_faces*/,
155 : std::vector<Xfem::CutFace> & /*cut_faces*/) const
156 : {
157 0 : mooseError("Invalid method for 2D mesh fragment cutting.");
158 : return false;
159 : }
160 :
161 : unsigned int
162 289 : MeshCut2DUserObjectBase::getNumberOfCrackFrontPoints() const
163 : {
164 289 : return _original_and_current_front_node_ids.size();
165 : }
166 :
167 : const std::vector<Point>
168 285 : MeshCut2DUserObjectBase::getCrackFrontPoints(unsigned int number_crack_front_points) const
169 : {
170 285 : std::vector<Point> crack_front_points(number_crack_front_points);
171 : // number_crack_front_points is updated via
172 : // _crack_front_definition->updateNumberOfCrackFrontPoints(_crack_front_points.size())
173 285 : if (number_crack_front_points != _original_and_current_front_node_ids.size())
174 0 : mooseError("Number of nodes in CrackFrontDefinition does not match the number of nodes in the "
175 0 : "cutter_mesh.\nCrackFrontDefinition nodes = " +
176 0 : Moose::stringify(number_crack_front_points) + "\ncutter_mesh nodes = " +
177 0 : Moose::stringify(_original_and_current_front_node_ids.size()));
178 :
179 1033 : for (unsigned int i = 0; i < number_crack_front_points; ++i)
180 : {
181 748 : dof_id_type id = _original_and_current_front_node_ids[i].second;
182 748 : Node * this_node = _cutter_mesh->node_ptr(id);
183 : mooseAssert(this_node, "Node is NULL");
184 : Point & this_point = *this_node;
185 748 : crack_front_points[i] = this_point;
186 : }
187 285 : return crack_front_points;
188 0 : }
189 :
190 : const std::vector<RealVectorValue>
191 285 : MeshCut2DUserObjectBase::getCrackPlaneNormals(unsigned int number_crack_front_points) const
192 : {
193 285 : if (number_crack_front_points != _original_and_current_front_node_ids.size())
194 0 : mooseError("MeshCut2DFractureUserObject::getCrackPlaneNormals: number_crack_front_points=" +
195 0 : Moose::stringify(number_crack_front_points) +
196 : " does not match the number of nodes given in "
197 0 : "_original_and_current_front_node_ids=" +
198 : Moose::stringify(_original_and_current_front_node_ids.size()) +
199 : ". This will happen if a crack front exits the boundary because the number of "
200 : "points in the CrackFrontDefinition is never updated.");
201 :
202 : std::vector<std::pair<dof_id_type, RealVectorValue>> crack_plane_normals;
203 5504 : for (const auto & elem : _cutter_mesh->element_ptr_range())
204 : {
205 2467 : dof_id_type id0 = elem->node_id(0);
206 2467 : dof_id_type id1 = elem->node_id(1);
207 : dof_id_type id;
208 :
209 2467 : auto it0 = std::find_if(_original_and_current_front_node_ids.begin(),
210 : _original_and_current_front_node_ids.end(),
211 : [&id0](const std::pair<dof_id_type, dof_id_type> & element)
212 5320 : { return element.second == id0; });
213 2467 : auto it1 = std::find_if(_original_and_current_front_node_ids.begin(),
214 : _original_and_current_front_node_ids.end(),
215 : [&id1](const std::pair<dof_id_type, dof_id_type> & element)
216 4868 : { return element.second == id1; });
217 :
218 : bool found_it0 = (it0 != _original_and_current_front_node_ids.end());
219 : bool found_it1 = (it1 != _original_and_current_front_node_ids.end());
220 :
221 : // Newly nucleated crack elements can have one normal if they are on the edge OR
222 : // two normals if they are in the bulk.
223 2467 : if (found_it0)
224 : {
225 : Point end_pt, connecting_pt;
226 :
227 152 : end_pt = elem->node_ref(0);
228 152 : connecting_pt = elem->node_ref(1);
229 152 : id = it0->first; // sort by original crack front node ids
230 :
231 : Point fracture_dir = end_pt - connecting_pt;
232 : // The crack normal is orthogonal to the crack extension direction (fracture_dir),
233 : // and is defined in this implementation as the cross product of the direction of crack
234 : // extension with the tangent direction, which is always (0, 0, 1) in 2D.
235 152 : RealVectorValue normal_dir{fracture_dir(1), -fracture_dir(0), 0};
236 152 : normal_dir /= normal_dir.norm();
237 152 : crack_plane_normals.push_back(std::make_pair(id, normal_dir));
238 : }
239 :
240 2467 : if (found_it1)
241 : {
242 : Point end_pt, connecting_pt;
243 :
244 596 : end_pt = elem->node_ref(1);
245 596 : connecting_pt = elem->node_ref(0);
246 596 : id = it1->first; // sort by original crack front node ids
247 :
248 : Point fracture_dir = end_pt - connecting_pt;
249 : // The crack normal is orthogonal to the crack extension direction (fracture_dir),
250 : // and is defined in this implementation as the cross product of the direction of crack
251 : // extension with the tangent direction, which is always (0, 0, 1) in 2D.
252 596 : RealVectorValue normal_dir{fracture_dir(1), -fracture_dir(0), 0};
253 596 : normal_dir /= normal_dir.norm();
254 596 : crack_plane_normals.push_back(std::make_pair(id, normal_dir));
255 : }
256 285 : }
257 : mooseAssert(
258 : _original_and_current_front_node_ids.size() == crack_plane_normals.size(),
259 : "Boundary nodes are attached to more than one element. This should not happen for a 1D "
260 : "cutter mesh."
261 : "\n Number of _original_and_current_front_node_ids=" +
262 : Moose::stringify(_original_and_current_front_node_ids.size()) +
263 : "\n Number of crack_plane_normals=" + Moose::stringify(crack_plane_normals.size()));
264 :
265 : // the crack_plane_normals are now sorted by the ORIGINAL crack front ids
266 285 : std::sort(crack_plane_normals.begin(), crack_plane_normals.end());
267 : std::vector<RealVectorValue> sorted_crack_plane_normals;
268 1033 : for (auto & crack : crack_plane_normals)
269 748 : sorted_crack_plane_normals.push_back(crack.second);
270 :
271 285 : return sorted_crack_plane_normals;
272 285 : }
273 :
274 : void
275 46 : MeshCut2DUserObjectBase::findOriginalCrackFrontNodes()
276 : {
277 46 : std::unique_ptr<PointLocatorBase> pl = _mesh.getPointLocator();
278 46 : pl->enable_out_of_mesh_mode();
279 46 : std::unordered_set boundary_nodes = MeshTools::find_boundary_nodes(*_cutter_mesh);
280 138 : for (const auto & node : boundary_nodes)
281 : {
282 92 : auto node_id = node;
283 92 : Node * this_node = _cutter_mesh->node_ptr(node_id);
284 : mooseAssert(this_node, "Node is NULL");
285 92 : Point & this_point = *this_node;
286 :
287 92 : const Elem * elem = (*pl)(this_point);
288 92 : if (elem != NULL)
289 54 : _original_and_current_front_node_ids.push_back(std::make_pair(node, node));
290 : }
291 46 : std::sort(_original_and_current_front_node_ids.begin(),
292 : _original_and_current_front_node_ids.end());
293 46 : }
294 :
295 : void
296 370 : MeshCut2DUserObjectBase::growFront()
297 : {
298 : dof_id_type current_front_node_id;
299 1320 : for (std::size_t i = 0; i < _original_and_current_front_node_ids.size(); ++i)
300 : {
301 950 : current_front_node_id = _original_and_current_front_node_ids[i].second;
302 : // check if node front node id is active
303 : auto direction_iter =
304 950 : std::find_if(_active_front_node_growth_vectors.begin(),
305 : _active_front_node_growth_vectors.end(),
306 : [¤t_front_node_id](const std::pair<dof_id_type, Point> & element)
307 540 : { return element.first == current_front_node_id; });
308 : // only add an element for active node front ids
309 950 : if (direction_iter != _active_front_node_growth_vectors.end())
310 : {
311 276 : Node * this_node = _cutter_mesh->node_ptr(current_front_node_id);
312 : mooseAssert(this_node, "Node is NULL");
313 : Point & this_point = *this_node;
314 :
315 276 : Point new_node_offset = direction_iter->second;
316 : Point x = this_point + new_node_offset;
317 :
318 : // TODO: Should check if cut line segment created between "this_point" and "x" crosses
319 : // another line element in the cutter mesh or solid mesh boundary.
320 : // Crossing another line element would be a special case that still needs to be handled,
321 : // however, it doesnot cause an error, it will just ignore the other line segment and recut
322 : // the solid mesh element.
323 : // Crossing a solid mesh boundary would be for aesthetics reasons so
324 : // that element was trimmed close to the boundary but would have not effect on the simulation.
325 : // Crossing a solid mesh boundary should be handled by something like
326 : // MeshCut2DRankTwoTensorNucleation::lineLineIntersect2D
327 :
328 : // add node to front
329 552 : this_node = Node::build(x, _cutter_mesh->n_nodes()).release();
330 276 : _cutter_mesh->add_node(this_node);
331 276 : dof_id_type new_front_node_id = _cutter_mesh->n_nodes() - 1;
332 :
333 : // add element to front
334 : std::vector<dof_id_type> elem;
335 276 : elem.push_back(current_front_node_id);
336 276 : elem.push_back(new_front_node_id);
337 276 : Elem * new_elem = Elem::build(EDGE2).release();
338 828 : for (unsigned int i = 0; i < new_elem->n_nodes(); ++i)
339 : {
340 : mooseAssert(_cutter_mesh->node_ptr(elem[i]) != nullptr, "Node is NULL");
341 552 : new_elem->set_node(i, _cutter_mesh->node_ptr(elem[i]));
342 : }
343 276 : _cutter_mesh->add_elem(new_elem);
344 : // now push to the end of _original_and_current_front_node_ids for tracking and fracture
345 : // integrals
346 276 : _original_and_current_front_node_ids[i].second = new_front_node_id;
347 276 : _is_mesh_modified = true;
348 276 : }
349 : }
350 370 : _cutter_mesh->prepare_for_use();
351 370 : }
352 :
353 : void
354 370 : MeshCut2DUserObjectBase::addNucleatedCracksToMesh()
355 : {
356 370 : if (_nucleate_uo)
357 : {
358 : std::map<unsigned int, std::pair<RealVectorValue, RealVectorValue>> nucleated_elems_map =
359 : _nucleate_uo->getNucleatedElemsMap();
360 184 : const Real nucleationRadius = _nucleate_uo->getNucleationRadius();
361 :
362 184 : removeNucleatedCracksTooCloseToEachOther(nucleated_elems_map, nucleationRadius);
363 184 : removeNucleatedCracksTooCloseToExistingCracks(nucleated_elems_map, nucleationRadius);
364 :
365 184 : std::unique_ptr<PointLocatorBase> pl = _mesh.getPointLocator();
366 184 : pl->enable_out_of_mesh_mode();
367 256 : for (const auto & elem_nodes : nucleated_elems_map)
368 : {
369 72 : std::pair<RealVectorValue, RealVectorValue> nodes = elem_nodes.second;
370 : // add nodes for the elements that define the nucleated cracks
371 144 : Node * node_0 = Node::build(nodes.first, _cutter_mesh->n_nodes()).release();
372 72 : _cutter_mesh->add_node(node_0);
373 72 : dof_id_type node_id_0 = _cutter_mesh->n_nodes() - 1;
374 144 : Node * node_1 = Node::build(nodes.second, _cutter_mesh->n_nodes()).release();
375 72 : _cutter_mesh->add_node(node_1);
376 72 : dof_id_type node_id_1 = _cutter_mesh->n_nodes() - 1;
377 : // add elements that define nucleated cracks
378 : std::vector<dof_id_type> elem;
379 72 : elem.push_back(node_id_0);
380 72 : elem.push_back(node_id_1);
381 72 : Elem * new_elem = Elem::build(EDGE2).release();
382 216 : for (unsigned int i = 0; i < new_elem->n_nodes(); ++i)
383 : {
384 : mooseAssert(_cutter_mesh->node_ptr(elem[i]) != nullptr, "Node is NULL");
385 144 : new_elem->set_node(i, _cutter_mesh->node_ptr(elem[i]));
386 : }
387 72 : _cutter_mesh->add_elem(new_elem);
388 : // now add the nucleated nodes to the crack id data struct
389 : // edge nucleated cracks will add one node to _original_and_current_front_node_ids
390 : // bulk nucleated cracks will add two nodes to _original_and_current_front_node_ids
391 72 : Point & point_0 = *node_0;
392 72 : const Elem * crack_front_elem_0 = (*pl)(point_0);
393 72 : if (crack_front_elem_0 != NULL)
394 36 : _original_and_current_front_node_ids.push_back(std::make_pair(node_id_0, node_id_0));
395 :
396 72 : Point & point_1 = *node_1;
397 72 : const Elem * crack_front_elem_1 = (*pl)(point_1);
398 72 : if (crack_front_elem_1 != NULL)
399 60 : _original_and_current_front_node_ids.push_back(std::make_pair(node_id_1, node_id_1));
400 :
401 72 : _is_mesh_modified = true;
402 72 : }
403 184 : _cutter_mesh->prepare_for_use();
404 184 : }
405 370 : }
406 :
407 : void
408 184 : MeshCut2DUserObjectBase::removeNucleatedCracksTooCloseToEachOther(
409 : std::map<unsigned int, std::pair<RealVectorValue, RealVectorValue>> & nucleated_elems_map,
410 : const Real nucleationRadius)
411 : {
412 : // remove nucleated elements that are too close too each other. Lowest key wins
413 658 : for (auto it1 = nucleated_elems_map.begin(); it1 != nucleated_elems_map.end(); ++it1)
414 : {
415 474 : std::pair<RealVectorValue, RealVectorValue> nodes = it1->second;
416 : Point p2 = nodes.first;
417 : Point p1 = nodes.second;
418 : Point p = p1 + (p2 - p1) / 2;
419 4132 : for (auto it2 = nucleated_elems_map.begin(); it2 != nucleated_elems_map.end();)
420 : {
421 3658 : if (it1 == it2)
422 : {
423 : ++it2;
424 474 : continue;
425 : }
426 :
427 : nodes = it2->second;
428 : p2 = nodes.first;
429 : p1 = nodes.second;
430 : Point q = p1 + (p2 - p1) / 2;
431 : Point pq = q - p;
432 3184 : if (pq.norm() <= nucleationRadius)
433 486 : it2 = nucleated_elems_map.erase(it2);
434 : else
435 : ++it2;
436 : }
437 : }
438 184 : }
439 :
440 : void
441 184 : MeshCut2DUserObjectBase::removeNucleatedCracksTooCloseToExistingCracks(
442 : std::map<unsigned int, std::pair<RealVectorValue, RealVectorValue>> & nucleated_elems_map,
443 : const Real nucleationRadius)
444 : {
445 658 : for (auto it = nucleated_elems_map.begin(); it != nucleated_elems_map.end();)
446 : {
447 474 : std::pair<RealVectorValue, RealVectorValue> nodes = it->second;
448 : Point p2 = nodes.first;
449 : Point p1 = nodes.second;
450 : Point p = p1 + (p2 - p1) / 2;
451 : bool removeNucleatedElem = false;
452 : for (const auto & cutter_elem :
453 5960 : as_range(_cutter_mesh->active_elements_begin(), _cutter_mesh->active_elements_end()))
454 : {
455 2671 : const Node * const * cutter_elem_nodes = cutter_elem->get_nodes();
456 2671 : Point m = *cutter_elem_nodes[1] - *cutter_elem_nodes[0];
457 2671 : Real t = m * (p - *cutter_elem_nodes[0]) / m.norm_sq();
458 : Real d = std::numeric_limits<Real>::max();
459 2671 : if (t <= 0)
460 : {
461 : Point j = p - *cutter_elem_nodes[0];
462 1070 : d = j.norm();
463 : }
464 1601 : else if (t >= 0)
465 : {
466 : Point j = p - *cutter_elem_nodes[1];
467 1601 : d = j.norm();
468 : }
469 : else
470 : {
471 : Point j = p - (*cutter_elem_nodes[0] + t * m);
472 0 : d = j.norm();
473 : }
474 2671 : if (d <= nucleationRadius)
475 : {
476 : removeNucleatedElem = true;
477 : break;
478 : }
479 474 : }
480 474 : if (removeNucleatedElem)
481 402 : it = nucleated_elems_map.erase(it);
482 : else
483 : ++it;
484 : }
485 184 : }
|
