www.mooseframework.org
SmoothSuperellipsoidBaseIC.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 
9 
10 // MOOSE includes
11 #include "MooseMesh.h"
12 #include "MooseVariable.h"
13 
14 template <>
15 InputParameters
17 {
18  InputParameters params = validParams<InitialCondition>();
19  params.addRequiredParam<Real>("invalue", "The variable value inside the superellipsoid");
20  params.addRequiredParam<Real>("outvalue", "The variable value outside the superellipsoid");
21  params.addParam<Real>(
22  "nestedvalue",
23  "The variable value for nested particles inside the superellipsoid in inverse configuration");
24  params.addParam<Real>(
25  "int_width", 0.0, "The interfacial width of the void surface. Defaults to sharp interface");
26  params.addParam<bool>("zero_gradient",
27  false,
28  "Set the gradient DOFs to zero. This can avoid "
29  "numerical problems with higher order shape "
30  "functions.");
31  params.addParam<unsigned int>("rand_seed", 12345, "Seed value for the random number generator");
32  return params;
33 }
34 
35 SmoothSuperellipsoidBaseIC::SmoothSuperellipsoidBaseIC(const InputParameters & parameters)
36  : InitialCondition(parameters),
37  _mesh(_fe_problem.mesh()),
38  _invalue(parameters.get<Real>("invalue")),
39  _outvalue(parameters.get<Real>("outvalue")),
40  _nestedvalue(isParamValid("nestedvalue") ? parameters.get<Real>("nestedvalue")
41  : parameters.get<Real>("outvalue")),
42  _int_width(parameters.get<Real>("int_width")),
43  _zero_gradient(parameters.get<bool>("zero_gradient"))
44 {
45  _random.seed(_tid, getParam<unsigned int>("rand_seed"));
46 }
47 
48 void
50 {
51  // Compute centers, semiaxes, exponents, and initialize vector sizes
55 
56  if (_centers.size() != _as.size())
57  mooseError("_center and semiaxis _as vectors are not the same size in the Superellipsoid IC");
58  if (_centers.size() != _bs.size())
59  mooseError("_center and semiaxis _bs vectors are not the same size in the Superellipsoid IC");
60  if (_centers.size() != _cs.size())
61  mooseError("_center and semiaxis _cs vectors are not the same size in the Superellipsoid IC");
62  if (_centers.size() != _ns.size())
63  mooseError("_center and exponent _ns vectors are not the same size in the Superellipsoid IC");
64 
65  if (_centers.size() < 1)
66  mooseError("_centers, _as, _bs, _cs, and _ns were not initialized in the Superellipsoid IC");
67 }
68 
69 Real
71 {
72  Real value = _outvalue;
73  Real val2 = 0.0;
74 
75  for (unsigned int ellip = 0; ellip < _centers.size() && value != _invalue; ++ellip)
76  {
78  p, _centers[ellip], _as[ellip], _bs[ellip], _cs[ellip], _ns[ellip]);
79  if ((val2 > value && _invalue > _outvalue) || (val2 < value && _outvalue > _invalue))
80  value = val2;
81  }
82 
83  return value;
84 }
85 
86 RealGradient
88 {
89  if (_zero_gradient)
90  return 0.0;
91 
92  RealGradient gradient = 0.0;
93  Real value = _outvalue;
94  Real val2 = 0.0;
95 
96  for (unsigned int ellip = 0; ellip < _centers.size(); ++ellip)
97  {
99  p, _centers[ellip], _as[ellip], _bs[ellip], _cs[ellip], _ns[ellip]);
100  if ((val2 > value && _invalue > _outvalue) || (val2 < value && _outvalue > _invalue))
101  {
102  value = val2;
104  p, _centers[ellip], _as[ellip], _bs[ellip], _cs[ellip], _ns[ellip]);
105  }
106  }
107 
108  return gradient;
109 }
110 
111 Real
113  const Point & p, const Point & center, Real a, Real b, Real c, Real n)
114 {
115  Point l_center = center;
116  Point l_p = p;
117  // Compute the distance between the current point and the center
118  Real dist = _mesh.minPeriodicDistance(_var.number(), l_p, l_center);
119 
120  // When dist is 0 we are exactly at the center of the superellipsoid so return _invalue
121  // Handle this case independently because we cannot calculate polar angles at this point
122  if (dist == 0.0)
123  return _invalue;
124 
125  // Compute the distance r from the center of the superellipsoid to its outside edge
126  // along the vector from the center to the current point
127  // This uses the equation for a superellipse in polar coordinates and substitutes
128  // distances for sin, cos functions
129  Point dist_vec = _mesh.minPeriodicVector(_var.number(), center, p);
130 
131  // First calculate rmn = r^(-n), replacing sin, cos functions with distances
132  Real rmn = (std::pow(std::abs(dist_vec(0) / dist / a), n) +
133  std::pow(std::abs(dist_vec(1) / dist / b), n) +
134  std::pow(std::abs(dist_vec(2) / dist / c), n));
135  // Then calculate r from rmn
136  Real r = std::pow(rmn, (-1.0 / n));
137 
138  Real value = _outvalue; // Outside superellipsoid
139 
140  if (dist <= r - _int_width / 2.0) // Inside superellipsoid
141  value = _invalue;
142  else if (dist < r + _int_width / 2.0) // Smooth interface
143  {
144  Real int_pos = (dist - r + _int_width / 2.0) / _int_width;
145  value = _outvalue + (_invalue - _outvalue) * (1.0 + std::cos(int_pos * libMesh::pi)) / 2.0;
146  }
147 
148  return value;
149 }
150 
151 // Following function does the same as computeSuperellipsoidValue but reverses invalue and outvalue
152 Real
154  const Point & p, const Point & center, Real a, Real b, Real c, Real n)
155 {
156  Point l_center = center;
157  Point l_p = p;
158  // Compute the distance between the current point and the center
159  Real dist = _mesh.minPeriodicDistance(_var.number(), l_p, l_center);
160 
161  // When dist is 0 we are exactly at the center of the superellipsoid so return _invalue
162  // Handle this case independently because we cannot calculate polar angles at this point
163  if (dist == 0.0)
164  return _nestedvalue;
165 
166  // Compute the distance r from the center of the superellipsoid to its outside edge
167  // along the vector from the center to the current point
168  // This uses the equation for a superellipse in polar coordinates and substitutes
169  // distances for sin, cos functions
170  Point dist_vec = _mesh.minPeriodicVector(_var.number(), center, p);
171 
172  // First calculate rmn = r^(-n), replacing sin, cos functions with distances
173  Real rmn = (std::pow(std::abs(dist_vec(0) / dist / a), n) +
174  std::pow(std::abs(dist_vec(1) / dist / b), n) +
175  std::pow(std::abs(dist_vec(2) / dist / c), n));
176  // Then calculate r from rmn
177  Real r = std::pow(rmn, (-1.0 / n));
178 
179  Real value = _invalue;
180 
181  if (dist <= r - _int_width / 2.0) // Reversing inside and outside values
182  value = _nestedvalue;
183  else if (dist < r + _int_width / 2.0) // Smooth interface
184  {
185  Real int_pos = (dist - r + _int_width / 2.0) / _int_width;
186  value = _invalue + (_nestedvalue - _invalue) * (1.0 + std::cos(int_pos * libMesh::pi)) / 2.0;
187  }
188 
189  return value;
190 }
191 
192 RealGradient
194  const Point & p, const Point & center, Real a, Real b, Real c, Real n)
195 {
196  Point l_center = center;
197  Point l_p = p;
198  // Compute the distance between the current point and the center
199  Real dist = _mesh.minPeriodicDistance(_var.number(), l_p, l_center);
200 
201  // When dist is 0 we are exactly at the center of the superellipsoid so return 0
202  // Handle this case independently because we cannot calculate polar angles at this point
203  if (dist == 0.0)
204  return 0.0;
205 
206  // Compute the distance r from the center of the superellipsoid to its outside edge
207  // along the vector from the center to the current point
208  // This uses the equation for a superellipse in polar coordinates and substitutes
209  // distances for sin, cos functions
210  Point dist_vec = _mesh.minPeriodicVector(_var.number(), center, p);
211  // First calculate rmn = r^(-n)
212  Real rmn = (std::pow(std::abs(dist_vec(0) / dist / a), n) +
213  std::pow(std::abs(dist_vec(1) / dist / b), n) +
214  std::pow(std::abs(dist_vec(2) / dist / c), n));
215  // Then calculate r from rmn
216  Real r = std::pow(rmn, (-1.0 / n));
217 
218  Real DvalueDr = 0.0;
219 
220  if (dist < r + _int_width / 2.0 && dist > r - _int_width / 2.0) // in interfacial region
221  {
222  Real int_pos = (dist - r + _int_width / 2.0) / _int_width;
223  Real Dint_posDr = 1.0 / _int_width;
224  DvalueDr = Dint_posDr * (_invalue - _outvalue) *
225  (-std::sin(int_pos * libMesh::pi) * libMesh::pi) / 2.0;
226  }
227 
228  return dist_vec * (DvalueDr / dist);
229 }
virtual RealGradient gradient(const Point &p)
virtual void computeSuperellipsoidSemiaxes()=0
virtual Real value(const Point &p)
InputParameters validParams< SmoothSuperellipsoidBaseIC >()
RealGradient computeSuperellipsoidGradient(const Point &p, const Point &center, Real a, Real b, Real c, Real n)
virtual Real computeSuperellipsoidValue(const Point &p, const Point &center, Real a, Real b, Real c, Real n)
virtual Real computeSuperellipsoidInverseValue(const Point &p, const Point &center, Real a, Real b, Real c, Real n)
virtual void computeSuperellipsoidCenters()=0
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
virtual void computeSuperellipsoidExponents()=0
SmoothSuperellipsoidBaseIC(const InputParameters &parameters)