14 #include "libmesh/face_quad4.h" 15 #include "libmesh/face_tri3.h" 24 "inner_radius",
"inner_radius>0.",
"The size of the inner circle.");
27 "The size of the outer circle." 28 " Logically, it has to be greater than inner_radius");
30 "nodes_per_ring",
"nodes_per_ring>5",
"Number of nodes on each ring.");
32 "use_tri6",
false,
"Generate mesh of TRI6 elements instead of TRI3 elements.");
34 "num_rings",
"num_rings>1",
"The number of rings.");
36 "cylinder_bid", 1,
"The boundary id to use for the cylinder (inner circle)");
38 "exterior_bid", 2,
"The boundary id to use for the exterior (outer circle)");
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");
45 " (it can also be TRI6 elements) on several rings.");
52 _inner_radius(getParam<
Real>(
"inner_radius")),
53 _outer_radius(getParam<
Real>(
"outer_radius")),
55 _nodes_per_ring(getParam<unsigned
int>(
"nodes_per_ring")),
56 _use_tri6(getParam<bool>(
"use_tri6")),
57 _num_rings(getParam<unsigned
int>(
"num_rings")),
60 _initial_delta_r(2 *
libMesh::
pi * _inner_radius / _nodes_per_ring)
64 mooseError(
"SpiralAnnularMesh: outer_radius must be greater than inner_radius");
67 std::unique_ptr<MooseMesh>
91 auto newton = [
this, n](
Real & f,
Real & df,
const Real & alpha)
93 f = (1. -
std::pow(alpha, n + 1)) / (1. - alpha) -
95 df = (-(n + 1) * (1 - alpha) *
std::pow(alpha, n) + (1. -
std::pow(alpha, n + 1))) /
96 (1. - alpha) / (1. - alpha);
101 newton(f, df, alpha);
103 while (
std::abs(f) > 1.e-9 && num_iter <= 25)
108 newton(f, df, alpha);
114 mooseError(
"Newton iteration failed to converge (more than 25 iterations).");
126 BoundaryInfo & boundary_info =
mesh.get_boundary_info();
129 std::vector<std::vector<Node *>> ring_nodes(
_num_rings);
136 unsigned int current_node_id = 0;
158 for (std::size_t r = 0; r <
_num_rings - 1; ++r)
169 Elem *
elem =
mesh.add_elem(
new Tri3);
170 elem->set_node(0) = ring_nodes[r][n];
171 elem->set_node(1) = ring_nodes[r + 1][n];
172 elem->set_node(2) = ring_nodes[r][np1];
184 Elem *
elem =
mesh.add_elem(
new Tri3);
185 elem->set_node(0) = ring_nodes[r + 1][n];
186 elem->set_node(1) = ring_nodes[r + 1][np1];
187 elem->set_node(2) = ring_nodes[r][np1];
203 Elem *
elem =
mesh.add_elem(
new Tri3);
204 elem->set_node(0) = ring_nodes[r][n];
205 elem->set_node(1) = ring_nodes[r + 1][np1];
206 elem->set_node(2) = ring_nodes[r][np1];
214 Elem *
elem =
mesh.add_elem(
new Tri3);
215 elem->set_node(0) = ring_nodes[r + 1][n];
216 elem->set_node(1) = ring_nodes[r + 1][np1];
217 elem->set_node(2) = ring_nodes[r][n];
230 for (
const auto &
elem :
mesh.element_ptr_range())
232 Point cp = (
elem->point(1) -
elem->point(0)).cross(
elem->point(2) -
elem->point(0));
234 mooseError(
"Invalid elem found with negative area");
242 mesh.prepare_for_use();
246 mesh.all_second_order(
true);
247 std::vector<unsigned int> nos;
253 for (
const auto &
elem :
mesh.element_ptr_range())
259 Real radii[3] = {
elem->point(0).norm(),
elem->point(1).norm(),
elem->point(2).norm()};
268 auto index = std::distance(std::begin(dr), std::min_element(std::begin(dr), std::end(dr)));
271 if (dr[index] > TOLERANCE)
272 mooseError(
"Error: element had no sides with nodes on same radius.");
277 nos =
elem->nodes_on_side(index);
280 Real theta0 = std::atan2(
elem->point(nos[0])(1),
elem->point(nos[0])(0)),
281 theta1 = std::atan2(
elem->point(nos[1])(1),
elem->point(nos[1])(0));
291 Real new_theta = 0.5 * (theta0 + theta1);
299 new_theta = 0.5 * (theta0 + theta1 + 2 *
libMesh::pi);
302 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...
std::unique_ptr< T > copyConstruct(const T &object)
Copy constructs the object object.
Factory & getFactory()
Retrieve a writable reference to the Factory associated with this App.
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
MooseApp & _app
The MOOSE application this is associated with.
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