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