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