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