www.mooseframework.org
AnnularMeshGenerator.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 "AnnularMeshGenerator.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 
18 
19 template <>
22 {
24 
25  params.addRangeCheckedParam<unsigned int>(
26  "nr", 1, "nr>0", "Number of elements in the radial direction");
27  params.addRequiredRangeCheckedParam<unsigned int>(
28  "nt", "nt>0", "Number of elements in the angular direction");
29  params.addRequiredRangeCheckedParam<Real>(
30  "rmin",
31  "rmin>=0.0",
32  "Inner radius. If rmin=0 then a disc mesh (with no central hole) will be created.");
33  params.addRequiredParam<Real>("rmax", "Outer radius");
34  params.addParam<Real>("tmin", 0.0, "Minimum angle, measured anticlockwise from x axis");
35  params.addParam<Real>("tmax",
36  2 * M_PI,
37  "Maximum angle, measured anticlockwise from x axis. If "
38  "tmin=0 and tmax=2Pi an annular mesh is created. "
39  "Otherwise, only a sector of an annulus is created");
40  params.addRangeCheckedParam<Real>(
41  "growth_r", 1.0, "growth_r>0.0", "The ratio of radial sizes of successive rings of elements");
42  params.addParam<SubdomainID>(
43  "quad_subdomain_id", 0, "The subdomain ID given to the QUAD4 elements");
44  params.addParam<SubdomainID>("tri_subdomain_id",
45  1,
46  "The subdomain ID given to the TRI3 elements "
47  "(these exist only if rmin=0, and they exist "
48  "at the center of the disc");
49  params.addClassDescription("For rmin>0: creates an annular mesh of QUAD4 elements. For rmin=0: "
50  "creates a disc mesh of QUAD4 and TRI3 elements. Boundary sidesets "
51  "are created at rmax and rmin, and given these names. If tmin!=0 and "
52  "tmax!=2Pi, a sector of an annulus or disc is created. In this case "
53  "boundary sidesets are also created a tmin and tmax, and "
54  "given these names");
55 
56  return params;
57 }
58 
60  : MeshGenerator(parameters),
61  _nr(getParam<unsigned int>("nr")),
62  _nt(getParam<unsigned int>("nt")),
63  _rmin(getParam<Real>("rmin")),
64  _rmax(getParam<Real>("rmax")),
65  _tmin(getParam<Real>("tmin")),
66  _tmax(getParam<Real>("tmax")),
67  _growth_r(getParam<Real>("growth_r")),
68  _len(_growth_r == 1.0 ? (_rmax - _rmin) / _nr
69  : (_rmax - _rmin) * (1.0 - _growth_r) / (1.0 - std::pow(_growth_r, _nr))),
70  _full_annulus(_tmin == 0.0 && _tmax == 2 * M_PI),
71  _quad_subdomain_id(getParam<SubdomainID>("quad_subdomain_id")),
72  _tri_subdomain_id(getParam<SubdomainID>("tri_subdomain_id"))
73 {
74  // catch likely user errors
75  if (_rmax <= _rmin)
76  mooseError("AnnularMesh: rmax must be greater than rmin");
77  if (_tmax <= _tmin)
78  mooseError("AnnularMesh: tmax must be greater than tmin");
79  if (_tmax - _tmin > 2 * M_PI)
80  mooseError("AnnularMesh: tmax - tmin must be <= 2 Pi");
81  if (_nt <= (_tmax - _tmin) / M_PI)
82  mooseError("AnnularMesh: nt must be greater than (tmax - tmin) / Pi in order to avoid inverted "
83  "elements");
85  mooseError("AnnularMesh: quad_subdomain_id must not equal tri_subdomain_id");
86 }
87 
88 std::unique_ptr<MeshBase>
90 {
91  // Have MOOSE construct the correct libMesh::Mesh object using Mesh block and CLI parameters.
92  auto mesh = _mesh->buildMeshBaseObject();
93 
94  const Real dt = (_tmax - _tmin) / _nt;
95 
96  mesh->set_mesh_dimension(2);
97  mesh->set_spatial_dimension(2);
98  BoundaryInfo & boundary_info = mesh->get_boundary_info();
99 
100  const unsigned num_angular_nodes = (_full_annulus ? _nt : _nt + 1);
101  const unsigned num_nodes =
102  (_rmin > 0.0 ? (_nr + 1) * num_angular_nodes : _nr * num_angular_nodes + 1);
103  const unsigned min_nonzero_layer_num = (_rmin > 0.0 ? 0 : 1);
104  std::vector<Node *> nodes(num_nodes);
105  unsigned node_id = 0;
106 
107  // add nodes at rmax that aren't yet connected to any elements
108  Real current_r = _rmax;
109  for (unsigned angle_num = 0; angle_num < num_angular_nodes; ++angle_num)
110  {
111  const Real angle = _tmin + angle_num * dt;
112  const Real x = current_r * std::cos(angle);
113  const Real y = current_r * std::sin(angle);
114  nodes[node_id] = mesh->add_point(Point(x, y, 0.0), node_id);
115  node_id++;
116  }
117 
118  // add nodes at smaller radii, and connect them with elements
119  for (unsigned layer_num = _nr; layer_num > min_nonzero_layer_num; --layer_num)
120  {
121  if (layer_num == 1)
122  current_r = _rmin; // account for precision loss
123  else
124  current_r -= _len * std::pow(_growth_r, layer_num - 1);
125 
126  // add node at angle = _tmin
127  nodes[node_id] = mesh->add_point(
128  Point(current_r * std::cos(_tmin), current_r * std::sin(_tmin), 0.0), node_id);
129  node_id++;
130  for (unsigned angle_num = 1; angle_num < num_angular_nodes; ++angle_num)
131  {
132  const Real angle = _tmin + angle_num * dt;
133  const Real x = current_r * std::cos(angle);
134  const Real y = current_r * std::sin(angle);
135  nodes[node_id] = mesh->add_point(Point(x, y, 0.0), node_id);
136  Elem * elem = mesh->add_elem(new Quad4);
137  elem->set_node(0) = nodes[node_id];
138  elem->set_node(1) = nodes[node_id - 1];
139  elem->set_node(2) = nodes[node_id - num_angular_nodes - 1];
140  elem->set_node(3) = nodes[node_id - num_angular_nodes];
141  elem->subdomain_id() = _quad_subdomain_id;
142  node_id++;
143 
144  if (layer_num == _nr)
145  // add outer boundary (boundary_id = 1)
146  boundary_info.add_side(elem, 2, 1);
147  if (layer_num == 1)
148  // add inner boundary (boundary_id = 0)
149  boundary_info.add_side(elem, 0, 0);
150  if (!_full_annulus && angle_num == 1)
151  // add tmin boundary (boundary_id = 2)
152  boundary_info.add_side(elem, 1, 2);
153  if (!_full_annulus && angle_num == num_angular_nodes - 1)
154  // add tmin boundary (boundary_id = 3)
155  boundary_info.add_side(elem, 3, 3);
156  }
157  if (_full_annulus)
158  {
159  // add element connecting to node at angle=0
160  Elem * elem = mesh->add_elem(new Quad4);
161  elem->set_node(0) = nodes[node_id - num_angular_nodes];
162  elem->set_node(1) = nodes[node_id - 1];
163  elem->set_node(2) = nodes[node_id - num_angular_nodes - 1];
164  elem->set_node(3) = nodes[node_id - 2 * num_angular_nodes];
165  elem->subdomain_id() = _quad_subdomain_id;
166 
167  if (layer_num == _nr)
168  // add outer boundary (boundary_id = 1)
169  boundary_info.add_side(elem, 2, 1);
170  if (layer_num == 1)
171  // add inner boundary (boundary_id = 0)
172  boundary_info.add_side(elem, 0, 0);
173  }
174  }
175 
176  // add single node at origin, if relevant
177  if (_rmin == 0.0)
178  {
179  nodes[node_id] = mesh->add_point(Point(0.0, 0.0, 0.0), node_id);
180  boundary_info.add_node(node_id, 0); // boundary_id=0 is centre
181  for (unsigned angle_num = 0; angle_num < num_angular_nodes - 1; ++angle_num)
182  {
183  Elem * elem = mesh->add_elem(new Tri3);
184  elem->set_node(0) = nodes[node_id];
185  elem->set_node(1) = nodes[node_id - num_angular_nodes + angle_num];
186  elem->set_node(2) = nodes[node_id - num_angular_nodes + angle_num + 1];
187  elem->subdomain_id() = _tri_subdomain_id;
188  }
189  if (_full_annulus)
190  {
191  Elem * elem = mesh->add_elem(new Tri3);
192  elem->set_node(0) = nodes[node_id];
193  elem->set_node(1) = nodes[node_id - 1];
194  elem->set_node(2) = nodes[node_id - num_angular_nodes];
195  elem->subdomain_id() = _tri_subdomain_id;
196  }
197  }
198 
199  boundary_info.sideset_name(0) = "rmin";
200  boundary_info.sideset_name(1) = "rmax";
201  boundary_info.nodeset_name(0) = "rmin";
202  boundary_info.nodeset_name(1) = "rmax";
203  if (!_full_annulus)
204  {
205  boundary_info.sideset_name(2) = "tmin";
206  boundary_info.sideset_name(3) = "tmax";
207  boundary_info.nodeset_name(2) = "tmin";
208  boundary_info.nodeset_name(3) = "tmax";
209  }
210 
211  mesh->prepare_for_use(false);
212 
213  return dynamic_pointer_cast<MeshBase>(mesh);
214 }
void addRequiredRangeCheckedParam(const std::string &name, const std::string &parsed_function, const std::string &doc_string)
These methods add an range checked parameters.
const unsigned _nt
Number of elements in angular direction.
const Real _tmin
Minimum angle.
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
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.
void mooseError(Args &&... args) const
Definition: MooseObject.h:147
InputParameters validParams< AnnularMeshGenerator >()
static PetscErrorCode Vec x
const SubdomainID _quad_subdomain_id
Subdomain ID of created quad elements.
void addRequiredParam(const std::string &name, const std::string &doc_string)
This method adds a parameter and documentation string to the InputParameters object that will be extr...
const SubdomainID _tri_subdomain_id
Subdomain ID of created tri elements (that only exist if rmin=0)
const bool _full_annulus
Whether a full annulus (as opposed to a sector) will needs to generate.
Real pow(Real x, int e)
Definition: MathUtils.C:211
const Real _growth_r
Bias on radial meshing.
const Real _rmin
Minimum radius.
std::unique_ptr< MeshBase > generate() override
Generate / modify the mesh.
subdomain_id_type SubdomainID
registerMooseObject("MooseApp", AnnularMeshGenerator)
const Real _tmax
Maximum angle.
AnnularMeshGenerator(const InputParameters &parameters)
const unsigned _nr
Number of elements in radial direction.
std::shared_ptr< MooseMesh > & _mesh
References to the mesh and displaced mesh (currently in the ActionWarehouse)
Definition: MeshGenerator.h:72
void addClassDescription(const std::string &doc_string)
This method adds a description of the class that will be displayed in the input file syntax dump...
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an option parameter and a documentation string to the InputParameters object...
void addRangeCheckedParam(const std::string &name, const T &value, const std::string &parsed_function, const std::string &doc_string)
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) ...
Generates an annular mesh given all the parameters.
InputParameters validParams< MeshGenerator >()
Definition: MeshGenerator.C:16
MeshGenerators are objects that can modify or add to an existing mesh.
Definition: MeshGenerator.h:30