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 14441 : TransfiniteMeshGenerator::validParams()
25 : {
26 14441 : InputParameters params = MeshGenerator::validParams();
27 14441 : params += FunctionParserUtils<false>::validParams();
28 :
29 14441 : MooseEnum edge_type("LINE=1 CIRCARC=2 DISCRETE=3 PARSED=4");
30 :
31 14441 : params.addRequiredParam<std::vector<Point>>("corners", "The x,y,z positions of the nodes");
32 :
33 : // Define edge types
34 14441 : params.addRequiredParam<MooseEnum>("bottom_type", edge_type, "type of the bottom (y) boundary");
35 14441 : params.addRequiredParam<MooseEnum>("top_type", edge_type, "type of the top (y) boundary");
36 14441 : params.addRequiredParam<MooseEnum>("left_type", edge_type, "type of the left (x) boundary");
37 14441 : 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 14441 : params.addRequiredParam<unsigned int>("nx",
42 : "Number of nodes on horizontal edges, including corners");
43 14441 : 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 14441 : params.addParam<std::string>("bottom_parameter", "", "Bottom side support parameter");
48 14441 : params.addParam<std::string>("top_parameter", "", "Top side support parameter");
49 14441 : params.addParam<std::string>("left_parameter", "", "Left side support parameter");
50 14441 : params.addParam<std::string>("right_parameter", "", "Right side support parameter");
51 :
52 43323 : params.addRangeCheckedParam<Real>(
53 : "bias_x",
54 28882 : 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 43323 : params.addRangeCheckedParam<Real>(
58 : "bias_y",
59 28882 : 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 14441 : 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 14441 : params.addParamNamesToGroup("bottom_type left_type top_type right_type", "Edge type");
78 14441 : params.addParamNamesToGroup("bottom_parameter left_parameter top_parameter right_parameter",
79 : "Edge");
80 14441 : params.addParamNamesToGroup("nx ny bias_x bias_y", "Number and distribution of points");
81 :
82 28882 : return params;
83 14441 : }
84 :
85 88 : TransfiniteMeshGenerator::TransfiniteMeshGenerator(const InputParameters & parameters)
86 : : MeshGenerator(parameters),
87 : FunctionParserUtils<false>(parameters),
88 88 : _corners(getParam<std::vector<Point>>("corners")),
89 88 : _nx(getParam<unsigned int>("nx")),
90 88 : _ny(getParam<unsigned int>("ny")),
91 88 : _bottom_type(getParam<MooseEnum>("bottom_type")),
92 88 : _top_type(getParam<MooseEnum>("top_type")),
93 88 : _left_type(getParam<MooseEnum>("left_type")),
94 88 : _right_type(getParam<MooseEnum>("right_type")),
95 88 : _bottom_parameter(getParam<std::string>("bottom_parameter")),
96 88 : _top_parameter(getParam<std::string>("top_parameter")),
97 88 : _left_parameter(getParam<std::string>("left_parameter")),
98 88 : _right_parameter(getParam<std::string>("right_parameter")),
99 88 : _bias_x(getParam<Real>("bias_x")),
100 176 : _bias_y(getParam<Real>("bias_y"))
101 : {
102 : // initialize parsed function
103 88 : _parsed_func = std::make_shared<SymFunction>();
104 88 : setParserFeatureFlags(_parsed_func);
105 88 : _parsed_func->AddConstant("pi", libMesh::pi);
106 88 : _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 88 : }
112 : std::unique_ptr<MeshBase>
113 88 : TransfiniteMeshGenerator::generate()
114 : {
115 88 : auto mesh = buildMeshBaseObject();
116 :
117 88 : mesh->set_mesh_dimension(2);
118 88 : mesh->set_spatial_dimension(2);
119 88 : 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 88 : const Point V00 = _corners[0];
123 88 : const Point V10 = _corners[1];
124 88 : const Point V11 = _corners[2];
125 88 : 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 88 : std::vector<Point> outward_vec(4);
131 88 : outward_vec[0] = Point(0.0, -1.0, 0.0);
132 88 : outward_vec[1] = Point(0.0, 1.0, 0.0);
133 88 : outward_vec[2] = Point(-1.0, 0.0, 0.0);
134 88 : outward_vec[3] = Point(1.0, 0.0, 0.0);
135 :
136 88 : 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 88 : Real edge_length = 1.0;
140 88 : std::vector<Real> param_x_dir = getPointsDistribution(edge_length, _nx, _bias_x);
141 88 : std::vector<Real> param_y_dir = getPointsDistribution(edge_length, _ny, _bias_y);
142 :
143 : std::vector<Point> edge_bottom =
144 88 : getEdge(V00, V10, _nx, _bottom_type, _bottom_parameter, outward_vec[0], param_x_dir);
145 : std::vector<Point> edge_top =
146 88 : getEdge(V01, V11, _nx, _top_type, _top_parameter, outward_vec[1], param_x_dir);
147 : std::vector<Point> edge_left =
148 88 : getEdge(V00, V01, _ny, _left_type, _left_parameter, outward_vec[2], param_y_dir);
149 : std::vector<Point> edge_right =
150 88 : 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 88 : std::vector<Node *> nodes(total_nodes);
157 88 : unsigned int node_id = 0;
158 :
159 88 : Point newPt;
160 : Real r1_basis, r2_basis, s1_basis, s2_basis;
161 :
162 957 : for (unsigned int idx = 0; idx < _nx; idx++)
163 : {
164 869 : rx_coord = param_x_dir[idx];
165 869 : r1_basis = 1 - rx_coord;
166 869 : r2_basis = rx_coord;
167 :
168 6435 : for (unsigned int idy = 0; idy < _ny; idy++)
169 : {
170 5566 : sy_coord = param_y_dir[idy];
171 5566 : s1_basis = 1 - sy_coord;
172 5566 : s2_basis = sy_coord;
173 :
174 : // this is the core of the algorithm and generates every internal point
175 11132 : newPt = r2_basis * edge_right[idy] + r1_basis * edge_left[idy] + s1_basis * edge_bottom[idx] +
176 11132 : s2_basis * edge_top[idx] - r1_basis * s1_basis * V00 - r1_basis * s2_basis * V01 -
177 11132 : r2_basis * s1_basis * V10 - r2_basis * s2_basis * V11;
178 :
179 5566 : nodes[node_id] = mesh->add_point(newPt, node_id);
180 5566 : node_id++;
181 : }
182 : }
183 :
184 869 : for (unsigned int idx = 0; idx < _nx - 1; idx++)
185 : {
186 4983 : for (unsigned int idy = 0; idy < _ny - 1; idy++)
187 : {
188 4202 : Elem * elem = mesh->add_elem(new Quad4);
189 4202 : elem->set_node(0, nodes[idy + idx * _ny]);
190 4202 : elem->set_node(1, nodes[idy + (idx + 1) * _ny]);
191 4202 : elem->set_node(2, nodes[idy + 1 + (idx + 1) * _ny]);
192 4202 : elem->set_node(3, nodes[idy + 1 + idx * _ny]);
193 :
194 4202 : if (idy == 0) // add bottom boundary (boundary_id = 0)
195 781 : boundary_info.add_side(elem, 0, 0);
196 :
197 4202 : if (idx == _nx - 2) // add right boundary (boundary_id = 1)
198 495 : boundary_info.add_side(elem, 1, 1);
199 :
200 4202 : if (idy == _ny - 2) // add top boundary (boundary_id = 2)
201 781 : boundary_info.add_side(elem, 2, 2);
202 :
203 4202 : if (idx == 0) // add left boundary (boundary_id = 3)
204 495 : boundary_info.add_side(elem, 3, 3);
205 : }
206 : }
207 :
208 88 : boundary_info.sideset_name(0) = "bottom";
209 88 : boundary_info.nodeset_name(0) = "bottom";
210 :
211 88 : boundary_info.sideset_name(1) = "right";
212 88 : boundary_info.nodeset_name(1) = "right";
213 :
214 88 : boundary_info.sideset_name(2) = "top";
215 88 : boundary_info.nodeset_name(2) = "top";
216 :
217 88 : boundary_info.sideset_name(3) = "left";
218 88 : boundary_info.nodeset_name(3) = "left";
219 :
220 88 : mesh->prepare_for_use();
221 176 : return dynamic_pointer_cast<MeshBase>(mesh);
222 88 : }
223 :
224 : std::vector<Point>
225 352 : 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 352 : std::vector<Point> edge;
234 :
235 352 : switch (type)
236 : {
237 220 : case 1:
238 : {
239 220 : 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 220 : break;
246 55 : case 2:
247 : {
248 55 : 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 55 : break;
255 11 : case 3:
256 : {
257 11 : 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 11 : break;
266 66 : case 4:
267 : {
268 66 : 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 66 : break;
275 : }
276 352 : if (edge.size() != np)
277 0 : mooseError("The generated edge does not match the number of points on the"
278 : "opposite edge.");
279 352 : return edge;
280 0 : }
281 :
282 : std::vector<Point>
283 220 : TransfiniteMeshGenerator::getLineEdge(const Point & P1,
284 : const Point & P2,
285 : const std::vector<Real> & param_vec)
286 : {
287 220 : std::vector<Point> edge(param_vec.size());
288 220 : auto it = 0;
289 :
290 1859 : for (auto rx : param_vec)
291 : {
292 1639 : edge[it] = P1 * (1.0 - rx) + P2 * rx;
293 1639 : it++;
294 : };
295 220 : return edge;
296 : }
297 :
298 : std::vector<Point>
299 66 : TransfiniteMeshGenerator::getParsedEdge(const std::string & parameter,
300 : const std::vector<Real> & param_vec)
301 : {
302 66 : std::vector<Point> edge(param_vec.size());
303 : Real x_coord, y_coord;
304 :
305 66 : std::vector<std::string> param_coords;
306 66 : MooseUtils::tokenize(parameter, param_coords, 1, ";");
307 :
308 66 : auto it = 0;
309 638 : for (auto rx : param_vec)
310 : {
311 572 : _parsed_func->Parse(param_coords[0], "r");
312 572 : x_coord = _parsed_func->Eval(&rx);
313 572 : _parsed_func->Parse(param_coords[1], "r");
314 572 : y_coord = _parsed_func->Eval(&rx);
315 572 : edge[it] = Point(x_coord, y_coord, 0.0);
316 572 : it++;
317 : }
318 :
319 132 : return edge;
320 66 : }
321 :
322 : std::vector<Point>
323 11 : TransfiniteMeshGenerator::getDiscreteEdge(const unsigned int np, const std::string & parameter)
324 : {
325 11 : std::vector<Point> edge(np);
326 11 : std::vector<std::string> string_points;
327 11 : MooseUtils::tokenize(parameter, string_points, 1, "\n");
328 11 : 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 11 : auto it = 0;
333 88 : for (unsigned int iter = 0; iter < string_points.size(); iter++)
334 : {
335 77 : std::vector<Real> point_vals;
336 77 : MooseUtils::tokenizeAndConvert(string_points[iter], point_vals, " ");
337 77 : edge[it] = Point(point_vals[0], point_vals[1], point_vals[2]);
338 77 : it++;
339 77 : }
340 :
341 22 : return edge;
342 11 : }
343 :
344 : std::vector<Point>
345 55 : 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 55 : std::vector<Point> edge(param_vec.size());
352 :
353 55 : std::vector<Real> param_coords;
354 55 : MooseUtils::tokenizeAndConvert(parameter, param_coords, ";");
355 55 : Point P3;
356 55 : if (param_coords.size() == 1)
357 : {
358 44 : P3 = computeMidPoint(P1, P2, param_coords[0], outward);
359 : }
360 : else
361 : {
362 11 : P3 = Point(param_coords[0], param_coords[1], param_coords[2]);
363 : }
364 :
365 55 : Real rad = computeRadius(P1, P2, P3);
366 55 : Point P0 = computeOrigin(P1, P2, P3);
367 :
368 : // need to shift to center of coordinates to find the corresponding radians
369 55 : Point x0 = (P1 - P0);
370 55 : 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 55 : Real arclength = std::acos((x0 * x1) / x0.norm() / x1.norm());
379 :
380 55 : Real a = std::atan2(x0(1), x0(0));
381 55 : Real b = std::atan2(x1(1), x1(0));
382 :
383 55 : if (a < 0)
384 22 : a = a + 2 * M_PI;
385 55 : if (std::abs(b - a) > M_PI)
386 22 : b = a + arclength;
387 :
388 55 : auto it = 0;
389 671 : for (auto rx : param_vec)
390 : {
391 616 : Real interval = getMapInterval(rx, 0.0, 1.0, a, b);
392 616 : Real x = P0(0) + rad * std::cos(interval);
393 616 : Real y = P0(1) + rad * std::sin(interval);
394 616 : edge[it] = Point(x, y, 0.0);
395 616 : it++;
396 : };
397 :
398 110 : return edge;
399 55 : }
400 :
401 : Real
402 2035 : 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 2035 : Real xcd = c + (d - c) / (b - a) * (xab - a);
409 :
410 2035 : return xcd;
411 : }
412 :
413 : std::vector<Real>
414 176 : TransfiniteMeshGenerator::getPointsDistribution(const Real edge_length,
415 : const unsigned int np,
416 : const Real bias) const
417 : {
418 176 : std::vector<Real> param_vec;
419 176 : Real rx = 0.0;
420 :
421 176 : if (bias != 1.0)
422 : {
423 33 : Real step = 0.0;
424 33 : Real factor = edge_length * (1.0 - std::abs(bias)) / (1.0 - std::pow(std::abs(bias), np - 1));
425 33 : param_vec.push_back(rx);
426 242 : for (unsigned int iter = 1; iter < np; iter++)
427 : {
428 209 : step = step + factor * std::pow(bias, Real(iter - 1));
429 209 : rx = getMapInterval(step, 0.0, edge_length, 0.0, 1.0);
430 209 : param_vec.push_back(rx);
431 : }
432 : }
433 : else
434 : {
435 143 : const Real interval = edge_length / Real(np - 1);
436 1353 : for (unsigned int iter = 0; iter < np; iter++)
437 : {
438 1210 : rx = getMapInterval(Real(iter) * interval, 0.0, edge_length, 0.0, 1.0);
439 1210 : param_vec.push_back(rx);
440 : }
441 : }
442 352 : return param_vec;
443 0 : }
444 :
445 : Real
446 55 : TransfiniteMeshGenerator::computeRadius(const Point & P1, const Point & P2, const Point & P3) const
447 : {
448 55 : Point temp1 = P1 - P2;
449 55 : Real a2 = temp1.norm_sq(); // a is the distance from P1 to P2, but we only need it squared
450 55 : Point temp2 = P3 - P1;
451 55 : Real b2 = temp2.norm_sq();
452 55 : Real dr = std::sqrt(b2 - a2 / 4.0);
453 55 : const Real rad = a2 / 8.0 / dr + dr / 2.0;
454 55 : return rad;
455 : }
456 :
457 : Point
458 55 : TransfiniteMeshGenerator::computeOrigin(const Point & P1, const Point & P2, const Point & P3) const
459 : {
460 : // define interim quantities
461 55 : const Real a1 = (2 * P3(0) - 2 * P1(0));
462 55 : const Real a2 = (2 * P3(1) - 2 * P1(1));
463 55 : const Real A = P3(1) * P3(1) - P1(0) * P1(0) + P3(0) * P3(0) - P1(1) * P1(1);
464 55 : const Real b1 = (2 * P3(0) - 2 * P2(0));
465 55 : const Real b2 = (2 * P3(1) - 2 * P2(1));
466 55 : 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 55 : const Real y0 = (a1 * B - A * b1) / (a1 * b2 - a2 * b1);
472 55 : const Real x0 = (A - y0 * a2) / a1;
473 :
474 : // fill in and return the origin point
475 55 : Point P0 = Point(x0, y0, 0.0);
476 55 : return P0;
477 : }
478 :
479 : Point
480 44 : TransfiniteMeshGenerator::computeMidPoint(const Point & P1,
481 : const Point & P2,
482 : const Real height,
483 : const Point & outward) const
484 : {
485 44 : const Real xm = (P1(0) + P2(0)) / 2;
486 44 : 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 44 : const Real dist = std::abs(height);
492 44 : const int orient = (height >= 0) ? 1 : -1;
493 44 : Point MidPoint;
494 :
495 : // this accounts for lines aligned with the system of coordinates
496 44 : if ((std::abs(P2(0) - P1(0)) < 1e-15) || (std::abs(P2(1) - P1(1)) < 1e-15))
497 11 : 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 44 : 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 33 : const Real m = -(P2(0) - P1(0)) / (P2(1) - P1(1));
505 : // The equation to determine allows two solutions
506 33 : const Real x_temp = dist * sqrt((m * m + 1)) / (m * m + 1);
507 33 : const Real factor = orient * outward(0) + m * orient * outward(1);
508 33 : int direction = (factor >= 0) ? 1 : -1;
509 :
510 33 : MidPoint = Point(direction * x_temp + xm, direction * x_temp * m + ym, 0.0);
511 : }
512 44 : return MidPoint;
513 : }
|
