Itasca C++ Interface
symtensor.h
Go to the documentation of this file.
1 #pragma once
8 #include "axes.h"
9 #include <cassert>
10 #include <limits>
11 #include <algorithm>
12 
13 class SymTensorInfo;
14 
19 PUSHWARNING
20 VSWARNING(26495) // member init (not initialized on purpose)
21 VSWARNING(4702) // Weird unreachable code warnings in release
22 class BASE_EXPORT SymTensor {
23 public:
25 #ifdef DEBUG
26  SymTensor() { d11_ = d22_ = d33_ = d12_ = d13_ = d23_ = initVal<T>(); }
27 #else
28  SymTensor() { } // NOLINT
29 #endif
31  SymTensor(const SymTensor &s) { operator=(s); }
34  explicit SymTensor(const double &i11,const double &i22=0.0,const double &i33=0.0,
35  const double &i12=0.0,const double &i13=0.0,const double &i23=0.0)
36  : d11_(i11), d22_(i22), d33_(i33), d12_(i12), d13_(i13), d23_(i23) { }
38  const SymTensor &operator=(const SymTensor &s);
39 
41  bool operator==(const SymTensor &s) const;
42  bool operator<(const SymTensor &s) const;
43 
45  const double &s11() const { return d11_; }
46  const double &s22() const { return d22_; }
47  const double &s33() const { return d33_; }
48  const double &s12() const { return d12_; }
49  const double &s13() const { return d13_; }
50  const double &s23() const { return d23_; }
51  const double &s21() const { return d12_; }
52  const double &s31() const { return d13_; }
53  const double &s32() const { return d23_; }
55  double &rs11() { return d11_; }
56  double &rs22() { return d22_; }
57  double &rs33() { return d33_; }
58  double &rs12() { return d12_; }
59  double &rs13() { return d13_; }
60  double &rs23() { return d23_; }
61  double &rs21() { return d12_; }
62  double &rs31() { return d13_; }
63  double &rs32() { return d23_; }
68  inline const double &operator[](unsigned int i) const;
69  inline double &operator[](unsigned int i);
70 
72  const double &operator()(unsigned int i,unsigned int j) const;
74  double &operator()(unsigned int i,unsigned int j);
75 
81  DVect3 getEigenInfo(SymTensorInfo *si=0) const;
83  double getTrace() const { return d11_+d22_+d33_; }
85  SymTensor getDeviatoric() const { double p=getTrace()/3.0; return SymTensor(d11_-p, d22_-p, d33_-p, d12_, d13_, d23_); }
87  double getI1() const { return getTrace(); }
89  double getI2() const { return d11_*d22_ + d22_*d33_ + d11_*d33_ - d12_*d12_ - d23_*d23_ - d13_*d13_; }
91  double getI3() const { return getDeterminate(); }
93  double getJ2() const;
95  double getJ2(SymTensor *dev, double *I1=nullptr);
96  //Keep only getJ2(), Von Mises and Octahedral stress/strain are calculated from J2.
98  //double getVonMises() const { return sqrt(3.0*getJ2()); }
100  //double getVonMisesStrain() const { return sqrt(4.0*getJ2()/9.0); }
102  //double getOctahedralShear() const { return sqrt(2.0*getJ2()/3.0); }
104  //double getOctahedralShearStrain() const { return sqrt(8.0*getJ2()/3.0); }
106  double getJ3() const;
108  double getLode(double *I1 = nullptr, double *J2 = nullptr, double *J3 = nullptr);
110  double getDeterminate() const;
112  double getNorm2() const { return d11_*d11_ + d22_*d22_ + d33_*d33_ + 2.0*(d12_*d12_ + d13_*d13_ + d23_*d23_); }
116  double getTotalMeasure() const;
118  SymTensor toGlobal(const Axes3D &a) const;
123  inline DVect3 operator*(const DVect3 &input) const;
124  inline DVect2 operator*(const DVect2 &input) const;
125  inline SymTensor operator*(const double &mul) const;
126  inline const SymTensor &operator*=(const double &mul);
127  inline SymTensor operator/(const double &mul) const;
128  inline const SymTensor &operator/=(const double &mul);
130  SymTensor mul(const double &d) const { SymTensor ret(d11_*d,d22_*d,d33_*d,d12_*d,d13_*d,d23_*d); return ret; }
132  const SymTensor &operator+=(const SymTensor &s) { d11_+=s.d11_; d22_+=s.d22_; d33_+=s.d33_; d12_+=s.d12_; d13_+=s.d13_; d23_+=s.d23_; return *this; }
133  const SymTensor &operator-=(const SymTensor &s) { d11_-=s.d11_; d22_-=s.d22_; d33_-=s.d33_; d12_-=s.d12_; d13_-=s.d13_; d23_-=s.d23_; return *this; }
134  inline SymTensor operator+(const SymTensor &s) const;
135  inline SymTensor operator-(const SymTensor &s) const;
137  static SymTensor fromPrincipal(const DVect3 &prin,const Axes3D &axes);
141  static SymTensor fromForceNormal(const DVect3 &normal, const DVect3 &force);
142  static inline uint32 doubleToSingleComponent(uint32 dof1,uint32 dof2);
144  bool isDiagonal(double tol = std::numeric_limits<double>::epsilon()*1000.0) const { return (abs(d12_) > tol || abs(d13_) > tol || abs(d23_) > tol) ? false : true; }
145  inline void adjustTrace(const double &newTrace);
146  inline void incrementDiagonal(const double &increment) { d11_ += increment; d22_ += increment; d33_ += increment; }
147  inline void rotate(const DVect3 &rot);
148  inline double maxAbs() const { return std::max(std::abs(d11_),std::max(std::abs(d22_),std::max(std::abs(d33_),std::max(std::abs(d12_),std::max(std::abs(d13_),std::abs(d23_)))))); }
149 private:
150  double d11_;
151  double d22_;
152  double d33_;
153  double d12_;
154  double d13_;
155  double d23_;
156 };
157 POPWARNING
158 
166 class SymTensorInfo {
167 public:
168  friend class SymTensor;
170  BASE_EXPORT SymTensorInfo(): type_(Type::ThreeDCube) { }
172  BASE_EXPORT SymTensorInfo(const SymTensorInfo &si);
174  BASE_EXPORT const SymTensorInfo &operator=(const SymTensorInfo &si);
176  BASE_EXPORT Axes3D getAxes() const;
178  BASE_EXPORT SymTensor resolve(const DVect3 &prin) const;
179 private:
180  enum class Type { ThreeDCube, ThreeDJacobi, ZMax, ZMid, ZMin };
181  Type type_;
182  Axes3D axes_;
183 };
184 
185 inline const SymTensor &SymTensor::operator=(const SymTensor &s) {
186  d11_ = s.d11_;
187  d22_ = s.d22_;
188  d33_ = s.d33_;
189  d12_ = s.d12_;
190  d13_ = s.d13_;
191  d23_ = s.d23_;
192  return *this;
193 }
194 
195 inline bool SymTensor::operator==(const SymTensor &s) const {
196  return (d11_ == s.d11_ && d22_ == s.d22_ &&
197  d33_ == s.d33_ && d12_ == s.d12_ && d13_ == s.d13_ &&
198  d23_ == s.d23_);
199 }
200 
201 inline const double &SymTensor::operator[](unsigned int i) const {
202  assert(i<6);
203  switch (i) {
204  case 0: return d11_;
205  case 1: return d22_;
206  case 2: return d33_;
207  case 3: return d12_;
208  case 4: return d13_;
209  case 5: return d23_;
210  }
211  return d11_;
212 }
213 
214 inline double &SymTensor::operator[](unsigned int i) {
215  assert(i<6);
216  switch (i) {
217  case 0: return d11_;
218  case 1: return d22_;
219  case 2: return d33_;
220  case 3: return d12_;
221  case 4: return d13_;
222  case 5: return d23_;
223  }
224  return d11_;
225 }
226 
227 inline double SymTensor::getTotalMeasure() const {
228  double I1 = getTrace();
229  double J2 = getJ2();
230  return sqrt(I1*I1/3.0 + 2.0*J2);
231 }
232 
233 inline DVect3 SymTensor::operator*(const DVect3 &normal) const {
234  return DVect3(normal.x()*s11() + normal.y()*s12() + normal.z()*s13(),
235  normal.x()*s21() + normal.y()*s22() + normal.z()*s23(),
236  normal.x()*s31() + normal.y()*s32() + normal.z()*s33());
237 }
238 
239 inline DVect2 SymTensor::operator*(const DVect2 &normal) const {
240  return DVect2(normal.x()*s11() + normal.y()*s12(),
241  normal.x()*s21() + normal.y()*s22());
242 }
243 
244 inline SymTensor SymTensor::operator*(const double &mul) const {
245  SymTensor ret;
246  ret.d11_ = d11_ * mul;
247  ret.d22_ = d22_ * mul;
248  ret.d33_ = d33_ * mul;
249  ret.d12_ = d12_ * mul;
250  ret.d13_ = d13_ * mul;
251  ret.d23_ = d23_ * mul;
252  return ret;
253 }
254 
255 inline const SymTensor &SymTensor::operator*=(const double &mul) {
256  d11_ *= mul;
257  d22_ *= mul;
258  d33_ *= mul;
259  d12_ *= mul;
260  d13_ *= mul;
261  d23_ *= mul;
262  return *this;
263 }
264 
265 inline SymTensor SymTensor::operator/(const double &mul) const {
266  SymTensor ret;
267  ret.d11_ = d11_ / mul;
268  ret.d22_ = d22_ / mul;
269  ret.d33_ = d33_ / mul;
270  ret.d12_ = d12_ / mul;
271  ret.d13_ = d13_ / mul;
272  ret.d23_ = d23_ / mul;
273  return ret;
274 }
275 
276 inline const SymTensor &SymTensor::operator/=(const double &mul) {
277  d11_ /= mul;
278  d22_ /= mul;
279  d33_ /= mul;
280  d12_ /= mul;
281  d13_ /= mul;
282  d23_ /= mul;
283  return *this;
284 }
285 
286 inline SymTensor SymTensor::operator+(const SymTensor &s) const {
287  SymTensor ret;
288  ret.d11_ = d11_ + s.d11_;
289  ret.d22_ = d22_ + s.d22_;
290  ret.d33_ = d33_ + s.d33_;
291  ret.d12_ = d12_ + s.d12_;
292  ret.d13_ = d13_ + s.d13_;
293  ret.d23_ = d23_ + s.d23_;
294  return ret;
295 }
296 
297 inline SymTensor SymTensor::operator-(const SymTensor &s) const {
298  SymTensor ret;
299  ret.d11_ = d11_ - s.d11_;
300  ret.d22_ = d22_ - s.d22_;
301  ret.d33_ = d33_ - s.d33_;
302  ret.d12_ = d12_ - s.d12_;
303  ret.d13_ = d13_ - s.d13_;
304  ret.d23_ = d23_ - s.d23_;
305  return ret;
306 }
307 
308 inline void SymTensor::adjustTrace(const double &newTrace) {
309  static constexpr double d1d3 = 1.0 / 3.0;
310  /* --- strain differences --- */
311  double dx = s11() - s22();
312  double dy = s22() - s33();
313  double dz = s33() - s11();
314  rs11() = (newTrace + dx - dz) * d1d3;
315  rs22() = (newTrace + dy - dx) * d1d3;
316  rs33() = (newTrace + dz - dy) * d1d3;
317 }
318 
319 inline void SymTensor::rotate(const DVect3 &rot) {
320  SymTensor copy(*this);
321  d11_ += 2.0*( copy.s12()*rot.x() + copy.s13()*rot.y());
322  d22_ += 2.0*(-copy.s12()*rot.x() + copy.s23()*rot.z());
323  d33_ += -2.0*( copy.s13()*rot.y() + copy.s23()*rot.z());
324  d12_ += ((copy.s22() - copy.s11())*rot.x() + copy.s23()*rot.x() + copy.s13()*rot.z());
325  d13_ += ( copy.s23()*rot.x() + (copy.s33() - copy.s11())*rot.x() - copy.s12()*rot.z());
326  d23_ += ( -copy.s13()*rot.x() - copy.s12()*rot.x() + (copy.s33() - copy.s22())*rot.z());
327 }
328 
329 // 11 0
330 // 22 1
331 // 33 2
332 // 12 3
333 // 13 4
334 // 23 5
335 uint32 SymTensor::doubleToSingleComponent(uint32 dof1,uint32 dof2) {
336  dof1 = std::clamp<uint32>(dof1,0,2);
337  dof2 = std::clamp<uint32>(dof2,0,2);
338  switch (dof1) {
339  case 0:
340  switch (dof2) {
341  case 0: return 0;
342  case 1: return 3;
343  case 2: return 4;
344  }
345  break;
346  case 1:
347  switch (dof2) {
348  case 0: return 3;
349  case 1: return 1;
350  case 2: return 5;
351  }
352  break;
353  case 2:
354  switch (dof2) {
355  case 0: return 4;
356  case 1: return 5;
357  case 2: return 3;
358  }
359  break;
360  }
361  return 0;
362 }
363 
364 namespace base {
365  BASE_EXPORT string ts(const SymTensor &s, int width=0, char notation = '\0', int precision = -1, char fill = ' ');
366 }
367 
369 // EoF
axes.h
2D and 3D cartesian Axes systems.
Axes3D
Class for specifying a particular 3D cartesian axes system, and converting to and from it.
Definition: axes.h:121
SymTensorInfo
SymTensor eigenvalue and direction helper class.
Definition: symtensor.h:166
SymTensorInfo::SymTensorInfo
BASE_EXPORT SymTensorInfo()
Default constructor.
Definition: symtensor.h:170
SymTensorInfo::getAxes
BASE_EXPORT Axes3D getAxes() const
Returns eigen directions (minimum, intermediate, maximum).
Definition: symtensor.cpp:567
SymTensorInfo::resolve
BASE_EXPORT SymTensor resolve(const DVect3 &prin) const
Regenerates full tensor from info + principal directions.
Definition: symtensor.cpp:574
SymTensorInfo::operator=
BASE_EXPORT const SymTensorInfo & operator=(const SymTensorInfo &si)
Equality operator.
Definition: symtensor.cpp:561
Matrix::operator*
constexpr Vector2< T > operator*(const Matrix< T, 2, 2 > &m, const Vector2< T > &v)
Definition: matrix.h:987
BASE_EXPORT
#define BASE_EXPORT
Definition: basedef.h:24
VSWARNING
PUSHWARNING VSWARNING(26495) VSWARNING(4702) class BASE_EXPORT SymTensor
A symmetric 2nd order tensor.
Definition: symtensor.h:20
Generated by doxygen 1.9.1