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 : // Creates multiple superellipsoids that are positioned randomly throughout the domain
11 : // each semiaxis can be varied by a uniform or normal distribution
12 :
13 : #include "MultiSmoothSuperellipsoidIC.h"
14 :
15 : // MOOSE includes
16 : #include "MooseMesh.h"
17 : #include "MooseUtils.h"
18 : #include "MooseVariable.h"
19 :
20 : registerMooseObject("PhaseFieldApp", MultiSmoothSuperellipsoidIC);
21 :
22 : InputParameters
23 46 : MultiSmoothSuperellipsoidIC::validParams()
24 : {
25 46 : InputParameters params = SmoothSuperellipsoidBaseIC::validParams();
26 46 : params.addClassDescription("Random distribution of smooth ellipse with given minimum spacing");
27 92 : params.addRequiredParam<std::vector<unsigned int>>("numbub",
28 : "Vector of the number of bubbles to place");
29 92 : params.addRequiredParam<std::vector<Real>>(
30 : "bubspac",
31 : "Vector of the minimum spacing of bubbles of one type, measured from center to center");
32 92 : params.addParam<unsigned int>("max_num_tries", 1000, "The number of tries");
33 92 : params.addRequiredParam<std::vector<Real>>(
34 : "semiaxis_a", "Vector of mean semiaxis values in the x direction for the ellipse");
35 92 : params.addRequiredParam<std::vector<Real>>(
36 : "semiaxis_b", "Vector of mean semiaxis values in the y direction for the ellipse");
37 92 : params.addRequiredParam<std::vector<Real>>("semiaxis_c",
38 : "Vector of mean semiaxis values in the z direction "
39 : "for the ellipse, must be greater than 0 even if 2D.");
40 46 : params.addParam<std::vector<Real>>(
41 : "exponent",
42 46 : std::vector<Real>(),
43 : "Vector of exponents for each superellipsoid, n=2 is a normal ellipse");
44 46 : params.addParam<std::vector<Real>>("semiaxis_a_variation",
45 46 : std::vector<Real>(),
46 : "Vector of plus or minus fractions of random variation in the "
47 : "bubble semiaxis in the x direction for uniform, standard "
48 : "deviation for normal");
49 46 : params.addParam<std::vector<Real>>("semiaxis_b_variation",
50 46 : std::vector<Real>(),
51 : "Vector of plus or minus fractions of random variation in the "
52 : "bubble semiaxis in the y direction for uniform, standard "
53 : "deviation for normal");
54 46 : params.addParam<std::vector<Real>>("semiaxis_c_variation",
55 46 : std::vector<Real>(),
56 : "Vector of plus or minus fractions of random variation in the "
57 : "bubble semiaxis in the z direction for uniform, standard "
58 : "deviation for normal. Must be set to 0 if 2D.");
59 92 : params.addParam<bool>("check_extremes",
60 92 : false,
61 : "Check all Superellipsoid extremes (center +- "
62 : "each semiaxis) for overlap, must have "
63 : "prevent_overlap set to True.");
64 92 : params.addParam<bool>("prevent_overlap",
65 92 : false,
66 : "Check all Superellipsoid centers for overlap with other superellipsoids.");
67 92 : params.addParam<bool>("vary_axes_independently",
68 92 : true,
69 : "If true the length of each semiaxis is randomly chosen "
70 : "within the provided parameters, if false then one random "
71 : "number is generated and applied to all semiaxes.");
72 92 : MooseEnum rand_options("uniform normal none", "none");
73 92 : params.addParam<MooseEnum>(
74 : "semiaxis_variation_type",
75 : rand_options,
76 : "Type of distribution that random superellipsoid semiaxes will follow");
77 46 : return params;
78 46 : }
79 :
80 24 : MultiSmoothSuperellipsoidIC::MultiSmoothSuperellipsoidIC(const InputParameters & parameters)
81 : : SmoothSuperellipsoidBaseIC(parameters),
82 24 : _max_num_tries(getParam<unsigned int>("max_num_tries")),
83 48 : _semiaxis_variation_type(getParam<MooseEnum>("semiaxis_variation_type")),
84 48 : _prevent_overlap(getParam<bool>("prevent_overlap")),
85 48 : _check_extremes(getParam<bool>("check_extremes")),
86 48 : _vary_axes_independently(getParam<bool>("vary_axes_independently")),
87 24 : _numbub(parameters.get<std::vector<unsigned int>>("numbub")),
88 24 : _bubspac(parameters.get<std::vector<Real>>("bubspac")),
89 24 : _exponent(parameters.get<std::vector<Real>>("exponent")),
90 24 : _semiaxis_a(parameters.get<std::vector<Real>>("semiaxis_a")),
91 24 : _semiaxis_b(parameters.get<std::vector<Real>>("semiaxis_b")),
92 24 : _semiaxis_c(parameters.get<std::vector<Real>>("semiaxis_c")),
93 24 : _semiaxis_a_variation(parameters.get<std::vector<Real>>("semiaxis_a_variation")),
94 24 : _semiaxis_b_variation(parameters.get<std::vector<Real>>("semiaxis_b_variation")),
95 48 : _semiaxis_c_variation(parameters.get<std::vector<Real>>("semiaxis_c_variation"))
96 : {
97 24 : }
98 :
99 : void
100 18 : MultiSmoothSuperellipsoidIC::initialSetup()
101 : {
102 18 : unsigned int nv = _numbub.size();
103 :
104 18 : if (nv != _bubspac.size() || nv != _exponent.size() || nv != _semiaxis_a.size() ||
105 36 : nv != _semiaxis_b.size() || nv != _semiaxis_c.size())
106 0 : mooseError("Vectors for numbub, bubspac, exponent, semiaxis_a, semiaxis_b, and semiaxis_c must "
107 : "be the same size.");
108 :
109 18 : if (_semiaxis_variation_type != 2 &&
110 18 : (nv != _semiaxis_a_variation.size() || nv != _semiaxis_b_variation.size() ||
111 : nv != _semiaxis_c_variation.size()))
112 0 : mooseError("Vectors for numbub, semiaxis_a_variation, semiaxis_b_variation, and "
113 : "semiaxis_c_variation must be the same size.");
114 :
115 18 : if (_semiaxis_variation_type == 2 &&
116 0 : (_semiaxis_a_variation.size() > 0 || _semiaxis_b_variation.size() > 0 ||
117 : _semiaxis_c_variation.size() > 0))
118 0 : mooseWarning(
119 : "Values were provided for semiaxis_a/b/c_variation but semiaxis_variation_type is set "
120 : "to 'none' in 'MultiSmoothSuperellipsoidIC'.");
121 :
122 54 : for (_gk = 0; _gk < nv; ++_gk)
123 : {
124 : // Set up domain bounds with mesh tools
125 144 : for (const auto i : make_range(Moose::dim))
126 : {
127 108 : _bottom_left(i) = _mesh.getMinInDimension(i);
128 108 : _top_right(i) = _mesh.getMaxInDimension(i);
129 : }
130 36 : _range = _top_right - _bottom_left;
131 :
132 36 : SmoothSuperellipsoidBaseIC::initialSetup();
133 : }
134 18 : }
135 :
136 : void
137 36 : MultiSmoothSuperellipsoidIC::computeSuperellipsoidSemiaxes()
138 : {
139 : Real randnum;
140 36 : unsigned int start = _as.size();
141 36 : _as.resize(start + _numbub[_gk]);
142 36 : _bs.resize(start + _numbub[_gk]);
143 36 : _cs.resize(start + _numbub[_gk]);
144 :
145 216 : for (unsigned int i = start; i < _as.size(); i++)
146 : {
147 180 : switch (_semiaxis_variation_type)
148 : {
149 180 : case 0: // Uniform distrubtion
150 180 : randnum = _random.rand(_tid);
151 180 : _as[i] = _semiaxis_a[_gk] * (1.0 + (1.0 - 2.0 * randnum) * _semiaxis_a_variation[_gk]);
152 360 : _bs[i] = _semiaxis_b[_gk] *
153 180 : (1.0 + (1.0 - 2.0 * (_vary_axes_independently ? _random.rand(_tid) : randnum)) *
154 180 : _semiaxis_b_variation[_gk]);
155 360 : _cs[i] = _semiaxis_c[_gk] *
156 180 : (1.0 + (1.0 - 2.0 * (_vary_axes_independently ? _random.rand(_tid) : randnum)) *
157 180 : _semiaxis_c_variation[_gk]);
158 180 : break;
159 :
160 0 : case 1: // Normal distribution
161 0 : randnum = _random.randNormal(_tid, 0, 1);
162 0 : _as[i] = _semiaxis_a[_gk] + (randnum * _semiaxis_a_variation[_gk]);
163 0 : _bs[i] = _semiaxis_b[_gk] +
164 0 : ((_vary_axes_independently ? _random.randNormal(_tid, 0, 1) : randnum) *
165 0 : _semiaxis_b_variation[_gk]);
166 0 : _cs[i] = _semiaxis_c[_gk] +
167 0 : ((_vary_axes_independently ? _random.randNormal(_tid, 0, 1) : randnum) *
168 0 : _semiaxis_c_variation[_gk]);
169 0 : break;
170 :
171 0 : case 2: // No variation
172 0 : _as[i] = _semiaxis_a[_gk];
173 0 : _bs[i] = _semiaxis_b[_gk];
174 0 : _cs[i] = _semiaxis_c[_gk];
175 : }
176 :
177 360 : _as[i] = _as[i] < 0.0 ? 0.0 : _as[i];
178 360 : _bs[i] = _bs[i] < 0.0 ? 0.0 : _bs[i];
179 360 : _cs[i] = _cs[i] < 0.0 ? 0.0 : _cs[i];
180 : }
181 36 : }
182 :
183 : void
184 36 : MultiSmoothSuperellipsoidIC::computeSuperellipsoidCenters()
185 : {
186 36 : unsigned int start = _centers.size();
187 36 : _centers.resize(start + _numbub[_gk]);
188 :
189 216 : for (unsigned int i = start; i < _centers.size(); i++)
190 : {
191 : // Vary circle center positions
192 : unsigned int num_tries = 0;
193 252 : while (num_tries < _max_num_tries)
194 : {
195 252 : num_tries++;
196 :
197 : RealTensorValue ran;
198 252 : ran(0, 0) = _random.rand(_tid);
199 252 : ran(1, 1) = _random.rand(_tid);
200 252 : ran(2, 2) = _random.rand(_tid);
201 :
202 252 : _centers[i] = _bottom_left + ran * _range;
203 :
204 1305 : for (unsigned int j = 0; j < i; ++j)
205 2241 : if (_mesh.minPeriodicDistance(_var.number(), _centers[i], _centers[j]) < _bubspac[_gk] ||
206 1116 : ellipsoidsOverlap(i, j))
207 72 : goto fail;
208 :
209 : // accept the position of the new center
210 180 : goto accept;
211 :
212 : // retry a new position until tries are exhausted
213 : fail:
214 : continue;
215 : }
216 :
217 0 : if (num_tries == _max_num_tries)
218 0 : mooseError("max_num_tries reached in 'MultiSmoothSuperellipsoidIC'.");
219 :
220 0 : accept:
221 : continue;
222 180 : }
223 36 : }
224 :
225 : void
226 36 : MultiSmoothSuperellipsoidIC::computeSuperellipsoidExponents()
227 : {
228 36 : unsigned int start = _ns.size();
229 36 : _ns.resize(start + _numbub[_gk]);
230 :
231 216 : for (unsigned int i = start; i < _ns.size(); ++i)
232 180 : _ns[i] = _exponent[_gk];
233 36 : }
234 :
235 : bool
236 1116 : MultiSmoothSuperellipsoidIC::ellipsoidsOverlap(unsigned int i, unsigned int j)
237 : {
238 : // Check for overlap between centers
239 1116 : const Point dist_vec = _mesh.minPeriodicVector(_var.number(), _centers[i], _centers[j]);
240 1116 : const Real dist = dist_vec.norm();
241 :
242 : // Handle this case independently because we cannot calculate polar angles at this point
243 1116 : if (MooseUtils::absoluteFuzzyEqual(dist, 0.0))
244 : return true;
245 :
246 : // Compute the distance r from the center of the superellipsoid to its outside edge
247 : // along the vector from the center to the current point
248 : // This uses the equation for a superellipse in polar coordinates and substitutes
249 : // distances for sin, cos functions
250 : Real rmn;
251 :
252 : // calculate rmn = r^(-n), replacing sin, cos functions with distances
253 1116 : rmn = (std::pow(std::abs(dist_vec(0) / dist / _as[i]), _ns[i]) +
254 1116 : std::pow(std::abs(dist_vec(1) / dist / _bs[i]), _ns[i]) +
255 1116 : std::pow(std::abs(dist_vec(2) / dist / _cs[i]), _ns[i]));
256 : // calculate r2 from rmn
257 1116 : const Real r1 = std::pow(rmn, (-1.0 / _ns[i]));
258 :
259 : // calculate rmn = r^(-n), replacing sin, cos functions with distances
260 1116 : rmn = (std::pow(std::abs(dist_vec(0) / dist / _as[j]), _ns[j]) +
261 1116 : std::pow(std::abs(dist_vec(1) / dist / _bs[j]), _ns[j]) +
262 1116 : std::pow(std::abs(dist_vec(2) / dist / _cs[j]), _ns[j]));
263 1116 : const Real r2 = std::pow(rmn, (-1.0 / _ns[j]));
264 :
265 1116 : if (dist < r1 + r2)
266 : return true;
267 :
268 : // Check for overlap between extremes of new ellipsoid candidate and the center
269 : // of accepted ellisoids if _check_extremes enabled
270 1053 : if (_check_extremes)
271 459 : return checkExtremes(i, j) || checkExtremes(j, i);
272 :
273 : // otherwise no overlap has been detected
274 : return false;
275 : }
276 :
277 : bool
278 918 : MultiSmoothSuperellipsoidIC::checkExtremes(unsigned int i, unsigned int j)
279 : {
280 : Point tmp_p;
281 6426 : for (unsigned int pc = 0; pc < 6; pc++)
282 : {
283 5508 : tmp_p = _centers[j];
284 : // Find extremes along semiaxis of candidate ellipsoids
285 : if (pc == 0)
286 918 : tmp_p(0) -= _as[j];
287 : else if (pc == 1)
288 918 : tmp_p(0) += _as[j];
289 : else if (pc == 2)
290 918 : tmp_p(1) -= _bs[j];
291 : else if (pc == 3)
292 918 : tmp_p(1) += _bs[j];
293 : else if (pc == 4)
294 918 : tmp_p(2) -= _cs[j];
295 : else
296 918 : tmp_p(2) += _cs[j];
297 :
298 5508 : const Point dist_vec = _mesh.minPeriodicVector(_var.number(), _centers[i], tmp_p);
299 5508 : const Real dist = dist_vec.norm();
300 :
301 : // Handle this case independently because we cannot calculate polar angles at this point
302 5508 : if (MooseUtils::absoluteFuzzyEqual(dist, 0.0))
303 0 : return true;
304 :
305 : // calculate rmn = r^(-n), replacing sin, cos functions with distances
306 5508 : Real rmn = (std::pow(std::abs(dist_vec(0) / dist / _as[i]), _ns[i]) +
307 5508 : std::pow(std::abs(dist_vec(1) / dist / _bs[i]), _ns[i]) +
308 5508 : std::pow(std::abs(dist_vec(2) / dist / _cs[i]), _ns[i]));
309 5508 : Real r = std::pow(rmn, (-1.0 / _ns[i]));
310 :
311 5508 : if (dist < r)
312 : return true;
313 : }
314 :
315 : return false;
316 : }
|
