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 "SpiralAnnularMesh.h"
11 :
12 : #include "MooseApp.h"
13 :
14 : #include "libmesh/face_quad4.h"
15 : #include "libmesh/face_tri3.h"
16 :
17 : registerMooseObject("MooseApp", SpiralAnnularMesh);
18 :
19 : InputParameters
20 14285 : SpiralAnnularMesh::validParams()
21 : {
22 14285 : InputParameters params = MooseMesh::validParams();
23 14285 : params.addRequiredRangeCheckedParam<Real>(
24 : "inner_radius", "inner_radius>0.", "The size of the inner circle.");
25 14285 : params.addRequiredRangeCheckedParam<Real>("outer_radius",
26 : "outer_radius>0.",
27 : "The size of the outer circle."
28 : " Logically, it has to be greater than inner_radius");
29 14285 : params.addRequiredRangeCheckedParam<unsigned int>(
30 : "nodes_per_ring", "nodes_per_ring>5", "Number of nodes on each ring.");
31 42855 : params.addParam<bool>(
32 28570 : "use_tri6", false, "Generate mesh of TRI6 elements instead of TRI3 elements.");
33 14285 : params.addRequiredRangeCheckedParam<unsigned int>(
34 : "num_rings", "num_rings>1", "The number of rings.");
35 42855 : params.addParam<boundary_id_type>(
36 28570 : "cylinder_bid", 1, "The boundary id to use for the cylinder (inner circle)");
37 42855 : params.addParam<boundary_id_type>(
38 28570 : "exterior_bid", 2, "The boundary id to use for the exterior (outer circle)");
39 14285 : params.addParam<Real>("initial_delta_r",
40 : "Width of the initial layer of elements around the cylinder."
41 : "This number should be approximately"
42 : " 2 * pi * inner_radius / nodes_per_ring to ensure that the"
43 : " initial layer of elements is almost equilateral");
44 14285 : params.addClassDescription("Creates an annual mesh based on TRI3 elements"
45 : " (it can also be TRI6 elements) on several rings.");
46 :
47 14285 : return params;
48 0 : }
49 :
50 10 : SpiralAnnularMesh::SpiralAnnularMesh(const InputParameters & parameters)
51 : : MooseMesh(parameters),
52 10 : _inner_radius(getParam<Real>("inner_radius")),
53 10 : _outer_radius(getParam<Real>("outer_radius")),
54 10 : _radial_bias(1.0),
55 10 : _nodes_per_ring(getParam<unsigned int>("nodes_per_ring")),
56 10 : _use_tri6(getParam<bool>("use_tri6")),
57 10 : _num_rings(getParam<unsigned int>("num_rings")),
58 10 : _cylinder_bid(getParam<boundary_id_type>("cylinder_bid")),
59 10 : _exterior_bid(getParam<boundary_id_type>("exterior_bid")),
60 10 : _initial_delta_r(2 * libMesh::pi * _inner_radius / _nodes_per_ring)
61 : {
62 : // catch likely user errors
63 10 : if (_outer_radius <= _inner_radius)
64 0 : mooseError("SpiralAnnularMesh: outer_radius must be greater than inner_radius");
65 10 : }
66 :
67 : std::unique_ptr<MooseMesh>
68 0 : SpiralAnnularMesh::safeClone() const
69 : {
70 0 : return _app.getFactory().copyConstruct(*this);
71 : }
72 :
73 : void
74 10 : SpiralAnnularMesh::buildMesh()
75 : {
76 : {
77 : // Compute the radial bias given:
78 : // .) the inner radius
79 : // .) the outer radius
80 : // .) the initial_delta_r
81 : // .) the desired number of intervals
82 : // Note: the exponent n used in the formula is one less than the
83 : // number of rings the user requests.
84 10 : Real alpha = 1.1;
85 10 : int n = _num_rings - 1;
86 :
87 : // lambda used to compute the residual and Jacobian for the Newton iterations.
88 : // We capture parameters which don't need to change from the current scope at
89 : // the time this lambda is declared. The values are not updated later, so we
90 : // can't use this for e.g. f, df, and alpha.
91 180 : auto newton = [this, n](Real & f, Real & df, const Real & alpha)
92 : {
93 60 : f = (1. - std::pow(alpha, n + 1)) / (1. - alpha) -
94 60 : (_outer_radius - _inner_radius) / _initial_delta_r;
95 60 : df = (-(n + 1) * (1 - alpha) * std::pow(alpha, n) + (1. - std::pow(alpha, n + 1))) /
96 60 : (1. - alpha) / (1. - alpha);
97 70 : };
98 :
99 : Real f, df;
100 10 : int num_iter = 1;
101 10 : newton(f, df, alpha);
102 :
103 60 : while (std::abs(f) > 1.e-9 && num_iter <= 25)
104 : {
105 : // Compute and apply update.
106 50 : Real dx = -f / df;
107 50 : alpha += dx;
108 50 : newton(f, df, alpha);
109 50 : num_iter++;
110 : }
111 :
112 : // In case the Newton iteration fails to converge.
113 10 : if (num_iter > 25)
114 0 : mooseError("Newton iteration failed to converge (more than 25 iterations).");
115 :
116 : // Set radial basis to the value of alpha that we computed with Newton.
117 10 : _radial_bias = alpha;
118 : }
119 :
120 : // The number of rings specified by the user does not include the ring at
121 : // the surface of the cylinder itself, so we increment it by one now.
122 10 : _num_rings += 1;
123 :
124 : // Mesh we are eventually going to create.
125 10 : MeshBase & mesh = getMesh();
126 10 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
127 :
128 : // Data structure that holds pointers to the Nodes of each ring.
129 10 : std::vector<std::vector<Node *>> ring_nodes(_num_rings);
130 :
131 : // Initialize radius and delta_r variables.
132 10 : Real radius = _inner_radius;
133 10 : Real delta_r = _initial_delta_r;
134 :
135 : // Node id counter.
136 10 : unsigned int current_node_id = 0;
137 :
138 120 : for (std::size_t r = 0; r < _num_rings; ++r)
139 : {
140 110 : ring_nodes[r].resize(_nodes_per_ring);
141 :
142 : // Add nodes starting from either theta=0 or theta=pi/nodes_per_ring
143 110 : Real theta = r % 2 == 0 ? 0 : (libMesh::pi / _nodes_per_ring);
144 2090 : for (std::size_t n = 0; n < _nodes_per_ring; ++n)
145 : {
146 1980 : ring_nodes[r][n] = mesh.add_point(Point(radius * std::cos(theta), radius * std::sin(theta)),
147 1980 : current_node_id++);
148 : // Update angle
149 1980 : theta += 2 * libMesh::pi / _nodes_per_ring;
150 : }
151 :
152 : // Go to next ring
153 110 : radius += delta_r;
154 110 : delta_r *= _radial_bias;
155 : }
156 :
157 : // Add elements
158 110 : for (std::size_t r = 0; r < _num_rings - 1; ++r)
159 : {
160 : // even -> odd ring
161 100 : if (r % 2 == 0)
162 : {
163 : // Inner ring (n, n*, n+1)
164 : // Starred indices refer to nodes on the "outer" ring of this pair.
165 950 : for (std::size_t n = 0; n < _nodes_per_ring; ++n)
166 : {
167 : // Wrap around
168 900 : unsigned int np1 = (n == _nodes_per_ring - 1) ? 0 : n + 1;
169 900 : Elem * elem = mesh.add_elem(new Tri3);
170 900 : elem->set_node(0, ring_nodes[r][n]);
171 900 : elem->set_node(1, ring_nodes[r + 1][n]);
172 900 : elem->set_node(2, ring_nodes[r][np1]);
173 :
174 : // Add interior faces to 'cylinder' sideset if we are on ring 0.
175 900 : if (r == 0)
176 180 : boundary_info.add_side(elem->id(), /*side=*/2, _cylinder_bid);
177 : }
178 :
179 : // Outer ring (n*, n+1*, n+1)
180 950 : for (std::size_t n = 0; n < _nodes_per_ring; ++n)
181 : {
182 : // Wrap around
183 900 : unsigned int np1 = (n == _nodes_per_ring - 1) ? 0 : n + 1;
184 900 : Elem * elem = mesh.add_elem(new Tri3);
185 900 : elem->set_node(0, ring_nodes[r + 1][n]);
186 900 : elem->set_node(1, ring_nodes[r + 1][np1]);
187 900 : elem->set_node(2, ring_nodes[r][np1]);
188 :
189 : // Add exterior faces to 'exterior' sideset if we're on the last ring.
190 : // Note: this code appears in two places since we could end on either an even or odd ring.
191 900 : if (r == _num_rings - 2)
192 0 : boundary_info.add_side(elem->id(), /*side=*/0, _exterior_bid);
193 : }
194 : }
195 : else
196 : {
197 : // odd -> even ring
198 : // Inner ring (n, n+1*, n+1)
199 950 : for (std::size_t n = 0; n < _nodes_per_ring; ++n)
200 : {
201 : // Wrap around
202 900 : unsigned int np1 = (n == _nodes_per_ring - 1) ? 0 : n + 1;
203 900 : Elem * elem = mesh.add_elem(new Tri3);
204 900 : elem->set_node(0, ring_nodes[r][n]);
205 900 : elem->set_node(1, ring_nodes[r + 1][np1]);
206 900 : elem->set_node(2, ring_nodes[r][np1]);
207 : }
208 :
209 : // Outer ring (n*, n+1*, n)
210 950 : for (std::size_t n = 0; n < _nodes_per_ring; ++n)
211 : {
212 : // Wrap around
213 900 : unsigned int np1 = (n == _nodes_per_ring - 1) ? 0 : n + 1;
214 900 : Elem * elem = mesh.add_elem(new Tri3);
215 900 : elem->set_node(0, ring_nodes[r + 1][n]);
216 900 : elem->set_node(1, ring_nodes[r + 1][np1]);
217 900 : elem->set_node(2, ring_nodes[r][n]);
218 :
219 : // Add exterior faces to 'exterior' sideset if we're on the last ring.
220 900 : if (r == _num_rings - 2)
221 180 : boundary_info.add_side(elem->id(), /*side=*/0, _exterior_bid);
222 : }
223 : }
224 : }
225 :
226 : // Sanity check: make sure all elements have positive area. Note: we
227 : // can't use elem->volume() for this, as that always returns a
228 : // positive area regardless of the node ordering.
229 : // We compute (p1-p0) \cross (p2-p0) and check that the z-component is positive.
230 3610 : for (const auto & elem : mesh.element_ptr_range())
231 : {
232 3600 : Point cp = (elem->point(1) - elem->point(0)).cross(elem->point(2) - elem->point(0));
233 3600 : if (cp(2) < 0.)
234 0 : mooseError("Invalid elem found with negative area");
235 10 : }
236 :
237 : // Create sideset names.
238 10 : boundary_info.sideset_name(_cylinder_bid) = "cylinder";
239 10 : boundary_info.sideset_name(_exterior_bid) = "exterior";
240 :
241 : // Find neighbors, etc.
242 10 : mesh.prepare_for_use();
243 :
244 10 : if (_use_tri6)
245 : {
246 10 : mesh.all_second_order(/*full_ordered=*/true);
247 10 : std::vector<unsigned int> nos;
248 :
249 : // Loop over the elements, moving mid-edge nodes onto the
250 : // nearest radius as applicable. For each element, exactly one
251 : // edge should lie on the same radius, so we move only that
252 : // mid-edge node.
253 3334 : for (const auto & elem : mesh.element_ptr_range())
254 : {
255 : // Make sure we are dealing only with triangles
256 : libmesh_assert(elem->n_vertices() == 3);
257 :
258 : // Compute vertex radii
259 3324 : Real radii[3] = {elem->point(0).norm(), elem->point(1).norm(), elem->point(2).norm()};
260 :
261 : // Compute absolute differences between radii so we can determine which two are on the same
262 : // circular arc.
263 3324 : Real dr[3] = {std::abs(radii[0] - radii[1]),
264 3324 : std::abs(radii[1] - radii[2]),
265 3324 : std::abs(radii[2] - radii[0])};
266 :
267 : // Compute index of minimum dr.
268 3324 : auto index = std::distance(std::begin(dr), std::min_element(std::begin(dr), std::end(dr)));
269 :
270 : // Make sure that the minimum found is also (almost) zero.
271 3324 : if (dr[index] > TOLERANCE)
272 0 : mooseError("Error: element had no sides with nodes on same radius.");
273 :
274 : // Get list of all local node ids on this side. The first
275 : // two entries in nos correspond to the vertices, the last
276 : // entry corresponds to the mid-edge node.
277 3324 : nos = elem->nodes_on_side(index);
278 :
279 : // Compute the angles associated with nodes nos[0] and nos[1].
280 3324 : Real theta0 = std::atan2(elem->point(nos[0])(1), elem->point(nos[0])(0)),
281 3324 : theta1 = std::atan2(elem->point(nos[1])(1), elem->point(nos[1])(0));
282 :
283 : // atan2 returns values in the range (-pi, pi). If theta0
284 : // and theta1 have the same sign, we can simply average them
285 : // to get half of the acute angle between them. On the other
286 : // hand, if theta0 and theta1 are of opposite sign _and_ both
287 : // are larger than pi/2, we need to add 2*pi when averaging,
288 : // otherwise we will get half of the _obtuse_ angle between
289 : // them, and the point will flip to the other side of the
290 : // circle (see below).
291 3324 : Real new_theta = 0.5 * (theta0 + theta1);
292 :
293 : // It should not be possible for both:
294 : // 1.) |theta0| > pi/2, and
295 : // 2.) |theta1| < pi/2
296 : // as this would not be a well-formed element.
297 3504 : if ((theta0 * theta1 < 0) && (std::abs(theta0) > 0.5 * libMesh::pi) &&
298 180 : (std::abs(theta1) > 0.5 * libMesh::pi))
299 180 : new_theta = 0.5 * (theta0 + theta1 + 2 * libMesh::pi);
300 :
301 : // The new radius will be the radius of point nos[0] or nos[1] (they are the same!).
302 3324 : Real new_r = elem->point(nos[0]).norm();
303 :
304 : // Finally, move the point to its new location.
305 3324 : elem->point(nos[2]) = Point(new_r * std::cos(new_theta), new_r * std::sin(new_theta), 0.);
306 10 : }
307 10 : }
308 10 : }
|
