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 "AnnularMeshGenerator.h"
11 : #include "CastUniquePointer.h"
12 :
13 : #include "libmesh/replicated_mesh.h"
14 : #include "libmesh/face_quad4.h"
15 : #include "libmesh/face_tri3.h"
16 : #include "libmesh/mesh_modification.h"
17 :
18 : registerMooseObject("MooseApp", AnnularMeshGenerator);
19 :
20 : InputParameters
21 14415 : AnnularMeshGenerator::validParams()
22 : {
23 14415 : InputParameters params = MeshGenerator::validParams();
24 :
25 43245 : params.addRangeCheckedParam<unsigned int>(
26 28830 : "nr", 1, "nr>0", "Number of elements in the radial direction");
27 14415 : params.addRequiredRangeCheckedParam<unsigned int>(
28 : "nt", "nt>0", "Number of elements in the angular direction");
29 14415 : params.addRequiredRangeCheckedParam<Real>(
30 : "rmin",
31 : "rmin>=0.0",
32 : "Inner radius. If rmin=0 then a disc mesh (with no central hole) will be created.");
33 14415 : params.addRequiredParam<Real>("rmax", "Outer radius");
34 14415 : params.addParam<std::vector<Real>>(
35 : "radial_positions", {}, "Directly prescribed positions of intermediate radial nodes");
36 43245 : params.addParam<bool>("equal_area",
37 28830 : false,
38 : "Whether to select the radial discretization "
39 : "to achieve equal areas of each ring");
40 43245 : params.addDeprecatedParam<Real>("tmin",
41 28830 : 0.0,
42 : "Minimum angle, measured in radians anticlockwise from x axis",
43 : "Use dmin instead");
44 43245 : params.addDeprecatedParam<Real>(
45 : "tmax",
46 28830 : 2 * M_PI,
47 : "Maximum angle, measured in radians anticlockwise from x axis. If "
48 : "tmin=0 and tmax=2Pi an annular mesh is created. "
49 : "Otherwise, only a sector of an annulus is created",
50 : "Use dmin instead");
51 43245 : params.addParam<Real>(
52 28830 : "dmin", 0.0, "Minimum degree, measured in degrees anticlockwise from x axis");
53 43245 : params.addParam<Real>("dmax",
54 28830 : 360.0,
55 : "Maximum angle, measured in degrees anticlockwise from x axis. If "
56 : "dmin=0 and dmax=360 an annular mesh is created. "
57 : "Otherwise, only a sector of an annulus is created");
58 43245 : params.addRangeCheckedParam<Real>("growth_r",
59 28830 : 1.0,
60 : "growth_r!=0.0",
61 : "The ratio of radial sizes of successive rings of elements");
62 43245 : params.addParam<SubdomainID>(
63 28830 : "quad_subdomain_id", 0, "The subdomain ID given to the QUAD4 elements");
64 43245 : params.addParam<SubdomainID>("tri_subdomain_id",
65 28830 : 1,
66 : "The subdomain ID given to the TRI3 elements "
67 : "(these exist only if rmin=0, and they exist "
68 : "at the center of the disc");
69 14415 : params.addClassDescription("For rmin>0: creates an annular mesh of QUAD4 elements. For rmin=0: "
70 : "creates a disc mesh of QUAD4 and TRI3 elements. Boundary sidesets "
71 : "are created at rmax and rmin, and given these names. If dmin!=0 and "
72 : "dmax!=360, a sector of an annulus or disc is created. In this case "
73 : "boundary sidesets are also created at dmin and dmax, and "
74 : "given these names");
75 14415 : params.addParam<BoundaryName>("boundary_name_prefix",
76 : "If provided, prefix the built in boundary names with this string");
77 43245 : params.addParam<boundary_id_type>(
78 28830 : "boundary_id_offset", 0, "This offset is added to the generated boundary IDs");
79 :
80 14415 : return params;
81 0 : }
82 :
83 85 : AnnularMeshGenerator::AnnularMeshGenerator(const InputParameters & parameters)
84 : : MeshGenerator(parameters),
85 85 : _nt(getParam<unsigned int>("nt")),
86 85 : _rmin(getParam<Real>("rmin")),
87 85 : _rmax(getParam<Real>("rmax")),
88 85 : _radial_positions(getParam<std::vector<Real>>("radial_positions")),
89 170 : _nr(parameters.isParamSetByUser("radial_positions") ? _radial_positions.size() + 1
90 85 : : getParam<unsigned int>("nr")),
91 170 : _dmin(parameters.isParamSetByUser("tmin") ? getParam<Real>("tmin") / M_PI * 180.0
92 85 : : getParam<Real>("dmin")),
93 170 : _dmax(parameters.isParamSetByUser("tmax") ? getParam<Real>("tmax") / M_PI * 180.0
94 85 : : getParam<Real>("dmax")),
95 85 : _radians((parameters.isParamSetByUser("tmin") || parameters.isParamSetByUser("tmax")) ? true
96 : : false),
97 85 : _growth_r(getParam<Real>("growth_r")),
98 85 : _len(_growth_r == 1.0 ? (_rmax - _rmin) / _nr
99 31 : : (_rmax - _rmin) * (1.0 - std::abs(_growth_r)) /
100 31 : (1.0 - std::pow(std::abs(_growth_r), _nr))),
101 85 : _full_annulus(_dmin == 0.0 && _dmax == 360),
102 85 : _quad_subdomain_id(getParam<SubdomainID>("quad_subdomain_id")),
103 85 : _tri_subdomain_id(getParam<SubdomainID>("tri_subdomain_id")),
104 85 : _equal_area(getParam<bool>("equal_area")),
105 255 : _boundary_name_prefix(isParamValid("boundary_name_prefix")
106 170 : ? getParam<BoundaryName>("boundary_name_prefix") + "_"
107 : : ""),
108 170 : _boundary_id_offset(getParam<boundary_id_type>("boundary_id_offset"))
109 : {
110 179 : if ((parameters.isParamSetByUser("tmin") || parameters.isParamSetByUser("tmax")) &&
111 94 : (parameters.isParamSetByUser("dmin") || parameters.isParamSetByUser("dmax")))
112 0 : paramError("tmin",
113 : "You specified the angles using both degrees and radians. Please use degrees.");
114 :
115 85 : if (_radial_positions.size() != 0)
116 : {
117 25 : if (parameters.isParamSetByUser("nr"))
118 4 : paramError("nr", "The 'nr' parameter cannot be specified together with 'radial_positions'");
119 21 : if (parameters.isParamSetByUser("equal_area"))
120 4 : paramError("equal_area",
121 : "The 'equal_area' parameter cannot be specified together with 'radial_positions'");
122 17 : if (parameters.isParamSetByUser("growth_r"))
123 4 : paramError("growth_r",
124 : "The 'growth_r' parameter cannot be specified together with 'radial_positions'");
125 31 : for (auto rpos : _radial_positions)
126 22 : if (rpos <= _rmin || rpos >= _rmax)
127 4 : paramError(
128 : "radial_positions",
129 : "The following provided value is not within the bounds between 'rmin' and 'rmax': ",
130 : rpos);
131 : }
132 :
133 69 : if (_equal_area && parameters.isParamSetByUser("growth_r"))
134 4 : paramError("growth_r", "The 'growth_r' parameter cannot be combined with 'equal_area'");
135 :
136 65 : if (_rmax <= _rmin)
137 0 : paramError("rmax", "rmax must be greater than rmin");
138 65 : if (_dmax <= _dmin)
139 0 : paramError("dmax", "dmax must be greater than dmin");
140 65 : if (_dmax - _dmin > 360)
141 0 : paramError("dmax", "dmax - dmin must be <= 360");
142 65 : if (_nt <= (_dmax - _dmin) / 180.0)
143 0 : paramError("nt",
144 : "nt must be greater than (dmax - dmin) / 180 in order to avoid inverted "
145 : "elements");
146 65 : if (_quad_subdomain_id == _tri_subdomain_id)
147 0 : paramError("quad_subdomain_id", "quad_subdomain_id must not equal tri_subdomain_id");
148 65 : }
149 :
150 : std::unique_ptr<MeshBase>
151 65 : AnnularMeshGenerator::generate()
152 : {
153 : // Have MOOSE construct the correct libMesh::Mesh object using Mesh block and CLI parameters.
154 65 : auto mesh = buildMeshBaseObject();
155 :
156 65 : const Real dt = (_dmax - _dmin) / _nt;
157 :
158 65 : mesh->set_mesh_dimension(2);
159 65 : mesh->set_spatial_dimension(2);
160 65 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
161 :
162 65 : const unsigned num_angular_nodes = (_full_annulus ? _nt : _nt + 1);
163 65 : const unsigned num_nodes =
164 65 : (_rmin > 0.0 ? (_nr + 1) * num_angular_nodes : _nr * num_angular_nodes + 1);
165 65 : const unsigned min_nonzero_layer_num = (_rmin > 0.0 ? 0 : 1);
166 65 : std::vector<Node *> nodes(num_nodes);
167 65 : unsigned node_id = 0;
168 :
169 : // add nodes at rmax that aren't yet connected to any elements
170 65 : Real current_r = _rmax;
171 987 : for (unsigned angle_num = 0; angle_num < num_angular_nodes; ++angle_num)
172 : {
173 922 : const Real angle = _dmin + angle_num * dt;
174 922 : const Real x = current_r * std::cos(angle * M_PI / 180.0);
175 922 : const Real y = current_r * std::sin(angle * M_PI / 180.0);
176 922 : nodes[node_id] = mesh->add_point(Point(x, y, 0.0), node_id);
177 922 : node_id++;
178 : }
179 :
180 : // first value for the ring outer radius, only used for _equal_area
181 65 : Real outer_r = _rmax;
182 :
183 : // add nodes at smaller radii, and connect them with elements
184 472 : for (unsigned layer_num = _nr; layer_num > min_nonzero_layer_num; --layer_num)
185 : {
186 407 : if (layer_num == 1)
187 65 : current_r = _rmin; // account for precision loss
188 342 : else if (_radial_positions.size() > 0)
189 18 : current_r = _radial_positions[layer_num - 2];
190 : else
191 : {
192 324 : if (_equal_area)
193 : {
194 81 : const Real ring_area = (_rmax * _rmax - _rmin * _rmin) / _nr;
195 81 : current_r = std::sqrt(outer_r * outer_r - ring_area);
196 81 : outer_r = current_r;
197 : }
198 : else
199 : {
200 243 : if (_growth_r > 0)
201 162 : current_r -= _len * std::pow(_growth_r, layer_num - 1);
202 : else
203 81 : current_r -= _len * std::pow(std::abs(_growth_r), _nr - layer_num);
204 : }
205 : }
206 :
207 : // add node at angle = _dmin
208 814 : nodes[node_id] = mesh->add_point(Point(current_r * std::cos(_dmin * M_PI / 180.0),
209 407 : current_r * std::sin(_dmin * M_PI / 180.0),
210 : 0.0),
211 : node_id);
212 407 : node_id++;
213 5368 : for (unsigned angle_num = 1; angle_num < num_angular_nodes; ++angle_num)
214 : {
215 4961 : const Real angle = _dmin + angle_num * dt;
216 4961 : const Real x = current_r * std::cos(angle * M_PI / 180.0);
217 4961 : const Real y = current_r * std::sin(angle * M_PI / 180.0);
218 4961 : nodes[node_id] = mesh->add_point(Point(x, y, 0.0), node_id);
219 4961 : Elem * elem = mesh->add_elem(new Quad4);
220 4961 : elem->set_node(0, nodes[node_id]);
221 4961 : elem->set_node(1, nodes[node_id - 1]);
222 4961 : elem->set_node(2, nodes[node_id - num_angular_nodes - 1]);
223 4961 : elem->set_node(3, nodes[node_id - num_angular_nodes]);
224 4961 : elem->subdomain_id() = _quad_subdomain_id;
225 4961 : node_id++;
226 :
227 4961 : if (layer_num == _nr)
228 : // add outer boundary (boundary_id = 1)
229 857 : boundary_info.add_side(elem, 2, 1);
230 4961 : if (layer_num == 1)
231 : // add inner boundary (boundary_id = 0)
232 857 : boundary_info.add_side(elem, 0, 0);
233 4961 : if (!_full_annulus && angle_num == 1)
234 : // add tmin boundary (boundary_id = 2)
235 396 : boundary_info.add_side(elem, 1, 2);
236 4961 : if (!_full_annulus && angle_num == num_angular_nodes - 1)
237 : // add tmin boundary (boundary_id = 3)
238 396 : boundary_info.add_side(elem, 3, 3);
239 : }
240 407 : if (_full_annulus)
241 : {
242 : // add element connecting to node at angle=0
243 11 : Elem * elem = mesh->add_elem(new Quad4);
244 11 : elem->set_node(0, nodes[node_id - num_angular_nodes]);
245 11 : elem->set_node(1, nodes[node_id - 1]);
246 11 : elem->set_node(2, nodes[node_id - num_angular_nodes - 1]);
247 11 : elem->set_node(3, nodes[node_id - 2 * num_angular_nodes]);
248 11 : elem->subdomain_id() = _quad_subdomain_id;
249 :
250 11 : if (layer_num == _nr)
251 : // add outer boundary (boundary_id = 1)
252 11 : boundary_info.add_side(elem, 2, 1);
253 11 : if (layer_num == 1)
254 : // add inner boundary (boundary_id = 0)
255 11 : boundary_info.add_side(elem, 0, 0);
256 : }
257 : }
258 :
259 : // add single node at origin, if relevant
260 65 : if (_rmin == 0.0)
261 : {
262 0 : nodes[node_id] = mesh->add_point(Point(0.0, 0.0, 0.0), node_id);
263 0 : boundary_info.add_node(node_id, 0); // boundary_id=0 is centre
264 0 : for (unsigned angle_num = 0; angle_num < num_angular_nodes - 1; ++angle_num)
265 : {
266 0 : Elem * elem = mesh->add_elem(new Tri3);
267 0 : elem->set_node(0, nodes[node_id]);
268 0 : elem->set_node(1, nodes[node_id - num_angular_nodes + angle_num]);
269 0 : elem->set_node(2, nodes[node_id - num_angular_nodes + angle_num + 1]);
270 0 : elem->subdomain_id() = _tri_subdomain_id;
271 : }
272 0 : if (_full_annulus)
273 : {
274 0 : Elem * elem = mesh->add_elem(new Tri3);
275 0 : elem->set_node(0, nodes[node_id]);
276 0 : elem->set_node(1, nodes[node_id - 1]);
277 0 : elem->set_node(2, nodes[node_id - num_angular_nodes]);
278 0 : elem->subdomain_id() = _tri_subdomain_id;
279 : }
280 : }
281 :
282 65 : boundary_info.sideset_name(0) = "rmin";
283 65 : boundary_info.sideset_name(1) = "rmax";
284 65 : boundary_info.nodeset_name(0) = "rmin";
285 65 : boundary_info.nodeset_name(1) = "rmax";
286 65 : if (!_full_annulus)
287 : {
288 54 : if (_radians)
289 : {
290 9 : boundary_info.sideset_name(2) = "tmin";
291 9 : boundary_info.sideset_name(3) = "tmax";
292 9 : boundary_info.nodeset_name(2) = "tmin";
293 9 : boundary_info.nodeset_name(3) = "tmax";
294 : }
295 : else
296 : {
297 45 : boundary_info.sideset_name(2) = "dmin";
298 45 : boundary_info.sideset_name(3) = "dmax";
299 45 : boundary_info.nodeset_name(2) = "dmin";
300 45 : boundary_info.nodeset_name(3) = "dmax";
301 : }
302 : }
303 65 : if (_boundary_id_offset != 0 || !_boundary_name_prefix.empty())
304 : {
305 : // apply boundary id offset and name prefix
306 9 : const auto mesh_boundary_ids = boundary_info.get_boundary_ids();
307 45 : for (auto rit = mesh_boundary_ids.rbegin(); rit != mesh_boundary_ids.rend(); ++rit)
308 : {
309 :
310 36 : const std::string old_sideset_name = boundary_info.sideset_name(*rit);
311 36 : const std::string old_nodeset_name = boundary_info.nodeset_name(*rit);
312 :
313 36 : if (_boundary_id_offset != 0)
314 36 : MeshTools::Modification::change_boundary_id(*mesh, *rit, *rit + _boundary_id_offset);
315 :
316 36 : if (!_boundary_name_prefix.empty())
317 : {
318 36 : boundary_info.sideset_name(*rit + _boundary_id_offset) =
319 72 : _boundary_name_prefix + old_sideset_name;
320 36 : boundary_info.nodeset_name(*rit + _boundary_id_offset) =
321 72 : _boundary_name_prefix + old_nodeset_name;
322 : }
323 36 : }
324 9 : }
325 :
326 65 : mesh->prepare_for_use();
327 130 : return dynamic_pointer_cast<MeshBase>(mesh);
328 65 : }
|
