Local axis follows origin node rotation

Question

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.

Answer 1

Let me start with this statement--I did not look at your code. I would still like to try to help you understand the math behind what I believe is your question.

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}

I do not know what you mean by nodes in your question, and I used that your "local axis" was a known point. You do not need to use trig functions to solve for orthogonality. Hope this helps.