13 #include "libmesh/replicated_mesh.h" 14 #include "libmesh/face_quad4.h" 15 #include "libmesh/face_tri3.h" 16 #include "libmesh/mesh_modification.h" 26 "nr", 1,
"nr>0",
"Number of elements in the radial direction");
28 "nt",
"nt>0",
"Number of elements in the angular direction");
32 "Inner radius. If rmin=0 then a disc mesh (with no central hole) will be created.");
35 "radial_positions", {},
"Directly prescribed positions of intermediate radial nodes");
38 "Whether to select the radial discretization " 39 "to achieve equal areas of each ring");
42 "Minimum angle, measured in radians anticlockwise from x axis",
47 "Maximum angle, measured in radians anticlockwise from x axis. If " 48 "tmin=0 and tmax=2Pi an annular mesh is created. " 49 "Otherwise, only a sector of an annulus is created",
52 "dmin", 0.0,
"Minimum degree, measured in degrees anticlockwise from x axis");
55 "Maximum angle, measured in degrees anticlockwise from x axis. If " 56 "dmin=0 and dmax=360 an annular mesh is created. " 57 "Otherwise, only a sector of an annulus is created");
61 "The ratio of radial sizes of successive rings of elements");
63 "quad_subdomain_id", 0,
"The subdomain ID given to the QUAD4 elements");
66 "The subdomain ID given to the TRI3 elements " 67 "(these exist only if rmin=0, and they exist " 68 "at the center of the disc");
69 params.
addClassDescription(
"For rmin>0: creates an annular mesh of QUAD4 elements. For rmin=0: " 70 "creates a disc mesh of QUAD4 and TRI3 elements. Boundary sidesets " 71 "are created at rmax and rmin, and given these names. If dmin!=0 and " 72 "dmax!=360, a sector of an annulus or disc is created. In this case " 73 "boundary sidesets are also created at dmin and dmax, and " 75 params.
addParam<BoundaryName>(
"boundary_name_prefix",
76 "If provided, prefix the built in boundary names with this string");
78 "boundary_id_offset", 0,
"This offset is added to the generated boundary IDs");
85 _nt(getParam<unsigned
int>(
"nt")),
86 _rmin(getParam<
Real>(
"rmin")),
87 _rmax(getParam<
Real>(
"rmax")),
88 _radial_positions(getParam<
std::vector<
Real>>(
"radial_positions")),
89 _nr(parameters.isParamSetByUser(
"radial_positions") ? _radial_positions.size() + 1
90 : getParam<unsigned
int>(
"nr")),
91 _dmin(parameters.isParamSetByUser(
"tmin") ? getParam<
Real>(
"tmin") / M_PI * 180.0
92 : getParam<
Real>(
"dmin")),
93 _dmax(parameters.isParamSetByUser(
"tmax") ? getParam<
Real>(
"tmax") / M_PI * 180.0
94 : getParam<
Real>(
"dmax")),
95 _radians((parameters.isParamSetByUser(
"tmin") || parameters.isParamSetByUser(
"tmax")) ? true
97 _growth_r(getParam<
Real>(
"growth_r")),
98 _len(_growth_r == 1.0 ? (_rmax - _rmin) / _nr
99 : (_rmax - _rmin) * (1.0 -
std::
abs(_growth_r)) /
101 _full_annulus(_dmin == 0.0 && _dmax == 360),
102 _quad_subdomain_id(getParam<
SubdomainID>(
"quad_subdomain_id")),
103 _tri_subdomain_id(getParam<
SubdomainID>(
"tri_subdomain_id")),
104 _equal_area(getParam<bool>(
"equal_area")),
105 _boundary_name_prefix(isParamValid(
"boundary_name_prefix")
106 ? getParam<BoundaryName>(
"boundary_name_prefix") +
"_" 113 "You specified the angles using both degrees and radians. Please use degrees.");
118 paramError(
"nr",
"The 'nr' parameter cannot be specified together with 'radial_positions'");
121 "The 'equal_area' parameter cannot be specified together with 'radial_positions'");
124 "The 'growth_r' parameter cannot be specified together with 'radial_positions'");
126 if (rpos <= _rmin || rpos >=
_rmax)
129 "The following provided value is not within the bounds between 'rmin' and 'rmax': ",
134 paramError(
"growth_r",
"The 'growth_r' parameter cannot be combined with 'equal_area'");
137 paramError(
"rmax",
"rmax must be greater than rmin");
139 paramError(
"dmax",
"dmax must be greater than dmin");
141 paramError(
"dmax",
"dmax - dmin must be <= 360");
144 "nt must be greater than (dmax - dmin) / 180 in order to avoid inverted " 147 paramError(
"quad_subdomain_id",
"quad_subdomain_id must not equal tri_subdomain_id");
150 std::unique_ptr<MeshBase>
158 mesh->set_mesh_dimension(2);
159 mesh->set_spatial_dimension(2);
160 BoundaryInfo & boundary_info =
mesh->get_boundary_info();
163 const unsigned num_nodes =
164 (
_rmin > 0.0 ? (
_nr + 1) * num_angular_nodes :
_nr * num_angular_nodes + 1);
165 const unsigned min_nonzero_layer_num = (
_rmin > 0.0 ? 0 : 1);
166 std::vector<Node *> nodes(num_nodes);
167 unsigned node_id = 0;
171 for (
unsigned angle_num = 0; angle_num < num_angular_nodes; ++angle_num)
173 const Real angle =
_dmin + angle_num * dt;
174 const Real x = current_r *
std::cos(angle * M_PI / 180.0);
175 const Real y = current_r *
std::sin(angle * M_PI / 180.0);
176 nodes[node_id] =
mesh->add_point(Point(x, y, 0.0), node_id);
184 for (
unsigned layer_num =
_nr; layer_num > min_nonzero_layer_num; --layer_num)
195 current_r =
std::sqrt(outer_r * outer_r - ring_area);
208 nodes[node_id] =
mesh->add_point(Point(current_r *
std::cos(
_dmin * M_PI / 180.0),
213 for (
unsigned angle_num = 1; angle_num < num_angular_nodes; ++angle_num)
215 const Real angle =
_dmin + angle_num * dt;
216 const Real x = current_r *
std::cos(angle * M_PI / 180.0);
217 const Real y = current_r *
std::sin(angle * M_PI / 180.0);
218 nodes[node_id] =
mesh->add_point(Point(x, y, 0.0), node_id);
219 Elem * elem =
mesh->add_elem(
new Quad4);
220 elem->set_node(0) = nodes[node_id];
221 elem->set_node(1) = nodes[node_id - 1];
222 elem->set_node(2) = nodes[node_id - num_angular_nodes - 1];
223 elem->set_node(3) = nodes[node_id - num_angular_nodes];
227 if (layer_num ==
_nr)
229 boundary_info.add_side(elem, 2, 1);
232 boundary_info.add_side(elem, 0, 0);
235 boundary_info.add_side(elem, 1, 2);
238 boundary_info.add_side(elem, 3, 3);
243 Elem * elem =
mesh->add_elem(
new Quad4);
244 elem->set_node(0) = nodes[node_id - num_angular_nodes];
245 elem->set_node(1) = nodes[node_id - 1];
246 elem->set_node(2) = nodes[node_id - num_angular_nodes - 1];
247 elem->set_node(3) = nodes[node_id - 2 * num_angular_nodes];
250 if (layer_num ==
_nr)
252 boundary_info.add_side(elem, 2, 1);
255 boundary_info.add_side(elem, 0, 0);
262 nodes[node_id] =
mesh->add_point(Point(0.0, 0.0, 0.0), node_id);
263 boundary_info.add_node(node_id, 0);
264 for (
unsigned angle_num = 0; angle_num < num_angular_nodes - 1; ++angle_num)
266 Elem * elem =
mesh->add_elem(
new Tri3);
267 elem->set_node(0) = nodes[node_id];
268 elem->set_node(1) = nodes[node_id - num_angular_nodes + angle_num];
269 elem->set_node(2) = nodes[node_id - num_angular_nodes + angle_num + 1];
274 Elem * elem =
mesh->add_elem(
new Tri3);
275 elem->set_node(0) = nodes[node_id];
276 elem->set_node(1) = nodes[node_id - 1];
277 elem->set_node(2) = nodes[node_id - num_angular_nodes];
282 boundary_info.sideset_name(0) =
"rmin";
283 boundary_info.sideset_name(1) =
"rmax";
284 boundary_info.nodeset_name(0) =
"rmin";
285 boundary_info.nodeset_name(1) =
"rmax";
290 boundary_info.sideset_name(2) =
"tmin";
291 boundary_info.sideset_name(3) =
"tmax";
292 boundary_info.nodeset_name(2) =
"tmin";
293 boundary_info.nodeset_name(3) =
"tmax";
297 boundary_info.sideset_name(2) =
"dmin";
298 boundary_info.sideset_name(3) =
"dmax";
299 boundary_info.nodeset_name(2) =
"dmin";
300 boundary_info.nodeset_name(3) =
"dmax";
306 const auto mesh_boundary_ids = boundary_info.get_boundary_ids();
307 for (
auto rit = mesh_boundary_ids.rbegin(); rit != mesh_boundary_ids.rend(); ++rit)
310 const std::string old_sideset_name = boundary_info.sideset_name(*rit);
311 const std::string old_nodeset_name = boundary_info.nodeset_name(*rit);
326 mesh->prepare_for_use();
const BoundaryName _boundary_name_prefix
prefix string for the boundary names
MetaPhysicL::DualNumber< V, D, asd > abs(const MetaPhysicL::DualNumber< V, D, asd > &a)
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template * sin(_arg) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(tan
const std::vector< Real > _radial_positions
Radial positions of intermediate rings of nodes (optional)
const unsigned _nt
Number of elements in angular direction.
const Real _dmax
Maximum angle in degrees.
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.
const Real _rmax
Maximum radius.
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
const SubdomainID _quad_subdomain_id
Subdomain ID of created quad elements.
const SubdomainID _tri_subdomain_id
Subdomain ID of created tri elements (that only exist if rmin=0)
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
const bool _full_annulus
Whether a full annulus (as opposed to a sector) will needs to generate.
CTSub CT_OPERATOR_BINARY CTMul CTCompareLess CTCompareGreater CTCompareEqual _arg template cos(_arg) *_arg.template D< dtag >()) CT_SIMPLE_UNARY_FUNCTION(cos
const Real _growth_r
Bias on radial meshing.
const Real _rmin
Minimum radius.
std::unique_ptr< MeshBase > generate() override
Generate / modify the mesh.
void paramError(const std::string ¶m, Args... args) const
Emits an error prefixed with the file and line number of the given param (from the input file) along ...
const boundary_id_type _boundary_id_offset
offset that is added to the boundary IDs
static InputParameters validParams()
registerMooseObject("MooseApp", AnnularMeshGenerator)
static InputParameters validParams()
const bool & _equal_area
Whether to construct rings to have equal areas.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
AnnularMeshGenerator(const InputParameters ¶meters)
const unsigned _nr
Number of elements in radial direction.
const InputParameters & parameters() const
Get the parameters of the object.
const bool _radians
Bool to check if radians are given in the input file.
const Real _len
rmax = rmin + len + len*g + len*g^2 + len*g^3 + ... + len*g^(nr-1) = rmin + len*(1 - g^nr)/(1 - g) ...
std::unique_ptr< MeshBase > buildMeshBaseObject(unsigned int dim=libMesh::invalid_uint)
Build a MeshBase object whose underlying type will be determined by the Mesh input file block...
Generates an annular mesh given all the parameters.
const Real _dmin
Minimum angle in degrees.
MooseUnits pow(const MooseUnits &, int)
MeshGenerators are objects that can modify or add to an existing mesh.
void ErrorVector unsigned int