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 : #include "SmoothSuperellipsoidBaseIC.h"
11 :
12 : // MOOSE includes
13 : #include "MooseMesh.h"
14 : #include "MooseVariable.h"
15 : #include "FEProblemBase.h"
16 :
17 : InputParameters
18 113 : SmoothSuperellipsoidBaseIC::validParams()
19 : {
20 113 : InputParameters params = InitialCondition::validParams();
21 226 : params.addRequiredParam<Real>("invalue", "The variable value inside the superellipsoid");
22 226 : params.addRequiredParam<Real>("outvalue", "The variable value outside the superellipsoid");
23 226 : params.addParam<Real>(
24 : "nestedvalue",
25 : "The variable value for nested particles inside the superellipsoid in inverse configuration");
26 226 : params.addParam<Real>(
27 226 : "int_width", 0.0, "The interfacial width of the void surface. Defaults to sharp interface");
28 226 : params.addParam<bool>("zero_gradient",
29 226 : false,
30 : "Set the gradient DOFs to zero. This can avoid "
31 : "numerical problems with higher order shape "
32 : "functions.");
33 226 : params.addParam<unsigned int>("rand_seed", 12345, "Seed value for the random number generator");
34 113 : return params;
35 0 : }
36 :
37 60 : SmoothSuperellipsoidBaseIC::SmoothSuperellipsoidBaseIC(const InputParameters & parameters)
38 : : InitialCondition(parameters),
39 60 : _mesh(_fe_problem.mesh()),
40 60 : _invalue(parameters.get<Real>("invalue")),
41 60 : _outvalue(parameters.get<Real>("outvalue")),
42 234 : _nestedvalue(isParamValid("nestedvalue") ? parameters.get<Real>("nestedvalue")
43 114 : : parameters.get<Real>("outvalue")),
44 60 : _int_width(parameters.get<Real>("int_width")),
45 120 : _zero_gradient(parameters.get<bool>("zero_gradient"))
46 : {
47 120 : _random.seed(_tid, getParam<unsigned int>("rand_seed"));
48 60 : }
49 :
50 : void
51 69 : SmoothSuperellipsoidBaseIC::initialSetup()
52 : {
53 : // Compute centers, semiaxes, exponents, and initialize vector sizes
54 69 : computeSuperellipsoidSemiaxes();
55 69 : computeSuperellipsoidExponents();
56 69 : computeSuperellipsoidCenters();
57 :
58 69 : if (_centers.size() != _as.size())
59 0 : mooseError("_center and semiaxis _as vectors are not the same size in the Superellipsoid IC");
60 69 : if (_centers.size() != _bs.size())
61 0 : mooseError("_center and semiaxis _bs vectors are not the same size in the Superellipsoid IC");
62 69 : if (_centers.size() != _cs.size())
63 0 : mooseError("_center and semiaxis _cs vectors are not the same size in the Superellipsoid IC");
64 69 : if (_centers.size() != _ns.size())
65 0 : mooseError("_center and exponent _ns vectors are not the same size in the Superellipsoid IC");
66 :
67 69 : if (_centers.size() < 1)
68 0 : mooseError("_centers, _as, _bs, _cs, and _ns were not initialized in the Superellipsoid IC");
69 69 : }
70 :
71 : Real
72 763104 : SmoothSuperellipsoidBaseIC::value(const Point & p)
73 : {
74 763104 : Real value = _outvalue;
75 : Real val2 = 0.0;
76 :
77 6608628 : for (unsigned int ellip = 0; ellip < _centers.size() && value != _invalue; ++ellip)
78 : {
79 5845524 : val2 = computeSuperellipsoidValue(
80 : p, _centers[ellip], _as[ellip], _bs[ellip], _cs[ellip], _ns[ellip]);
81 5845524 : if ((val2 > value && _invalue > _outvalue) || (val2 < value && _outvalue > _invalue))
82 : value = val2;
83 : }
84 :
85 763104 : return value;
86 : }
87 :
88 : RealGradient
89 120000 : SmoothSuperellipsoidBaseIC::gradient(const Point & p)
90 : {
91 120000 : if (_zero_gradient)
92 : return 0.0;
93 :
94 : RealGradient gradient = 0.0;
95 120000 : Real value = _outvalue;
96 : Real val2 = 0.0;
97 :
98 690000 : for (unsigned int ellip = 0; ellip < _centers.size(); ++ellip)
99 : {
100 570000 : val2 = computeSuperellipsoidValue(
101 : p, _centers[ellip], _as[ellip], _bs[ellip], _cs[ellip], _ns[ellip]);
102 570000 : if ((val2 > value && _invalue > _outvalue) || (val2 < value && _outvalue > _invalue))
103 : {
104 : value = val2;
105 51456 : gradient = computeSuperellipsoidGradient(
106 : p, _centers[ellip], _as[ellip], _bs[ellip], _cs[ellip], _ns[ellip]);
107 : }
108 : }
109 :
110 120000 : return gradient;
111 : }
112 :
113 : Real
114 6425524 : SmoothSuperellipsoidBaseIC::computeSuperellipsoidValue(
115 : const Point & p, const Point & center, Real a, Real b, Real c, Real n)
116 : {
117 6425524 : Point l_center = center;
118 6425524 : Point l_p = p;
119 : // Compute the distance between the current point and the center
120 6425524 : 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 6425524 : if (dist == 0.0)
125 0 : 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 6425524 : Point dist_vec = _mesh.minPeriodicVector(_var.number(), center, p);
132 :
133 : // First calculate rmn = r^(-n), replacing sin, cos functions with distances
134 6425524 : Real rmn = (std::pow(std::abs(dist_vec(0) / dist / a), n) +
135 6425524 : std::pow(std::abs(dist_vec(1) / dist / b), n) +
136 6425524 : std::pow(std::abs(dist_vec(2) / dist / c), n));
137 : // Then calculate r from rmn
138 6425524 : Real r = std::pow(rmn, (-1.0 / n));
139 :
140 6425524 : Real value = _outvalue; // Outside superellipsoid
141 :
142 6425524 : if (dist <= r - _int_width / 2.0) // Inside superellipsoid
143 74256 : value = _invalue;
144 6351268 : else if (dist < r + _int_width / 2.0) // Smooth interface
145 : {
146 38336 : Real int_pos = (dist - r + _int_width / 2.0) / _int_width;
147 38336 : 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
155 80000 : SmoothSuperellipsoidBaseIC::computeSuperellipsoidInverseValue(
156 : const Point & p, const Point & center, Real a, Real b, Real c, Real n)
157 : {
158 80000 : Point l_center = center;
159 80000 : Point l_p = p;
160 : // Compute the distance between the current point and the center
161 80000 : 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 80000 : if (dist == 0.0)
166 0 : 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 80000 : Point dist_vec = _mesh.minPeriodicVector(_var.number(), center, p);
173 :
174 : // First calculate rmn = r^(-n), replacing sin, cos functions with distances
175 80000 : Real rmn = (std::pow(std::abs(dist_vec(0) / dist / a), n) +
176 80000 : std::pow(std::abs(dist_vec(1) / dist / b), n) +
177 80000 : std::pow(std::abs(dist_vec(2) / dist / c), n));
178 : // Then calculate r from rmn
179 80000 : Real r = std::pow(rmn, (-1.0 / n));
180 :
181 80000 : Real value = _invalue;
182 :
183 80000 : if (dist <= r - _int_width / 2.0) // Reversing inside and outside values
184 896 : value = _nestedvalue;
185 79104 : else if (dist < r + _int_width / 2.0) // Smooth interface
186 : {
187 0 : Real int_pos = (dist - r + _int_width / 2.0) / _int_width;
188 0 : value = _invalue + (_nestedvalue - _invalue) * (1.0 + std::cos(int_pos * libMesh::pi)) / 2.0;
189 : }
190 :
191 : return value;
192 : }
193 :
194 : RealGradient
195 51456 : SmoothSuperellipsoidBaseIC::computeSuperellipsoidGradient(
196 : const Point & p, const Point & center, Real a, Real b, Real c, Real n)
197 : {
198 51456 : Point l_center = center;
199 51456 : Point l_p = p;
200 : // Compute the distance between the current point and the center
201 51456 : 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 51456 : 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 51456 : Point dist_vec = _mesh.minPeriodicVector(_var.number(), center, p);
213 : // First calculate rmn = r^(-n)
214 51456 : Real rmn = (std::pow(std::abs(dist_vec(0) / dist / a), n) +
215 51456 : std::pow(std::abs(dist_vec(1) / dist / b), n) +
216 51456 : std::pow(std::abs(dist_vec(2) / dist / c), n));
217 : // Then calculate r from rmn
218 51456 : Real r = std::pow(rmn, (-1.0 / n));
219 :
220 : Real DvalueDr = 0.0;
221 :
222 51456 : if (dist < r + _int_width / 2.0 && dist > r - _int_width / 2.0) // in interfacial region
223 : {
224 21216 : Real int_pos = (dist - r + _int_width / 2.0) / _int_width;
225 21216 : Real Dint_posDr = 1.0 / _int_width;
226 21216 : DvalueDr = Dint_posDr * (_invalue - _outvalue) *
227 21216 : (-std::sin(int_pos * libMesh::pi) * libMesh::pi) / 2.0;
228 : }
229 :
230 51456 : return dist_vec * (DvalueDr / dist);
231 : }
|
