Itasca C++ Interface
Loading...
Searching...
No Matches
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
13class SymTensorInfo;
14
19PUSHWARNING
20VSWARNING(26495) // member init (not initialized on purpose)
21VSWARNING(4702) // Weird unreachable code warnings in release
23public:
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 bool isIsotropic(double tol = std::numeric_limits<double>::epsilon()*1000.0) const;
146 inline void adjustTrace(const double &newTrace);
147 inline void incrementDiagonal(const double &increment) { d11_ += increment; d22_ += increment; d33_ += increment; }
148 inline void rotate(const DVect3 &rot);
149 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_)))))); }
150 inline bool zero() const { return (d11_==0.0) and (d22_==0.0) and (d33_==0.0) and (d12_==0.0) and (d13_==0.0) and (d23_==0.0); }
151private:
152 double d11_;
153 double d22_;
154 double d33_;
155 double d12_;
156 double d13_;
157 double d23_;
158};
159POPWARNING
160
169public:
170 friend class SymTensor;
172 BASE_EXPORT SymTensorInfo(): type_(Type::ThreeDCube) { }
178 BASE_EXPORT Axes3D getAxes() const;
180 BASE_EXPORT SymTensor resolve(const DVect3 &prin) const;
181private:
182 enum class Type { ThreeDCube, ThreeDJacobi, ZMax, ZMid, ZMin };
183 Type type_;
184 Axes3D axes_;
185};
186
187inline const SymTensor &SymTensor::operator=(const SymTensor &s) {
188 d11_ = s.d11_;
189 d22_ = s.d22_;
190 d33_ = s.d33_;
191 d12_ = s.d12_;
192 d13_ = s.d13_;
193 d23_ = s.d23_;
194 return *this;
195}
196
197inline bool SymTensor::operator==(const SymTensor &s) const {
198 return (d11_ == s.d11_ && d22_ == s.d22_ &&
199 d33_ == s.d33_ && d12_ == s.d12_ && d13_ == s.d13_ &&
200 d23_ == s.d23_);
201}
202
203inline const double &SymTensor::operator[](unsigned int i) const {
204 assert(i<6);
205 switch (i) {
206 case 0: return d11_;
207 case 1: return d22_;
208 case 2: return d33_;
209 case 3: return d12_;
210 case 4: return d13_;
211 case 5: return d23_;
212 }
213 return d11_;
214}
215
216inline double &SymTensor::operator[](unsigned int i) {
217 assert(i<6);
218 switch (i) {
219 case 0: return d11_;
220 case 1: return d22_;
221 case 2: return d33_;
222 case 3: return d12_;
223 case 4: return d13_;
224 case 5: return d23_;
225 }
226 return d11_;
227}
228
229inline double SymTensor::getTotalMeasure() const {
230 double I1 = getTrace();
231 double J2 = getJ2();
232 return sqrt(I1*I1/3.0 + 2.0*J2);
233}
234
235inline DVect3 SymTensor::operator*(const DVect3 &normal) const {
236 return DVect3(normal.x()*s11() + normal.y()*s12() + normal.z()*s13(),
237 normal.x()*s21() + normal.y()*s22() + normal.z()*s23(),
238 normal.x()*s31() + normal.y()*s32() + normal.z()*s33());
239}
240
241inline DVect2 SymTensor::operator*(const DVect2 &normal) const {
242 return DVect2(normal.x()*s11() + normal.y()*s12(),
243 normal.x()*s21() + normal.y()*s22());
244}
245
246inline SymTensor SymTensor::operator*(const double &mul) const {
247 SymTensor ret;
248 ret.d11_ = d11_ * mul;
249 ret.d22_ = d22_ * mul;
250 ret.d33_ = d33_ * mul;
251 ret.d12_ = d12_ * mul;
252 ret.d13_ = d13_ * mul;
253 ret.d23_ = d23_ * mul;
254 return ret;
255}
256
257inline const SymTensor &SymTensor::operator*=(const double &mul) {
258 d11_ *= mul;
259 d22_ *= mul;
260 d33_ *= mul;
261 d12_ *= mul;
262 d13_ *= mul;
263 d23_ *= mul;
264 return *this;
265}
266
267inline SymTensor SymTensor::operator/(const double &mul) const {
268 SymTensor ret;
269 ret.d11_ = d11_ / mul;
270 ret.d22_ = d22_ / mul;
271 ret.d33_ = d33_ / mul;
272 ret.d12_ = d12_ / mul;
273 ret.d13_ = d13_ / mul;
274 ret.d23_ = d23_ / mul;
275 return ret;
276}
277
278inline const SymTensor &SymTensor::operator/=(const double &mul) {
279 d11_ /= mul;
280 d22_ /= mul;
281 d33_ /= mul;
282 d12_ /= mul;
283 d13_ /= mul;
284 d23_ /= mul;
285 return *this;
286}
287
288inline SymTensor SymTensor::operator+(const SymTensor &s) const {
289 SymTensor ret;
290 ret.d11_ = d11_ + s.d11_;
291 ret.d22_ = d22_ + s.d22_;
292 ret.d33_ = d33_ + s.d33_;
293 ret.d12_ = d12_ + s.d12_;
294 ret.d13_ = d13_ + s.d13_;
295 ret.d23_ = d23_ + s.d23_;
296 return ret;
297}
298
299inline SymTensor SymTensor::operator-(const SymTensor &s) const {
300 SymTensor ret;
301 ret.d11_ = d11_ - s.d11_;
302 ret.d22_ = d22_ - s.d22_;
303 ret.d33_ = d33_ - s.d33_;
304 ret.d12_ = d12_ - s.d12_;
305 ret.d13_ = d13_ - s.d13_;
306 ret.d23_ = d23_ - s.d23_;
307 return ret;
308}
309
310inline bool SymTensor::isIsotropic(double tol) const {
311 double dtol = std::max(std::max(std::abs(s11()),std::abs(s22())),std::abs(s33())) * tol;
312 if (std::abs(s11()-s22()) > dtol) return false;
313 if (std::abs(s11()-s33()) > dtol) return false;
314 if (std::abs(s12()) > dtol) return false;
315 if (std::abs(s13()) > dtol) return false;
316 if (std::abs(s23()) > dtol) return false;
317 return true;
318}
319
320inline void SymTensor::adjustTrace(const double &newTrace) {
321 static constexpr double d1d3 = 1.0 / 3.0;
322 /* --- strain differences --- */
323 double dx = s11() - s22();
324 double dy = s22() - s33();
325 double dz = s33() - s11();
326 rs11() = (newTrace + dx - dz) * d1d3;
327 rs22() = (newTrace + dy - dx) * d1d3;
328 rs33() = (newTrace + dz - dy) * d1d3;
329}
330
331inline void SymTensor::rotate(const DVect3 &rot) {
332 SymTensor copy(*this);
333 d11_ += 2.0*( copy.s12()*rot.x() + copy.s13()*rot.y());
334 d22_ += 2.0*(-copy.s12()*rot.x() + copy.s23()*rot.z());
335 d33_ += -2.0*( copy.s13()*rot.y() + copy.s23()*rot.z());
336 d12_ += ((copy.s22() - copy.s11())*rot.x() + copy.s23()*rot.x() + copy.s13()*rot.z());
337 d13_ += ( copy.s23()*rot.x() + (copy.s33() - copy.s11())*rot.x() - copy.s12()*rot.z());
338 d23_ += ( -copy.s13()*rot.x() - copy.s12()*rot.x() + (copy.s33() - copy.s22())*rot.z());
339}
340
341// 11 0
342// 22 1
343// 33 2
344// 12 3
345// 13 4
346// 23 5
347uint32 SymTensor::doubleToSingleComponent(uint32 dof1,uint32 dof2) {
348 dof1 = std::clamp<uint32>(dof1,0,2);
349 dof2 = std::clamp<uint32>(dof2,0,2);
350 switch (dof1) {
351 case 0:
352 switch (dof2) {
353 case 0: return 0;
354 case 1: return 3;
355 case 2: return 4;
356 }
357 break;
358 case 1:
359 switch (dof2) {
360 case 0: return 3;
361 case 1: return 1;
362 case 2: return 5;
363 }
364 break;
365 case 2:
366 switch (dof2) {
367 case 0: return 4;
368 case 1: return 5;
369 case 2: return 3;
370 }
371 break;
372 }
373 return 0;
374}
375
376namespace base {
377 BASE_EXPORT string ts(const SymTensor &s, int width=0, char notation = '\0', int precision = -1, char fill = ' ');
378}
379
381// EoF
2D and 3D cartesian Axes systems.
Class for specifying a particular 3D cartesian axes system, and converting to and from it.
Definition axes.h:121
A symmetric 2nd order tensor.
Definition symtensor.h:22
double & rs31()
Reference access to components, note that s12 is equivalent to s21 (for instance).
Definition symtensor.h:62
double & rs22()
Reference access to components, note that s12 is equivalent to s21 (for instance).
Definition symtensor.h:56
double & rs13()
Reference access to components, note that s12 is equivalent to s21 (for instance).
Definition symtensor.h:59
double & rs11()
Reference access to components, note that s12 is equivalent to s21 (for instance).
Definition symtensor.h:55
double getI3() const
Returns the third invariant, or I3.
Definition symtensor.h:91
SymTensor(const double &i11, const double &i22=0.0, const double &i33=0.0, const double &i12=0.0, const double &i13=0.0, const double &i23=0.0)
Definition symtensor.h:34
const SymTensor & operator+=(const SymTensor &s)
+= operator for a SymTensor
Definition symtensor.h:132
double & rs12()
Reference access to components, note that s12 is equivalent to s21 (for instance).
Definition symtensor.h:58
const double & s23() const
Component access, note that s12 is equivalent to s21 (for instance).
Definition symtensor.h:50
const double & s21() const
Component access, note that s12 is equivalent to s21 (for instance).
Definition symtensor.h:51
const double & s22() const
Component access, note that s12 is equivalent to s21 (for instance).
Definition symtensor.h:46
double getI1() const
Same as getTrace() - returns the first invariant.
Definition symtensor.h:87
const double & s33() const
Component access, note that s12 is equivalent to s21 (for instance).
Definition symtensor.h:47
SymTensor mul(const double &d) const
Returns a SymTensor with every component multiplied by a scalar value.
Definition symtensor.h:130
const double & s13() const
Component access, note that s12 is equivalent to s21 (for instance).
Definition symtensor.h:49
double & rs33()
Reference access to components, note that s12 is equivalent to s21 (for instance).
Definition symtensor.h:57
bool isDiagonal(double tol=std::numeric_limits< double >::epsilon() *1000.0) const
Determines whether or not the SymTensor is diagonal.
Definition symtensor.h:144
const double & s11() const
Component access, note that s12 is equivalent to s21 (for instance).
Definition symtensor.h:45
double getI2() const
Returns the second invariant.
Definition symtensor.h:89
double & rs21()
Reference access to components, note that s12 is equivalent to s21 (for instance).
Definition symtensor.h:61
double getTrace() const
Returns the trace of the tensor (11+22+33). I1.
Definition symtensor.h:83
const double & s31() const
Component access, note that s12 is equivalent to s21 (for instance).
Definition symtensor.h:52
double & rs23()
Reference access to components, note that s12 is equivalent to s21 (for instance).
Definition symtensor.h:60
SymTensor(const SymTensor &s)
Copy constructor.
Definition symtensor.h:31
double getJ2() const
Returns the second invariant of the deviatoric – J2.
Definition symtensor.cpp:156
SymTensor()
Default constructor, no data initialization.
Definition symtensor.h:28
double & rs32()
Reference access to components, note that s12 is equivalent to s21 (for instance).
Definition symtensor.h:63
double getNorm2() const
Returns a scalar norm (magnitude) value for the tensor, can be used for tolerance checking,...
Definition symtensor.h:112
SymTensor getDeviatoric() const
Returns the deviatoric tensor.
Definition symtensor.h:85
const double & s12() const
Component access, note that s12 is equivalent to s21 (for instance).
Definition symtensor.h:48
const double & s32() const
Component access, note that s12 is equivalent to s21 (for instance).
Definition symtensor.h:53
SymTensor eigenvalue and direction helper class.
Definition symtensor.h:168
BASE_EXPORT SymTensorInfo()
Default constructor.
Definition symtensor.h:172
BASE_EXPORT const SymTensorInfo & operator=(const SymTensorInfo &si)
Equality operator.
Definition symtensor.cpp:559
const SymTensor & operator=(const SymTensor &s)
Assignment operator.
Definition symtensor.h:187
double getTotalMeasure() const
Definition symtensor.h:229
bool operator==(const SymTensor &s) const
Equality operator.
Definition symtensor.h:197
#define BASE_EXPORT
Definition basedef.h:24
const double & operator[](unsigned int i) const
Allows Index access to tensor components.
Definition symtensor.h:203
DVect3 operator*(const DVect3 &input) const
Performs the linear mapping represented by the tensor on the vector input.
Definition symtensor.h:235
constexpr Vector2< T > max(const Vector2< T > &v1, const Vector2< T > &v2)
Template specialization for max, min.
Definition vect.h:433