www.mooseframework.org
RinglebMeshGenerator.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 "RinglebMeshGenerator.h"
11 #include "CastUniquePointer.h"
12 
13 #include "libmesh/replicated_mesh.h"
14 #include "libmesh/mesh_modification.h"
15 #include "libmesh/face_quad4.h"
16 #include "libmesh/face_tri3.h"
17 
19 
20 template <>
23 {
25 
26  params.addRequiredParam<Real>("gamma", "Gamma parameter");
27  params.addRequiredParam<Real>("kmax", "Value of k on the inner wall.");
28  params.addRequiredParam<Real>("kmin", "Value of k on the outer wall.");
29  params.addRequiredParam<int>("num_q_pts",
30  "How many points to discretize the range q = (0.5, k) into.");
31  params.addRequiredParam<int>("n_extra_q_pts",
32  "How many 'extra' points should be inserted in the final element"
33  " *in addition to* the equispaced q points.");
34  params.addRequiredParam<int>("num_k_pts", "How many points in the range k=(kmin, kmax).");
35  params.addParam<boundary_id_type>("inflow_bid", 1, "The boundary id to use for the inflow");
36  params.addParam<boundary_id_type>(
37  "inner_wall_bid", 2, "The boundary id to use for the inner wall");
38  params.addParam<boundary_id_type>("outflow_bid", 3, "The boundary id to use for the outflow");
39  params.addParam<boundary_id_type>(
40  "outer_wall_bid", 4, "The boundary id to use for the outer wall");
41  params.addParam<bool>(
42  "triangles", false, "If true, all the quadrilateral elements will be split into triangles");
43  params.addClassDescription("Creates a mesh for the Ringleb problem.");
44 
45  return params;
46 }
47 
49  : MeshGenerator(parameters),
50  _gamma(getParam<Real>("gamma")),
51  _kmax(getParam<Real>("kmax")),
52  _kmin(getParam<Real>("kmin")),
53  _num_q_pts(getParam<int>("num_q_pts")),
54  _n_extra_q_pts(getParam<int>("n_extra_q_pts")),
55  _num_k_pts(getParam<int>("num_k_pts")),
56  _inflow_bid(getParam<boundary_id_type>("inflow_bid")),
57  _outflow_bid(getParam<boundary_id_type>("outflow_bid")),
58  _inner_wall_bid(getParam<boundary_id_type>("inner_wall_bid")),
59  _outer_wall_bid(getParam<boundary_id_type>("outer_wall_bid")),
60  _triangles(getParam<bool>("triangles"))
61 {
62  // catch likely user errors
63  if (_kmax <= _kmin)
64  mooseError("RinglebMesh: kmax must be greater than kmin");
65 }
66 
67 std::vector<Real>
68 RinglebMeshGenerator::arhopj(const Real & gamma, const std::vector<Real> & q, const int & index)
69 {
70  std::vector<Real> values(4);
71  Real a = std::sqrt(1 - ((gamma - 1) / 2.) * std::pow(q[index], 2));
72  Real rho = std::pow(a, 2. / (gamma - 1));
73  Real p = (1. / gamma) * std::pow(a, 2 * gamma / (gamma - 1));
74  Real J = 1. / a + 1. / (3. * std::pow(a, 3)) + 1. / (5. * std::pow(a, 5)) -
75  0.5 * std::log((1 + a) / (1 - a));
76  values = {a, rho, p, J};
77  return values;
78 }
79 
80 std::vector<Real>
81 RinglebMeshGenerator::computexy(const std::vector<Real> values,
82  const int & i,
83  const int & index,
84  const std::vector<Real> & ks,
85  const std::vector<Real> & q)
86 {
87  std::vector<Real> xy(2);
88 
89  // Compute x(q,k)
90  xy[0] = 0.5 / values[1] * (2. / ks[i] / ks[i] - 1. / q[index] / q[index]) - 0.5 * values[3];
91 
92  // Compute the term that goes under the sqrt sign
93  // If 1 - (q/k)^2 is slightly negative, we make it zero.
94  Real sqrt_term = 1. - q[index] * q[index] / ks[i] / ks[i];
95  sqrt_term = std::max(sqrt_term, 0.);
96 
97  // Compute y(q,k)
98  xy[1] = 1. / (ks[i] * values[1] * q[index]) * std::sqrt(sqrt_term);
99 
100  return xy;
101 }
102 
103 std::unique_ptr<MeshBase>
105 {
106  std::unique_ptr<ReplicatedMesh> mesh = libmesh_make_unique<ReplicatedMesh>(comm(), 2);
107 
109  std::vector<std::vector<Node *>> stream_nodes(_num_k_pts);
110 
112  int current_node_id = 0;
113 
115  std::vector<Real> ks(_num_k_pts);
116  Real diff = (_kmax - _kmin) / (_num_k_pts - 1);
117  for (int i = 0; i < _num_k_pts; i++)
118  ks[i] = _kmin + i * diff;
119 
120  for (int i = 0; i < _num_k_pts; i++)
121  {
122  stream_nodes[i].resize(2 * (_num_q_pts + _n_extra_q_pts));
123 
125  std::vector<Real> q(_num_q_pts);
126  Real diffq = (ks[i] - 0.5) / (_num_q_pts - 1);
127  for (int j = 0; j < _num_q_pts; j++)
128  q[j] = 0.5 + j * diffq;
129 
131  for (int j = _num_q_pts; j < _num_q_pts + _n_extra_q_pts; j++)
132  {
133  std::vector<Real>::iterator it = q.end();
134  q.insert(--it, 0.3 * q[j - 2] + 0.7 * q[j - 1]);
135  }
136 
137  std::vector<Real> vals(4);
138  std::vector<Real> xy(2);
140  for (int j = 0; j < _num_q_pts + _n_extra_q_pts; j++)
141  {
142  // Compute the different parameters
143  vals = arhopj(_gamma, q, j);
144 
145  // Compute x and y
146  xy = computexy(vals, i, j, ks, q);
147 
148  // Create a node with (x,y) coordinates as it's on the upper part of the mesh
149  if (j != _num_q_pts + _n_extra_q_pts - 1)
150  stream_nodes[i][j] = mesh->add_point(Point(xy[0], xy[1]), current_node_id++);
151  }
152 
154  for (int j = _num_q_pts + _n_extra_q_pts; j < 2 * (_num_q_pts + _n_extra_q_pts); j++)
155  {
156  int index = 2 * (_num_q_pts + _n_extra_q_pts) - 1 - j;
157  // Compute the different parameters
158  vals = arhopj(_gamma, q, index);
159 
160  // Compute x and y
161  xy = computexy(vals, i, index, ks, q);
162 
163  // Create a node with (x,-y) coordinates as it's on the lower part of the mesh
164  stream_nodes[i][j] = mesh->add_point(Point(xy[0], -xy[1]), current_node_id++);
165  }
166  }
167 
169  for (int i = 0; i < _num_k_pts - 1; i++)
170  {
171  for (int j = 0; j < 2 * (_num_q_pts + _n_extra_q_pts) - 1; j++)
172  {
174  if (j != _num_q_pts + _n_extra_q_pts - 1 and j != _num_q_pts + _n_extra_q_pts - 2)
175  {
176  Elem * elem = mesh->add_elem(new Quad4);
177  elem->set_node(0) = stream_nodes[i][j];
178  elem->set_node(1) = stream_nodes[i][j + 1];
179  elem->set_node(2) = stream_nodes[i + 1][j + 1];
180  elem->set_node(3) = stream_nodes[i + 1][j];
181 
182  if (i == 0)
183  mesh->boundary_info->add_side(elem->id(), /*side=*/0, _outer_wall_bid);
184  if (j == 0)
185  mesh->boundary_info->add_side(elem->id(), /*side=*/3, _inflow_bid);
186  if (j == 2 * (_num_q_pts + _n_extra_q_pts) - 2)
187  mesh->boundary_info->add_side(elem->id(), /*side=*/1, _outflow_bid);
188  if (i == _num_k_pts - 2)
189  mesh->boundary_info->add_side(elem->id(), /*side=*/2, _inner_wall_bid);
190  }
191  else if (j == _num_q_pts + _n_extra_q_pts - 2)
192  {
193  Elem * elem = mesh->add_elem(new Quad4);
194  elem->set_node(0) = stream_nodes[i][j];
195  elem->set_node(1) = stream_nodes[i][j + 2];
196  elem->set_node(2) = stream_nodes[i + 1][j + 2];
197  elem->set_node(3) = stream_nodes[i + 1][j];
198 
199  if (i == 0)
200  mesh->boundary_info->add_side(elem->id(), /*side=*/0, _outer_wall_bid);
201  if (i == _num_k_pts - 2)
202  mesh->boundary_info->add_side(elem->id(), /*side=*/2, _inner_wall_bid);
203  }
204  }
205  }
206 
208  mesh->prepare_for_use();
209 
211  if (_triangles)
212  MeshTools::Modification::all_tri(*mesh);
213 
215  mesh->boundary_info->sideset_name(_inflow_bid) = "inflow";
216  mesh->boundary_info->sideset_name(_outflow_bid) = "outflow";
217  mesh->boundary_info->sideset_name(_inner_wall_bid) = "inner_wall";
218  mesh->boundary_info->sideset_name(_outer_wall_bid) = "outer_wall";
219 
220  return dynamic_pointer_cast<MeshBase>(mesh);
221 }
const Real & _kmax
k is a streamline parameter, i.e.
const Real & _kmin
kmin corresponds to the outer wall
std::vector< Real > computexy(const std::vector< Real > values, const int &i, const int &index, const std::vector< Real > &ks, const std::vector< Real > &q)
InputParameters validParams< RinglebMeshGenerator >()
const int & _n_extra_q_pts
how many "extra" points should be inserted in the nearest element from the horizontal in additi /// o...
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.
void mooseError(Args &&... args) const
Definition: MooseObject.h:147
const boundary_id_type _outflow_bid
const boundary_id_type _inflow_bid
The boundary ids to use for the ringleb mesh.
Generates a mesh given all the parameters.
const boundary_id_type _inner_wall_bid
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...
std::unique_ptr< MeshBase > generate() override
Generate / modify the mesh.
std::vector< Real > arhopj(const Real &gamma, const std::vector< Real > &q, const int &index)
const Real & _gamma
Gamma.
Real pow(Real x, int e)
Definition: MathUtils.C:211
const boundary_id_type _outer_wall_bid
const bool & _triangles
This parameter, if true, allows to split the quadrilateral elements into triangular elements...
registerMooseObject("MooseApp", RinglebMeshGenerator)
const int & _num_q_pts
How many points to discretize the range q = (0.5, k) into.
MPI_Comm comm
const int & _num_k_pts
how many points in the range k=(kmin, kmax).
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...
RinglebMeshGenerator(const InputParameters &parameters)
InputParameters validParams< MeshGenerator >()
Definition: MeshGenerator.C:16
MeshGenerators are objects that can modify or add to an existing mesh.
Definition: MeshGenerator.h:30