13 #include "libmesh/replicated_mesh.h" 14 #include "libmesh/face_quad4.h" 15 #include "libmesh/face_tri3.h" 25 "inner_radius",
"inner_radius>0.",
"The size of the inner circle.");
28 "The size of the outer circle." 29 " Logically, it has to be greater than inner_radius");
31 "nodes_per_ring",
"nodes_per_ring>5",
"Number of nodes on each ring.");
33 "use_tri6",
false,
"Generate mesh of TRI6 elements instead of TRI3 elements.");
35 "num_rings",
"num_rings>1",
"The number of rings.");
37 "cylinder_bid", 1,
"The boundary id to use for the cylinder (inner circle)");
39 "exterior_bid", 2,
"The boundary id to use for the exterior (outer circle)");
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");
46 "Creates an annular mesh based on TRI3 or TRI6 elements on several rings.");
53 _inner_radius(getParam<
Real>(
"inner_radius")),
54 _outer_radius(getParam<
Real>(
"outer_radius")),
56 _nodes_per_ring(getParam<unsigned
int>(
"nodes_per_ring")),
57 _use_tri6(getParam<bool>(
"use_tri6")),
58 _num_rings(getParam<unsigned
int>(
"num_rings")),
61 _initial_delta_r(2 *
libMesh::
pi * _inner_radius / _nodes_per_ring)
67 mooseError(
"SpiralAnnularMesh: outer_radius must be greater than inner_radius");
70 std::unique_ptr<MeshBase>
74 BoundaryInfo & boundary_info =
mesh->get_boundary_info();
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).");
125 std::vector<std::vector<Node *>> ring_nodes(
_num_rings);
132 unsigned int current_node_id = 0;
154 for (std::size_t r = 0; r <
_num_rings - 1; ++r)
165 Elem * elem =
mesh->add_elem(
new Tri3);
166 elem->set_node(0) = ring_nodes[r][n];
167 elem->set_node(1) = ring_nodes[r + 1][n];
168 elem->set_node(2) = ring_nodes[r][np1];
180 Elem * elem =
mesh->add_elem(
new Tri3);
181 elem->set_node(0) = ring_nodes[r + 1][n];
182 elem->set_node(1) = ring_nodes[r + 1][np1];
183 elem->set_node(2) = ring_nodes[r][np1];
199 Elem * elem =
mesh->add_elem(
new Tri3);
200 elem->set_node(0) = ring_nodes[r][n];
201 elem->set_node(1) = ring_nodes[r + 1][np1];
202 elem->set_node(2) = ring_nodes[r][np1];
210 Elem * elem =
mesh->add_elem(
new Tri3);
211 elem->set_node(0) = ring_nodes[r + 1][n];
212 elem->set_node(1) = ring_nodes[r + 1][np1];
213 elem->set_node(2) = ring_nodes[r][n];
226 for (
const auto & elem :
mesh->element_ptr_range())
228 Point cp = (elem->point(1) - elem->point(0)).cross(elem->point(2) - elem->point(0));
230 mooseError(
"Invalid elem found with negative area");
238 mesh->prepare_for_use();
242 mesh->all_second_order(
true);
243 std::vector<unsigned int> nos;
249 for (
const auto & elem :
mesh->element_ptr_range())
255 Real radii[3] = {elem->point(0).norm(), elem->point(1).norm(), elem->point(2).norm()};
264 auto index = std::distance(std::begin(dr), std::min_element(std::begin(dr), std::end(dr)));
267 if (dr[index] > TOLERANCE)
268 mooseError(
"Error: element had no sides with nodes on same radius.");
273 nos = elem->nodes_on_side(index);
276 Real theta0 = std::atan2(elem->point(nos[0])(1), elem->point(nos[0])(0)),
277 theta1 = std::atan2(elem->point(nos[1])(1), elem->point(nos[1])(0));
287 Real new_theta = 0.5 * (theta0 + theta1);
295 new_theta = 0.5 * (theta0 + theta1 + 2 *
libMesh::pi);
298 Real new_r = elem->point(nos[0]).norm();
301 elem->point(nos[2]) = Point(new_r *
std::cos(new_theta), new_r *
std::sin(new_theta), 0.);
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template * sin(_arg) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(tan
std::unique_ptr< MeshBase > generate() override
Generate / modify the mesh.
std::unique_ptr< ReplicatedMesh > buildReplicatedMesh(unsigned int dim=libMesh::invalid_uint)
Build a replicated mesh.
const Real _outer_radius
Radius of the outer circle. Logically, it's bigger that inner_radius.
const boundary_id_type _cylinder_bid
The boundary id to use for the cylinder.
Real _radial_bias
Factor to increase initial_delta_r for each ring.
const bool _use_tri6
Generate mesh of TRI6 elements instead of TRI3 elements.
std::unique_ptr< T_DEST, T_DELETER > dynamic_pointer_cast(std::unique_ptr< T_SRC, T_DELETER > &src)
These are reworked from https://stackoverflow.com/a/11003103.
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
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...
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
SpiralAnnularMeshGenerator(const InputParameters ¶meters)
registerMooseObject("MooseApp", SpiralAnnularMeshGenerator)
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template cos(_arg) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(cos
Generates a spiral annular mesh given all the parameters.
const Real _inner_radius
Radius of the inner circle.
static InputParameters validParams()
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static InputParameters validParams()
const Real _initial_delta_r
Width of the initial layer of elements around the cylinder.
void mooseError(Args &&... args) const
Emits an error prefixed with object name and type.
T & declareMeshProperty(const std::string &data_name, Args &&... args)
Methods for writing out attributes to the mesh meta-data store, which can be retrieved from most othe...
const unsigned int _nodes_per_ring
Number of nodes on each ring.
const boundary_id_type _exterior_bid
MooseUnits pow(const MooseUnits &, int)
MeshGenerators are objects that can modify or add to an existing mesh.
void ErrorVector unsigned int