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