www.mooseframework.org
AnnularMesh.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
15 #include "AnnularMesh.h"
16 
17 #include "libmesh/face_quad4.h"
18 #include "libmesh/face_tri3.h"
19 
20 template <>
23 {
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  return params;
56 }
57 
59  : MooseMesh(parameters),
60  _nr(getParam<unsigned int>("nr")),
61  _nt(getParam<unsigned int>("nt")),
62  _rmin(getParam<Real>("rmin")),
63  _rmax(getParam<Real>("rmax")),
64  _tmin(getParam<Real>("tmin")),
65  _tmax(getParam<Real>("tmax")),
66  _growth_r(getParam<Real>("growth_r")),
67  _len(_growth_r == 1.0 ? (_rmax - _rmin) / _nr
68  : (_rmax - _rmin) * (1.0 - _growth_r) / (1.0 - std::pow(_growth_r, _nr))),
69  _full_annulus(_tmin == 0.0 && _tmax == 2 * M_PI),
70  _quad_subdomain_id(getParam<SubdomainID>("quad_subdomain_id")),
71  _tri_subdomain_id(getParam<SubdomainID>("tri_subdomain_id"))
72 {
73  // catch likely user errors
74  if (_rmax <= _rmin)
75  mooseError("AnnularMesh: rmax must be greater than rmin");
76  if (_tmax <= _tmin)
77  mooseError("AnnularMesh: tmax must be greater than tmin");
78  if (_tmax - _tmin > 2 * M_PI)
79  mooseError("AnnularMesh: tmax - tmin must be <= 2 Pi");
80  if (_nt <= (_tmax - _tmin) / M_PI)
81  mooseError("AnnularMesh: nt must be greater than (tmax - tmin) / Pi in order to avoid inverted "
82  "elements");
84  mooseError("AnnularMesh: quad_subdomain_id must not equal tri_subdomain_id");
85 }
86 
87 Real
88 AnnularMesh::getMinInDimension(unsigned int component) const
89 {
90  switch (component)
91  {
92  case 0:
93  return -_rmax;
94  case 1:
95  return -_rmax;
96  case 2:
97  return 0.0;
98  default:
99  mooseError("Invalid component");
100  }
101 }
102 
103 Real
104 AnnularMesh::getMaxInDimension(unsigned int component) const
105 {
106  switch (component)
107  {
108  case 0:
109  return _rmax;
110  case 1:
111  return _rmax;
112  case 2:
113  return 0.0;
114  default:
115  mooseError("Invalid component");
116  }
117 }
118 
119 MooseMesh &
121 {
122  return *(new AnnularMesh(*this));
123 }
124 
125 void
127 {
128  const Real dt = (_tmax - _tmin) / _nt;
129 
130  MeshBase & mesh = getMesh();
131  mesh.clear();
132  mesh.set_mesh_dimension(2);
133  mesh.set_spatial_dimension(2);
134  BoundaryInfo & boundary_info = mesh.get_boundary_info();
135 
136  const unsigned num_angular_nodes = (_full_annulus ? _nt : _nt + 1);
137  const unsigned num_nodes =
138  (_rmin > 0.0 ? (_nr + 1) * num_angular_nodes : _nr * num_angular_nodes + 1);
139  const unsigned min_nonzero_layer_num = (_rmin > 0.0 ? 0 : 1);
140  std::vector<Node *> nodes(num_nodes);
141  unsigned node_id = 0;
142 
143  // add nodes at rmax that aren't yet connected to any elements
144  Real current_r = _rmax;
145  for (unsigned angle_num = 0; angle_num < num_angular_nodes; ++angle_num)
146  {
147  const Real angle = _tmin + angle_num * dt;
148  const Real x = current_r * std::cos(angle);
149  const Real y = current_r * std::sin(angle);
150  nodes[node_id] = mesh.add_point(Point(x, y, 0.0), node_id);
151  node_id++;
152  }
153 
154  // add nodes at smaller radii, and connect them with elements
155  for (unsigned layer_num = _nr; layer_num > min_nonzero_layer_num; --layer_num)
156  {
157  if (layer_num == 1)
158  current_r = _rmin; // account for precision loss
159  else
160  current_r -= _len * std::pow(_growth_r, layer_num - 1);
161 
162  // add node at angle = _tmin
163  nodes[node_id] = mesh.add_point(
164  Point(current_r * std::cos(_tmin), current_r * std::sin(_tmin), 0.0), node_id);
165  node_id++;
166  for (unsigned angle_num = 1; angle_num < num_angular_nodes; ++angle_num)
167  {
168  const Real angle = _tmin + angle_num * dt;
169  const Real x = current_r * std::cos(angle);
170  const Real y = current_r * std::sin(angle);
171  nodes[node_id] = mesh.add_point(Point(x, y, 0.0), node_id);
172  Elem * elem = mesh.add_elem(new Quad4);
173  elem->set_node(0) = nodes[node_id];
174  elem->set_node(1) = nodes[node_id - 1];
175  elem->set_node(2) = nodes[node_id - num_angular_nodes - 1];
176  elem->set_node(3) = nodes[node_id - num_angular_nodes];
177  elem->subdomain_id() = _quad_subdomain_id;
178  node_id++;
179 
180  if (layer_num == _nr)
181  // add outer boundary (boundary_id = 1)
182  boundary_info.add_side(elem, 2, 1);
183  if (layer_num == 1)
184  // add inner boundary (boundary_id = 0)
185  boundary_info.add_side(elem, 0, 0);
186  if (!_full_annulus && angle_num == 1)
187  // add tmin boundary (boundary_id = 2)
188  boundary_info.add_side(elem, 1, 2);
189  if (!_full_annulus && angle_num == num_angular_nodes - 1)
190  // add tmin boundary (boundary_id = 3)
191  boundary_info.add_side(elem, 3, 3);
192  }
193  if (_full_annulus)
194  {
195  // add element connecting to node at angle=0
196  Elem * elem = mesh.add_elem(new Quad4);
197  elem->set_node(0) = nodes[node_id - num_angular_nodes];
198  elem->set_node(1) = nodes[node_id - 1];
199  elem->set_node(2) = nodes[node_id - num_angular_nodes - 1];
200  elem->set_node(3) = nodes[node_id - 2 * num_angular_nodes];
201  elem->subdomain_id() = _quad_subdomain_id;
202 
203  if (layer_num == _nr)
204  // add outer boundary (boundary_id = 1)
205  boundary_info.add_side(elem, 2, 1);
206  if (layer_num == 1)
207  // add inner boundary (boundary_id = 0)
208  boundary_info.add_side(elem, 0, 0);
209  }
210  }
211 
212  // add single node at origin, if relevant
213  if (_rmin == 0.0)
214  {
215  nodes[node_id] = mesh.add_point(Point(0.0, 0.0, 0.0), node_id);
216  boundary_info.add_node(node_id, 0); // boundary_id=0 is centre
217  for (unsigned angle_num = 0; angle_num < num_angular_nodes - 1; ++angle_num)
218  {
219  Elem * elem = mesh.add_elem(new Tri3);
220  elem->set_node(0) = nodes[node_id];
221  elem->set_node(1) = nodes[node_id - num_angular_nodes + angle_num];
222  elem->set_node(2) = nodes[node_id - num_angular_nodes + angle_num + 1];
223  elem->subdomain_id() = _tri_subdomain_id;
224  }
225  if (_full_annulus)
226  {
227  Elem * elem = mesh.add_elem(new Tri3);
228  elem->set_node(0) = nodes[node_id];
229  elem->set_node(1) = nodes[node_id - 1];
230  elem->set_node(2) = nodes[node_id - num_angular_nodes];
231  elem->subdomain_id() = _tri_subdomain_id;
232  }
233  }
234 
235  boundary_info.sideset_name(0) = "rmin";
236  boundary_info.sideset_name(1) = "rmax";
237  boundary_info.nodeset_name(0) = "rmin";
238  boundary_info.nodeset_name(1) = "rmax";
239  if (!_full_annulus)
240  {
241  boundary_info.sideset_name(2) = "tmin";
242  boundary_info.sideset_name(3) = "tmax";
243  boundary_info.nodeset_name(2) = "tmin";
244  boundary_info.nodeset_name(3) = "tmax";
245  }
246 
247  mesh.prepare_for_use(false);
248 }
virtual void buildMesh() override
Must be overridden by child classes.
Definition: AnnularMesh.C:126
const Real _growth_r
Bias on radial meshing.
Definition: AnnularMesh.h:62
void addRequiredRangeCheckedParam(const std::string &name, const std::string &parsed_function, const std::string &doc_string)
These methods add an range checked parameters.
const SubdomainID _tri_subdomain_id
Subdomain ID of created tri elements (that only exist if rmin=0)
Definition: AnnularMesh.h:74
const Real _rmin
Minimum radius.
Definition: AnnularMesh.h:50
subdomain_id_type SubdomainID
Definition: MooseTypes.h:77
virtual Real getMinInDimension(unsigned int component) const override
Returns the min or max of the requested dimension respectively.
Definition: AnnularMesh.C:88
const Real _rmax
Maximum radius.
Definition: AnnularMesh.h:53
const unsigned _nr
Number of elements in radial direction.
Definition: AnnularMesh.h:44
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
InputParameters validParams< AnnularMesh >()
Definition: AnnularMesh.C:22
virtual Real getMaxInDimension(unsigned int component) const override
Definition: AnnularMesh.C:104
AnnularMesh(const InputParameters &parameters)
Definition: AnnularMesh.C:58
static PetscErrorCode Vec x
virtual MooseMesh & clone() const override
Clone method.
Definition: AnnularMesh.C:120
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 bool _full_annulus
Whether a full annulus (as opposed to a sector) will needs to generate.
Definition: AnnularMesh.h:68
InputParameters validParams< MooseMesh >()
Definition: MooseMesh.C:65
MeshBase & getMesh()
Accessor for the underlying libMesh Mesh object.
Definition: MooseMesh.C:2355
const SubdomainID _quad_subdomain_id
Subdomain ID of created quad elements.
Definition: AnnularMesh.h:71
MooseMesh wraps a libMesh::Mesh object and enhances its capabilities by caching additional data and s...
Definition: MooseMesh.h:74
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) ...
Definition: AnnularMesh.h:65
const Real _tmin
Minimum angle.
Definition: AnnularMesh.h:56
const Real _tmax
Maximum angle.
Definition: AnnularMesh.h:59
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 mooseError(Args &&...args) const
Definition: MooseObject.h:80
void addRangeCheckedParam(const std::string &name, const T &value, const std::string &parsed_function, const std::string &doc_string)
virtual Elem * elem(const dof_id_type i)
Various accessors (pointers/references) for Elem "i".
Definition: MooseMesh.C:2048
const unsigned _nt
Number of elements in angular direction.
Definition: AnnularMesh.h:47