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 "MeshCut2DRankTwoTensorNucleation.h"
11 :
12 : #include "MooseError.h"
13 : #include "MooseTypes.h"
14 : #include "libmesh/quadrature.h"
15 : #include "RankTwoTensor.h"
16 : #include "RankTwoScalarTools.h"
17 : #include "Assembly.h"
18 :
19 : registerMooseObject("XFEMApp", MeshCut2DRankTwoTensorNucleation);
20 :
21 : InputParameters
22 34 : MeshCut2DRankTwoTensorNucleation::validParams()
23 : {
24 34 : InputParameters params = MeshCut2DNucleationBase::validParams();
25 34 : params.addClassDescription(
26 : "Nucleate a crack in MeshCut2D UO based on a scalar extracted from a RankTwoTensor");
27 102 : params.addRangeCheckedParam<Real>(
28 : "edge_extension_factor",
29 68 : 1e-5,
30 : "edge_extension_factor >= 0",
31 : "Factor by which the length of cracks that have been adjusted to coincide with element edges "
32 : "are increased to avoid missing an intersection with the edge due to numerical tolerance "
33 : "issues.");
34 68 : params.addRangeCheckedParam<Real>(
35 : "nucleation_length",
36 : "nucleation_length >= 0",
37 : "Size of crack to nucleate. Must be larger than the length of the element in which the "
38 : "crack is nucleated, unless 'nucleate_across_full_element' is set to 'true'.");
39 68 : params.addParam<MooseEnum>(
40 : "scalar_type",
41 68 : RankTwoScalarTools::scalarOptions(),
42 : "Scalar quantity to be computed from tensor and used as a failure criterion");
43 68 : params.addRequiredParam<std::string>("tensor", "The material tensor name.");
44 68 : params.addRequiredCoupledVar(
45 : "nucleation_threshold",
46 : "Threshold for the scalar quantity of the RankTwoTensor to nucleate new cracks");
47 34 : params.addParam<Point>(
48 : "point1",
49 34 : Point(0, 0, 0),
50 : "Start point for axis used to calculate some cylindrical material tensor quantities");
51 34 : params.addParam<Point>(
52 : "point2",
53 34 : Point(0, 1, 0),
54 : "End point for axis used to calculate some cylindrical material tensor quantities");
55 68 : params.addParam<bool>(
56 : "nucleate_across_full_element",
57 68 : false,
58 : "Controls the behavior of nucleating cracks if 'nucleation_length' is smaller than the "
59 : "length required to travserse an element with a nucleating crack. If this is set to 'false', "
60 : "that condition will result in an error, but if set to 'true', a crack with the length "
61 : "needed "
62 : "to traverse the element will be inserted.");
63 34 : return params;
64 0 : }
65 :
66 17 : MeshCut2DRankTwoTensorNucleation::MeshCut2DRankTwoTensorNucleation(
67 17 : const InputParameters & parameters)
68 : : MeshCut2DNucleationBase(parameters),
69 17 : _edge_extension_factor(getParam<Real>("edge_extension_factor")),
70 34 : _nucleate_across_full_element(getParam<bool>("nucleate_across_full_element")),
71 34 : _is_nucleation_length_provided(isParamValid("nucleation_length")),
72 30 : _nucleation_length(_is_nucleation_length_provided ? getParam<Real>("nucleation_length") : 0),
73 34 : _tensor(getMaterialProperty<RankTwoTensor>(getParam<std::string>("tensor"))),
74 17 : _nucleation_threshold(coupledValue("nucleation_threshold")),
75 51 : _scalar_type(getParam<MooseEnum>("scalar_type")),
76 17 : _point1(parameters.get<Point>("point1")),
77 17 : _point2(parameters.get<Point>("point2")),
78 17 : _JxW(_assembly.JxW()),
79 17 : _coord(_assembly.coordTransformation())
80 : {
81 17 : }
82 :
83 : bool
84 48972 : MeshCut2DRankTwoTensorNucleation::doesElementCrack(
85 : std::pair<RealVectorValue, RealVectorValue> & cutterElemNodes)
86 : {
87 : bool does_it_crack = false;
88 48972 : unsigned int numqp = _qrule->n_points();
89 : Point zero; // Used for checking whether cutter Element is zero length
90 :
91 : Real average_threshold = 0.0;
92 48972 : RankTwoTensor average_tensor;
93 : Point average_point;
94 244860 : for (unsigned int qp = 0; qp < numqp; ++qp)
95 : {
96 195888 : if (_nucleation_threshold[qp] <= 0.0)
97 0 : mooseError("Threshold must be strictly positive in MeshCut2DRankTwoTensorNucleation");
98 195888 : average_threshold += _JxW[qp] * _coord[qp] * _nucleation_threshold[qp];
99 195888 : average_tensor += _JxW[qp] * _coord[qp] * _tensor[qp];
100 195888 : average_point += _JxW[qp] * _coord[qp] * _q_point[qp];
101 : }
102 : Point point_dir;
103 48972 : Real tensor_quantity = RankTwoScalarTools::getQuantity(
104 48972 : average_tensor, _scalar_type, _point1, _point2, average_point, point_dir);
105 48972 : if (point_dir.absolute_fuzzy_equals(zero))
106 0 : mooseError("Cutter element has zero length in MeshCut2DRankTwoTensorNucleation");
107 :
108 48972 : if (tensor_quantity > average_threshold)
109 : {
110 : does_it_crack = true;
111 :
112 749 : const Point elem_center(_current_elem->vertex_average());
113 749 : Point crack_dir = point_dir.cross(Point(0, 0, 1));
114 :
115 : // get max element length to use for ray length
116 749 : std::unique_ptr<const Elem> edge;
117 : Real circumference = 0;
118 3745 : for (const auto e : _current_elem->edge_index_range())
119 : {
120 2996 : _current_elem->build_edge_ptr(edge, e);
121 2996 : circumference += edge->length(0, 1);
122 : }
123 :
124 : // Temporaries for doing things below
125 : Real intersection_distance;
126 :
127 : Point point_0;
128 : Point point_1;
129 : bool is_point_0_on_external_boundary = false;
130 : bool is_point_1_on_external_boundary = false;
131 3745 : for (const auto e : _current_elem->edge_index_range())
132 : {
133 2996 : _current_elem->build_edge_ptr(edge, e);
134 : const auto & edge0 = edge->point(0);
135 : const auto & edge1 = edge->point(1);
136 2996 : if (lineLineIntersect2D(
137 : elem_center, -crack_dir, circumference, edge0, edge1, point_0, intersection_distance))
138 749 : is_point_0_on_external_boundary = (_current_elem->neighbor_ptr(e) == nullptr);
139 2247 : else if (lineLineIntersect2D(elem_center,
140 : crack_dir,
141 : circumference,
142 : edge0,
143 : edge1,
144 : point_1,
145 : intersection_distance))
146 749 : is_point_1_on_external_boundary = (_current_elem->neighbor_ptr(e) == nullptr);
147 : }
148 :
149 : // traverse_length is the length of a cut that goes across a single elment
150 : // extend_length is the amount that needs to be added to each side of the traverse_length cut to
151 : // make its length equal to the nucleation length
152 : // We also add a small amount to the crack length to make sure it fully cuts the element edge.
153 : // This factor is based on the element length.
154 749 : Real traverse_length = (point_0 - point_1).norm();
155 749 : Real extend_length = (_nucleation_length - traverse_length) / 2.0;
156 749 : if (_nucleation_length < traverse_length)
157 : {
158 25 : if (_is_nucleation_length_provided && !_nucleate_across_full_element)
159 1 : mooseError(
160 : "Trying to nucleate crack smaller than element length, increase nucleation_length.\n "
161 : "Location of crack being nucleated: ",
162 : point_0,
163 : "\n nucleation_length: ",
164 1 : _nucleation_length,
165 : "\n Length to traverse element: ",
166 : traverse_length,
167 : "\n This error can be suppressed by setting nucleate_across_full_element=True which "
168 : "will nucleate a crack the size of the element.");
169 : else
170 : //_edge_extension_factor is used to make sure cut will cut both edges of the element.
171 24 : extend_length = (traverse_length * _edge_extension_factor) / 2.0;
172 : }
173 :
174 : // First create a cut assuming bulk nucleation
175 748 : point_0 = point_0 - extend_length * crack_dir.unit();
176 748 : point_1 = point_1 + extend_length * crack_dir.unit();
177 :
178 : // modify edges of cut to put them on the edge of the domain and have a length
179 : // equal to nucleation_length
180 748 : if (is_point_0_on_external_boundary)
181 : {
182 : // first move point back to edge
183 198 : point_0 = point_0 + extend_length * crack_dir.unit();
184 : // second move point a little bit over edge
185 198 : point_0 = point_0 - (traverse_length * _edge_extension_factor) / 2.0 * crack_dir.unit();
186 : // third, move point 1 so that crack equals nucleation length
187 198 : point_1 = point_1 + extend_length * crack_dir.unit();
188 : }
189 550 : else if (is_point_1_on_external_boundary)
190 : {
191 9 : point_1 = point_1 - extend_length * crack_dir.unit();
192 9 : point_1 = point_1 + (traverse_length * _edge_extension_factor) / 2.0 * crack_dir.unit();
193 9 : point_0 = point_0 - extend_length * crack_dir.unit();
194 : }
195 748 : cutterElemNodes = {point_0, point_1};
196 748 : }
197 :
198 48971 : return does_it_crack;
199 : }
200 :
201 : bool
202 5243 : MeshCut2DRankTwoTensorNucleation::lineLineIntersect2D(const Point & start,
203 : const Point & direction,
204 : const Real length,
205 : const Point & v0,
206 : const Point & v1,
207 : Point & intersection_point,
208 : Real & intersection_distance)
209 : {
210 : /// The standard "looser" tolerance to use in ray tracing when having difficulty finding intersection
211 : const Real TRACE_TOLERANCE = 1e-5;
212 :
213 : const auto r = direction * length;
214 : const auto s = v1 - v0;
215 5243 : const auto rxs = r(0) * s(1) - r(1) * s(0);
216 :
217 : // Lines are parallel or colinear
218 5243 : if (std::abs(rxs) < TRACE_TOLERANCE)
219 : return false;
220 :
221 : const auto v0mu0 = v0 - start;
222 5131 : const auto t = (v0mu0(0) * s(1) - v0mu0(1) * s(0)) / rxs;
223 5131 : if (0 >= t + TRACE_TOLERANCE || t - TRACE_TOLERANCE > 1.0)
224 : {
225 : return false;
226 : }
227 :
228 2666 : const auto u = (v0mu0(0) * r(1) - v0mu0(1) * r(0)) / rxs;
229 2666 : if (0 < u + TRACE_TOLERANCE && u - TRACE_TOLERANCE <= 1.0)
230 : {
231 1498 : intersection_point = start + r * t;
232 1498 : intersection_distance = t * length;
233 1498 : return true;
234 : }
235 :
236 : // Not parallel, but don't intersect
237 : return false;
238 : }
|
