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 "TransfiniteMeshGenerator.h"
11 : #include "CastUniquePointer.h"
12 : #include "Conversion.h"
13 : #include "MooseMeshUtils.h"
14 :
15 : #include "libmesh/parsed_function.h"
16 : #include "libmesh/replicated_mesh.h"
17 : #include "libmesh/face_quad4.h"
18 : #include "libmesh/fparser_ad.hh"
19 : #include "libmesh/elem.h"
20 :
21 : registerMooseObject("MooseApp", TransfiniteMeshGenerator);
22 :
23 : InputParameters
24 14425 : TransfiniteMeshGenerator::validParams()
25 : {
26 14425 : InputParameters params = MeshGenerator::validParams();
27 14425 : params += FunctionParserUtils<false>::validParams();
28 :
29 14425 : MooseEnum edge_type("LINE=1 CIRCARC=2 DISCRETE=3 PARSED=4");
30 :
31 14425 : params.addRequiredParam<std::vector<Point>>("corners", "The x,y,z positions of the nodes");
32 :
33 : // Define edge types
34 14425 : params.addRequiredParam<MooseEnum>("bottom_type", edge_type, "type of the bottom (y) boundary");
35 14425 : params.addRequiredParam<MooseEnum>("top_type", edge_type, "type of the top (y) boundary");
36 14425 : params.addRequiredParam<MooseEnum>("left_type", edge_type, "type of the left (x) boundary");
37 14425 : params.addRequiredParam<MooseEnum>("right_type", edge_type, "type of the right (x) boundary");
38 :
39 : // We need to know the number of points on opposite sides, and if no other parameter is available
40 : // we shall assume them to be equally distributed
41 14425 : params.addRequiredParam<unsigned int>("nx",
42 : "Number of nodes on horizontal edges, including corners");
43 14425 : params.addRequiredParam<unsigned int>("ny",
44 : "Number of Nodes on vertical edges, including corners");
45 :
46 : // each edge has a different parameter according to its type
47 14425 : params.addParam<std::string>("bottom_parameter", "", "Bottom side support parameter");
48 14425 : params.addParam<std::string>("top_parameter", "", "Top side support parameter");
49 14425 : params.addParam<std::string>("left_parameter", "", "Left side support parameter");
50 14425 : params.addParam<std::string>("right_parameter", "", "Right side support parameter");
51 :
52 43275 : params.addRangeCheckedParam<Real>(
53 : "bias_x",
54 28850 : 1.,
55 : "bias_x>=1.0 & bias_x<=2",
56 : "The amount by which to grow (or shrink) the cells in the x-direction.");
57 43275 : params.addRangeCheckedParam<Real>(
58 : "bias_y",
59 28850 : 1.,
60 : "bias_y>=1.0 & bias_y<=2",
61 : "The amount by which to grow (or shrink) the cells in the y-direction.");
62 :
63 14425 : params.addClassDescription(
64 : "Creates a QUAD4 mesh given a set of corner vertices and edge types. "
65 : "The edge type can be either LINE, CIRCARC, DISCRETE or PARSED, with "
66 : "LINE as the default option. "
67 : "For the non-default options the user needs to specify additional "
68 : "parameters via the edge_parameter option "
69 : "as follows: for CIRCARC the deviation of the midpoint from an "
70 : "arccircle, for DISCRETE a set of points, or "
71 : "a paramterization via the PARSED option. Opposite edges may have "
72 : "different distributions s long as the "
73 : "number of points is identical. Along opposite edges a different point "
74 : "distribution can be prescribed "
75 : "via the options bias_x or bias_y for opposing edges.");
76 :
77 14425 : params.addParamNamesToGroup("bottom_type left_type top_type right_type", "Edge type");
78 14425 : params.addParamNamesToGroup("bottom_parameter left_parameter top_parameter right_parameter",
79 : "Edge");
80 14425 : params.addParamNamesToGroup("nx ny bias_x bias_y", "Number and distribution of points");
81 :
82 28850 : return params;
83 14425 : }
84 :
85 80 : TransfiniteMeshGenerator::TransfiniteMeshGenerator(const InputParameters & parameters)
86 : : MeshGenerator(parameters),
87 : FunctionParserUtils<false>(parameters),
88 80 : _corners(getParam<std::vector<Point>>("corners")),
89 80 : _nx(getParam<unsigned int>("nx")),
90 80 : _ny(getParam<unsigned int>("ny")),
91 80 : _bottom_type(getParam<MooseEnum>("bottom_type")),
92 80 : _top_type(getParam<MooseEnum>("top_type")),
93 80 : _left_type(getParam<MooseEnum>("left_type")),
94 80 : _right_type(getParam<MooseEnum>("right_type")),
95 80 : _bottom_parameter(getParam<std::string>("bottom_parameter")),
96 80 : _top_parameter(getParam<std::string>("top_parameter")),
97 80 : _left_parameter(getParam<std::string>("left_parameter")),
98 80 : _right_parameter(getParam<std::string>("right_parameter")),
99 80 : _bias_x(getParam<Real>("bias_x")),
100 160 : _bias_y(getParam<Real>("bias_y"))
101 : {
102 : // initialize parsed function
103 80 : _parsed_func = std::make_shared<SymFunction>();
104 80 : setParserFeatureFlags(_parsed_func);
105 80 : _parsed_func->AddConstant("pi", libMesh::pi);
106 80 : _func_params.resize(1);
107 :
108 : mooseAssert((_nx > 1) && (_ny > 1),
109 : "A minimum of 2 points is needed on each edge, i.e. the user needs to consider edge "
110 : "vertices as well.");
111 80 : }
112 : std::unique_ptr<MeshBase>
113 80 : TransfiniteMeshGenerator::generate()
114 : {
115 80 : auto mesh = buildMeshBaseObject();
116 :
117 80 : mesh->set_mesh_dimension(2);
118 80 : mesh->set_spatial_dimension(2);
119 80 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
120 :
121 : // explicitly assign corners since they will be used extensively and the ordering may be confusing
122 80 : const Point V00 = _corners[0];
123 80 : const Point V10 = _corners[1];
124 80 : const Point V11 = _corners[2];
125 80 : const Point V01 = _corners[3];
126 :
127 : // we construct a vector that mimics the normal only to account for inward and outward arcircles
128 : // note we do not need to define inward directions since they would multiply the sign of the
129 : // arc circle parameter, which if negative would give inward vectors
130 80 : std::vector<Point> outward_vec(4);
131 80 : outward_vec[0] = Point(0.0, -1.0, 0.0);
132 80 : outward_vec[1] = Point(0.0, 1.0, 0.0);
133 80 : outward_vec[2] = Point(-1.0, 0.0, 0.0);
134 80 : outward_vec[3] = Point(1.0, 0.0, 0.0);
135 :
136 80 : const unsigned long int total_nodes = _nx * _ny;
137 :
138 : // we take [0,1] as the reference interval and we need to set the biases upfront
139 80 : Real edge_length = 1.0;
140 80 : std::vector<Real> param_x_dir = getPointsDistribution(edge_length, _nx, _bias_x);
141 80 : std::vector<Real> param_y_dir = getPointsDistribution(edge_length, _ny, _bias_y);
142 :
143 : std::vector<Point> edge_bottom =
144 80 : getEdge(V00, V10, _nx, _bottom_type, _bottom_parameter, outward_vec[0], param_x_dir);
145 : std::vector<Point> edge_top =
146 80 : getEdge(V01, V11, _nx, _top_type, _top_parameter, outward_vec[1], param_x_dir);
147 : std::vector<Point> edge_left =
148 80 : getEdge(V00, V01, _ny, _left_type, _left_parameter, outward_vec[2], param_y_dir);
149 : std::vector<Point> edge_right =
150 80 : getEdge(V10, V11, _ny, _right_type, _right_parameter, outward_vec[3], param_y_dir);
151 :
152 : // Used for the parametrization on edge pairs, provided by the point distribution according to
153 : // biases
154 : Real rx_coord, sy_coord;
155 :
156 80 : std::vector<Node *> nodes(total_nodes);
157 80 : unsigned int node_id = 0;
158 :
159 80 : Point newPt;
160 : Real r1_basis, r2_basis, s1_basis, s2_basis;
161 :
162 870 : for (unsigned int idx = 0; idx < _nx; idx++)
163 : {
164 790 : rx_coord = param_x_dir[idx];
165 790 : r1_basis = 1 - rx_coord;
166 790 : r2_basis = rx_coord;
167 :
168 5850 : for (unsigned int idy = 0; idy < _ny; idy++)
169 : {
170 5060 : sy_coord = param_y_dir[idy];
171 5060 : s1_basis = 1 - sy_coord;
172 5060 : s2_basis = sy_coord;
173 :
174 : // this is the core of the algorithm and generates every internal point
175 10120 : newPt = r2_basis * edge_right[idy] + r1_basis * edge_left[idy] + s1_basis * edge_bottom[idx] +
176 10120 : s2_basis * edge_top[idx] - r1_basis * s1_basis * V00 - r1_basis * s2_basis * V01 -
177 10120 : r2_basis * s1_basis * V10 - r2_basis * s2_basis * V11;
178 :
179 5060 : nodes[node_id] = mesh->add_point(newPt, node_id);
180 5060 : node_id++;
181 : }
182 : }
183 :
184 790 : for (unsigned int idx = 0; idx < _nx - 1; idx++)
185 : {
186 4530 : for (unsigned int idy = 0; idy < _ny - 1; idy++)
187 : {
188 3820 : Elem * elem = mesh->add_elem(new Quad4);
189 3820 : elem->set_node(0, nodes[idy + idx * _ny]);
190 3820 : elem->set_node(1, nodes[idy + (idx + 1) * _ny]);
191 3820 : elem->set_node(2, nodes[idy + 1 + (idx + 1) * _ny]);
192 3820 : elem->set_node(3, nodes[idy + 1 + idx * _ny]);
193 :
194 3820 : if (idy == 0) // add bottom boundary (boundary_id = 0)
195 710 : boundary_info.add_side(elem, 0, 0);
196 :
197 3820 : if (idx == _nx - 2) // add right boundary (boundary_id = 1)
198 450 : boundary_info.add_side(elem, 1, 1);
199 :
200 3820 : if (idy == _ny - 2) // add top boundary (boundary_id = 2)
201 710 : boundary_info.add_side(elem, 2, 2);
202 :
203 3820 : if (idx == 0) // add left boundary (boundary_id = 3)
204 450 : boundary_info.add_side(elem, 3, 3);
205 : }
206 : }
207 :
208 80 : boundary_info.sideset_name(0) = "bottom";
209 80 : boundary_info.nodeset_name(0) = "bottom";
210 :
211 80 : boundary_info.sideset_name(1) = "right";
212 80 : boundary_info.nodeset_name(1) = "right";
213 :
214 80 : boundary_info.sideset_name(2) = "top";
215 80 : boundary_info.nodeset_name(2) = "top";
216 :
217 80 : boundary_info.sideset_name(3) = "left";
218 80 : boundary_info.nodeset_name(3) = "left";
219 :
220 80 : mesh->prepare_for_use();
221 160 : return dynamic_pointer_cast<MeshBase>(mesh);
222 80 : }
223 :
224 : std::vector<Point>
225 320 : TransfiniteMeshGenerator::getEdge(const Point & P1,
226 : const Point & P2,
227 : const unsigned int np,
228 : const MooseEnum & type,
229 : const std::string & parameter,
230 : const Point & outward,
231 : const std::vector<Real> & param_vec)
232 : {
233 320 : std::vector<Point> edge;
234 :
235 320 : switch (type)
236 : {
237 200 : case 1:
238 : {
239 200 : edge = getLineEdge(P1, P2, param_vec);
240 : mooseAssert(MooseUtils::absoluteFuzzyEqual((edge[0] - P1).norm(), 0.0),
241 : "The line does not fit the first vertex on the edge.");
242 : mooseAssert(MooseUtils::absoluteFuzzyEqual((edge[np - 1] - P2).norm(), 0.0),
243 : "The line does not fit the end vertex on the edge.");
244 : }
245 200 : break;
246 50 : case 2:
247 : {
248 50 : edge = getCircarcEdge(P1, P2, parameter, outward, param_vec);
249 : mooseAssert(MooseUtils::absoluteFuzzyEqual((edge[0] - P1).norm(), 0.0),
250 : "No arccircle parametrization can be found to fit the first vertex on the edge.");
251 : mooseAssert(MooseUtils::absoluteFuzzyEqual((edge[np - 1] - P2).norm(), 0.0),
252 : "No arccircle parametrization can be found to fit the end vertex on the edge.");
253 : }
254 50 : break;
255 10 : case 3:
256 : {
257 10 : edge = getDiscreteEdge(np, parameter);
258 : mooseAssert(MooseUtils::absoluteFuzzyEqual((edge[0] - P1).norm(), 0.0),
259 : "The first discrete point does not fit the corresponding edge vertex."
260 : "Note: discrete points need to replicate the edge corners.");
261 : mooseAssert(MooseUtils::absoluteFuzzyEqual((edge[np - 1] - P2).norm(), 0.0),
262 : "The last discrete point does not fit the corresponding edge vertex."
263 : "Note: discrete points need to replicate the edge corners.");
264 : }
265 10 : break;
266 60 : case 4:
267 : {
268 60 : edge = getParsedEdge(parameter, param_vec);
269 : mooseAssert(MooseUtils::absoluteFuzzyEqual((edge[0] - P1).norm(), 0.0),
270 : "The parametrization does not fit the first vertex on the edge.");
271 : mooseAssert(MooseUtils::absoluteFuzzyEqual((edge[np - 1] - P2).norm(), 0.0),
272 : "The parametrization does not fit the end vertex on the edge.");
273 : }
274 60 : break;
275 : }
276 320 : if (edge.size() != np)
277 0 : mooseError("The generated edge does not match the number of points on the"
278 : "opposite edge.");
279 320 : return edge;
280 0 : }
281 :
282 : std::vector<Point>
283 200 : TransfiniteMeshGenerator::getLineEdge(const Point & P1,
284 : const Point & P2,
285 : const std::vector<Real> & param_vec)
286 : {
287 200 : std::vector<Point> edge(param_vec.size());
288 200 : auto it = 0;
289 :
290 1690 : for (auto rx : param_vec)
291 : {
292 1490 : edge[it] = P1 * (1.0 - rx) + P2 * rx;
293 1490 : it++;
294 : };
295 200 : return edge;
296 : }
297 :
298 : std::vector<Point>
299 60 : TransfiniteMeshGenerator::getParsedEdge(const std::string & parameter,
300 : const std::vector<Real> & param_vec)
301 : {
302 60 : std::vector<Point> edge(param_vec.size());
303 : Real x_coord, y_coord;
304 :
305 60 : std::vector<std::string> param_coords;
306 60 : MooseUtils::tokenize(parameter, param_coords, 1, ";");
307 :
308 60 : auto it = 0;
309 580 : for (auto rx : param_vec)
310 : {
311 520 : _parsed_func->Parse(param_coords[0], "r");
312 520 : x_coord = _parsed_func->Eval(&rx);
313 520 : _parsed_func->Parse(param_coords[1], "r");
314 520 : y_coord = _parsed_func->Eval(&rx);
315 520 : edge[it] = Point(x_coord, y_coord, 0.0);
316 520 : it++;
317 : }
318 :
319 120 : return edge;
320 60 : }
321 :
322 : std::vector<Point>
323 10 : TransfiniteMeshGenerator::getDiscreteEdge(const unsigned int np, const std::string & parameter)
324 : {
325 10 : std::vector<Point> edge(np);
326 10 : std::vector<std::string> string_points;
327 10 : MooseUtils::tokenize(parameter, string_points, 1, "\n");
328 10 : if (string_points.size() != np)
329 0 : mooseError("DISCRETE: the number of discrete points does not match the number of points on the"
330 : "opposite edge.");
331 :
332 10 : auto it = 0;
333 80 : for (unsigned int iter = 0; iter < string_points.size(); iter++)
334 : {
335 70 : std::vector<Real> point_vals;
336 70 : MooseUtils::tokenizeAndConvert(string_points[iter], point_vals, " ");
337 70 : edge[it] = Point(point_vals[0], point_vals[1], point_vals[2]);
338 70 : it++;
339 70 : }
340 :
341 20 : return edge;
342 10 : }
343 :
344 : std::vector<Point>
345 50 : TransfiniteMeshGenerator::getCircarcEdge(const Point & P1,
346 : const Point & P2,
347 : const std::string & parameter,
348 : const Point & outward,
349 : const std::vector<Real> & param_vec)
350 : {
351 50 : std::vector<Point> edge(param_vec.size());
352 :
353 50 : std::vector<Real> param_coords;
354 50 : MooseUtils::tokenizeAndConvert(parameter, param_coords, ";");
355 50 : Point P3;
356 50 : if (param_coords.size() == 1)
357 : {
358 40 : P3 = computeMidPoint(P1, P2, param_coords[0], outward);
359 : }
360 : else
361 : {
362 10 : P3 = Point(param_coords[0], param_coords[1], param_coords[2]);
363 : }
364 :
365 50 : Real rad = computeRadius(P1, P2, P3);
366 50 : Point P0 = computeOrigin(P1, P2, P3);
367 :
368 : // need to shift to center of coordinates to find the corresponding radians
369 50 : Point x0 = (P1 - P0);
370 50 : Point x1 = (P2 - P0);
371 :
372 : // The case when the edge spans quadrants 1 and 4 requires special treatment
373 : // to periodically switch we compute the angle that needs added to one edge
374 : // to identify the entire edge span
375 : mooseAssert(x0.norm() > 0.0 && x1.norm() > 0.0,
376 : "The point provided cannot generate an arc circle on the edge specified");
377 :
378 50 : Real arclength = std::acos((x0 * x1) / x0.norm() / x1.norm());
379 :
380 50 : Real a = std::atan2(x0(1), x0(0));
381 50 : Real b = std::atan2(x1(1), x1(0));
382 :
383 50 : if (a < 0)
384 20 : a = a + 2 * M_PI;
385 50 : if (std::abs(b - a) > M_PI)
386 20 : b = a + arclength;
387 :
388 50 : auto it = 0;
389 610 : for (auto rx : param_vec)
390 : {
391 560 : Real interval = getMapInterval(rx, 0.0, 1.0, a, b);
392 560 : Real x = P0(0) + rad * std::cos(interval);
393 560 : Real y = P0(1) + rad * std::sin(interval);
394 560 : edge[it] = Point(x, y, 0.0);
395 560 : it++;
396 : };
397 :
398 100 : return edge;
399 50 : }
400 :
401 : Real
402 1850 : TransfiniteMeshGenerator::getMapInterval(
403 : const Real xab, const Real a, const Real b, const Real c, const Real d) const
404 : // this routine maps a point x\in[a, b] to the corresponding point in the interval [c, d]
405 : {
406 : mooseAssert(std::abs(b - a) > 0.0,
407 : "The input interval [a, b] is empty, check that a and b are not identical.");
408 1850 : Real xcd = c + (d - c) / (b - a) * (xab - a);
409 :
410 1850 : return xcd;
411 : }
412 :
413 : std::vector<Real>
414 160 : TransfiniteMeshGenerator::getPointsDistribution(const Real edge_length,
415 : const unsigned int np,
416 : const Real bias) const
417 : {
418 160 : std::vector<Real> param_vec;
419 160 : Real rx = 0.0;
420 :
421 160 : if (bias != 1.0)
422 : {
423 30 : Real step = 0.0;
424 30 : Real factor = edge_length * (1.0 - std::abs(bias)) / (1.0 - std::pow(std::abs(bias), np - 1));
425 30 : param_vec.push_back(rx);
426 220 : for (unsigned int iter = 1; iter < np; iter++)
427 : {
428 190 : step = step + factor * std::pow(bias, Real(iter - 1));
429 190 : rx = getMapInterval(step, 0.0, edge_length, 0.0, 1.0);
430 190 : param_vec.push_back(rx);
431 : }
432 : }
433 : else
434 : {
435 130 : const Real interval = edge_length / Real(np - 1);
436 1230 : for (unsigned int iter = 0; iter < np; iter++)
437 : {
438 1100 : rx = getMapInterval(Real(iter) * interval, 0.0, edge_length, 0.0, 1.0);
439 1100 : param_vec.push_back(rx);
440 : }
441 : }
442 320 : return param_vec;
443 0 : }
444 :
445 : Real
446 50 : TransfiniteMeshGenerator::computeRadius(const Point & P1, const Point & P2, const Point & P3) const
447 : {
448 50 : Point temp1 = P1 - P2;
449 50 : Real a2 = temp1.norm_sq(); // a is the distance from P1 to P2, but we only need it squared
450 50 : Point temp2 = P3 - P1;
451 50 : Real b2 = temp2.norm_sq();
452 50 : Real dr = std::sqrt(b2 - a2 / 4.0);
453 50 : const Real rad = a2 / 8.0 / dr + dr / 2.0;
454 50 : return rad;
455 : }
456 :
457 : Point
458 50 : TransfiniteMeshGenerator::computeOrigin(const Point & P1, const Point & P2, const Point & P3) const
459 : {
460 : // define interim quantities
461 50 : const Real a1 = (2 * P3(0) - 2 * P1(0));
462 50 : const Real a2 = (2 * P3(1) - 2 * P1(1));
463 50 : const Real A = P3(1) * P3(1) - P1(0) * P1(0) + P3(0) * P3(0) - P1(1) * P1(1);
464 50 : const Real b1 = (2 * P3(0) - 2 * P2(0));
465 50 : const Real b2 = (2 * P3(1) - 2 * P2(1));
466 50 : const Real B = P3(1) * P3(1) - P2(0) * P2(0) + P3(0) * P3(0) - P2(1) * P2(1);
467 : mooseAssert(std::abs(b2) > 0.0 && std::abs(a1) > 0.0 && std::abs(a1 * b2 - a2 * b1) > 0.0,
468 : "The point provided cannot generate an arc circle on the edge specified."
469 : "The origin of the corresponding circle cannot be computed.");
470 :
471 50 : const Real y0 = (a1 * B - A * b1) / (a1 * b2 - a2 * b1);
472 50 : const Real x0 = (A - y0 * a2) / a1;
473 :
474 : // fill in and return the origin point
475 50 : Point P0 = Point(x0, y0, 0.0);
476 50 : return P0;
477 : }
478 :
479 : Point
480 40 : TransfiniteMeshGenerator::computeMidPoint(const Point & P1,
481 : const Point & P2,
482 : const Real height,
483 : const Point & outward) const
484 : {
485 40 : const Real xm = (P1(0) + P2(0)) / 2;
486 40 : const Real ym = (P1(1) + P2(1)) / 2;
487 : // The arc can be inverted into the domain (concave) or outward (convex)
488 : // we use the convention that if the height is given as negative we seek a concave arc
489 : // if "height" is positive then we seek a convex arc.
490 : // The outward vector is sufficient to determine the direction of the arc.
491 40 : const Real dist = std::abs(height);
492 40 : const int orient = (height >= 0) ? 1 : -1;
493 40 : Point MidPoint;
494 :
495 : // this accounts for lines aligned with the system of coordinates
496 40 : if ((std::abs(P2(0) - P1(0)) < 1e-15) || (std::abs(P2(1) - P1(1)) < 1e-15))
497 10 : MidPoint = Point(xm, ym, 0.0) + dist * orient * outward;
498 :
499 : // some of the cases are not covered by this strategy and we need to use normals info
500 : // this is implemented and tested in the prototype and is soon to be updated
501 40 : if (std::abs(P2(0) - P1(0)) > 1e-15 && std::abs(P2(1) - P1(1)) > 1e-15)
502 : {
503 : // m is the slope of the line orthogonal to the midpoint
504 30 : const Real m = -(P2(0) - P1(0)) / (P2(1) - P1(1));
505 : // The equation to determine allows two solutions
506 30 : const Real x_temp = dist * sqrt((m * m + 1)) / (m * m + 1);
507 30 : const Real factor = orient * outward(0) + m * orient * outward(1);
508 30 : int direction = (factor >= 0) ? 1 : -1;
509 :
510 30 : MidPoint = Point(direction * x_temp + xm, direction * x_temp * m + ym, 0.0);
511 : }
512 40 : return MidPoint;
513 : }
|
