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 14461 : AnnularMesh::validParams()
21 : {
22 14461 : InputParameters params = MooseMesh::validParams();
23 43383 : params.addRangeCheckedParam<unsigned int>(
24 28922 : "nr", 1, "nr>0", "Number of elements in the radial direction");
25 14461 : params.addRequiredRangeCheckedParam<unsigned int>(
26 : "nt", "nt>0", "Number of elements in the angular direction");
27 14461 : 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 14461 : params.addRequiredParam<Real>("rmax", "Outer radius");
32 43383 : params.addDeprecatedParam<Real>("tmin",
33 28922 : 0.0,
34 : "Minimum angle, measured in radians anticlockwise from x axis",
35 : "Use dmin instead");
36 43383 : params.addDeprecatedParam<Real>(
37 : "tmax",
38 28922 : 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 43383 : params.addParam<Real>(
44 28922 : "dmin", 0.0, "Minimum degree, measured in degrees anticlockwise from x axis");
45 43383 : params.addParam<Real>("dmax",
46 28922 : 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 43383 : params.addRangeCheckedParam<Real>(
51 28922 : "growth_r", 1.0, "growth_r>0.0", "The ratio of radial sizes of successive rings of elements");
52 43383 : params.addParam<SubdomainID>(
53 28922 : "quad_subdomain_id", 0, "The subdomain ID given to the QUAD4 elements");
54 43383 : params.addParam<SubdomainID>("tri_subdomain_id",
55 28922 : 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 14461 : 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 14461 : return params;
66 0 : }
67 :
68 118 : AnnularMesh::AnnularMesh(const InputParameters & parameters)
69 : : MooseMesh(parameters),
70 118 : _nr(getParam<unsigned int>("nr")),
71 118 : _nt(getParam<unsigned int>("nt")),
72 118 : _rmin(getParam<Real>("rmin")),
73 118 : _rmax(getParam<Real>("rmax")),
74 236 : _dmin(parameters.isParamSetByUser("tmin") ? getParam<Real>("tmin") / M_PI * 180.0
75 118 : : getParam<Real>("dmin")),
76 236 : _dmax(parameters.isParamSetByUser("tmax") ? getParam<Real>("tmax") / M_PI * 180.0
77 118 : : getParam<Real>("dmax")),
78 118 : _radians((parameters.isParamSetByUser("tmin") || parameters.isParamSetByUser("tmax")) ? true
79 : : false),
80 118 : _growth_r(getParam<Real>("growth_r")),
81 118 : _len(_growth_r == 1.0 ? (_rmax - _rmin) / _nr
82 118 : : (_rmax - _rmin) * (1.0 - _growth_r) / (1.0 - std::pow(_growth_r, _nr))),
83 118 : _full_annulus(_dmin == 0.0 && _dmax == 360),
84 118 : _quad_subdomain_id(getParam<SubdomainID>("quad_subdomain_id")),
85 118 : _tri_subdomain_id(getParam<SubdomainID>("tri_subdomain_id")),
86 118 : _dims_may_have_changed(false)
87 : {
88 278 : if ((parameters.isParamSetByUser("tmin") || parameters.isParamSetByUser("tmax")) &&
89 160 : (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 118 : if (_rmax <= _rmin)
94 8 : paramError("rmax", "rmax must be greater than rmin");
95 110 : if (_dmax <= _dmin)
96 8 : paramError("dmax", "dmax must be greater than dmin");
97 102 : if (_dmax - _dmin > 360)
98 8 : paramError("dmax", "dmax - dmin must be <= 360");
99 94 : 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 82 : if (_quad_subdomain_id == _tri_subdomain_id)
104 4 : paramError("quad_subdomain_id", "quad_subdomain_id must not equal tri_subdomain_id");
105 78 : }
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 144 : AnnularMesh::getMinInDimension(unsigned int component) const
119 : {
120 144 : if (_dims_may_have_changed)
121 0 : return MooseMesh::getMinInDimension(component);
122 :
123 144 : switch (component)
124 : {
125 0 : case 0:
126 0 : return -_rmax;
127 72 : case 1:
128 72 : return -_rmax;
129 72 : case 2:
130 72 : return 0.0;
131 0 : default:
132 0 : mooseError("Invalid component");
133 : }
134 : }
135 :
136 : Real
137 144 : AnnularMesh::getMaxInDimension(unsigned int component) const
138 : {
139 144 : if (_dims_may_have_changed)
140 0 : return MooseMesh::getMaxInDimension(component);
141 :
142 144 : switch (component)
143 : {
144 0 : case 0:
145 0 : return _rmax;
146 72 : case 1:
147 72 : return _rmax;
148 72 : case 2:
149 72 : 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 72 : AnnularMesh::buildMesh()
163 : {
164 72 : const Real dt = (_dmax - _dmin) / _nt;
165 :
166 72 : MeshBase & mesh = getMesh();
167 72 : mesh.clear();
168 72 : mesh.set_mesh_dimension(2);
169 72 : mesh.set_spatial_dimension(2);
170 72 : libMesh::BoundaryInfo & boundary_info = mesh.get_boundary_info();
171 :
172 72 : const unsigned num_angular_nodes = (_full_annulus ? _nt : _nt + 1);
173 72 : const unsigned num_nodes =
174 72 : (_rmin > 0.0 ? (_nr + 1) * num_angular_nodes : _nr * num_angular_nodes + 1);
175 72 : const unsigned min_nonzero_layer_num = (_rmin > 0.0 ? 0 : 1);
176 72 : std::vector<Node *> nodes(num_nodes);
177 72 : unsigned node_id = 0;
178 :
179 : // add nodes at rmax that aren't yet connected to any elements
180 72 : Real current_r = _rmax;
181 984 : for (unsigned angle_num = 0; angle_num < num_angular_nodes; ++angle_num)
182 : {
183 912 : const Real angle = _dmin + angle_num * dt;
184 912 : const Real x = current_r * std::cos(angle * M_PI / 180.0);
185 912 : const Real y = current_r * std::sin(angle * M_PI / 180.0);
186 912 : nodes[node_id] = mesh.add_point(Point(x, y, 0.0), node_id);
187 912 : node_id++;
188 : }
189 :
190 : // add nodes at smaller radii, and connect them with elements
191 756 : for (unsigned layer_num = _nr; layer_num > min_nonzero_layer_num; --layer_num)
192 : {
193 684 : if (layer_num == 1)
194 36 : current_r = _rmin; // account for precision loss
195 : else
196 648 : current_r -= _len * std::pow(_growth_r, layer_num - 1);
197 :
198 : // add node at angle = _dmin
199 684 : nodes[node_id] = mesh.add_point(Point(current_r * std::cos(_dmin * M_PI / 180.0),
200 684 : current_r * std::sin(_dmin * M_PI / 180.0),
201 : 0.0),
202 : node_id);
203 684 : node_id++;
204 8664 : for (unsigned angle_num = 1; angle_num < num_angular_nodes; ++angle_num)
205 : {
206 7980 : const Real angle = _dmin + angle_num * dt;
207 7980 : const Real x = current_r * std::cos(angle * M_PI / 180.0);
208 7980 : const Real y = current_r * std::sin(angle * M_PI / 180.0);
209 7980 : nodes[node_id] = mesh.add_point(Point(x, y, 0.0), node_id);
210 7980 : Elem * elem = mesh.add_elem(new libMesh::Quad4);
211 7980 : elem->set_node(0, nodes[node_id]);
212 7980 : elem->set_node(1, nodes[node_id - 1]);
213 7980 : elem->set_node(2, nodes[node_id - num_angular_nodes - 1]);
214 7980 : elem->set_node(3, nodes[node_id - num_angular_nodes]);
215 7980 : elem->subdomain_id() = _quad_subdomain_id;
216 7980 : node_id++;
217 :
218 7980 : if (layer_num == _nr)
219 : // add outer boundary (boundary_id = 1)
220 840 : boundary_info.add_side(elem, 2, 1);
221 7980 : if (layer_num == 1)
222 : // add inner boundary (boundary_id = 0)
223 420 : boundary_info.add_side(elem, 0, 0);
224 7980 : if (!_full_annulus && angle_num == 1)
225 : // add tmin boundary (boundary_id = 2)
226 456 : boundary_info.add_side(elem, 1, 2);
227 7980 : if (!_full_annulus && angle_num == num_angular_nodes - 1)
228 : // add tmin boundary (boundary_id = 3)
229 456 : boundary_info.add_side(elem, 3, 3);
230 : }
231 684 : if (_full_annulus)
232 : {
233 : // add element connecting to node at angle=0
234 228 : Elem * elem = mesh.add_elem(new libMesh::Quad4);
235 228 : elem->set_node(0, nodes[node_id - num_angular_nodes]);
236 228 : elem->set_node(1, nodes[node_id - 1]);
237 228 : elem->set_node(2, nodes[node_id - num_angular_nodes - 1]);
238 228 : elem->set_node(3, nodes[node_id - 2 * num_angular_nodes]);
239 228 : elem->subdomain_id() = _quad_subdomain_id;
240 :
241 228 : if (layer_num == _nr)
242 : // add outer boundary (boundary_id = 1)
243 24 : boundary_info.add_side(elem, 2, 1);
244 228 : if (layer_num == 1)
245 : // add inner boundary (boundary_id = 0)
246 12 : boundary_info.add_side(elem, 0, 0);
247 : }
248 : }
249 :
250 : // add single node at origin, if relevant
251 72 : if (_rmin == 0.0)
252 : {
253 36 : nodes[node_id] = mesh.add_point(Point(0.0, 0.0, 0.0), node_id);
254 36 : boundary_info.add_node(node_id, 0); // boundary_id=0 is centre
255 456 : for (unsigned angle_num = 0; angle_num < num_angular_nodes - 1; ++angle_num)
256 : {
257 420 : Elem * elem = mesh.add_elem(new libMesh::Tri3);
258 420 : elem->set_node(0, nodes[node_id]);
259 420 : elem->set_node(1, nodes[node_id - num_angular_nodes + angle_num]);
260 420 : elem->set_node(2, nodes[node_id - num_angular_nodes + angle_num + 1]);
261 420 : elem->subdomain_id() = _tri_subdomain_id;
262 : }
263 36 : if (_full_annulus)
264 : {
265 12 : Elem * elem = mesh.add_elem(new libMesh::Tri3);
266 12 : elem->set_node(0, nodes[node_id]);
267 12 : elem->set_node(1, nodes[node_id - 1]);
268 12 : elem->set_node(2, nodes[node_id - num_angular_nodes]);
269 12 : elem->subdomain_id() = _tri_subdomain_id;
270 : }
271 : }
272 :
273 72 : boundary_info.sideset_name(0) = "rmin";
274 72 : boundary_info.sideset_name(1) = "rmax";
275 72 : boundary_info.nodeset_name(0) = "rmin";
276 72 : boundary_info.nodeset_name(1) = "rmax";
277 72 : if (!_full_annulus)
278 : {
279 48 : if (_radians)
280 : {
281 24 : boundary_info.sideset_name(2) = "tmin";
282 24 : boundary_info.sideset_name(3) = "tmax";
283 24 : boundary_info.nodeset_name(2) = "tmin";
284 24 : boundary_info.nodeset_name(3) = "tmax";
285 : }
286 : else
287 : {
288 24 : boundary_info.sideset_name(2) = "dmin";
289 24 : boundary_info.sideset_name(3) = "dmax";
290 24 : boundary_info.nodeset_name(2) = "dmin";
291 24 : boundary_info.nodeset_name(3) = "dmax";
292 : }
293 : }
294 :
295 72 : mesh.prepare_for_use();
296 72 : }
|
