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