returned grid Position incorrect when certian coordinate system combinations are used

[GoogleCodeExporter]

What is the expected output? What do you see instead?

Expected output from GRID.Position() is the grid coordinates in the basic 
rectangular coordinate system. However, when the grid is defined in:
 * A spherical coordinate system that is defined in CS 0
 * A CS that is defined in a cylindrical or spherical CS
the results are incorrect.

Please see the attached "Position_bug.pptx" for further info.

Looking through the code it seems like something is missing from the 
cylindrical and spherical transformations. I will try to look at it more deeply.

What steps will reproduce the problem?

Define a grid in a coordinate system as shown above, then use the 
GRID.Position() method to get the position.

An example is attached. Unzip "Position_bug.zip" and run the script 
"test_position.py". This reads in the bulk data and calculates the grid 
positions. All of the positions should be [30., 40., 50.]; in most cases they 
are not.

What version of the product are you using? On what operating system? 

Python 2.7.5 on Windows 7
pyNastran version 0.7 dev, svn rev 2049

Original issue reported on code.google.com by jeffrey....@gmail.com on 23 Jun 2014 at 2:44

Attachments:

• Position_bug.pptx
• Position_bug.zip

[GoogleCodeExporter]

After looking at the code, I think the problem could be in 
coordinateSystems.Coord.setup

        try:
            assert len(self.e1) == 3, self.e1
            assert len(self.e2) == 3, self.e2
            assert len(self.e3) == 3, self.e3
            # e_{13}
   ---->    e13 = self.e3 - self.e1
            # e_{12}
   ---->    e12 = self.e2 - self.e1
            #print("e13 = %s" % e13)
            #print("e12 = %s" % e12)

This is not correct for cylindrical and spherical coordinate systems where the 
one component (for cylindrical) or two components (for spherical) coordinates 
are angles.

For example, with a rectangular coordinate system 10 defined in a cylindrical 
coordinate system 1:

    CORD2C*                1               0              0.              0.
    *                     0.              0.              0.              1.*       
    *                     1.              0.              1.
    $ rectangular defined in cylindrical
    CORD2R        10       1     10.      0.      5.     10.     90.      5.   
                 10.      0.      6.

Using the definitions in the code:

     e1 = [10.,  0., 5.]
     e2 = [10., 90., 5.]
     e3 = [10.,  0., 6.]

     e12 = e2 - e1 = [0., 90., 0.]
     e13 = e3 - e1 = [0., 0., 1.]

     k = e12 / norm(e12) = [0, 1, 0]
     j  = cross(e13, k) / norm(cross(e13, k)) = [-1, 0, 0]
     i = cross(j, k) =  [0, 0, -1]

The correct unit vectors for csys 10 are:

    i = [0, 0, 1]
    j = [.7071, .7071, 0]
    k = [-.7071, .7071, 0]

Which do not match the vectors produced by the code.

Original comment by jeffrey....@gmail.com on 23 Jun 2014 at 8:10

[GoogleCodeExporter]

I'll need to take a look at the definition of the CORD2C and CORD2S cards and 
the code in more detail, but it looks like e1, e2, and e3 need to be resolved 
into coordinate system 0 before being used (I thought they were).

The spherical coordinate defined with a rectangular system should have the 
correct i, j, and k vectors, so there is probably an issue with the transform 
to global or from local.

Original comment by mesheb82 on 23 Jun 2014 at 8:46

[GoogleCodeExporter]

The bugs have been fixed.

1.  There was a bug in the spherical transform (theta & phi were flipped).

2.  e1, e2, and e3 were not resolved into cid=0 before being used (fixes i,j,k 
bug).

3.  The origin (in global xyz) should be used as opposed to e1 in the transform 
step.  The origin has been added as self.origin.

Additionally, the argument resolveAltCoord in the coordinate system method 
transformToGlobal has been removed to reduce confusion.

Also, thanks for the very thorough job on the problem.

Original comment by mesheb82 on 6 Jul 2014 at 1:52

• Changed state: Fixed