Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
HomXform.cpp
Go to the documentation of this file.
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 
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 
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