12 #include "libmesh/face_quad4.h" 13 #include "libmesh/face_tri3.h" 22 "inner_radius",
"inner_radius>0.",
"The size of the inner circle.");
25 "The size of the outer circle." 26 " Logically, it has to be greater than inner_radius");
28 "nodes_per_ring",
"nodes_per_ring>5",
"Number of nodes on each ring.");
30 "use_tri6",
false,
"Generate mesh of TRI6 elements instead of TRI3 elements.");
32 "num_rings",
"num_rings>1",
"The number of rings.");
34 "cylinder_bid", 1,
"The boundary id to use for the cylinder (inner circle)");
36 "exterior_bid", 2,
"The boundary id to use for the exterior (outer circle)");
38 "Width of the initial layer of elements around the cylinder." 39 "This number should be approximately" 40 " 2 * pi * inner_radius / nodes_per_ring to ensure that the" 41 " initial layer of elements is almost equilateral");
43 " (it can also be TRI6 elements) on several rings.");
50 _inner_radius(getParam<
Real>(
"inner_radius")),
51 _outer_radius(getParam<
Real>(
"outer_radius")),
53 _nodes_per_ring(getParam<unsigned
int>(
"nodes_per_ring")),
54 _use_tri6(getParam<bool>(
"use_tri6")),
55 _num_rings(getParam<unsigned
int>(
"num_rings")),
58 _initial_delta_r(2 *
libMesh::
pi * _inner_radius / _nodes_per_ring)
62 mooseError(
"SpiralAnnularMesh: outer_radius must be greater than inner_radius");
65 std::unique_ptr<MooseMesh>
68 return std::make_unique<SpiralAnnularMesh>(*this);
89 auto newton = [
this, n](
Real & f,
Real & df,
const Real & alpha)
91 f = (1. -
std::pow(alpha, n + 1)) / (1. - alpha) -
93 df = (-(n + 1) * (1 - alpha) *
std::pow(alpha, n) + (1. -
std::pow(alpha, n + 1))) /
94 (1. - alpha) / (1. - alpha);
101 while (
std::abs(f) > 1.e-9 && num_iter <= 25)
106 newton(f, df, alpha);
112 mooseError(
"Newton iteration failed to converge (more than 25 iterations).");
124 BoundaryInfo & boundary_info =
mesh.get_boundary_info();
127 std::vector<std::vector<Node *>> ring_nodes(
_num_rings);
134 unsigned int current_node_id = 0;
156 for (std::size_t r = 0; r <
_num_rings - 1; ++r)
167 Elem *
elem =
mesh.add_elem(
new Tri3);
168 elem->set_node(0) = ring_nodes[r][n];
169 elem->set_node(1) = ring_nodes[r + 1][n];
170 elem->set_node(2) = ring_nodes[r][np1];
182 Elem *
elem =
mesh.add_elem(
new Tri3);
183 elem->set_node(0) = ring_nodes[r + 1][n];
184 elem->set_node(1) = ring_nodes[r + 1][np1];
185 elem->set_node(2) = ring_nodes[r][np1];
201 Elem *
elem =
mesh.add_elem(
new Tri3);
202 elem->set_node(0) = ring_nodes[r][n];
203 elem->set_node(1) = ring_nodes[r + 1][np1];
204 elem->set_node(2) = ring_nodes[r][np1];
212 Elem *
elem =
mesh.add_elem(
new Tri3);
213 elem->set_node(0) = ring_nodes[r + 1][n];
214 elem->set_node(1) = ring_nodes[r + 1][np1];
215 elem->set_node(2) = ring_nodes[r][n];
228 for (
const auto &
elem :
mesh.element_ptr_range())
230 Point cp = (
elem->point(1) -
elem->point(0)).cross(
elem->point(2) -
elem->point(0));
232 mooseError(
"Invalid elem found with negative area");
240 mesh.prepare_for_use();
244 mesh.all_second_order(
true);
245 std::vector<unsigned int> nos;
251 for (
const auto &
elem :
mesh.element_ptr_range())
257 Real radii[3] = {
elem->point(0).norm(),
elem->point(1).norm(),
elem->point(2).norm()};
266 auto index = std::distance(std::begin(dr), std::min_element(std::begin(dr), std::end(dr)));
269 if (dr[index] > TOLERANCE)
270 mooseError(
"Error: element had no sides with nodes on same radius.");
275 nos =
elem->nodes_on_side(index);
278 Real theta0 = std::atan2(
elem->point(nos[0])(1),
elem->point(nos[0])(0)),
279 theta1 = std::atan2(
elem->point(nos[1])(1),
elem->point(nos[1])(0));
289 Real new_theta = 0.5 * (theta0 + theta1);
297 new_theta = 0.5 * (theta0 + theta1 + 2 *
libMesh::pi);
300 Real new_r =
elem->point(nos[0]).norm();
static InputParameters validParams()
Typical "Moose-style" constructor and copy constructor.
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template * sin(_arg) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(tan
registerMooseObject("MooseApp", SpiralAnnularMesh)
virtual std::unique_ptr< MooseMesh > safeClone() const override
A safer version of the clone() method that hands back an allocated object wrapped in a smart pointer...
const Real _outer_radius
Radius of the outer circle. Logically, it's bigger that inner_radius.
const Real _inner_radius
Radius of the inner circle.
Real _radial_bias
Factor to increase initial_delta_r for each ring.
const boundary_id_type _exterior_bid
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
unsigned int _num_rings
Number of rings.You can't specify both the number of rings and the radial bias if you want to match a...
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template cos(_arg) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(cos
MeshBase & getMesh()
Accessor for the underlying libMesh Mesh object.
SpiralAnnularMesh(const InputParameters ¶meters)
MooseMesh wraps a libMesh::Mesh object and enhances its capabilities by caching additional data and s...
const Real _initial_delta_r
static InputParameters validParams()
virtual void buildMesh() override
Must be overridden by child classes.
const boundary_id_type _cylinder_bid
The boundary id to use for the cylinder.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void mooseError(Args &&... args) const
Emits an error prefixed with object name and type.
Mesh generated from parameters.
const bool _use_tri6
Generate mesh of TRI6 elements instead of TRI3 elements.
virtual Elem * elem(const dof_id_type i)
Various accessors (pointers/references) for Elem "i".
const unsigned int _nodes_per_ring
Number of nodes on each ring.
MooseUnits pow(const MooseUnits &, int)
void ErrorVector unsigned int