|
moab
|
Functions | |
| template<typename Matrix > | |
| Matrix | inverse (const Matrix &d, const double i) |
| template<typename Matrix > | |
| bool | positive_definite (const Matrix &d, double &det) |
| template<typename Matrix > | |
| Matrix | transpose (const Matrix &d) |
| template<typename Matrix > | |
| Matrix | mmult3 (const Matrix &a, const Matrix &b) |
| template<typename Vector , typename Matrix > | |
| Matrix | outer_product (const Vector &u, const Vector &v, Matrix &m) |
| template<typename Matrix > | |
| double | determinant3 (const Matrix &d) |
| template<typename Matrix > | |
| const Matrix | inverse (const Matrix &d) |
| template<typename Vector , typename Matrix > | |
| Vector | vector_matrix (const Vector &v, const Matrix &m) |
| template<typename Vector , typename Matrix > | |
| Vector | matrix_vector (const Matrix &m, const Vector &v) |
| template<typename Matrix , typename Vector > | |
| ErrorCode | EigenDecomp (const Matrix &_a, double w[3], Vector v[3]) |
| double moab::Matrix::determinant3 | ( | const Matrix & | d | ) | [inline] |
Definition at line 95 of file Matrix3.hpp.
{
return (d(0) * d(4) * d(8)
+ d(1) * d(5) * d(6)
+ d(2) * d(3) * d(7)
- d(0) * d(5) * d(7)
- d(1) * d(3) * d(8)
- d(2) * d(4) * d(6));
}
| ErrorCode moab::Matrix::EigenDecomp | ( | const Matrix & | _a, |
| double | w[3], | ||
| Vector | v[3] | ||
| ) |
Definition at line 143 of file Matrix3.hpp.
{
const int MAX_ROTATIONS = 20;
const double one_ninth = 1./9;
int i, j, k, iq, ip, numPos;
double tresh, theta, tau, t, sm, s, h, g, c, tmp;
double b[3], z[3];
Matrix a( _a);
// initialize
for (ip=0; ip<3; ip++) {
for (iq=0; iq<3; iq++){
v[ip][iq] = 0.0;
}
v[ip][ip] = 1.0;
}
for (ip=0; ip<3; ip++) {
b[ip] = w[ip] = a[ip][ip];
z[ip] = 0.0;
}
// begin rotation sequence
for (i=0; i<MAX_ROTATIONS; i++){
sm = 0.0;
for (ip=0; ip<2; ip++){
for (iq=ip+1; iq<3; iq++){ sm += fabs(a[ip][iq]); }
}
if ( sm == 0.0 ){ break; }
// first 3 sweeps
tresh = (i < 3)? 0.2*sm*one_ninth : 0.0;
for (ip=0; ip<2; ip++) {
for (iq=ip+1; iq<3; iq++) {
g = 100.0*fabs(a[ip][iq]);
// after 4 sweeps
if ( i > 3 && (fabs(w[ip])+g) == fabs(w[ip])
&& (fabs(w[iq])+g) == fabs(w[iq])) {
a[ip][iq] = 0.0;
}
else if ( fabs(a[ip][iq]) > tresh) {
h = w[iq] - w[ip];
if ( (fabs(h)+g) == fabs(h)){ t = (a[ip][iq]) / h; }
else {
theta = 0.5*h / (a[ip][iq]);
t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
if (theta < 0.0) { t = -t;}
}
c = 1.0 / sqrt(1+t*t);
s = t*c;
tau = s/(1.0+c);
h = t*a[ip][iq];
z[ip] -= h;
z[iq] += h;
w[ip] -= h;
w[iq] += h;
a[ip][iq]=0.0;
// ip already shifted left by 1 unit
for (j = 0;j <= ip-1;j++) { VTK_ROTATE(a,j,ip,j,iq); }
// ip and iq already shifted left by 1 unit
for (j = ip+1;j <= iq-1;j++) { VTK_ROTATE(a,ip,j,j,iq); }
// iq already shifted left by 1 unit
for (j=iq+1; j<3; j++) { VTK_ROTATE(a,ip,j,iq,j); }
for (j=0; j<3; j++) { VTK_ROTATE(v,j,ip,j,iq); }
}
}
}
for (ip=0; ip<3; ip++) {
b[ip] += z[ip];
w[ip] = b[ip];
z[ip] = 0.0;
}
}
if ( i >= MAX_ROTATIONS ) {
std::cerr << "Matrix3D: Error extracting eigenfunctions" << std::endl;
return MB_FAILURE;
}
// sort eigenfunctions these changes do not affect accuracy
for (j=0; j<2; j++){ // boundary incorrect
k = j;
tmp = w[k];
for (i=j+1; i<3; i++){ // boundary incorrect, shifted already
if (w[i] >= tmp){ // why exchage if same?
k = i;
tmp = w[k];
}
}
if (k != j){
w[k] = w[j];
w[j] = tmp;
for (i=0; i<3; i++){
tmp = v[i][j];
v[i][j] = v[i][k];
v[i][k] = tmp;
}
}
}
// insure eigenvector consistency (i.e., Jacobi can compute vectors that
// are negative of one another (.707,.707,0) and (-.707,-.707,0). This can
// reek havoc in hyperstreamline/other stuff. We will select the most
// positive eigenvector.
int ceil_half_n = (3 >> 1) + (3 & 1);
for (j=0; j<3; j++) {
for (numPos=0, i=0; i<3; i++) {
if ( v[i][j] >= 0.0 ) { numPos++; }
}
// if ( numPos < ceil(double(n)/double(2.0)) )
if ( numPos < ceil_half_n) {
for(i=0; i<3; i++) { v[i][j] *= -1.0; }
}
}
return MB_SUCCESS;
}
| Matrix moab::Matrix::inverse | ( | const Matrix & | d, |
| const double | i | ||
| ) |
Definition at line 43 of file Matrix3.hpp.
{
Matrix m( d);
m( 0) = i * (d(4) * d(8) - d(5) * d(7));
m( 1) = i * (d(2) * d(7) - d(8) * d(1));
m( 2) = i * (d(1) * d(5) - d(4) * d(2));
m( 3) = i * (d(5) * d(6) - d(8) * d(3));
m( 4) = i * (d(0) * d(8) - d(6) * d(2));
m( 5) = i * (d(2) * d(3) - d(5) * d(0));
m( 6) = i * (d(3) * d(7) - d(6) * d(4));
m( 7) = i * (d(1) * d(6) - d(7) * d(0));
m( 8) = i * (d(0) * d(4) - d(3) * d(1));
return m;
}
| const Matrix moab::Matrix::inverse | ( | const Matrix & | d | ) | [inline] |
Definition at line 105 of file Matrix3.hpp.
{
const double det = 1.0/determinant3( d);
return inverse( d, det);
}
| Vector moab::Matrix::matrix_vector | ( | const Matrix & | m, |
| const Vector & | v | ||
| ) | [inline] |
Definition at line 118 of file Matrix3.hpp.
{
Vector res = v;
res[ 0] = v[0] * m(0,0) + v[1] * m(0,1) + v[2] * m(0,2);
res[ 1] = v[0] * m(1,0) + v[1] * m(1,1) + v[2] * m(1,2);
res[ 2] = v[0] * m(2,0) + v[1] * m(2,1) + v[2] * m(2,2);
return res;
}
| Matrix moab::Matrix::mmult3 | ( | const Matrix & | a, |
| const Matrix & | b | ||
| ) | [inline] |
Definition at line 73 of file Matrix3.hpp.
{
return Matrix( a(0,0) * b(0,0) + a(0,1) * b(1,0) + a(0,2) * b(2,0),
a(0,0) * b(0,1) + a(0,1) * b(1,1) + a(0,2) * b(2,1),
a(0,0) * b(0,2) + a(0,1) * b(1,2) + a(0,2) * b(2,2),
a(1,0) * b(0,0) + a(1,1) * b(1,0) + a(1,2) * b(2,0),
a(1,0) * b(0,1) + a(1,1) * b(1,1) + a(1,2) * b(2,1),
a(1,0) * b(0,2) + a(1,1) * b(1,2) + a(1,2) * b(2,2),
a(2,0) * b(0,0) + a(2,1) * b(1,0) + a(2,2) * b(2,0),
a(2,0) * b(0,1) + a(2,1) * b(1,1) + a(2,2) * b(2,1),
a(2,0) * b(0,2) + a(2,1) * b(1,2) + a(2,2) * b(2,2) );
}
| Matrix moab::Matrix::outer_product | ( | const Vector & | u, |
| const Vector & | v, | ||
| Matrix & | m | ||
| ) | [inline] |
Definition at line 86 of file Matrix3.hpp.
{
m = Matrix( u[0] * v[0], u[0] * v[1], u[0] * v[2],
u[1] * v[0], u[1] * v[1], u[1] * v[2],
u[2] * v[0], u[2] * v[1], u[2] * v[2] );
return m;
}
| bool moab::Matrix::positive_definite | ( | const Matrix & | d, |
| double & | det | ||
| ) | [inline] |
Definition at line 58 of file Matrix3.hpp.
{
double subdet6 = d(1)*d(5)-d(2)*d(4);
double subdet7 = d(2)*d(3)-d(0)*d(5);
double subdet8 = d(0)*d(4)-d(1)*d(3);
det = d(6)*subdet6 + d(7)*subdet7 + d(8)*subdet8;
return d(0) > 0 && subdet8 > 0 && det > 0;
}
| Matrix moab::Matrix::transpose | ( | const Matrix & | d | ) | [inline] |
Definition at line 67 of file Matrix3.hpp.
{
return Matrix( d(0), d(3), d(6),
d(1), d(4), d(7),
d(2), d(5), d(8) );
}
| Vector moab::Matrix::vector_matrix | ( | const Vector & | v, |
| const Matrix & | m | ||
| ) | [inline] |
Definition at line 111 of file Matrix3.hpp.
{
return Vector( v[0] * m(0,0) + v[1] * m(1,0) + v[2] * m(2,0),
v[0] * m(0,1) + v[1] * m(1,1) + v[2] * m(2,1),
v[0] * m(0,2) + v[1] * m(1,2) + v[2] * m(2,2) );
}
1.7.6.1