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 14965 : ParsedCurveGenerator::validParams()
28 : {
29 14965 : InputParameters params = MeshGenerator::validParams();
30 14965 : params += FunctionParserUtils<false>::validParams();
31 :
32 14965 : MooseEnum edge_elem_type("EDGE2 EDGE3 EDGE4", "EDGE2");
33 :
34 14965 : params.addRequiredParam<std::string>("x_formula", "Function expression of x(t)");
35 14965 : params.addRequiredParam<std::string>("y_formula", "Function expression of y(t)");
36 14965 : params.addParam<std::string>("z_formula", "0", "Function expression of z(t)");
37 14965 : 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 14965 : 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 14965 : params.addParam<std::vector<std::string>>(
49 : "constant_names", "Vector of constants used in the parsed function (use this for kB etc.)");
50 14965 : 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 14965 : params.addParam<bool>("is_closed_loop", false, "Whether the curve is closed or not.");
54 14965 : 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 14965 : 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 44895 : params.addRangeCheckedParam<Real>("oversample_factor",
65 29930 : 10.0,
66 : "oversample_factor>2.0",
67 : "Oversample factor to help make node distance nearly uniform.");
68 44895 : params.addRangeCheckedParam<unsigned int>(
69 : "max_oversample_number_factor",
70 29930 : 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 44895 : params.addParam<std::vector<BoundaryName>>(
75 29930 : "edge_nodesets", std::vector<BoundaryName>(), "Nodesets on both edges of the parsed curves");
76 14965 : params.addParam<MooseEnum>(
77 : "edge_element_type", edge_elem_type, "Type of the EDGE elements to be generated.");
78 14965 : 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 29930 : return params;
82 14965 : }
83 :
84 364 : ParsedCurveGenerator::ParsedCurveGenerator(const InputParameters & parameters)
85 : : MeshGenerator(parameters),
86 : FunctionParserUtils<false>(parameters),
87 364 : _function_x(getParam<std::string>("x_formula")),
88 364 : _function_y(getParam<std::string>("y_formula")),
89 364 : _function_z(getParam<std::string>("z_formula")),
90 364 : _nums_segments(getParam<std::vector<unsigned int>>("nums_segments")),
91 364 : _section_bounding_t_values(getParam<std::vector<Real>>("section_bounding_t_values")),
92 364 : _is_closed_loop(getParam<bool>("is_closed_loop")),
93 364 : _point_overlapping_tolerance(getParam<Real>("point_overlapping_tolerance")),
94 728 : _forced_closing_num_segments(isParamValid("forced_closing_num_segments")
95 364 : ? getParam<unsigned int>("forced_closing_num_segments")
96 : : 0),
97 364 : _oversample_factor(getParam<Real>("oversample_factor")),
98 364 : _max_oversample_number_factor(getParam<unsigned int>("max_oversample_number_factor")),
99 364 : _node_set_boundaries(getParam<std::vector<BoundaryName>>("edge_nodesets")),
100 728 : _order(getParam<MooseEnum>("edge_element_type") == "EDGE2"
101 430 : ? 1
102 1158 : : (getParam<MooseEnum>("edge_element_type") == "EDGE3" ? 2 : 3))
103 : {
104 364 : if (std::adjacent_find(_section_bounding_t_values.begin(), _section_bounding_t_values.end()) !=
105 728 : _section_bounding_t_values.end())
106 4 : paramError("section_bounding_t_values", "elements must be unique.");
107 364 : if (!std::is_sorted(_section_bounding_t_values.begin(), _section_bounding_t_values.end()) &&
108 364 : !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 356 : 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 352 : _func_Fx = std::make_shared<SymFunction>();
115 352 : _func_Fy = std::make_shared<SymFunction>();
116 352 : _func_Fz = std::make_shared<SymFunction>();
117 :
118 : // set FParser internal feature flags
119 352 : setParserFeatureFlags(_func_Fx);
120 352 : setParserFeatureFlags(_func_Fy);
121 352 : setParserFeatureFlags(_func_Fz);
122 :
123 : // add the constant expressions; note that the three functions share one set of constants
124 352 : if (isParamValid("constant_names") && isParamValid("constant_expressions"))
125 : {
126 249 : addFParserConstants(_func_Fx,
127 : getParam<std::vector<std::string>>("constant_names"),
128 : getParam<std::vector<std::string>>("constant_expressions"));
129 249 : addFParserConstants(_func_Fy,
130 : getParam<std::vector<std::string>>("constant_names"),
131 : getParam<std::vector<std::string>>("constant_expressions"));
132 249 : 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 352 : 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 348 : 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 344 : 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 340 : _func_params.resize(1);
161 :
162 340 : 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 336 : 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 336 : 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 336 : }
171 :
172 : std::unique_ptr<MeshBase>
173 320 : ParsedCurveGenerator::generate()
174 : {
175 320 : auto mesh = buildReplicatedMesh(2);
176 :
177 : // Do oversampling for each section of the curve as defined by "section_bounding_t_values"
178 833 : for (unsigned int i = 0; i < _nums_segments.size(); i++)
179 : {
180 521 : std::vector<Real> t_sect_space;
181 521 : std::vector<Real> dis_sect_space;
182 521 : tSectionSpaceDefiner(_section_bounding_t_values[i],
183 521 : _section_bounding_t_values[i + 1],
184 : t_sect_space,
185 : dis_sect_space,
186 521 : _nums_segments[i] * _order,
187 521 : _is_closed_loop && _nums_segments.size() == 1,
188 521 : _oversample_factor);
189 513 : if (i > 0)
190 : {
191 : // Remove the last t value and distance value of previous section to avert overlapping
192 201 : _t_space.pop_back();
193 : // Add the last distance value as this is a cumulative distance
194 39929 : for (auto & dis_elem : dis_sect_space)
195 39728 : dis_elem += _dis_space.back();
196 201 : _dis_space.pop_back();
197 : }
198 513 : _t_space.insert(_t_space.end(), t_sect_space.begin(), t_sect_space.end());
199 513 : _dis_space.insert(_dis_space.end(), dis_sect_space.begin(), dis_sect_space.end());
200 513 : }
201 :
202 : // Use linear interpolation to help achieve nearly uniform intervals for each section
203 312 : std::unique_ptr<LinearInterpolation> linear_t;
204 312 : linear_t = std::make_unique<LinearInterpolation>(_dis_space, _t_space);
205 :
206 : std::vector<Node *> nodes(
207 312 : std::accumulate(_nums_segments.begin(), _nums_segments.end(), 0) * _order + 1);
208 6667 : for (unsigned int i = 0; i < nodes.size(); i++)
209 : {
210 6355 : _func_params[0] = linear_t->sample((Real)i);
211 6355 : Point side_p = Point(evaluate(_func_Fx), evaluate(_func_Fy), evaluate(_func_Fz));
212 6355 : 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 312 : if (_is_closed_loop)
217 : {
218 426 : if (MooseUtils::absoluteFuzzyEqual(
219 213 : (*nodes.back() - *nodes.front()).norm(), 0.0, _point_overlapping_tolerance))
220 : {
221 : // Remove the overlapped nodes for a closed loop
222 169 : mesh->delete_node(nodes.back());
223 169 : nodes.resize(nodes.size() - 1);
224 169 : 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 44 : else if (_forced_closing_num_segments * _order > 1)
230 : {
231 : // Add extra nodes on the curve section used to forcefully close the loop
232 33 : const Point start_pt(*nodes.back());
233 33 : const Point end_pt(*nodes.front());
234 33 : const unsigned int num_nodes_0(nodes.size());
235 330 : for (unsigned int i = 1; i < _forced_closing_num_segments * _order; i++)
236 : {
237 297 : Point side_p = start_pt + (end_pt - start_pt) * (Real)i /
238 594 : (Real)(_forced_closing_num_segments * _order);
239 297 : nodes.push_back(mesh->add_point(side_p, num_nodes_0 + i - 1));
240 : }
241 : }
242 : }
243 :
244 308 : const unsigned int num_elems = (nodes.size() - !_is_closed_loop) / _order;
245 5193 : for (unsigned int i = 0; i < num_elems; i++)
246 : {
247 4885 : std::unique_ptr<Elem> new_elem;
248 4885 : new_elem = std::make_unique<Edge2>();
249 4885 : if (_order > 1)
250 : {
251 946 : new_elem = std::make_unique<Edge3>();
252 946 : if (_order == 3)
253 : {
254 473 : new_elem = std::make_unique<Edge4>();
255 473 : new_elem->set_node(3, nodes[i * _order + 2]);
256 : }
257 946 : new_elem->set_node(2, nodes[i * _order + 1]);
258 : }
259 4885 : new_elem->set_node(0, nodes[i * _order]);
260 4885 : new_elem->set_node(1, nodes[((i + 1) * _order) % nodes.size()]);
261 :
262 4885 : new_elem->subdomain_id() = 1;
263 4885 : mesh->add_elem(std::move(new_elem));
264 4885 : }
265 :
266 308 : if (_node_set_boundaries.size())
267 : {
268 : // Add boundary nodesets to boundary info
269 11 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
270 11 : int i = 0;
271 33 : for (auto & side_name : _node_set_boundaries)
272 22 : boundary_info.nodeset_name(i++) = side_name;
273 :
274 11 : boundary_info.add_node(*nodes.begin(), boundary_info.get_id_by_name(_node_set_boundaries[0]));
275 11 : if (_node_set_boundaries.size() > 1)
276 11 : boundary_info.add_node(*(nodes.end() - 1),
277 11 : boundary_info.get_id_by_name(_node_set_boundaries[1]));
278 : }
279 :
280 616 : return dynamic_pointer_cast<MeshBase>(mesh);
281 308 : }
282 :
283 : void
284 521 : 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 521 : std::vector<Point> pt_sect_space;
293 521 : t_sect_space.push_back(t_start);
294 521 : t_sect_space.push_back(t_end);
295 521 : pt_sect_space.push_back(pointCalculator(t_start));
296 521 : pt_sect_space.push_back(pointCalculator(t_end));
297 : const Real total_distance =
298 : is_closed_loop
299 521 : ? euclideanDistance(pt_sect_space.front(), pointCalculator((t_start + t_end) / 2.0))
300 521 : : euclideanDistance(pt_sect_space.front(), pt_sect_space.back());
301 521 : 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 517 : const Real target_fine_interval = total_distance / (Real)num_segments / oversample_factor;
305 517 : unsigned int t_space_id = 0;
306 517 : dis_sect_space.push_back(0.0);
307 : // Use a maximum oversampling points number and a counter to avert infinite loop
308 517 : const unsigned int max_binary_search = _max_oversample_number_factor * (num_segments + 1);
309 517 : unsigned int binary_search_counter(0);
310 : // Use binary algorithm to make all the intervals smaller than the target value
311 337978 : while (t_space_id + 1 < t_sect_space.size() && binary_search_counter <= max_binary_search)
312 : {
313 337461 : binary_search_counter++;
314 337461 : if (euclideanDistance(pt_sect_space[t_space_id], pt_sect_space[t_space_id + 1]) >
315 : target_fine_interval)
316 : {
317 168492 : const Real new_t_value = (t_sect_space[t_space_id] + t_sect_space[t_space_id + 1]) / 2.0;
318 168492 : t_sect_space.insert(t_sect_space.begin() + t_space_id + 1, new_t_value);
319 168492 : pt_sect_space.insert(pt_sect_space.begin() + t_space_id + 1, pointCalculator(new_t_value));
320 : }
321 : else
322 : {
323 168969 : dis_sect_space.push_back(
324 168969 : dis_sect_space.back() +
325 168969 : euclideanDistance(pt_sect_space[t_space_id], pt_sect_space[t_space_id + 1]));
326 168969 : t_space_id++;
327 : }
328 : }
329 517 : 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 513 : 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 120011 : for (auto & dist : dis_sect_space)
337 119498 : dist = dist / total_dis_seg * (Real)num_segments;
338 513 : }
339 :
340 : Point
341 169647 : ParsedCurveGenerator::pointCalculator(const Real t_param)
342 : {
343 169647 : _func_params[0] = t_param;
344 169647 : return Point(evaluate(_func_Fx), evaluate(_func_Fy), evaluate(_func_Fz));
345 : }
346 :
347 : Real
348 506951 : ParsedCurveGenerator::euclideanDistance(const Point p1, const Point p2)
349 : {
350 506951 : return (p1 - p2).norm();
351 : }
|
