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 14401 : AnnularMeshGenerator::validParams()
22 : {
23 14401 : InputParameters params = MeshGenerator::validParams();
24 :
25 43203 : params.addRangeCheckedParam<unsigned int>(
26 28802 : "nr", 1, "nr>0", "Number of elements in the radial direction");
27 14401 : params.addRequiredRangeCheckedParam<unsigned int>(
28 : "nt", "nt>0", "Number of elements in the angular direction");
29 14401 : 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 14401 : params.addRequiredParam<Real>("rmax", "Outer radius");
34 14401 : params.addParam<std::vector<Real>>(
35 : "radial_positions", {}, "Directly prescribed positions of intermediate radial nodes");
36 43203 : params.addParam<bool>("equal_area",
37 28802 : false,
38 : "Whether to select the radial discretization "
39 : "to achieve equal areas of each ring");
40 43203 : params.addDeprecatedParam<Real>("tmin",
41 28802 : 0.0,
42 : "Minimum angle, measured in radians anticlockwise from x axis",
43 : "Use dmin instead");
44 43203 : params.addDeprecatedParam<Real>(
45 : "tmax",
46 28802 : 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 43203 : params.addParam<Real>(
52 28802 : "dmin", 0.0, "Minimum degree, measured in degrees anticlockwise from x axis");
53 43203 : params.addParam<Real>("dmax",
54 28802 : 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 43203 : params.addRangeCheckedParam<Real>("growth_r",
59 28802 : 1.0,
60 : "growth_r!=0.0",
61 : "The ratio of radial sizes of successive rings of elements");
62 43203 : params.addParam<SubdomainID>(
63 28802 : "quad_subdomain_id", 0, "The subdomain ID given to the QUAD4 elements");
64 43203 : params.addParam<SubdomainID>("tri_subdomain_id",
65 28802 : 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 14401 : 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 14401 : params.addParam<BoundaryName>("boundary_name_prefix",
76 : "If provided, prefix the built in boundary names with this string");
77 43203 : params.addParam<boundary_id_type>(
78 28802 : "boundary_id_offset", 0, "This offset is added to the generated boundary IDs");
79 :
80 14401 : return params;
81 0 : }
82 :
83 78 : AnnularMeshGenerator::AnnularMeshGenerator(const InputParameters & parameters)
84 : : MeshGenerator(parameters),
85 78 : _nt(getParam<unsigned int>("nt")),
86 78 : _rmin(getParam<Real>("rmin")),
87 78 : _rmax(getParam<Real>("rmax")),
88 78 : _radial_positions(getParam<std::vector<Real>>("radial_positions")),
89 156 : _nr(parameters.isParamSetByUser("radial_positions") ? _radial_positions.size() + 1
90 78 : : getParam<unsigned int>("nr")),
91 156 : _dmin(parameters.isParamSetByUser("tmin") ? getParam<Real>("tmin") / M_PI * 180.0
92 78 : : getParam<Real>("dmin")),
93 156 : _dmax(parameters.isParamSetByUser("tmax") ? getParam<Real>("tmax") / M_PI * 180.0
94 78 : : getParam<Real>("dmax")),
95 78 : _radians((parameters.isParamSetByUser("tmin") || parameters.isParamSetByUser("tmax")) ? true
96 : : false),
97 78 : _growth_r(getParam<Real>("growth_r")),
98 78 : _len(_growth_r == 1.0 ? (_rmax - _rmin) / _nr
99 28 : : (_rmax - _rmin) * (1.0 - std::abs(_growth_r)) /
100 28 : (1.0 - std::pow(std::abs(_growth_r), _nr))),
101 78 : _full_annulus(_dmin == 0.0 && _dmax == 360),
102 78 : _quad_subdomain_id(getParam<SubdomainID>("quad_subdomain_id")),
103 78 : _tri_subdomain_id(getParam<SubdomainID>("tri_subdomain_id")),
104 78 : _equal_area(getParam<bool>("equal_area")),
105 234 : _boundary_name_prefix(isParamValid("boundary_name_prefix")
106 156 : ? getParam<BoundaryName>("boundary_name_prefix") + "_"
107 : : ""),
108 156 : _boundary_id_offset(getParam<boundary_id_type>("boundary_id_offset"))
109 : {
110 164 : if ((parameters.isParamSetByUser("tmin") || parameters.isParamSetByUser("tmax")) &&
111 86 : (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 78 : if (_radial_positions.size() != 0)
116 : {
117 24 : if (parameters.isParamSetByUser("nr"))
118 4 : paramError("nr", "The 'nr' parameter cannot be specified together with 'radial_positions'");
119 20 : if (parameters.isParamSetByUser("equal_area"))
120 4 : paramError("equal_area",
121 : "The 'equal_area' parameter cannot be specified together with 'radial_positions'");
122 16 : if (parameters.isParamSetByUser("growth_r"))
123 4 : paramError("growth_r",
124 : "The 'growth_r' parameter cannot be specified together with 'radial_positions'");
125 28 : for (auto rpos : _radial_positions)
126 20 : 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 62 : if (_equal_area && parameters.isParamSetByUser("growth_r"))
134 4 : paramError("growth_r", "The 'growth_r' parameter cannot be combined with 'equal_area'");
135 :
136 58 : if (_rmax <= _rmin)
137 0 : paramError("rmax", "rmax must be greater than rmin");
138 58 : if (_dmax <= _dmin)
139 0 : paramError("dmax", "dmax must be greater than dmin");
140 58 : if (_dmax - _dmin > 360)
141 0 : paramError("dmax", "dmax - dmin must be <= 360");
142 58 : 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 58 : if (_quad_subdomain_id == _tri_subdomain_id)
147 0 : paramError("quad_subdomain_id", "quad_subdomain_id must not equal tri_subdomain_id");
148 58 : }
149 :
150 : std::unique_ptr<MeshBase>
151 58 : AnnularMeshGenerator::generate()
152 : {
153 : // Have MOOSE construct the correct libMesh::Mesh object using Mesh block and CLI parameters.
154 58 : auto mesh = buildMeshBaseObject();
155 :
156 58 : const Real dt = (_dmax - _dmin) / _nt;
157 :
158 58 : mesh->set_mesh_dimension(2);
159 58 : mesh->set_spatial_dimension(2);
160 58 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
161 :
162 58 : const unsigned num_angular_nodes = (_full_annulus ? _nt : _nt + 1);
163 58 : const unsigned num_nodes =
164 58 : (_rmin > 0.0 ? (_nr + 1) * num_angular_nodes : _nr * num_angular_nodes + 1);
165 58 : const unsigned min_nonzero_layer_num = (_rmin > 0.0 ? 0 : 1);
166 58 : std::vector<Node *> nodes(num_nodes);
167 58 : unsigned node_id = 0;
168 :
169 : // add nodes at rmax that aren't yet connected to any elements
170 58 : Real current_r = _rmax;
171 882 : for (unsigned angle_num = 0; angle_num < num_angular_nodes; ++angle_num)
172 : {
173 824 : const Real angle = _dmin + angle_num * dt;
174 824 : const Real x = current_r * std::cos(angle * M_PI / 180.0);
175 824 : const Real y = current_r * std::sin(angle * M_PI / 180.0);
176 824 : nodes[node_id] = mesh->add_point(Point(x, y, 0.0), node_id);
177 824 : node_id++;
178 : }
179 :
180 : // first value for the ring outer radius, only used for _equal_area
181 58 : Real outer_r = _rmax;
182 :
183 : // add nodes at smaller radii, and connect them with elements
184 420 : for (unsigned layer_num = _nr; layer_num > min_nonzero_layer_num; --layer_num)
185 : {
186 362 : if (layer_num == 1)
187 58 : current_r = _rmin; // account for precision loss
188 304 : else if (_radial_positions.size() > 0)
189 16 : current_r = _radial_positions[layer_num - 2];
190 : else
191 : {
192 288 : if (_equal_area)
193 : {
194 72 : const Real ring_area = (_rmax * _rmax - _rmin * _rmin) / _nr;
195 72 : current_r = std::sqrt(outer_r * outer_r - ring_area);
196 72 : outer_r = current_r;
197 : }
198 : else
199 : {
200 216 : if (_growth_r > 0)
201 144 : current_r -= _len * std::pow(_growth_r, layer_num - 1);
202 : else
203 72 : current_r -= _len * std::pow(std::abs(_growth_r), _nr - layer_num);
204 : }
205 : }
206 :
207 : // add node at angle = _dmin
208 724 : nodes[node_id] = mesh->add_point(Point(current_r * std::cos(_dmin * M_PI / 180.0),
209 362 : current_r * std::sin(_dmin * M_PI / 180.0),
210 : 0.0),
211 : node_id);
212 362 : node_id++;
213 4776 : for (unsigned angle_num = 1; angle_num < num_angular_nodes; ++angle_num)
214 : {
215 4414 : const Real angle = _dmin + angle_num * dt;
216 4414 : const Real x = current_r * std::cos(angle * M_PI / 180.0);
217 4414 : const Real y = current_r * std::sin(angle * M_PI / 180.0);
218 4414 : nodes[node_id] = mesh->add_point(Point(x, y, 0.0), node_id);
219 4414 : Elem * elem = mesh->add_elem(new Quad4);
220 4414 : elem->set_node(0, nodes[node_id]);
221 4414 : elem->set_node(1, nodes[node_id - 1]);
222 4414 : elem->set_node(2, nodes[node_id - num_angular_nodes - 1]);
223 4414 : elem->set_node(3, nodes[node_id - num_angular_nodes]);
224 4414 : elem->subdomain_id() = _quad_subdomain_id;
225 4414 : node_id++;
226 :
227 4414 : if (layer_num == _nr)
228 : // add outer boundary (boundary_id = 1)
229 766 : boundary_info.add_side(elem, 2, 1);
230 4414 : if (layer_num == 1)
231 : // add inner boundary (boundary_id = 0)
232 766 : boundary_info.add_side(elem, 0, 0);
233 4414 : if (!_full_annulus && angle_num == 1)
234 : // add tmin boundary (boundary_id = 2)
235 352 : boundary_info.add_side(elem, 1, 2);
236 4414 : if (!_full_annulus && angle_num == num_angular_nodes - 1)
237 : // add tmin boundary (boundary_id = 3)
238 352 : boundary_info.add_side(elem, 3, 3);
239 : }
240 362 : if (_full_annulus)
241 : {
242 : // add element connecting to node at angle=0
243 10 : Elem * elem = mesh->add_elem(new Quad4);
244 10 : elem->set_node(0, nodes[node_id - num_angular_nodes]);
245 10 : elem->set_node(1, nodes[node_id - 1]);
246 10 : elem->set_node(2, nodes[node_id - num_angular_nodes - 1]);
247 10 : elem->set_node(3, nodes[node_id - 2 * num_angular_nodes]);
248 10 : elem->subdomain_id() = _quad_subdomain_id;
249 :
250 10 : if (layer_num == _nr)
251 : // add outer boundary (boundary_id = 1)
252 10 : boundary_info.add_side(elem, 2, 1);
253 10 : if (layer_num == 1)
254 : // add inner boundary (boundary_id = 0)
255 10 : boundary_info.add_side(elem, 0, 0);
256 : }
257 : }
258 :
259 : // add single node at origin, if relevant
260 58 : 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 58 : boundary_info.sideset_name(0) = "rmin";
283 58 : boundary_info.sideset_name(1) = "rmax";
284 58 : boundary_info.nodeset_name(0) = "rmin";
285 58 : boundary_info.nodeset_name(1) = "rmax";
286 58 : if (!_full_annulus)
287 : {
288 48 : if (_radians)
289 : {
290 8 : boundary_info.sideset_name(2) = "tmin";
291 8 : boundary_info.sideset_name(3) = "tmax";
292 8 : boundary_info.nodeset_name(2) = "tmin";
293 8 : boundary_info.nodeset_name(3) = "tmax";
294 : }
295 : else
296 : {
297 40 : boundary_info.sideset_name(2) = "dmin";
298 40 : boundary_info.sideset_name(3) = "dmax";
299 40 : boundary_info.nodeset_name(2) = "dmin";
300 40 : boundary_info.nodeset_name(3) = "dmax";
301 : }
302 : }
303 58 : if (_boundary_id_offset != 0 || !_boundary_name_prefix.empty())
304 : {
305 : // apply boundary id offset and name prefix
306 8 : const auto mesh_boundary_ids = boundary_info.get_boundary_ids();
307 40 : for (auto rit = mesh_boundary_ids.rbegin(); rit != mesh_boundary_ids.rend(); ++rit)
308 : {
309 :
310 32 : const std::string old_sideset_name = boundary_info.sideset_name(*rit);
311 32 : const std::string old_nodeset_name = boundary_info.nodeset_name(*rit);
312 :
313 32 : if (_boundary_id_offset != 0)
314 32 : MeshTools::Modification::change_boundary_id(*mesh, *rit, *rit + _boundary_id_offset);
315 :
316 32 : if (!_boundary_name_prefix.empty())
317 : {
318 32 : boundary_info.sideset_name(*rit + _boundary_id_offset) =
319 64 : _boundary_name_prefix + old_sideset_name;
320 32 : boundary_info.nodeset_name(*rit + _boundary_id_offset) =
321 64 : _boundary_name_prefix + old_nodeset_name;
322 : }
323 32 : }
324 8 : }
325 :
326 58 : mesh->prepare_for_use();
327 116 : return dynamic_pointer_cast<MeshBase>(mesh);
328 58 : }
|
