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