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 "ParsedCurveGenerator.h"
11 : #include "LinearInterpolation.h"
12 : #include "MooseUtils.h"
13 : #include "CastUniquePointer.h"
14 : #include "MooseMeshUtils.h"
15 :
16 : #include "libmesh/fparser_ad.hh"
17 : #include "libmesh/edge_edge2.h"
18 : #include "libmesh/edge_edge3.h"
19 : #include "libmesh/edge_edge4.h"
20 :
21 : // C++ includes
22 : #include <cmath>
23 :
24 : registerMooseObject("MooseApp", ParsedCurveGenerator);
25 :
26 : InputParameters
27 14917 : ParsedCurveGenerator::validParams()
28 : {
29 14917 : InputParameters params = MeshGenerator::validParams();
30 14917 : params += FunctionParserUtils<false>::validParams();
31 :
32 14917 : MooseEnum edge_elem_type("EDGE2 EDGE3 EDGE4", "EDGE2");
33 :
34 14917 : params.addRequiredParam<std::string>("x_formula", "Function expression of x(t)");
35 14917 : params.addRequiredParam<std::string>("y_formula", "Function expression of y(t)");
36 14917 : params.addParam<std::string>("z_formula", "0", "Function expression of z(t)");
37 14917 : params.addRequiredRangeCheckedParam<std::vector<unsigned int>>(
38 : "nums_segments",
39 : "nums_segments>=1",
40 : "Numbers of segments (EDGE elements) of each section of the curve to be generated. The "
41 : "number of entries in this parameter should be equal to one less than the number of entries "
42 : "in 'section_bounding_t_values'");
43 14917 : params.addRequiredParam<std::vector<Real>>(
44 : "section_bounding_t_values",
45 : "The 't' values that bound the sections of the curve. Start and end points must be included. "
46 : "The number of entries in 'nums_segments' should be equal to one less than the number of "
47 : "entries in this parameter.");
48 14917 : params.addParam<std::vector<std::string>>(
49 : "constant_names", "Vector of constants used in the parsed function (use this for kB etc.)");
50 14917 : params.addParam<std::vector<std::string>>(
51 : "constant_expressions",
52 : "Vector of values for the constants in constant_names (can be an FParser expression)");
53 14917 : params.addParam<bool>("is_closed_loop", false, "Whether the curve is closed or not.");
54 14917 : params.addRangeCheckedParam<Real>("point_overlapping_tolerance",
55 : libMesh::TOLERANCE,
56 : "point_overlapping_tolerance>0.0",
57 : "The point-to-point distance tolerance that is used to "
58 : "determine whether the two points are overlapped.");
59 14917 : params.addRangeCheckedParam<unsigned int>(
60 : "forced_closing_num_segments",
61 : "forced_closing_num_segments>1",
62 : "Number of segments (EDGE elements) of the curve section that is generated to forcefully "
63 : "close the loop.");
64 44751 : params.addRangeCheckedParam<Real>("oversample_factor",
65 29834 : 10.0,
66 : "oversample_factor>2.0",
67 : "Oversample factor to help make node distance nearly uniform.");
68 44751 : params.addRangeCheckedParam<unsigned int>(
69 : "max_oversample_number_factor",
70 29834 : 1000,
71 : "max_oversample_number_factor>100",
72 : "For each section of the curve, the maximum number of oversampling points is the product of "
73 : "this factor and the number of nodes on the curve.");
74 44751 : params.addParam<std::vector<BoundaryName>>(
75 29834 : "edge_nodesets", std::vector<BoundaryName>(), "Nodesets on both edges of the parsed curves");
76 14917 : params.addParam<MooseEnum>(
77 : "edge_element_type", edge_elem_type, "Type of the EDGE elements to be generated.");
78 14917 : params.addClassDescription("This ParsedCurveGenerator object is designed to generate a mesh of a "
79 : "curve that consists of EDGE2, EDGE3, or EDGE4 elements.");
80 :
81 29834 : return params;
82 14917 : }
83 :
84 340 : ParsedCurveGenerator::ParsedCurveGenerator(const InputParameters & parameters)
85 : : MeshGenerator(parameters),
86 : FunctionParserUtils<false>(parameters),
87 340 : _function_x(getParam<std::string>("x_formula")),
88 340 : _function_y(getParam<std::string>("y_formula")),
89 340 : _function_z(getParam<std::string>("z_formula")),
90 340 : _nums_segments(getParam<std::vector<unsigned int>>("nums_segments")),
91 340 : _section_bounding_t_values(getParam<std::vector<Real>>("section_bounding_t_values")),
92 340 : _is_closed_loop(getParam<bool>("is_closed_loop")),
93 340 : _point_overlapping_tolerance(getParam<Real>("point_overlapping_tolerance")),
94 680 : _forced_closing_num_segments(isParamValid("forced_closing_num_segments")
95 340 : ? getParam<unsigned int>("forced_closing_num_segments")
96 : : 0),
97 340 : _oversample_factor(getParam<Real>("oversample_factor")),
98 340 : _max_oversample_number_factor(getParam<unsigned int>("max_oversample_number_factor")),
99 340 : _node_set_boundaries(getParam<std::vector<BoundaryName>>("edge_nodesets")),
100 680 : _order(getParam<MooseEnum>("edge_element_type") == "EDGE2"
101 400 : ? 1
102 1080 : : (getParam<MooseEnum>("edge_element_type") == "EDGE3" ? 2 : 3))
103 : {
104 340 : if (std::adjacent_find(_section_bounding_t_values.begin(), _section_bounding_t_values.end()) !=
105 680 : _section_bounding_t_values.end())
106 4 : paramError("section_bounding_t_values", "elements must be unique.");
107 340 : if (!std::is_sorted(_section_bounding_t_values.begin(), _section_bounding_t_values.end()) &&
108 340 : !std::is_sorted(_section_bounding_t_values.rbegin(), _section_bounding_t_values.rend()))
109 4 : paramError("section_bounding_t_values", "elements must change monotonically.");
110 332 : if (_nums_segments.size() != _section_bounding_t_values.size() - 1)
111 4 : paramError(
112 : "nums_segments",
113 : "The size of this parameter must be one less than size of section_bounding_t_values.");
114 328 : _func_Fx = std::make_shared<SymFunction>();
115 328 : _func_Fy = std::make_shared<SymFunction>();
116 328 : _func_Fz = std::make_shared<SymFunction>();
117 :
118 : // set FParser internal feature flags
119 328 : setParserFeatureFlags(_func_Fx);
120 328 : setParserFeatureFlags(_func_Fy);
121 328 : setParserFeatureFlags(_func_Fz);
122 :
123 : // add the constant expressions; note that the three functions share one set of constants
124 328 : if (isParamValid("constant_names") && isParamValid("constant_expressions"))
125 : {
126 234 : addFParserConstants(_func_Fx,
127 : getParam<std::vector<std::string>>("constant_names"),
128 : getParam<std::vector<std::string>>("constant_expressions"));
129 234 : addFParserConstants(_func_Fy,
130 : getParam<std::vector<std::string>>("constant_names"),
131 : getParam<std::vector<std::string>>("constant_expressions"));
132 234 : addFParserConstants(_func_Fz,
133 : getParam<std::vector<std::string>>("constant_names"),
134 : getParam<std::vector<std::string>>("constant_expressions"));
135 : }
136 :
137 : // parse functions
138 328 : if (_func_Fx->Parse(_function_x, "t") >= 0)
139 8 : mooseError("Invalid function x(t)\n",
140 4 : _function_x,
141 : "\nin ParsedCurveGenerator ",
142 4 : name(),
143 : ".\n",
144 4 : _func_Fx->ErrorMsg());
145 324 : if (_func_Fy->Parse(_function_y, "t") >= 0)
146 8 : mooseError("Invalid function y(t)\n",
147 4 : _function_y,
148 : "\nin ParsedCurveGenerator ",
149 4 : name(),
150 : ".\n",
151 4 : _func_Fy->ErrorMsg());
152 320 : if (_func_Fz->Parse(_function_z, "t") >= 0)
153 8 : mooseError("Invalid function z(t)\n",
154 4 : _function_z,
155 : "\nin ParsedCurveGenerator ",
156 4 : name(),
157 : ".\n",
158 4 : _func_Fz->ErrorMsg());
159 :
160 316 : _func_params.resize(1);
161 :
162 316 : if (!_is_closed_loop && _forced_closing_num_segments > 0)
163 4 : paramError("forced_closing_num_segments",
164 : "this parameter is not needed if the curve to be generated is not a closed loop.");
165 :
166 312 : if (_node_set_boundaries.size() > 1 && _is_closed_loop)
167 0 : paramError("edge_nodesets", "Cannot add more than one edge nodeset on a closed loop");
168 312 : if (_node_set_boundaries.size() && _node_set_boundaries.size() != 2 && !_is_closed_loop)
169 0 : paramError("edge_nodesets", "Exactly two edges in an open loop");
170 312 : }
171 :
172 : std::unique_ptr<MeshBase>
173 296 : ParsedCurveGenerator::generate()
174 : {
175 296 : auto mesh = buildReplicatedMesh(2);
176 :
177 : // Do oversampling for each section of the curve as defined by "section_bounding_t_values"
178 770 : for (unsigned int i = 0; i < _nums_segments.size(); i++)
179 : {
180 482 : std::vector<Real> t_sect_space;
181 482 : std::vector<Real> dis_sect_space;
182 482 : tSectionSpaceDefiner(_section_bounding_t_values[i],
183 482 : _section_bounding_t_values[i + 1],
184 : t_sect_space,
185 : dis_sect_space,
186 482 : _nums_segments[i] * _order,
187 482 : _is_closed_loop && _nums_segments.size() == 1,
188 482 : _oversample_factor);
189 474 : if (i > 0)
190 : {
191 : // Remove the last t value and distance value of previous section to avert overlapping
192 186 : _t_space.pop_back();
193 : // Add the last distance value as this is a cumulative distance
194 36678 : for (auto & dis_elem : dis_sect_space)
195 36492 : dis_elem += _dis_space.back();
196 186 : _dis_space.pop_back();
197 : }
198 474 : _t_space.insert(_t_space.end(), t_sect_space.begin(), t_sect_space.end());
199 474 : _dis_space.insert(_dis_space.end(), dis_sect_space.begin(), dis_sect_space.end());
200 474 : }
201 :
202 : // Use linear interpolation to help achieve nearly uniform intervals for each section
203 288 : std::unique_ptr<LinearInterpolation> linear_t;
204 288 : linear_t = std::make_unique<LinearInterpolation>(_dis_space, _t_space);
205 :
206 : std::vector<Node *> nodes(
207 288 : std::accumulate(_nums_segments.begin(), _nums_segments.end(), 0) * _order + 1);
208 6146 : for (unsigned int i = 0; i < nodes.size(); i++)
209 : {
210 5858 : _func_params[0] = linear_t->sample((Real)i);
211 5858 : Point side_p = Point(evaluate(_func_Fx), evaluate(_func_Fy), evaluate(_func_Fz));
212 5858 : nodes[i] = mesh->add_point(side_p, i);
213 : }
214 :
215 : // For a closed loop, need to check if the first and last Points are overlapped
216 288 : if (_is_closed_loop)
217 : {
218 396 : if (MooseUtils::absoluteFuzzyEqual(
219 198 : (*nodes.back() - *nodes.front()).norm(), 0.0, _point_overlapping_tolerance))
220 : {
221 : // Remove the overlapped nodes for a closed loop
222 158 : mesh->delete_node(nodes.back());
223 158 : nodes.resize(nodes.size() - 1);
224 158 : if (_forced_closing_num_segments > 0)
225 4 : paramError("forced_closing_num_segments",
226 : "this parameter is not needed if the first and last points of the curve to be "
227 : "generated are overlapped.");
228 : }
229 40 : else if (_forced_closing_num_segments * _order > 1)
230 : {
231 : // Add extra nodes on the curve section used to forcefully close the loop
232 30 : const Point start_pt(*nodes.back());
233 30 : const Point end_pt(*nodes.front());
234 30 : const unsigned int num_nodes_0(nodes.size());
235 300 : for (unsigned int i = 1; i < _forced_closing_num_segments * _order; i++)
236 : {
237 270 : Point side_p = start_pt + (end_pt - start_pt) * (Real)i /
238 540 : (Real)(_forced_closing_num_segments * _order);
239 270 : nodes.push_back(mesh->add_point(side_p, num_nodes_0 + i - 1));
240 : }
241 : }
242 : }
243 :
244 284 : const unsigned int num_elems = (nodes.size() - !_is_closed_loop) / _order;
245 4794 : for (unsigned int i = 0; i < num_elems; i++)
246 : {
247 4510 : std::unique_ptr<Elem> new_elem;
248 4510 : new_elem = std::make_unique<Edge2>();
249 4510 : if (_order > 1)
250 : {
251 860 : new_elem = std::make_unique<Edge3>();
252 860 : if (_order == 3)
253 : {
254 430 : new_elem = std::make_unique<Edge4>();
255 430 : new_elem->set_node(3, nodes[i * _order + 2]);
256 : }
257 860 : new_elem->set_node(2, nodes[i * _order + 1]);
258 : }
259 4510 : new_elem->set_node(0, nodes[i * _order]);
260 4510 : new_elem->set_node(1, nodes[((i + 1) * _order) % nodes.size()]);
261 :
262 4510 : new_elem->subdomain_id() = 1;
263 4510 : mesh->add_elem(std::move(new_elem));
264 4510 : }
265 :
266 284 : if (_node_set_boundaries.size())
267 : {
268 : // Add boundary nodesets to boundary info
269 10 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
270 10 : int i = 0;
271 30 : for (auto & side_name : _node_set_boundaries)
272 20 : boundary_info.nodeset_name(i++) = side_name;
273 :
274 10 : boundary_info.add_node(*nodes.begin(), boundary_info.get_id_by_name(_node_set_boundaries[0]));
275 10 : if (_node_set_boundaries.size() > 1)
276 10 : boundary_info.add_node(*(nodes.end() - 1),
277 10 : boundary_info.get_id_by_name(_node_set_boundaries[1]));
278 : }
279 :
280 568 : return dynamic_pointer_cast<MeshBase>(mesh);
281 284 : }
282 :
283 : void
284 482 : ParsedCurveGenerator::tSectionSpaceDefiner(const Real t_start,
285 : const Real t_end,
286 : std::vector<Real> & t_sect_space,
287 : std::vector<Real> & dis_sect_space,
288 : unsigned int num_segments,
289 : const bool is_closed_loop,
290 : const Real oversample_factor)
291 : {
292 482 : std::vector<Point> pt_sect_space;
293 482 : t_sect_space.push_back(t_start);
294 482 : t_sect_space.push_back(t_end);
295 482 : pt_sect_space.push_back(pointCalculator(t_start));
296 482 : pt_sect_space.push_back(pointCalculator(t_end));
297 : const Real total_distance =
298 : is_closed_loop
299 482 : ? euclideanDistance(pt_sect_space.front(), pointCalculator((t_start + t_end) / 2.0))
300 482 : : euclideanDistance(pt_sect_space.front(), pt_sect_space.back());
301 482 : if (MooseUtils::absoluteFuzzyEqual(total_distance, 0.0))
302 4 : paramError("section_bounding_t_values",
303 : "The curve has at least one cross, which is not supported.");
304 478 : const Real target_fine_interval = total_distance / (Real)num_segments / oversample_factor;
305 478 : unsigned int t_space_id = 0;
306 478 : dis_sect_space.push_back(0.0);
307 : // Use a maximum oversampling points number and a counter to avert infinite loop
308 478 : const unsigned int max_binary_search = _max_oversample_number_factor * (num_segments + 1);
309 478 : unsigned int binary_search_counter(0);
310 : // Use binary algorithm to make all the intervals smaller than the target value
311 317468 : while (t_space_id + 1 < t_sect_space.size() && binary_search_counter <= max_binary_search)
312 : {
313 316990 : binary_search_counter++;
314 316990 : if (euclideanDistance(pt_sect_space[t_space_id], pt_sect_space[t_space_id + 1]) >
315 : target_fine_interval)
316 : {
317 158276 : const Real new_t_value = (t_sect_space[t_space_id] + t_sect_space[t_space_id + 1]) / 2.0;
318 158276 : t_sect_space.insert(t_sect_space.begin() + t_space_id + 1, new_t_value);
319 158276 : pt_sect_space.insert(pt_sect_space.begin() + t_space_id + 1, pointCalculator(new_t_value));
320 : }
321 : else
322 : {
323 158714 : dis_sect_space.push_back(
324 158714 : dis_sect_space.back() +
325 158714 : euclideanDistance(pt_sect_space[t_space_id], pt_sect_space[t_space_id + 1]));
326 158714 : t_space_id++;
327 : }
328 : }
329 478 : if (binary_search_counter > max_binary_search)
330 4 : paramError(
331 : "max_oversample_number_factor",
332 : "Maximum oversampling points number has been exceeded. Please consider adding more t "
333 : "values into 'section_bounding_t_values' or increase 'max_oversample_number_factor'.");
334 474 : const Real total_dis_seg = dis_sect_space.back();
335 : // Normalization to make the normalized length of the curve equals to number of segments
336 109678 : for (auto & dist : dis_sect_space)
337 109204 : dist = dist / total_dis_seg * (Real)num_segments;
338 474 : }
339 :
340 : Point
341 159346 : ParsedCurveGenerator::pointCalculator(const Real t_param)
342 : {
343 159346 : _func_params[0] = t_param;
344 159346 : return Point(evaluate(_func_Fx), evaluate(_func_Fy), evaluate(_func_Fz));
345 : }
346 :
347 : Real
348 476186 : ParsedCurveGenerator::euclideanDistance(const Point p1, const Point p2)
349 : {
350 476186 : return (p1 - p2).norm();
351 : }
|
