MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MatrixSymmetric.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 "MatrixSymmetric.h"
27 #include "ExtendedMath.h"
28 
34 MatrixSymmetric3D::operator Matrix3D() const
35 {
36  return {XX, XY, XZ, XY, YY, YZ, XZ, YZ, ZZ};
37 }
38 
43 {
44  setZero();
45 }
46 
56 MatrixSymmetric3D::MatrixSymmetric3D(const Mdouble xx, const Mdouble xy, const Mdouble xz, const Mdouble yy,
57  const Mdouble yz, const Mdouble zz)
58 {
59  XX = xx;
60  XY = xy;
61  XZ = xz;
62  YY = yy;
63  YZ = yz;
64  ZZ = zz;
65 }
66 
71 {
72  XX = XY = XZ = YY = YZ = ZZ = 0.0;
73 }
74 
81 {
82  return (XX + YY + ZZ) / 3;
83 }
84 
91 {
92  return MatrixSymmetric3D(XX + A.XX, XY + A.XY, XZ + A.XZ, YY + A.YY, YZ + A.YZ, ZZ + A.ZZ);
93 }
94 
101 {
102  return MatrixSymmetric3D(XX - A.XX, XY - A.XY, XZ - A.XZ, YY - A.YY, YZ - A.YZ, ZZ - A.ZZ);
103 }
104 
111 {
112  return MatrixSymmetric3D(XX + a, XY + a, XZ + a, YY + a, YZ + a, ZZ + a);
113 }
114 
121 {
122  return MatrixSymmetric3D(XX - a, XY - a, XZ - a, YY - a, YZ - a, ZZ - a);
123 }
124 
133 {
134  return Vec3D(A.XX * b.X + A.XY * b.Y + A.XZ * b.Z,
135  A.XY * b.X + A.YY * b.Y + A.YZ * b.Z,
136  A.XZ * b.X + A.YZ * b.Y + A.ZZ * b.Z);
137 }
138 
145 {
146  return MatrixSymmetric3D(XX * a, XY * a, XZ * a,
147  YY * a, YZ * a, ZZ * a);
148 }
149 
156 {
157  return MatrixSymmetric3D(XX / a, XY / a, XZ / a, YY / a, YZ / a, ZZ / a);
158 }
159 
166 std::ostream& operator<<(std::ostream& os, const MatrixSymmetric3D& A)
167 {
168  os << A.XX << ' ' << A.XY << ' ' << A.XZ << " " << A.YY << ' ' << A.YZ << " " << A.ZZ;
169  return os;
170 }
171 
178 std::istream& operator>>(std::istream& is, MatrixSymmetric3D& A)
179 {
180  is >> A.XX >> A.XY >> A.XZ >> A.YY >> A.YZ >> A.ZZ;
181  return is;
182 }
183 
190 {
191  XX += A.XX;
192  XY += A.XY;
193  XZ += A.XZ;
194  YY += A.YY;
195  YZ += A.YZ;
196  ZZ += A.ZZ;
197  return *this;
198 }
199 
206 {
207  XX -= A.XX;
208  XY -= A.XY;
209  XZ -= A.XZ;
210  YY -= A.YY;
211  YZ -= A.YZ;
212  ZZ -= A.ZZ;
213  return *this;
214 }
215 
222 {
223  XX /= a;
224  XY /= a;
225  XZ /= a;
226  YY /= a;
227  YZ /= a;
228  ZZ /= a;
229  return *this;
230 }
231 
239 {
242 }
243 
251 {
252  return MatrixSymmetric3D(std::sqrt(A.XX), std::sqrt(A.XY), std::sqrt(A.XZ), std::sqrt(A.YY), std::sqrt(A.YZ),
253  std::sqrt(A.ZZ));
254 }
255 
264 {
265  return MatrixSymmetric3D(a.X * a.X, a.X * a.Y, a.X * a.Z, a.Y * a.Y, a.Y * a.Z, a.Z * a.Z);
266 }
267 
278 {
279  return MatrixSymmetric3D(a.X * b.X, 0.5 * (a.X * b.Y + b.X * a.Y), 0.5 * (a.X * b.Z + b.X * a.Z), a.Y * b.Y,
280  0.5 * (a.Y * b.Z + b.Y * a.Z), a.Z * b.Z);
281 }
282 
289 {
290  return A.inverse();
291 }
292 
299 {
300  MatrixSymmetric3D result;
302  logger.assert(det != 0,
303  "determinant is not zero"); //should be replaced by the condition number or sth. like fabs(det)>1e-5*norm.
304  result.XX = (YY * ZZ - YZ * YZ) * det;
305  result.XY = -(XY * ZZ - XZ * YZ) * det;
306  result.XZ = (XY * YZ - XZ * YY) * det;
307  result.YY = (XX * ZZ - XZ * XZ) * det;
308  result.YZ = -(XX * YZ - XZ * XY) * det;
309  result.ZZ = (XX * YY - XY * XY) * det;
310  return result;
311 }
312 
314 {
315  return 1. / (A.XX * (A.YY * A.ZZ - A.YZ * A.YZ)
316  + A.XY * (A.YZ * A.XZ - A.XY * A.ZZ)
317  + A.XZ * (A.XY * A.YZ - A.YY * A.XZ));
318 }
319 
327 {
328 // //define sin(A)=y/r, cos(A)=x/r
329 // Mdouble r = std::sqrt(p.X*p.X+p.Y*p.Y);
330 // Mdouble s = p.Y/r;
331 // Mdouble c = p.X/r;
332 // if (r==0) {
333 // s=0;
334 // c=1;
335 // }
336  return *this;
337 }
static MatrixSymmetric3D square(const MatrixSymmetric3D &A)
Calculates the pointwise square.
Mdouble X
the vector components
Definition: Vector.h:65
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
MatrixSymmetric3D & operator-=(const MatrixSymmetric3D &A)
Matrix substraction.
std::istream & operator>>(std::istream &is, MatrixSymmetric3D &A)
std::ostream & operator<<(std::ostream &os, const MatrixSymmetric3D &A)
MatrixSymmetric3D & operator/=(Mdouble a)
Scalar division.
MatrixSymmetric3D operator+(const MatrixSymmetric3D &A) const
Matrix addition.
Vec3D operator*(const MatrixSymmetric3D &A, const Vec3D &b)
MatrixSymmetric3D inverse() const
Computes the inverse of a matrix; exits if the inverse doesn't exist.
void setZero()
Sets all elements to zero.
MatrixSymmetric3D operator-(const MatrixSymmetric3D &A) const
Matrix substraction.
MatrixSymmetric3D & operator+=(const MatrixSymmetric3D &A)
Matrix addition.
Mdouble Y
Definition: Vector.h:65
static Mdouble determinant(const MatrixSymmetric3D &A)
Computes the determinant of a matrix.
static MatrixSymmetric3D symmetrisedDyadic(const Vec3D &a, const Vec3D &b)
Calculates the symmetrised dyadic product of two Vec3D: .
static MatrixSymmetric3D selfDyadic(const Vec3D &a)
Calculates the dyadic product of a Vec3D with itself: .
MatrixSymmetric3D operator/(Mdouble a) const
Scalar division.
MatrixSymmetric3D()
Default constructor.
friend Vec3D operator*(const MatrixSymmetric3D &A, const Vec3D &b)
Vector multiplication.
static MatrixSymmetric3D inverse(const MatrixSymmetric3D &A)
Computes the inverse of a matrix; exits if the inverse doesn't exist.
MatrixSymmetric3D getCylindricalTensorField(const Vec3D &p) const
Returns the matrix in cylindrical coordinates.
Implementation of a 3D matrix.
Definition: Matrix.h:37
Definition: Vector.h:49
T square(const T val)
squares a number
Definition: ExtendedMath.h:104
Mdouble trace() const
Returns the MEAN of the diagonal elements (i.e. the trace divided by three).
static MatrixSymmetric3D sqrt(const MatrixSymmetric3D &A)
Calculates the pointwise square root.
Mdouble Z
Definition: Vector.h:65
Implementation of a 3D symmetric matrix.
Mdouble XX
The six distinctive matrix elements.