MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Matrix.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 #include "Matrix.h"
27 #include "ExtendedMath.h"
28 
29 using mathsFunc::square;
30 
35 {
36  setZero();
37 }
38 
45 Matrix3D::Matrix3D(const Mdouble xx, const Mdouble xy, const Mdouble xz, const Mdouble yx, const Mdouble yy,
46  const Mdouble yz, const Mdouble zx, const Mdouble zy, const Mdouble zz)
47 {
48  XX = xx;
49  XY = xy;
50  XZ = xz;
51  YX = yx;
52  YY = yy;
53  YZ = yz;
54  ZX = zx;
55  ZY = zy;
56  ZZ = zz;
57 }
58 
60 {
61  XX = matrix(0, 0);
62  XY = matrix(0, 1);
63  XZ = matrix(0, 2);
64  YX = matrix(1, 0);
65  YY = matrix(1, 1);
66  YZ = matrix(1, 2);
67  ZX = matrix(2, 0);
68  ZY = matrix(2, 1);
69  ZZ = matrix(2, 2);
70 }
71 
76 {
77  XX = XY = XZ = YX = YY = YZ = ZX = ZY = ZZ = 0.0;
78 }
79 
85 {
86  return XX + YY + ZZ;
87 }
88 
94 {
95  return Vec3D(XX, YY, ZZ);
96 }
97 
98 
104 {
105  const Mdouble P = trace() / 3.0;
106  return std::sqrt(0.5 * (mathsFunc::square(XX - P) + mathsFunc::square(YY - P) + mathsFunc::square(ZZ - P)) +
108 }
109 
116 {
117  return Matrix3D(XX + A.XX, XY + A.XY, XZ + A.XZ,
118  YX + A.YX, YY + A.YY, YZ + A.YZ,
119  ZX + A.ZX, ZY + A.ZY, ZZ + A.ZZ);
120 }
121 
129 {
130  return Matrix3D(XX - A.XX, XY - A.XY, XZ - A.XZ,
131  YX - A.YX, YY - A.YY, YZ - A.YZ,
132  ZX - A.ZX, ZY - A.ZY, ZZ - A.ZZ);
133 }
134 
141 {
142  return Matrix3D(XX + a, XY + a, XZ + a,
143  YX + a, YY + a, YZ + a,
144  ZX + a, ZY + a, ZZ + a);
145 }
146 
153 {
154  return Matrix3D(XX - a, XY - a, XZ - a,
155  YX - a, YY - a, YZ - a,
156  ZX - a, ZY - a, ZZ - a);
157 }
158 
165 {
166  return Matrix3D(XX * a, XY * a, XZ * a,
167  YX * a, YY * a, YZ * a,
168  ZX * a, ZY * a, ZZ * a);
169 }
170 
177 {
178  return Vec3D(XX * a.X + XY * a.Y + XZ * a.Z,
179  YX * a.X + YY * a.Y + YZ * a.Z,
180  ZX * a.X + ZY * a.Y + ZZ * a.Z);
181 }
182 
185 {
186  return Matrix3D(XX * a.XX + XY * a.YX + XZ * a.ZX,
187  XX * a.XY + XY * a.YY + XZ * a.ZY,
188  XX * a.XZ + XY * a.YZ + XZ * a.ZZ,
189  YX * a.XX + YY * a.YX + YZ * a.ZX,
190  YX * a.XY + YY * a.YY + YZ * a.ZY,
191  YX * a.XZ + YY * a.YZ + YZ * a.ZZ,
192  ZX * a.XX + ZY * a.YX + ZZ * a.ZX,
193  ZX * a.XY + ZY * a.YY + ZZ * a.ZY,
194  ZX * a.XZ + ZY * a.YZ + ZZ * a.ZZ
195  );
196 }
197 
204 {
205  return Matrix3D(XX / a, XY / a, XZ / a,
206  YX / a, YY / a, YZ / a,
207  ZX / a, ZY / a, ZZ / a);
208 }
209 
216 std::ostream& operator<<(std::ostream& os, const Matrix3D& A)
217 {
218  os << A.XX << ' ' << A.XY << ' ' << A.XZ << ' '
219  << A.YX << ' ' << A.YY << ' ' << A.YZ << ' '
220  << A.ZX << ' ' << A.ZY << ' ' << A.ZZ;
221  return os;
222 }
223 
230 std::istream& operator>>(std::istream& is, Matrix3D& A)
231 {
232  is >> A.XX >> A.XY >> A.XZ >> A.YX >> A.YY >> A.YZ >> A.ZX >> A.ZY >> A.ZZ;
233  return is;
234 }
235 
242 {
243  XX += A.XX;
244  XY += A.XY;
245  XZ += A.XZ;
246  YX += A.YX;
247  YY += A.YY;
248  YZ += A.YZ;
249  ZX += A.ZX;
250  ZY += A.ZY;
251  ZZ += A.ZZ;
252  return *this;
253 }
254 
261 {
262  XX -= A.XX;
263  XY -= A.XY;
264  XZ -= A.XZ;
265  YX -= A.YX;
266  YY -= A.YY;
267  YZ -= A.YZ;
268  ZX -= A.ZX;
269  ZY -= A.ZY;
270  ZZ -= A.ZZ;
271  return *this;
272 }
273 
280 {
281  XX /= a;
282  XY /= a;
283  XZ /= a;
284  YX /= a;
285  YY /= a;
286  YZ /= a;
287  ZX /= a;
288  ZY /= a;
289  ZZ /= a;
290  return *this;
291 }
292 
299 {
303 }
304 
311 {
312  return Matrix3D(std::sqrt(A.XX), std::sqrt(A.XY), std::sqrt(A.XZ),
313  std::sqrt(A.YX), std::sqrt(A.YY), std::sqrt(A.YZ),
314  std::sqrt(A.ZX), std::sqrt(A.ZY), std::sqrt(A.ZZ));
315 }
316 
324 {
325  return Matrix3D(a.X * b.X, a.X * b.Y, a.X * b.Z,
326  a.Y * b.X, a.Y * b.Y, a.Y * b.Z,
327  a.Z * b.X, a.Z * b.Y, a.Z * b.Z);
328 }
329 
338 {
339  return Matrix3D(
340  a.Y * B.ZX - a.Z * B.YX, a.Y * B.ZY - a.Z * B.YY, a.Y * B.ZZ - a.Z * B.YZ,
341  a.Z * B.XX - a.X * B.ZX, a.Z * B.XY - a.X * B.ZY, a.Z * B.XZ - a.X * B.ZZ,
342  a.X * B.YX - a.Y * B.XX, a.X * B.YY - a.Y * B.XY, a.X * B.YZ - a.Y * B.XZ);
343 }
344 
350 {
351  Matrix3D result;
352  Mdouble det;
353  det = 1. / (A.XX * (A.YY * A.ZZ - A.YZ * A.ZY)
354  + A.XY * (A.YZ * A.ZX - A.YX * A.ZZ)
355  + A.XZ * (A.YX * A.ZY - A.YY * A.ZX));
356  result.XX = (A.YY * A.ZZ - A.YZ * A.ZY) * det;
357  result.XY = -(A.XY * A.ZZ - A.XZ * A.ZY) * det;
358  result.XZ = (A.XY * A.YZ - A.XZ * A.YY) * det;
359  result.YX = -(A.YX * A.ZZ - A.YZ * A.ZX) * det;
360  result.YY = (A.XX * A.ZZ - A.XZ * A.ZX) * det;
361  result.YZ = -(A.XX * A.YZ - A.XZ * A.YX) * det;
362  result.ZX = (A.YX * A.ZY - A.YY * A.ZX) * det;
363  result.ZY = -(A.XX * A.ZY - A.XY * A.ZX) * det;
364  result.ZZ = (A.XX * A.YY - A.XY * A.YX) * det;
365  return result;
366 }
367 
374 {
375  Mdouble invdet = 1. / (XX * (YY * ZZ - YZ * ZY)
376  + XY * (YZ * ZX - YX * ZZ)
377  + XZ * (YX * ZY - YY * ZX));
378  return Vec3D(+(XY * YZ - YY * XZ) * b.Z
379  - (ZY * YZ - ZZ * YY) * b.X
380  + (ZY * XZ - ZZ * XY) * b.Y,
381  -(XX * YZ - XZ * YX) * b.Z
382  + (ZX * YZ - ZZ * YX) * b.X
383  - (ZX * XZ - XX * ZZ) * b.Y,
384  +(XX * YY - YX * XY) * b.Z
385  - (ZX * YY - YX * ZY) * b.X
386  + (ZX * XY - XX * ZY) * b.Y) * invdet;
387 }
388 
395 {
396  //define sin(A)=y/r, cos(A)=x/r
397  const Mdouble r = std::sqrt(p.X * p.X + p.Y * p.Y);
398  Mdouble s = p.Y / r;
399  Mdouble c = p.X / r;
400  if (r == 0)
401  {
402  s = 0;
403  c = 1;
404  }
405  const Matrix3D M = Matrix3D(XX * c + XY * s, -XX * s + XY * c, XZ,
406  YX * c + YY * s, -YX * s + YY * c, YZ,
407  ZX * c + ZY * s, -ZX * s + ZY * c, ZZ);
408  return Matrix3D(M.XX * c + M.YX * s, -M.XX * s + M.YX * c, M.ZX,
409  M.XY * c + M.YY * s, -M.XY * s + M.YY * c, M.ZY,
410  M.XZ * c + M.YZ * s, -M.XZ * s + M.YZ * c, M.ZZ);
411 }
Mdouble XY
Definition: Matrix.h:43
Matrix3D operator-(const Matrix3D &A) const
Matrix subtraction.
Definition: Matrix.cc:128
double trace() const
Sum of the diagonal elements.
Definition: Matrix.cc:84
Matrix3D & operator+=(const Matrix3D &A)
Matrix addition.
Definition: Matrix.cc:241
Mdouble X
the vector components
Definition: Vector.h:65
Mdouble ZY
Definition: Matrix.h:43
double Mdouble
Definition: GeneralDefine.h:34
Mdouble XZ
Definition: Matrix.h:43
double deviator() const
Deviator.
Definition: Matrix.cc:103
Mdouble ZX
Definition: Matrix.h:43
static Matrix3D inverse(const Matrix3D &A)
Computes the inverse of a matrix.
Definition: Matrix.cc:349
Mdouble ZZ
Definition: Matrix.h:43
Matrix3D operator/(Mdouble a) const
Scalar division.
Definition: Matrix.cc:203
Mdouble YZ
Definition: Matrix.h:43
static Matrix3D dyadic(const Vec3D &a, const Vec3D &b)
Calculates the dyadic product of a two Vec3D: .
Definition: Matrix.cc:323
Matrix3D operator*(Mdouble a) const
Scalar multiplication.
Definition: Matrix.cc:164
static Matrix3D cross(const Vec3D &a, const Matrix3D &b)
'Special' cross product; CP of vector with each column of a matrix
Definition: Matrix.cc:337
Mdouble YX
Definition: Matrix.h:43
Matrix3D operator+(const Matrix3D &A) const
Matrix addition.
Definition: Matrix.cc:115
Matrix3D & operator/=(Mdouble a)
Scalar division.
Definition: Matrix.cc:279
Mdouble Y
Definition: Vector.h:65
Matrix3D & operator-=(const Matrix3D &A)
Matrix substraction.
Definition: Matrix.cc:260
Matrix3D getCylindricalTensorField(const Vec3D &p) const
Returns the matrix in cylindrical coordinates.
Definition: Matrix.cc:394
Vec3D ldivide(const Vec3D &b)
A.ldivide(b) computes the solution x to A*x=b.
Definition: Matrix.cc:373
static Matrix3D sqrt(const Matrix3D &A)
Calculates the pointwise square root.
Definition: Matrix.cc:310
std::ostream & operator<<(std::ostream &os, const Matrix3D &A)
Definition: Matrix.cc:216
Mdouble XX
all nine matrix elements
Definition: Matrix.h:43
void setZero()
Sets all elements to zero.
Definition: Matrix.cc:75
Mdouble YY
Definition: Matrix.h:43
Implementation of a 3D matrix.
Definition: Matrix.h:37
Definition: Vector.h:49
std::istream & operator>>(std::istream &is, Matrix3D &A)
Definition: Matrix.cc:230
Data type for small dense matrix.
Definition: SmallMatrix.h:67
T square(const T val)
squares a number
Definition: ExtendedMath.h:104
Mdouble Z
Definition: Vector.h:65
Vec3D diag() const
The diagonal elements.
Definition: Matrix.cc:93
Matrix3D()
default constructor
Definition: Matrix.cc:34
static Matrix3D square(const Matrix3D &A)
Calculates the pointwise square.
Definition: Matrix.cc:298