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 #include "libmesh/mesh_modification.h"
17 
19 
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");
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<std::vector<Real>>(
35  "radial_positions", {}, "Directly prescribed positions of intermediate radial nodes");
36  params.addParam<bool>("equal_area",
37  false,
38  "Whether to select the radial discretization "
39  "to achieve equal areas of each ring");
40  params.addDeprecatedParam<Real>("tmin",
41  0.0,
42  "Minimum angle, measured in radians anticlockwise from x axis",
43  "Use dmin instead");
44  params.addDeprecatedParam<Real>(
45  "tmax",
46  2 * M_PI,
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",
50  "Use dmin instead");
51  params.addParam<Real>(
52  "dmin", 0.0, "Minimum degree, measured in degrees anticlockwise from x axis");
53  params.addParam<Real>("dmax",
54  360.0,
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");
58  params.addRangeCheckedParam<Real>("growth_r",
59  1.0,
60  "growth_r!=0.0",
61  "The ratio of radial sizes of successive rings of elements");
62  params.addParam<SubdomainID>(
63  "quad_subdomain_id", 0, "The subdomain ID given to the QUAD4 elements");
64  params.addParam<SubdomainID>("tri_subdomain_id",
65  1,
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 "
74  "given these names");
75  params.addParam<BoundaryName>("boundary_name_prefix",
76  "If provided, prefix the built in boundary names with this string");
77  params.addParam<boundary_id_type>(
78  "boundary_id_offset", 0, "This offset is added to the generated boundary IDs");
79 
80  return params;
81 }
82 
84  : MeshGenerator(parameters),
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
96  : false),
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)) /
100  (1.0 - std::pow(std::abs(_growth_r), _nr))),
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") + "_"
107  : ""),
108  _boundary_id_offset(getParam<boundary_id_type>("boundary_id_offset"))
109 {
110  if ((parameters.isParamSetByUser("tmin") || parameters.isParamSetByUser("tmax")) &&
112  paramError("tmin",
113  "You specified the angles using both degrees and radians. Please use degrees.");
114 
115  if (_radial_positions.size() != 0)
116  {
117  if (parameters.isParamSetByUser("nr"))
118  paramError("nr", "The 'nr' parameter cannot be specified together with 'radial_positions'");
119  if (parameters.isParamSetByUser("equal_area"))
120  paramError("equal_area",
121  "The 'equal_area' parameter cannot be specified together with 'radial_positions'");
122  if (parameters.isParamSetByUser("growth_r"))
123  paramError("growth_r",
124  "The 'growth_r' parameter cannot be specified together with 'radial_positions'");
125  for (auto rpos : _radial_positions)
126  if (rpos <= _rmin || rpos >= _rmax)
127  paramError(
128  "radial_positions",
129  "The following provided value is not within the bounds between 'rmin' and 'rmax': ",
130  rpos);
131  }
132 
133  if (_equal_area && parameters.isParamSetByUser("growth_r"))
134  paramError("growth_r", "The 'growth_r' parameter cannot be combined with 'equal_area'");
135 
136  if (_rmax <= _rmin)
137  paramError("rmax", "rmax must be greater than rmin");
138  if (_dmax <= _dmin)
139  paramError("dmax", "dmax must be greater than dmin");
140  if (_dmax - _dmin > 360)
141  paramError("dmax", "dmax - dmin must be <= 360");
142  if (_nt <= (_dmax - _dmin) / 180.0)
143  paramError("nt",
144  "nt must be greater than (dmax - dmin) / 180 in order to avoid inverted "
145  "elements");
147  paramError("quad_subdomain_id", "quad_subdomain_id must not equal tri_subdomain_id");
148 }
149 
150 std::unique_ptr<MeshBase>
152 {
153  // Have MOOSE construct the correct libMesh::Mesh object using Mesh block and CLI parameters.
154  auto mesh = buildMeshBaseObject();
155 
156  const Real dt = (_dmax - _dmin) / _nt;
157 
158  mesh->set_mesh_dimension(2);
159  mesh->set_spatial_dimension(2);
160  BoundaryInfo & boundary_info = mesh->get_boundary_info();
161 
162  const unsigned num_angular_nodes = (_full_annulus ? _nt : _nt + 1);
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;
168 
169  // add nodes at rmax that aren't yet connected to any elements
170  Real current_r = _rmax;
171  for (unsigned angle_num = 0; angle_num < num_angular_nodes; ++angle_num)
172  {
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);
177  node_id++;
178  }
179 
180  // first value for the ring outer radius, only used for _equal_area
181  Real outer_r = _rmax;
182 
183  // add nodes at smaller radii, and connect them with elements
184  for (unsigned layer_num = _nr; layer_num > min_nonzero_layer_num; --layer_num)
185  {
186  if (layer_num == 1)
187  current_r = _rmin; // account for precision loss
188  else if (_radial_positions.size() > 0)
189  current_r = _radial_positions[layer_num - 2];
190  else
191  {
192  if (_equal_area)
193  {
194  const Real ring_area = (_rmax * _rmax - _rmin * _rmin) / _nr;
195  current_r = std::sqrt(outer_r * outer_r - ring_area);
196  outer_r = current_r;
197  }
198  else
199  {
200  if (_growth_r > 0)
201  current_r -= _len * std::pow(_growth_r, layer_num - 1);
202  else
203  current_r -= _len * std::pow(std::abs(_growth_r), _nr - layer_num);
204  }
205  }
206 
207  // add node at angle = _dmin
208  nodes[node_id] = mesh->add_point(Point(current_r * std::cos(_dmin * M_PI / 180.0),
209  current_r * std::sin(_dmin * M_PI / 180.0),
210  0.0),
211  node_id);
212  node_id++;
213  for (unsigned angle_num = 1; angle_num < num_angular_nodes; ++angle_num)
214  {
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];
224  elem->subdomain_id() = _quad_subdomain_id;
225  node_id++;
226 
227  if (layer_num == _nr)
228  // add outer boundary (boundary_id = 1)
229  boundary_info.add_side(elem, 2, 1);
230  if (layer_num == 1)
231  // add inner boundary (boundary_id = 0)
232  boundary_info.add_side(elem, 0, 0);
233  if (!_full_annulus && angle_num == 1)
234  // add tmin boundary (boundary_id = 2)
235  boundary_info.add_side(elem, 1, 2);
236  if (!_full_annulus && angle_num == num_angular_nodes - 1)
237  // add tmin boundary (boundary_id = 3)
238  boundary_info.add_side(elem, 3, 3);
239  }
240  if (_full_annulus)
241  {
242  // add element connecting to node at angle=0
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];
248  elem->subdomain_id() = _quad_subdomain_id;
249 
250  if (layer_num == _nr)
251  // add outer boundary (boundary_id = 1)
252  boundary_info.add_side(elem, 2, 1);
253  if (layer_num == 1)
254  // add inner boundary (boundary_id = 0)
255  boundary_info.add_side(elem, 0, 0);
256  }
257  }
258 
259  // add single node at origin, if relevant
260  if (_rmin == 0.0)
261  {
262  nodes[node_id] = mesh->add_point(Point(0.0, 0.0, 0.0), node_id);
263  boundary_info.add_node(node_id, 0); // boundary_id=0 is centre
264  for (unsigned angle_num = 0; angle_num < num_angular_nodes - 1; ++angle_num)
265  {
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];
270  elem->subdomain_id() = _tri_subdomain_id;
271  }
272  if (_full_annulus)
273  {
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];
278  elem->subdomain_id() = _tri_subdomain_id;
279  }
280  }
281 
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";
286  if (!_full_annulus)
287  {
288  if (_radians)
289  {
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";
294  }
295  else
296  {
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";
301  }
302  }
303  if (_boundary_id_offset != 0 || !_boundary_name_prefix.empty())
304  {
305  // apply boundary id offset and name prefix
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)
308  {
309 
310  const std::string old_sideset_name = boundary_info.sideset_name(*rit);
311  const std::string old_nodeset_name = boundary_info.nodeset_name(*rit);
312 
313  if (_boundary_id_offset != 0)
314  MeshTools::Modification::change_boundary_id(*mesh, *rit, *rit + _boundary_id_offset);
315 
316  if (!_boundary_name_prefix.empty())
317  {
318  boundary_info.sideset_name(*rit + _boundary_id_offset) =
319  _boundary_name_prefix + old_sideset_name;
320  boundary_info.nodeset_name(*rit + _boundary_id_offset) =
321  _boundary_name_prefix + old_nodeset_name;
322  }
323  }
324  }
325 
326  mesh->prepare_for_use();
327  return dynamic_pointer_cast<MeshBase>(mesh);
328 }
const BoundaryName _boundary_name_prefix
prefix string for the boundary names
MetaPhysicL::DualNumber< V, D, asd > abs(const MetaPhysicL::DualNumber< V, D, asd > &a)
Definition: ADReal.h:42
void addRequiredRangeCheckedParam(const std::string &name, const std::string &parsed_function, const std::string &doc_string)
These methods add an range checked parameters.
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)
void addDeprecatedParam(const std::string &name, const T &value, const std::string &doc_string, const std::string &deprecation_message)
const unsigned _nt
Number of elements in angular direction.
MeshBase & mesh
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
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.
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)
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.
T pow(const T &x)
int8_t boundary_id_type
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 &param, 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()
Definition: MeshGenerator.C:23
registerMooseObject("MooseApp", AnnularMeshGenerator)
bool isParamSetByUser(const std::string &name) const
Method returns true if the parameter was by the user.
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 &parameters)
const unsigned _nr
Number of elements in radial direction.
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...
const InputParameters & parameters() const
Get the parameters of the object.
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 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)
Definition: Units.C:537
MeshGenerators are objects that can modify or add to an existing mesh.
Definition: MeshGenerator.h:32
void ErrorVector unsigned int