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 30 : MultiSmoothSuperellipsoidIC::validParams()
24 : {
25 30 : InputParameters params = SmoothSuperellipsoidBaseIC::validParams();
26 30 : params.addClassDescription("Random distribution of smooth ellipse with given minimum spacing");
27 60 : params.addRequiredParam<std::vector<unsigned int>>("numbub",
28 : "Vector of the number of bubbles to place");
29 60 : params.addRequiredParam<std::vector<Real>>(
30 : "bubspac",
31 : "Vector of the minimum spacing of bubbles of one type, measured from center to center");
32 60 : params.addParam<unsigned int>("max_num_tries", 1000, "The number of tries");
33 60 : params.addRequiredParam<std::vector<Real>>(
34 : "semiaxis_a", "Vector of mean semiaxis values in the x direction for the ellipse");
35 60 : params.addRequiredParam<std::vector<Real>>(
36 : "semiaxis_b", "Vector of mean semiaxis values in the y direction for the ellipse");
37 60 : 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 30 : params.addParam<std::vector<Real>>(
41 : "exponent",
42 30 : std::vector<Real>(),
43 : "Vector of exponents for each superellipsoid, n=2 is a normal ellipse");
44 30 : params.addParam<std::vector<Real>>("semiaxis_a_variation",
45 30 : 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 30 : params.addParam<std::vector<Real>>("semiaxis_b_variation",
50 30 : 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 30 : params.addParam<std::vector<Real>>("semiaxis_c_variation",
55 30 : 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 60 : params.addParam<bool>("check_extremes",
60 60 : false,
61 : "Check all Superellipsoid extremes (center +- "
62 : "each semiaxis) for overlap, must have "
63 : "prevent_overlap set to True.");
64 60 : params.addParam<bool>("prevent_overlap",
65 60 : false,
66 : "Check all Superellipsoid centers for overlap with other superellipsoids.");
67 60 : params.addParam<bool>("vary_axes_independently",
68 60 : 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 60 : MooseEnum rand_options("uniform normal none", "none");
73 60 : params.addParam<MooseEnum>(
74 : "semiaxis_variation_type",
75 : rand_options,
76 : "Type of distribution that random superellipsoid semiaxes will follow");
77 30 : return params;
78 30 : }
79 :
80 16 : MultiSmoothSuperellipsoidIC::MultiSmoothSuperellipsoidIC(const InputParameters & parameters)
81 : : SmoothSuperellipsoidBaseIC(parameters),
82 16 : _max_num_tries(getParam<unsigned int>("max_num_tries")),
83 32 : _semiaxis_variation_type(getParam<MooseEnum>("semiaxis_variation_type")),
84 32 : _prevent_overlap(getParam<bool>("prevent_overlap")),
85 32 : _check_extremes(getParam<bool>("check_extremes")),
86 32 : _vary_axes_independently(getParam<bool>("vary_axes_independently")),
87 16 : _numbub(parameters.get<std::vector<unsigned int>>("numbub")),
88 16 : _bubspac(parameters.get<std::vector<Real>>("bubspac")),
89 16 : _exponent(parameters.get<std::vector<Real>>("exponent")),
90 16 : _semiaxis_a(parameters.get<std::vector<Real>>("semiaxis_a")),
91 16 : _semiaxis_b(parameters.get<std::vector<Real>>("semiaxis_b")),
92 16 : _semiaxis_c(parameters.get<std::vector<Real>>("semiaxis_c")),
93 16 : _semiaxis_a_variation(parameters.get<std::vector<Real>>("semiaxis_a_variation")),
94 16 : _semiaxis_b_variation(parameters.get<std::vector<Real>>("semiaxis_b_variation")),
95 32 : _semiaxis_c_variation(parameters.get<std::vector<Real>>("semiaxis_c_variation"))
96 : {
97 16 : }
98 :
99 : void
100 14 : MultiSmoothSuperellipsoidIC::initialSetup()
101 : {
102 14 : unsigned int nv = _numbub.size();
103 :
104 14 : if (nv != _bubspac.size() || nv != _exponent.size() || nv != _semiaxis_a.size() ||
105 28 : 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 14 : if (_semiaxis_variation_type != 2 &&
110 14 : (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 14 : 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 42 : for (_gk = 0; _gk < nv; ++_gk)
123 : {
124 : // Set up domain bounds with mesh tools
125 112 : for (const auto i : make_range(Moose::dim))
126 : {
127 84 : _bottom_left(i) = _mesh.getMinInDimension(i);
128 84 : _top_right(i) = _mesh.getMaxInDimension(i);
129 : }
130 28 : _range = _top_right - _bottom_left;
131 :
132 28 : SmoothSuperellipsoidBaseIC::initialSetup();
133 : }
134 14 : }
135 :
136 : void
137 28 : MultiSmoothSuperellipsoidIC::computeSuperellipsoidSemiaxes()
138 : {
139 : Real randnum;
140 28 : unsigned int start = _as.size();
141 28 : _as.resize(start + _numbub[_gk]);
142 28 : _bs.resize(start + _numbub[_gk]);
143 28 : _cs.resize(start + _numbub[_gk]);
144 :
145 168 : for (unsigned int i = start; i < _as.size(); i++)
146 : {
147 140 : switch (_semiaxis_variation_type)
148 : {
149 140 : case 0: // Uniform distrubtion
150 140 : randnum = _random.rand(_tid);
151 140 : _as[i] = _semiaxis_a[_gk] * (1.0 + (1.0 - 2.0 * randnum) * _semiaxis_a_variation[_gk]);
152 280 : _bs[i] = _semiaxis_b[_gk] *
153 140 : (1.0 + (1.0 - 2.0 * (_vary_axes_independently ? _random.rand(_tid) : randnum)) *
154 140 : _semiaxis_b_variation[_gk]);
155 280 : _cs[i] = _semiaxis_c[_gk] *
156 140 : (1.0 + (1.0 - 2.0 * (_vary_axes_independently ? _random.rand(_tid) : randnum)) *
157 140 : _semiaxis_c_variation[_gk]);
158 140 : 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 280 : _as[i] = _as[i] < 0.0 ? 0.0 : _as[i];
178 280 : _bs[i] = _bs[i] < 0.0 ? 0.0 : _bs[i];
179 280 : _cs[i] = _cs[i] < 0.0 ? 0.0 : _cs[i];
180 : }
181 28 : }
182 :
183 : void
184 28 : MultiSmoothSuperellipsoidIC::computeSuperellipsoidCenters()
185 : {
186 28 : unsigned int start = _centers.size();
187 28 : _centers.resize(start + _numbub[_gk]);
188 :
189 168 : for (unsigned int i = start; i < _centers.size(); i++)
190 : {
191 : // Vary circle center positions
192 : unsigned int num_tries = 0;
193 196 : while (num_tries < _max_num_tries)
194 : {
195 196 : num_tries++;
196 :
197 : RealTensorValue ran;
198 196 : ran(0, 0) = _random.rand(_tid);
199 196 : ran(1, 1) = _random.rand(_tid);
200 196 : ran(2, 2) = _random.rand(_tid);
201 :
202 196 : _centers[i] = _bottom_left + ran * _range;
203 :
204 1015 : for (unsigned int j = 0; j < i; ++j)
205 1743 : if (_mesh.minPeriodicDistance(_var, _centers[i], _centers[j]) < _bubspac[_gk] ||
206 868 : ellipsoidsOverlap(i, j))
207 56 : goto fail;
208 :
209 : // accept the position of the new center
210 140 : 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 140 : }
223 28 : }
224 :
225 : void
226 28 : MultiSmoothSuperellipsoidIC::computeSuperellipsoidExponents()
227 : {
228 28 : unsigned int start = _ns.size();
229 28 : _ns.resize(start + _numbub[_gk]);
230 :
231 168 : for (unsigned int i = start; i < _ns.size(); ++i)
232 140 : _ns[i] = _exponent[_gk];
233 28 : }
234 :
235 : bool
236 868 : MultiSmoothSuperellipsoidIC::ellipsoidsOverlap(unsigned int i, unsigned int j)
237 : {
238 : // Check for overlap between centers
239 868 : const Point dist_vec = _mesh.minPeriodicVector(_var, _centers[i], _centers[j]);
240 868 : const Real dist = dist_vec.norm();
241 :
242 : // Handle this case independently because we cannot calculate polar angles at this point
243 868 : 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 868 : rmn = (std::pow(std::abs(dist_vec(0) / dist / _as[i]), _ns[i]) +
254 868 : std::pow(std::abs(dist_vec(1) / dist / _bs[i]), _ns[i]) +
255 868 : std::pow(std::abs(dist_vec(2) / dist / _cs[i]), _ns[i]));
256 : // calculate r2 from rmn
257 868 : const Real r1 = std::pow(rmn, (-1.0 / _ns[i]));
258 :
259 : // calculate rmn = r^(-n), replacing sin, cos functions with distances
260 868 : rmn = (std::pow(std::abs(dist_vec(0) / dist / _as[j]), _ns[j]) +
261 868 : std::pow(std::abs(dist_vec(1) / dist / _bs[j]), _ns[j]) +
262 868 : std::pow(std::abs(dist_vec(2) / dist / _cs[j]), _ns[j]));
263 868 : const Real r2 = std::pow(rmn, (-1.0 / _ns[j]));
264 :
265 868 : 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 819 : if (_check_extremes)
271 357 : return checkExtremes(i, j) || checkExtremes(j, i);
272 :
273 : // otherwise no overlap has been detected
274 : return false;
275 : }
276 :
277 : bool
278 714 : MultiSmoothSuperellipsoidIC::checkExtremes(unsigned int i, unsigned int j)
279 : {
280 : Point tmp_p;
281 4998 : for (unsigned int pc = 0; pc < 6; pc++)
282 : {
283 4284 : tmp_p = _centers[j];
284 : // Find extremes along semiaxis of candidate ellipsoids
285 : if (pc == 0)
286 714 : tmp_p(0) -= _as[j];
287 : else if (pc == 1)
288 714 : tmp_p(0) += _as[j];
289 : else if (pc == 2)
290 714 : tmp_p(1) -= _bs[j];
291 : else if (pc == 3)
292 714 : tmp_p(1) += _bs[j];
293 : else if (pc == 4)
294 714 : tmp_p(2) -= _cs[j];
295 : else
296 714 : tmp_p(2) += _cs[j];
297 :
298 4284 : const Point dist_vec = _mesh.minPeriodicVector(_var, _centers[i], tmp_p);
299 4284 : const Real dist = dist_vec.norm();
300 :
301 : // Handle this case independently because we cannot calculate polar angles at this point
302 4284 : if (MooseUtils::absoluteFuzzyEqual(dist, 0.0))
303 0 : return true;
304 :
305 : // calculate rmn = r^(-n), replacing sin, cos functions with distances
306 4284 : Real rmn = (std::pow(std::abs(dist_vec(0) / dist / _as[i]), _ns[i]) +
307 4284 : std::pow(std::abs(dist_vec(1) / dist / _bs[i]), _ns[i]) +
308 4284 : std::pow(std::abs(dist_vec(2) / dist / _cs[i]), _ns[i]));
309 4284 : Real r = std::pow(rmn, (-1.0 / _ns[i]));
310 :
311 4284 : if (dist < r)
312 : return true;
313 : }
314 :
315 : return false;
316 : }
|
