 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 /**
2 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3 * storing and accessing finite element mesh data.
4 *
5 * Copyright 2004 Sandia Corporation. Under the terms of Contract
6 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7 * retains certain rights in this software.
8 *
9 * This library is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public
11 * License as published by the Free Software Foundation; either
12 * version 2.1 of the License, or (at your option) any later version.
13 *
14 */
15
16 #include "moab/HomXform.hpp"
17 #include <cassert>
18
19 namespace moab
20 {
21
22 HomCoord HomCoord::unitv[3] = { HomCoord( 1, 0, 0 ), HomCoord( 0, 1, 0 ), HomCoord( 0, 0, 1 ) };
23 HomCoord HomCoord::IDENTITY( 1, 1, 1 );
24
25 HomCoord& HomCoord::getUnitv( int c )
26 {
27 return unitv[c];
28 }
29
30 int dum[] = { 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1 };
31 HomXform HomXform::IDENTITY( dum );
32
33 void HomXform::three_pt_xform( const HomCoord& p1,
34 const HomCoord& q1,
35 const HomCoord& p2,
36 const HomCoord& q2,
37 const HomCoord& p3,
38 const HomCoord& q3 )
39 {
40 // pmin and pmax are min and max bounding box corners which are mapped to
41 // qmin and qmax, resp. qmin and qmax are not necessarily min/max corners,
42 // since the mapping can change the orientation of the box in the q reference
43 // frame. Re-interpreting the min/max bounding box corners does not change
44 // the mapping.
45
46 // change that: base on three points for now (figure out whether we can
47 // just use two later); three points are assumed to define an orthogonal
48 // system such that (p2-p1)%(p3-p1) = 0
49
50 // use the three point rule to compute the mapping, from Mortensen,
51 // "Geometric Modeling". If p1, p2, p3 and q1, q2, q3 are three points in
52 // the two coordinate systems, the three pt rule is:
53 //
54 // v1 = p2 - p1
55 // v2 = p3 - p1
56 // v3 = v1 x v2
57 // w1-w3 similar, with q1-q3
58 // V = matrix with v1-v3 as rows
59 // W similar, with w1-w3
60 // R = V^-1 * W
61 // t = q1 - p1 * W
62 // Form transform matrix M from R, t
63
64 // check to see whether unity transform applies
65 if( p1 == q1 && p2 == q2 && p3 == q3 )
66 {
67 *this = HomXform::IDENTITY;
68 return;
69 }
70
71 // first, construct 3 pts from input
72 HomCoord v1 = p2 - p1;
73 assert( v1.i() != 0 || v1.j() != 0 || v1.k() != 0 );
74 HomCoord v2 = p3 - p1;
75 HomCoord v3 = v1 * v2;
76
77 if( v3.length_squared() == 0 )
78 {
79 // 1d coordinate system; set one of v2's coordinates such that
80 // it's orthogonal to v1
81 if( v1.i() == 0 )
82 v2.set( 1, 0, 0 );
83 else if( v1.j() == 0 )
84 v2.set( 0, 1, 0 );
85 else if( v1.k() == 0 )
86 v2.set( 0, 0, 1 );
87 else
88 assert( false );
89 v3 = v1 * v2;
90 assert( v3.length_squared() != 0 );
91 }
92 // assert to make sure they're each orthogonal
93 assert( v1 % v2 == 0 && v1 % v3 == 0 && v2 % v3 == 0 );
94 v1.normalize();
95 v2.normalize();
96 v3.normalize();
97 // Make sure h is set to zero here, since it'll mess up things if it's one
98 v1.homCoord[3] = v2.homCoord[3] = v3.homCoord[3] = 0;
99
100 HomCoord w1 = q2 - q1;
101 assert( w1.i() != 0 || w1.j() != 0 || w1.k() != 0 );
102 HomCoord w2 = q3 - q1;
103 HomCoord w3 = w1 * w2;
104
105 if( w3.length_squared() == 0 )
106 {
107 // 1d coordinate system; set one of w2's coordinates such that
108 // it's orthogonal to w1
109 if( w1.i() == 0 )
110 w2.set( 1, 0, 0 );
111 else if( w1.j() == 0 )
112 w2.set( 0, 1, 0 );
113 else if( w1.k() == 0 )
114 w2.set( 0, 0, 1 );
115 else
116 assert( false );
117 w3 = w1 * w2;
118 assert( w3.length_squared() != 0 );
119 }
120 // assert to make sure they're each orthogonal
121 assert( w1 % w2 == 0 && w1 % w3 == 0 && w2 % w3 == 0 );
122 w1.normalize();
123 w2.normalize();
124 w3.normalize();
125 // Make sure h is set to zero here, since it'll mess up things if it's one
126 w1.homCoord[3] = w2.homCoord[3] = w3.homCoord[3] = 0;
127
128 // form v^-1 as transpose (ok for orthogonal vectors); put directly into
129 // transform matrix, since it's eventually going to become R
130 *this = HomXform( v1.i(), v2.i(), v3.i(), 0, v1.j(), v2.j(), v3.j(), 0, v1.k(), v2.k(), v3.k(), 0, 0, 0, 0, 1 );
131
132 // multiply by w to get R
133 *this *= HomXform( w1.i(), w1.j(), w1.k(), 0, w2.i(), w2.j(), w2.k(), 0, w3.i(), w3.j(), w3.k(), 0, 0, 0, 0, 1 );
134
135 // compute t and put into last row
136 HomCoord t = q1 - p1 * *this;
137 ( *this ).XFORM( 3, 0 ) = t.i();
138 ( *this ).XFORM( 3, 1 ) = t.j();
139 ( *this ).XFORM( 3, 2 ) = t.k();
140
141 // M should transform p to q
142 assert( ( q1 == p1 * *this ) && ( q2 == p2 * *this ) && ( q3 == p3 * *this ) );
143 }
144
145 } // namespace moab
1.9.1.