Local axis follows origin node rotation

I'd like to define a local axis (unit vectors l, m and n) which once defined follow the rotation of the origin node, i.e. regardless of the deformation the local axis should be basically the same as the original one only rotated by some angles.

I know the coordinates of node 1 and node 2 which define the l (x) local axis of the element.

I have to define the m (y) local axis by specifying a vector normal to the x local axis. I want this vector to be on the global Z axis - local x axis plane and initially toward the positive part of the global Z axis (and not the negative part of the global Z axis).

Having defined the local x and y axis the local z axis is automatically defined.

My problem now is that with the code I've developed (see below) the local y axis will always be upward regardless of the rotation of the origin node, which should not happen. The local y axis should be initially upwards but then should follow the origin node rotation.

I've tried to define the local Z axis based on the initial local Y axis orientation but was not successful.

I have the following code in fortran to define the transformation matrix (direction cosine matrix):

C COORDINATES OF NODE 1 and NODE 2 (GIVEN), they define the local x (l) axis along the element's length.
c
X1=COORDS(1,2)+U(7) !NODE 1
Y1=COORDS(2,2)+U(8)
Z1=COORDS(3,2)+U(9)
C
X2=COORDS(1,3)+U(13) !NODE 2
Y2=COORDS(2,3)+U(14)
Z2=COORDS(3,3)+U(15)
C
X21=X2-X1
Y21=Y2-Y1
Z21=Z2-Z1
C
C DIRECTION COSINE MATRIX
C initialize matrix
DO K1 = 1, NDOFEL
DO K2 = 1, NDOFEL
T(K2,K1) = ZERO
END DO
END DO
c
lX=ZERO
lY=ZERO
lZ=ZERO
mX=ZERO
mY=ZERO
mZ=ZERO
nX=ZERO
nY=ZERO
nZ=ZERO
C
C THIRD NODE: auxiliary node in the local x-y (l-m) plane
P1X3=X1+X21*ONE
P1Y3=Y1+Y21*ONE
P1Z3=Z1+Z21*ONE
! Determine coordinates of fourth NODE
! Define local y axis always up initially
IF (X2-X1.EQ.ZERO) THEN
if (Z2-Z1.EQ.ZERO) then
write(6,*) 'vector along Y'
P2X3=ZERO+P1X3
P2Y3=-Z21*PI+P1Y3
P2Z3=+Y21*PI+P1Z3
if (P2Z3.LT.P1Z3) then
P2Y3=Z21*PI+P1Y3
P2Z3=-Y21*PI+P1Z3
endif
elseif (Y2-Y1.EQ.ZERO) then
write(6,*) 'vector along Z'
P2X3=ZERO+P1X3
P2Y3=-Z21*PI+P1Y3
P2Z3=+Y21*PI+P1Z3
if (P2Z3.LT.P1Z3) then
P2Y3=Z21*PI+P1Y3
P2Z3=-Y21*PI+P1Z3
endif
else
write(6,*) 'vector Y-Z plane'
P2X3=ZERO+P1X3
P2Y3=-Z21*PI+P1Y3
P2Z3=+Y21*PI+P1Z3
if (P2Z3.LT.P1Z3) then
P2Y3=Z21*PI+P1Y3
P2Z3=-Y21*PI+P1Z3
endif
endif
ELSEIF (Y2-Y1.EQ.ZERO) THEN
if (X2-X1.EQ.ZERO) then
write(6,*) 'vector along Z'
P2X3=ZERO+P1X3
P2Y3=-Z21*PI+P1Y3
P2Z3=+Y21*PI+P1Z3
if (P2Z3.LT.P1Z3) then
P2Y3=Z21*PI+P1Y3
P2Z3=-Y21*PI+P1Z3
endif
elseif (Z2-Z1.EQ.ZERO) then
write(6,*) 'vector along X'
P2X3=-Z21*PI+P1X3
P2Y3=ZERO+P1Y3
P2Z3=+X21*PI+P1Z3
if (P2Z3.LT.P1Z3) then
P2X3=Z21*PI+P1X3
P2Z3=-X21*PI+P1Z3
endif
else
write(6,*) 'vector along X-Z plane'
P2X3=-Z21*PI+P1X3
P2Y3=ZERO+P1Y3
P2Z3=+X21*PI+P1Z3
if (P2Z3.LT.P1Z3) then
P2X3=Z21*PI+P1X3
P2Z3=-X21*PI+P1Z3
endif
endif
ELSEIF (Z2-Z1.EQ.ZERO) THEN
if (X2-X1.EQ.ZERO) then
write(6,*) 'vector along Y'
P2X3=ZERO+P1X3
P2Y3=-Z21*PI+P1Y3
P2Z3=+Y21*PI+P1Z3
if (P2Z3.LT.P1Z3) then
P2Y3=Z21*PI+P1Y3
P2Z3=-Y21*PI+P1Z3
endif
elseif (Y2-Y1.EQ.ZERO) then
write(6,*) 'vector along X'
P2X3=-Z21*PI+P1X3
P2Y3=ZERO+P1Y3
P2Z3=+X21*PI+P1Z3
if (P2Z3.LT.P1Z3) then
P2X3=Z21*PI+P1X3
P2Z3=-X21*PI+P1Z3
endif
else
write(6,*) 'vector along X-Y plane'
P2X3=-Z21*PI+P1X3
P2Y3=ZERO+P1Y3
P2Z3=+X21*PI+P1Z3
if (P2Z3.LT.P1Z3) then
P2X3=Z21*PI+P1X3
P2Z3=-X21*PI+P1Z3
endif
endif
ELSE
write(6,*) '3D frame'
P2X3=X1+X21*TWO
P2Y3=Y1+Y21*TWO
c
P2Z3=(-(P2X3-P1X3)*(X2-X1)-(P2Y3-P1Y3)*(Y2-Y1))/(Z2-Z1)+P1Z3
if (P2Y3.LT.P1Y3) then
P2X3=X1-X21*TWO
P2Y3=Y1-Y21*TWO
P2Z3=(-(P2X3-P1X3)*(X2-X1)-(P2Y3-P1Y3)*(Y2-Y1))/(Z2-Z1)+P1Z3
endif
ENDIF
c
X3=P2X3
Y3=P2Y3
Z3=P2Z3
C
X31=X3-X1
Y31=Y3-Y1
Z31=Z3-Z1
C
X32=X3-X2
Y32=Y3-Y2
Z32=Z3-Z2
C
lX=X21/LCE
mX=Y21/LCE
nX=Z21/LCE
C
mA=((Y31*Z21-Y21*Z31)**2+(Z31*X21-Z21*X31)**2+
1   (X31*Y21-X21*Y31)**2)**0.5
C
lZ=1/mA*(Y21*Z31-Y31*Z21) !Right Handed local system, x along el. length
mZ=1/mA*(Z21*X31-Z31*X21)
nZ=1/mA*(X21*Y31-X31*Y21)
C
lY=mZ*nX-nZ*mX
mY=nZ*lX-lZ*nX
nY=lZ*mX-mZ*lX
c
C T(i,j): DIRECTION COSINE MATRIX
T(1,1)=lX
T(1,2)=mX
T(1,3)=nX
T(2,1)=lY
T(2,2)=mY
T(2,3)=nY
T(3,1)=lZ
T(3,2)=mZ
T(3,3)=nZ
C
T(4,4)=lX
T(4,5)=mX
T(4,6)=nX
T(5,4)=lY
T(5,5)=mY
T(5,6)=nY
T(6,4)=lZ
T(6,5)=mZ
T(6,6)=nZ
C
T(7,7)=lX
T(7,8)=mX
T(7,9)=nX
T(8,7)=lY
T(8,8)=mY
T(8,9)=nY
T(9,7)=lZ
T(9,8)=mZ
T(9,9)=nZ
c
T(10,10)=lX
T(10,11)=mX
T(10,12)=nX
T(11,10)=lY
T(11,11)=mY
T(11,12)=nY
T(12,10)=lZ
T(12,11)=mZ
T(12,12)=nZ


Any ideas? Many thanks

Joao

P.S.: To maximize the chances of solving the problem I'll post this also on stack overflow forum. I'll post the solution here if it is solved there.

• It might be helpful if you could provide more details using math rather than posting code as you'd likely get better and more responses. Regards – Amzoti Mar 18 '13 at 14:30

First, my version of your question: given one vector (your l(x) I will call $\mathbf{l}_x$) how to find another vector (your m(y) I will call $\mathbf{m}_y$) that is orthogonal to it, in the $y=0$ plane, and with positive z-coordinate.
This is easy. Let $$\mathbf{l}_x = \pmatrix{a & b & c}$$ Solve $$\mathbf{m}_ y\mathbf{l}_x^\top = \pmatrix{e & 0 & f}\pmatrix{a \\ b \\ c} = 0$$ with the constraints $$f > 0$$ $$e^2 + f^2 = 1$$ In other words, given $a$ and $c$ solve the equation: $$ea+fc=0$$ The solution to this equation is unique when taken along with these two others: \begin{align} e^2 + f^2 &= 1 \quad\text{the vector should be a unit vector} \\ f &> 0 \quad \text{the vector should point in the positive z-direction} \\ \end{align}