Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
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 :
18 : registerMooseObject("MooseApp", RinglebMeshGenerator);
19 :
20 : InputParameters
21 14281 : RinglebMeshGenerator::validParams()
22 : {
23 14281 : InputParameters params = MeshGenerator::validParams();
24 :
25 14281 : params.addRequiredParam<Real>("gamma", "Gamma parameter");
26 14281 : params.addRequiredParam<Real>("kmax", "Value of k on the inner wall.");
27 14281 : params.addRequiredParam<Real>("kmin", "Value of k on the outer wall.");
28 14281 : params.addRequiredParam<int>("num_q_pts",
29 : "How many points to discretize the range q = (0.5, k) into.");
30 14281 : params.addRequiredParam<int>("n_extra_q_pts",
31 : "How many 'extra' points should be inserted in the final element"
32 : " *in addition to* the equispaced q points.");
33 14281 : params.addRequiredParam<int>("num_k_pts", "How many points in the range k=(kmin, kmax).");
34 14281 : params.addParam<boundary_id_type>("inflow_bid", 1, "The boundary id to use for the inflow");
35 42843 : params.addParam<boundary_id_type>(
36 28562 : "inner_wall_bid", 2, "The boundary id to use for the inner wall");
37 14281 : params.addParam<boundary_id_type>("outflow_bid", 3, "The boundary id to use for the outflow");
38 42843 : params.addParam<boundary_id_type>(
39 28562 : "outer_wall_bid", 4, "The boundary id to use for the outer wall");
40 42843 : params.addParam<bool>(
41 28562 : "triangles", false, "If true, all the quadrilateral elements will be split into triangles");
42 14281 : params.addClassDescription("Creates a mesh for the Ringleb problem.");
43 :
44 14281 : return params;
45 0 : }
46 :
47 8 : RinglebMeshGenerator::RinglebMeshGenerator(const InputParameters & parameters)
48 : : MeshGenerator(parameters),
49 8 : _gamma(getParam<Real>("gamma")),
50 8 : _kmax(getParam<Real>("kmax")),
51 8 : _kmin(getParam<Real>("kmin")),
52 8 : _num_q_pts(getParam<int>("num_q_pts")),
53 8 : _n_extra_q_pts(getParam<int>("n_extra_q_pts")),
54 8 : _num_k_pts(getParam<int>("num_k_pts")),
55 8 : _inflow_bid(getParam<boundary_id_type>("inflow_bid")),
56 8 : _outflow_bid(getParam<boundary_id_type>("outflow_bid")),
57 8 : _inner_wall_bid(getParam<boundary_id_type>("inner_wall_bid")),
58 8 : _outer_wall_bid(getParam<boundary_id_type>("outer_wall_bid")),
59 16 : _triangles(getParam<bool>("triangles"))
60 : {
61 8 : declareMeshProperty("use_distributed_mesh", false);
62 :
63 : // catch likely user errors
64 8 : if (_kmax <= _kmin)
65 0 : mooseError("RinglebMesh: kmax must be greater than kmin");
66 8 : }
67 :
68 : std::vector<Real>
69 3168 : RinglebMeshGenerator::arhopj(const Real & gamma, const std::vector<Real> & q, const int & index)
70 : {
71 3168 : std::vector<Real> values(4);
72 3168 : Real a = std::sqrt(1 - ((gamma - 1) / 2.) * std::pow(q[index], 2));
73 3168 : Real rho = std::pow(a, 2. / (gamma - 1));
74 3168 : Real p = (1. / gamma) * std::pow(a, 2 * gamma / (gamma - 1));
75 3168 : Real J = 1. / a + 1. / (3. * std::pow(a, 3)) + 1. / (5. * std::pow(a, 5)) -
76 3168 : 0.5 * std::log((1 + a) / (1 - a));
77 3168 : values = {a, rho, p, J};
78 3168 : return values;
79 0 : }
80 :
81 : std::vector<Real>
82 3168 : RinglebMeshGenerator::computexy(const std::vector<Real> values,
83 : const int & i,
84 : const int & index,
85 : const std::vector<Real> & ks,
86 : const std::vector<Real> & q)
87 : {
88 3168 : std::vector<Real> xy(2);
89 :
90 : // Compute x(q,k)
91 3168 : xy[0] = 0.5 / values[1] * (2. / ks[i] / ks[i] - 1. / q[index] / q[index]) - 0.5 * values[3];
92 :
93 : // Compute the term that goes under the sqrt sign
94 : // If 1 - (q/k)^2 is slightly negative, we make it zero.
95 3168 : Real sqrt_term = 1. - q[index] * q[index] / ks[i] / ks[i];
96 3168 : sqrt_term = std::max(sqrt_term, 0.);
97 :
98 : // Compute y(q,k)
99 3168 : xy[1] = 1. / (ks[i] * values[1] * q[index]) * std::sqrt(sqrt_term);
100 :
101 6336 : return xy;
102 : }
103 :
104 : std::unique_ptr<MeshBase>
105 8 : RinglebMeshGenerator::generate()
106 : {
107 8 : std::unique_ptr<ReplicatedMesh> mesh = buildReplicatedMesh(2);
108 8 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
109 :
110 : /// Data structure that holds pointers to the Nodes of each streamline.
111 8 : std::vector<std::vector<Node *>> stream_nodes(_num_k_pts);
112 :
113 : /// Node id counter.
114 8 : int current_node_id = 0;
115 :
116 : /// Vector containing the regularly spaced k values
117 8 : std::vector<Real> ks(_num_k_pts);
118 8 : Real diff = (_kmax - _kmin) / (_num_k_pts - 1);
119 80 : for (int i = 0; i < _num_k_pts; i++)
120 72 : ks[i] = _kmin + i * diff;
121 :
122 80 : for (int i = 0; i < _num_k_pts; i++)
123 : {
124 72 : stream_nodes[i].resize(2 * (_num_q_pts + _n_extra_q_pts));
125 :
126 : /// Vector containing the regularly spaced (and the extra q points) q values
127 72 : std::vector<Real> q(_num_q_pts);
128 72 : Real diffq = (ks[i] - 0.5) / (_num_q_pts - 1);
129 1512 : for (int j = 0; j < _num_q_pts; j++)
130 1440 : q[j] = 0.5 + j * diffq;
131 :
132 : /// Add the extra q points
133 216 : for (int j = _num_q_pts; j < _num_q_pts + _n_extra_q_pts; j++)
134 : {
135 144 : std::vector<Real>::iterator it = q.end();
136 144 : q.insert(--it, 0.3 * q[j - 2] + 0.7 * q[j - 1]);
137 : }
138 :
139 72 : std::vector<Real> vals(4);
140 72 : std::vector<Real> xy(2);
141 : /// Create the nodes for the upper part
142 1656 : for (int j = 0; j < _num_q_pts + _n_extra_q_pts; j++)
143 : {
144 : // Compute the different parameters
145 1584 : vals = arhopj(_gamma, q, j);
146 :
147 : // Compute x and y
148 1584 : xy = computexy(vals, i, j, ks, q);
149 :
150 : // Create a node with (x,y) coordinates as it's on the upper part of the mesh
151 1584 : if (j != _num_q_pts + _n_extra_q_pts - 1)
152 1512 : stream_nodes[i][j] = mesh->add_point(Point(xy[0], xy[1]), current_node_id++);
153 : }
154 :
155 : /// Create the nodes for the lower part
156 1656 : for (int j = _num_q_pts + _n_extra_q_pts; j < 2 * (_num_q_pts + _n_extra_q_pts); j++)
157 : {
158 1584 : int index = 2 * (_num_q_pts + _n_extra_q_pts) - 1 - j;
159 : // Compute the different parameters
160 1584 : vals = arhopj(_gamma, q, index);
161 :
162 : // Compute x and y
163 1584 : xy = computexy(vals, i, index, ks, q);
164 :
165 : // Create a node with (x,-y) coordinates as it's on the lower part of the mesh
166 1584 : stream_nodes[i][j] = mesh->add_point(Point(xy[0], -xy[1]), current_node_id++);
167 : }
168 72 : }
169 :
170 : /// Add elements for the whole mesh
171 72 : for (int i = 0; i < _num_k_pts - 1; i++)
172 : {
173 2816 : for (int j = 0; j < 2 * (_num_q_pts + _n_extra_q_pts) - 1; j++)
174 : {
175 : /// This is done in order to avoid having two nodes at the exact same location (on the straight /// line y=0)
176 2752 : if (j != _num_q_pts + _n_extra_q_pts - 1 and j != _num_q_pts + _n_extra_q_pts - 2)
177 : {
178 2624 : Elem * elem = mesh->add_elem(new Quad4);
179 2624 : elem->set_node(0, stream_nodes[i][j]);
180 2624 : elem->set_node(1, stream_nodes[i][j + 1]);
181 2624 : elem->set_node(2, stream_nodes[i + 1][j + 1]);
182 2624 : elem->set_node(3, stream_nodes[i + 1][j]);
183 :
184 2624 : if (i == 0)
185 328 : boundary_info.add_side(elem->id(), /*side=*/0, _outer_wall_bid);
186 2624 : if (j == 0)
187 64 : boundary_info.add_side(elem->id(), /*side=*/3, _inflow_bid);
188 2624 : if (j == 2 * (_num_q_pts + _n_extra_q_pts) - 2)
189 64 : boundary_info.add_side(elem->id(), /*side=*/1, _outflow_bid);
190 2624 : if (i == _num_k_pts - 2)
191 328 : boundary_info.add_side(elem->id(), /*side=*/2, _inner_wall_bid);
192 2624 : }
193 128 : else if (j == _num_q_pts + _n_extra_q_pts - 2)
194 : {
195 64 : Elem * elem = mesh->add_elem(new Quad4);
196 64 : elem->set_node(0, stream_nodes[i][j]);
197 64 : elem->set_node(1, stream_nodes[i][j + 2]);
198 64 : elem->set_node(2, stream_nodes[i + 1][j + 2]);
199 64 : elem->set_node(3, stream_nodes[i + 1][j]);
200 :
201 64 : if (i == 0)
202 8 : boundary_info.add_side(elem->id(), /*side=*/0, _outer_wall_bid);
203 64 : if (i == _num_k_pts - 2)
204 8 : boundary_info.add_side(elem->id(), /*side=*/2, _inner_wall_bid);
205 : }
206 : }
207 : }
208 :
209 : /// Find neighbors, etc.
210 8 : mesh->prepare_for_use();
211 :
212 : /// Create the triangular elements if required by the user
213 8 : if (_triangles)
214 8 : MeshTools::Modification::all_tri(*mesh);
215 :
216 : /// Create sideset names.
217 8 : boundary_info.sideset_name(_inflow_bid) = "inflow";
218 8 : boundary_info.sideset_name(_outflow_bid) = "outflow";
219 8 : boundary_info.sideset_name(_inner_wall_bid) = "inner_wall";
220 8 : boundary_info.sideset_name(_outer_wall_bid) = "outer_wall";
221 :
222 16 : return dynamic_pointer_cast<MeshBase>(mesh);
223 8 : }
|
