www.mooseframework.org
RankThreeTensor.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 #include "RankThreeTensor.h"
8 #include "RankTwoTensor.h"
9 #include "RankFourTensor.h"
10 
11 // MOOSE includes
12 #include "MooseEnum.h"
13 #include "MooseException.h"
14 #include "MooseUtils.h"
15 #include "MatrixTools.h"
16 #include "MaterialProperty.h"
17 #include "PermutationTensor.h"
18 
19 #include "libmesh/utility.h"
20 
21 // C++ includes
22 #include <iomanip>
23 #include <ostream>
24 
25 template <>
26 void
28 {
29  v.zero();
30 }
31 
32 template <>
33 void
34 dataStore(std::ostream & stream, RankThreeTensor & rtht, void * context)
35 {
36  dataStore(stream, rtht._vals, context);
37 }
38 
39 template <>
40 void
41 dataLoad(std::istream & stream, RankThreeTensor & rtht, void * context)
42 {
43  dataLoad(stream, rtht._vals, context);
44 }
45 
46 MooseEnum RankThreeTensor::fillMethodEnum() // TODO: Need new fillMethodEnum() -- for now we will
47  // just use general (at most 27 components)
48 {
49  return MooseEnum("general");
50 }
51 
53 {
54  mooseAssert(N == 3, "RankThreeTensor is currently only tested for 3 dimensions.");
55 
56  for (unsigned int i = 0; i < N3; ++i)
57  _vals[i] = 0;
58 }
59 
61 {
62  switch (init)
63  {
64  case initNone:
65  break;
66 
67  default:
68  mooseError("Unknown RankThreeTensor initialization pattern.");
69  }
70 }
71 
72 RankThreeTensor::RankThreeTensor(const std::vector<Real> & input, FillMethod fill_method)
73 {
74  fillFromInputVector(input, fill_method);
75 }
76 
77 void
79 {
80  for (unsigned int i = 0; i < N3; ++i)
81  _vals[i] = 0;
82 }
83 
86 {
87  for (unsigned int i = 0; i < N3; ++i)
88  _vals[i] = a._vals[i];
89 
90  return *this;
91 }
92 
94 {
95  RealVectorValue result;
96 
97  for (unsigned int i = 0; i < N; ++i)
98  {
99  Real sum = 0;
100  unsigned int i1 = i * N2;
101  for (unsigned int j1 = 0; j1 < N2; j1 += N)
102  for (unsigned int k = 0; k < N; ++k)
103  sum += _vals[i1 + j1 + k] * a._vals[j1 + k];
104  result(i) = sum;
105  }
106 
107  return result;
108 }
109 
111 {
112  RankThreeTensor result;
113 
114  for (unsigned int i = 0; i < N3; ++i)
115  result._vals[i] = _vals[i] * b;
116 
117  return result;
118 }
119 
122 {
123  for (unsigned int i = 0; i < N3; ++i)
124  _vals[i] *= a;
125 
126  return *this;
127 }
128 
130 RankThreeTensor::operator/(const Real b) const
131 {
132  RankThreeTensor result;
133 
134  for (unsigned int i = 0; i < N3; ++i)
135  result._vals[i] = _vals[i] / b;
136 
137  return result;
138 }
139 
142 {
143  for (unsigned int i = 0; i < N3; ++i)
144  _vals[i] /= a;
145 
146  return *this;
147 }
148 
151 {
152  for (unsigned int i = 0; i < N3; ++i)
153  _vals[i] += a._vals[i];
154 
155  return *this;
156 }
157 
160 {
161  RankThreeTensor result;
162 
163  for (unsigned int i = 0; i < N3; ++i)
164  result._vals[i] = _vals[i] + b._vals[i];
165 
166  return result;
167 }
168 
171 {
172  for (unsigned int i = 0; i < N3; ++i)
173  _vals[i] -= a._vals[i];
174 
175  return *this;
176 }
177 
180 {
181  RankThreeTensor result;
182 
183  for (unsigned int i = 0; i < N3; ++i)
184  result._vals[i] = _vals[i] - b._vals[i];
185 
186  return result;
187 }
188 
191 {
192  RankThreeTensor result;
193 
194  for (unsigned int i = 0; i < N3; ++i)
195  result._vals[i] = -_vals[i];
196 
197  return result;
198 }
199 
200 Real
202 {
203  Real l2 = 0;
204 
205  for (unsigned int i = 0; i < N3; ++i)
206  l2 += Utility::pow<2>(_vals[i]);
207 
208  return std::sqrt(l2);
209 }
210 
211 void
212 RankThreeTensor::fillFromInputVector(const std::vector<Real> & input, FillMethod fill_method)
213 {
214  zero();
215 
216  switch (fill_method)
217  {
218  case general:
220  break;
221  default:
222  mooseError("fillFromInputVector called with unknown fill_method of ", fill_method);
223  }
224 }
225 
226 void
228 {
229  unsigned int index = 0;
230  for (unsigned int i = 0; i < N; ++i)
231  {
232  const Real a = input(i);
233  for (unsigned int j = 0; j < N; ++j)
234  {
235  const Real b = input(j);
236  for (unsigned int k = 0; k < N; ++k)
237  {
238  const Real c = input(k);
239  Real sum = 0;
240  sum = -2 * a * b * c;
241  if (i == j)
242  sum += c;
243  if (i == k)
244  sum += b;
245  _vals[index++] = sum / 2.0;
246  }
247  }
248  }
249 }
250 
253 {
254  RankFourTensor result;
255 
256  unsigned int index = 0;
257  for (unsigned int i = 0; i < N; ++i)
258  for (unsigned int j = 0; j < N; ++j)
259  for (unsigned int k = 0; k < N; ++k)
260  for (unsigned int l = 0; l < N; ++l)
261  {
262  for (unsigned int m = 0; m < N; ++m)
263  for (unsigned int n = 0; n < N; ++n)
264  result._vals[index] += (*this)(m, i, j) * a._vals[m * N + n] * (*this)(n, k, l);
265  index++;
266  }
267 
268  return result;
269 }
270 
271 void
273 {
274  RankThreeTensor old = *this;
275 
276  unsigned int index = 0;
277  for (unsigned int i = 0; i < N; ++i)
278  for (unsigned int j = 0; j < N; ++j)
279  for (unsigned int k = 0; k < N; ++k)
280  {
281  Real sum = 0.0;
282  unsigned int index2 = 0;
283  for (unsigned int m = 0; m < N; ++m)
284  {
285  Real a = R(i, m);
286  for (unsigned int n = 0; n < N; ++n)
287  {
288  Real ab = a * R(j, n);
289  for (unsigned int o = 0; o < N; ++o)
290  sum += ab * R(k, o) * old._vals[index2++];
291  }
292  }
293  _vals[index++] = sum;
294  }
295 }
296 
297 void
299 {
300  RankThreeTensor old = *this;
301 
302  unsigned int index = 0;
303  unsigned int i1 = 0;
304  for (unsigned int i = 0; i < N; ++i)
305  {
306  unsigned int j1 = 0;
307  for (unsigned int j = 0; j < N; ++j)
308  {
309  unsigned int k1 = 0;
310  for (unsigned int k = 0; k < N; ++k)
311  {
312  Real sum = 0.0;
313  unsigned int index2 = 0;
314  for (unsigned int m = 0; m < N; ++m)
315  {
316  Real a = R._vals[i1 + m];
317  for (unsigned int n = 0; n < N; ++n)
318  {
319  Real ab = a * R._vals[j1 + n];
320  for (unsigned int o = 0; o < N; ++o)
321  sum += ab * R._vals[k1 + o] * old._vals[index2++];
322  }
323  }
324  _vals[index++] = sum;
325  k1 += N;
326  }
327  j1 += N;
328  }
329  i1 += N;
330  }
331 }
332 
333 void
334 RankThreeTensor::fillGeneralFromInputVector(const std::vector<Real> & input)
335 {
336  if (input.size() != 27)
337  mooseError("To use fillGeneralFromInputVector, your input must have size 27. Yours has size ",
338  input.size());
339 
340  for (unsigned int i = 0; i < N3; ++i)
341  _vals[i] = input[i];
342 }
343 
345 {
346  RankTwoTensor result;
347 
348  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
349  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
350  for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
351  result(i, j) += p(k) * b(k, i, j);
352 
353  return result;
354 }
355 
358 {
359  RealVectorValue result;
360 
361  for (unsigned int i = 0; i < N; ++i)
362  for (unsigned int j = 0; j < N2; ++j)
363  result(i) += _vals[i * N2 + j] * b._vals[j];
364 
365  return result;
366 }
RankThreeTensor & operator+=(const RankThreeTensor &a)
r_ijk += a_ijk for all i, j, k
InitMethod
Initialization method.
RankFourTensor mixedProductRankFour(const RankTwoTensor &a) const
Creates fourth order tensor D=A*b*A where A is rank 3 and b is rank 2.
VectorValue< Real > RealVectorValue
Definition: Assembly.h:40
RankThreeTensor & operator=(const RankThreeTensor &a)
copies values from a into this tensor
void fillFromPlaneNormal(const RealVectorValue &input)
Fills RankThreeTensor from plane normal vectors ref.
Real L2norm() const
(r_ijk*r_ijk)
static constexpr unsigned int N2
Real _vals[N2]
The values of the rank-two tensor stored by index=(i * LIBMESH_DIM + j)
static MooseEnum fillMethodEnum()
Static method for use in validParams for getting the "fill_method".
void dataStore(std::ostream &stream, RankThreeTensor &rtht, void *context)
void mooseError(Args &&...args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:182
FillMethod
To fill up the 27 entries in the 3rd-order tensor, fillFromInputVector is called with one of the foll...
static constexpr unsigned int N
Dimensionality of rank-three tensor.
void fillFromInputVector(const std::vector< Real > &input, FillMethod fill_method)
fillFromInputVector takes some number of inputs to fill the Rank-3 tensor.
RankThreeTensor operator+(const RankThreeTensor &a) const
r_ijkl + a_ijk
PetscInt m
RankThreeTensor operator/(const Real a) const
r_ijk/a
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:37
RankTwoTensor is designed to handle the Stress or Strain Tensor for a fully anisotropic material...
Definition: RankTwoTensor.h:45
RankThreeTensor is designed to handle any N-dimensional third order tensor, r.
Real _vals[N4]
The values of the rank-four tensor stored by index=(((i * LIBMESH_DIM + j) * LIBMESH_DIM + k) * LIBME...
Real _vals[N3]
The values of the rank-three tensor stored by index=((i * LIBMESH_DIM + j) * LIBMESH_DIM + k) ...
void fillGeneralFromInputVector(const std::vector< Real > &input)
void zero()
Zeros out the tensor.
RealVectorValue operator*(const RankTwoTensor &a) const
r_ijk*a_kl
void mooseSetToZero< RankThreeTensor >(RankThreeTensor &v)
Helper function template specialization to set an object to zero.
PetscInt n
RealVectorValue doubleContraction(const RankTwoTensor &b) const
Creates a vector from the double contraction of a rank three and rank two tensor. ...
RankThreeTensor & operator/=(const Real a)
r_ijk /= a for all i, j, k
RankThreeTensor()
Default constructor; fills to zero.
void rotate(const T &R)
Rotate the tensor using r_ijk = R_im R_in R_ko r_mno.
static constexpr unsigned int N3
RankThreeTensor operator-() const
-r_ijk
TensorValue< Real > RealTensorValue
Definition: Assembly.h:45
void dataLoad(std::istream &stream, RankThreeTensor &rtht, void *context)
RankThreeTensor & operator-=(const RankThreeTensor &a)
r_ijk -= a_ijk
RankFourTensor is designed to handle any N-dimensional fourth order tensor, C.
RankThreeTensor & operator*=(const Real a)
r_ijk *= a