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