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 "AnnularMesh.h"
11 :
12 : #include "MooseApp.h"
13 :
14 : #include "libmesh/face_quad4.h"
15 : #include "libmesh/face_tri3.h"
16 :
17 : registerMooseObject("MooseApp", AnnularMesh);
18 :
19 : InputParameters
20 14449 : AnnularMesh::validParams()
21 : {
22 14449 : InputParameters params = MooseMesh::validParams();
23 43347 : params.addRangeCheckedParam<unsigned int>(
24 28898 : "nr", 1, "nr>0", "Number of elements in the radial direction");
25 14449 : params.addRequiredRangeCheckedParam<unsigned int>(
26 : "nt", "nt>0", "Number of elements in the angular direction");
27 14449 : params.addRequiredRangeCheckedParam<Real>(
28 : "rmin",
29 : "rmin>=0.0",
30 : "Inner radius. If rmin=0 then a disc mesh (with no central hole) will be created.");
31 14449 : params.addRequiredParam<Real>("rmax", "Outer radius");
32 43347 : params.addDeprecatedParam<Real>("tmin",
33 28898 : 0.0,
34 : "Minimum angle, measured in radians anticlockwise from x axis",
35 : "Use dmin instead");
36 43347 : params.addDeprecatedParam<Real>(
37 : "tmax",
38 28898 : 2 * M_PI,
39 : "Maximum angle, measured in radians anticlockwise from x axis. If "
40 : "tmin=0 and tmax=2Pi an annular mesh is created. "
41 : "Otherwise, only a sector of an annulus is created",
42 : "Use dmin instead");
43 43347 : params.addParam<Real>(
44 28898 : "dmin", 0.0, "Minimum degree, measured in degrees anticlockwise from x axis");
45 43347 : params.addParam<Real>("dmax",
46 28898 : 360.0,
47 : "Maximum angle, measured in degrees anticlockwise from x axis. If "
48 : "dmin=0 and dmax=360 an annular mesh is created. "
49 : "Otherwise, only a sector of an annulus is created");
50 43347 : params.addRangeCheckedParam<Real>(
51 28898 : "growth_r", 1.0, "growth_r>0.0", "The ratio of radial sizes of successive rings of elements");
52 43347 : params.addParam<SubdomainID>(
53 28898 : "quad_subdomain_id", 0, "The subdomain ID given to the QUAD4 elements");
54 43347 : params.addParam<SubdomainID>("tri_subdomain_id",
55 28898 : 1,
56 : "The subdomain ID given to the TRI3 elements "
57 : "(these exist only if rmin=0, and they exist "
58 : "at the center of the disc");
59 14449 : params.addClassDescription("For rmin>0: creates an annular mesh of QUAD4 elements. For rmin=0: "
60 : "creates a disc mesh of QUAD4 and TRI3 elements. Boundary sidesets "
61 : "are created at rmax and rmin, and given these names. If dmin!=0 and "
62 : "dmax!=360, a sector of an annulus or disc is created. In this case "
63 : "boundary sidesets are also created a dmin and dmax, and "
64 : "given these names");
65 14449 : return params;
66 0 : }
67 :
68 112 : AnnularMesh::AnnularMesh(const InputParameters & parameters)
69 : : MooseMesh(parameters),
70 112 : _nr(getParam<unsigned int>("nr")),
71 112 : _nt(getParam<unsigned int>("nt")),
72 112 : _rmin(getParam<Real>("rmin")),
73 112 : _rmax(getParam<Real>("rmax")),
74 224 : _dmin(parameters.isParamSetByUser("tmin") ? getParam<Real>("tmin") / M_PI * 180.0
75 112 : : getParam<Real>("dmin")),
76 224 : _dmax(parameters.isParamSetByUser("tmax") ? getParam<Real>("tmax") / M_PI * 180.0
77 112 : : getParam<Real>("dmax")),
78 112 : _radians((parameters.isParamSetByUser("tmin") || parameters.isParamSetByUser("tmax")) ? true
79 : : false),
80 112 : _growth_r(getParam<Real>("growth_r")),
81 112 : _len(_growth_r == 1.0 ? (_rmax - _rmin) / _nr
82 112 : : (_rmax - _rmin) * (1.0 - _growth_r) / (1.0 - std::pow(_growth_r, _nr))),
83 112 : _full_annulus(_dmin == 0.0 && _dmax == 360),
84 112 : _quad_subdomain_id(getParam<SubdomainID>("quad_subdomain_id")),
85 112 : _tri_subdomain_id(getParam<SubdomainID>("tri_subdomain_id")),
86 112 : _dims_may_have_changed(false)
87 : {
88 264 : if ((parameters.isParamSetByUser("tmin") || parameters.isParamSetByUser("tmax")) &&
89 152 : (parameters.isParamSetByUser("dmin") || parameters.isParamSetByUser("dmax")))
90 0 : paramError("tmin",
91 : "You specified the angles using both degrees and radians. Please use degrees.");
92 :
93 112 : if (_rmax <= _rmin)
94 8 : paramError("rmax", "rmax must be greater than rmin");
95 104 : if (_dmax <= _dmin)
96 8 : paramError("dmax", "dmax must be greater than dmin");
97 96 : if (_dmax - _dmin > 360)
98 8 : paramError("dmax", "dmax - dmin must be <= 360");
99 88 : if (_nt <= (_dmax - _dmin) / 180.0)
100 12 : paramError("nt",
101 : "nt must be greater than (dmax - dmin) / 180 in order to avoid inverted "
102 : "elements");
103 76 : if (_quad_subdomain_id == _tri_subdomain_id)
104 4 : paramError("quad_subdomain_id", "quad_subdomain_id must not equal tri_subdomain_id");
105 72 : }
106 :
107 : void
108 0 : AnnularMesh::prepared(bool state)
109 : {
110 0 : MooseMesh::prepared(state);
111 :
112 : // Fall back on scanning the mesh for coordinates instead of using input parameters for queries
113 0 : if (!state)
114 0 : _dims_may_have_changed = true;
115 0 : }
116 :
117 : Real
118 132 : AnnularMesh::getMinInDimension(unsigned int component) const
119 : {
120 132 : if (_dims_may_have_changed)
121 0 : return MooseMesh::getMinInDimension(component);
122 :
123 132 : switch (component)
124 : {
125 0 : case 0:
126 0 : return -_rmax;
127 66 : case 1:
128 66 : return -_rmax;
129 66 : case 2:
130 66 : return 0.0;
131 0 : default:
132 0 : mooseError("Invalid component");
133 : }
134 : }
135 :
136 : Real
137 132 : AnnularMesh::getMaxInDimension(unsigned int component) const
138 : {
139 132 : if (_dims_may_have_changed)
140 0 : return MooseMesh::getMaxInDimension(component);
141 :
142 132 : switch (component)
143 : {
144 0 : case 0:
145 0 : return _rmax;
146 66 : case 1:
147 66 : return _rmax;
148 66 : case 2:
149 66 : return 0.0;
150 0 : default:
151 0 : mooseError("Invalid component");
152 : }
153 : }
154 :
155 : std::unique_ptr<MooseMesh>
156 0 : AnnularMesh::safeClone() const
157 : {
158 0 : return _app.getFactory().copyConstruct(*this);
159 : }
160 :
161 : void
162 66 : AnnularMesh::buildMesh()
163 : {
164 66 : const Real dt = (_dmax - _dmin) / _nt;
165 :
166 66 : MeshBase & mesh = getMesh();
167 66 : mesh.clear();
168 66 : mesh.set_mesh_dimension(2);
169 66 : mesh.set_spatial_dimension(2);
170 66 : libMesh::BoundaryInfo & boundary_info = mesh.get_boundary_info();
171 :
172 66 : const unsigned num_angular_nodes = (_full_annulus ? _nt : _nt + 1);
173 66 : const unsigned num_nodes =
174 66 : (_rmin > 0.0 ? (_nr + 1) * num_angular_nodes : _nr * num_angular_nodes + 1);
175 66 : const unsigned min_nonzero_layer_num = (_rmin > 0.0 ? 0 : 1);
176 66 : std::vector<Node *> nodes(num_nodes);
177 66 : unsigned node_id = 0;
178 :
179 : // add nodes at rmax that aren't yet connected to any elements
180 66 : Real current_r = _rmax;
181 902 : for (unsigned angle_num = 0; angle_num < num_angular_nodes; ++angle_num)
182 : {
183 836 : const Real angle = _dmin + angle_num * dt;
184 836 : const Real x = current_r * std::cos(angle * M_PI / 180.0);
185 836 : const Real y = current_r * std::sin(angle * M_PI / 180.0);
186 836 : nodes[node_id] = mesh.add_point(Point(x, y, 0.0), node_id);
187 836 : node_id++;
188 : }
189 :
190 : // add nodes at smaller radii, and connect them with elements
191 693 : for (unsigned layer_num = _nr; layer_num > min_nonzero_layer_num; --layer_num)
192 : {
193 627 : if (layer_num == 1)
194 33 : current_r = _rmin; // account for precision loss
195 : else
196 594 : current_r -= _len * std::pow(_growth_r, layer_num - 1);
197 :
198 : // add node at angle = _dmin
199 627 : nodes[node_id] = mesh.add_point(Point(current_r * std::cos(_dmin * M_PI / 180.0),
200 627 : current_r * std::sin(_dmin * M_PI / 180.0),
201 : 0.0),
202 : node_id);
203 627 : node_id++;
204 7942 : for (unsigned angle_num = 1; angle_num < num_angular_nodes; ++angle_num)
205 : {
206 7315 : const Real angle = _dmin + angle_num * dt;
207 7315 : const Real x = current_r * std::cos(angle * M_PI / 180.0);
208 7315 : const Real y = current_r * std::sin(angle * M_PI / 180.0);
209 7315 : nodes[node_id] = mesh.add_point(Point(x, y, 0.0), node_id);
210 7315 : Elem * elem = mesh.add_elem(new libMesh::Quad4);
211 7315 : elem->set_node(0, nodes[node_id]);
212 7315 : elem->set_node(1, nodes[node_id - 1]);
213 7315 : elem->set_node(2, nodes[node_id - num_angular_nodes - 1]);
214 7315 : elem->set_node(3, nodes[node_id - num_angular_nodes]);
215 7315 : elem->subdomain_id() = _quad_subdomain_id;
216 7315 : node_id++;
217 :
218 7315 : if (layer_num == _nr)
219 : // add outer boundary (boundary_id = 1)
220 770 : boundary_info.add_side(elem, 2, 1);
221 7315 : if (layer_num == 1)
222 : // add inner boundary (boundary_id = 0)
223 385 : boundary_info.add_side(elem, 0, 0);
224 7315 : if (!_full_annulus && angle_num == 1)
225 : // add tmin boundary (boundary_id = 2)
226 418 : boundary_info.add_side(elem, 1, 2);
227 7315 : if (!_full_annulus && angle_num == num_angular_nodes - 1)
228 : // add tmin boundary (boundary_id = 3)
229 418 : boundary_info.add_side(elem, 3, 3);
230 : }
231 627 : if (_full_annulus)
232 : {
233 : // add element connecting to node at angle=0
234 209 : Elem * elem = mesh.add_elem(new libMesh::Quad4);
235 209 : elem->set_node(0, nodes[node_id - num_angular_nodes]);
236 209 : elem->set_node(1, nodes[node_id - 1]);
237 209 : elem->set_node(2, nodes[node_id - num_angular_nodes - 1]);
238 209 : elem->set_node(3, nodes[node_id - 2 * num_angular_nodes]);
239 209 : elem->subdomain_id() = _quad_subdomain_id;
240 :
241 209 : if (layer_num == _nr)
242 : // add outer boundary (boundary_id = 1)
243 22 : boundary_info.add_side(elem, 2, 1);
244 209 : if (layer_num == 1)
245 : // add inner boundary (boundary_id = 0)
246 11 : boundary_info.add_side(elem, 0, 0);
247 : }
248 : }
249 :
250 : // add single node at origin, if relevant
251 66 : if (_rmin == 0.0)
252 : {
253 33 : nodes[node_id] = mesh.add_point(Point(0.0, 0.0, 0.0), node_id);
254 33 : boundary_info.add_node(node_id, 0); // boundary_id=0 is centre
255 418 : for (unsigned angle_num = 0; angle_num < num_angular_nodes - 1; ++angle_num)
256 : {
257 385 : Elem * elem = mesh.add_elem(new libMesh::Tri3);
258 385 : elem->set_node(0, nodes[node_id]);
259 385 : elem->set_node(1, nodes[node_id - num_angular_nodes + angle_num]);
260 385 : elem->set_node(2, nodes[node_id - num_angular_nodes + angle_num + 1]);
261 385 : elem->subdomain_id() = _tri_subdomain_id;
262 : }
263 33 : if (_full_annulus)
264 : {
265 11 : Elem * elem = mesh.add_elem(new libMesh::Tri3);
266 11 : elem->set_node(0, nodes[node_id]);
267 11 : elem->set_node(1, nodes[node_id - 1]);
268 11 : elem->set_node(2, nodes[node_id - num_angular_nodes]);
269 11 : elem->subdomain_id() = _tri_subdomain_id;
270 : }
271 : }
272 :
273 66 : boundary_info.sideset_name(0) = "rmin";
274 66 : boundary_info.sideset_name(1) = "rmax";
275 66 : boundary_info.nodeset_name(0) = "rmin";
276 66 : boundary_info.nodeset_name(1) = "rmax";
277 66 : if (!_full_annulus)
278 : {
279 44 : if (_radians)
280 : {
281 22 : boundary_info.sideset_name(2) = "tmin";
282 22 : boundary_info.sideset_name(3) = "tmax";
283 22 : boundary_info.nodeset_name(2) = "tmin";
284 22 : boundary_info.nodeset_name(3) = "tmax";
285 : }
286 : else
287 : {
288 22 : boundary_info.sideset_name(2) = "dmin";
289 22 : boundary_info.sideset_name(3) = "dmax";
290 22 : boundary_info.nodeset_name(2) = "dmin";
291 22 : boundary_info.nodeset_name(3) = "dmax";
292 : }
293 : }
294 :
295 66 : mesh.prepare_for_use();
296 66 : }
|
