SD: Implementing directional cosine matrices and section properties f…

…or rectangular members (#1413)

.gitignore

@@ -41,6 +41,7 @@ vs-build/
 .vscode
 .atom
 .fortls
+.devcontainer
 # backup files
 *.asv
 ~$*.xlsx

docs/source/user/subdyn/input_files.rst

@@ -427,8 +427,9 @@ A member is one of the three following types (see :numref:`SD_FEM`):
 
 - Rigid link (*MType=3*)
 
-**COSMID** refers to the IDs of the members’ cosine matrices for
-noncircular members; the current release ignores this column.
+**COSMID** refers to the IDs of the members' cosine matrices for
+noncircular members; the current release uses SubDyn's default direction
+ cosine convention if it's not present or when COSMID values are -1.
 
 
 An example of member table is given below
@@ -543,17 +544,13 @@ An example of rigid link properties table is given below
 
 Member Cosine Matrices COSM (i,j)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-This table is not currently used by SubDyn, but in future releases it
-will need to be populated if members with cross-sections other than
-circular will be employed.
-
 **NCOSMs** rows, one for each unique member orientation set, will need
 to be provided. Each row of the table will list the nine entries of the
-direction cosine matrices (COSM11, COSM12,…COSM33) for matrix elements
-(1,1), (1,2),…(3,3) that establish the orientation of the local member
-axes (*x*,\ *y* principal axes in the cross-sectional plane, *z* along
-the member longitudinal axis) with respect to the SS coordinate system
-(local-to-global transformation matrices).
+direction cosine matrices (COSM11, COSM12,…COSM33) for matrix elements.
+Each row is a vector in the global coordinate system for principal axes 
+in the x, y and z directions respectively. These vectors need to be 
+specified with an extremely high level of precision for results to be
+equivalent to an internal calculation.
 
 Joint Additional Concentrated Masses
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

modules/subdyn/src/FEM.f90

@@ -26,8 +26,9 @@ MODULE FEM
   INTEGER, PARAMETER  :: LaKi = R8Ki  ! Define the kind to be used for LaPack
 
   INTERFACE FINDLOCI ! In the future, use FINDLOC from intrinsic
-     MODULE PROCEDURE FINDLOCI_ReKi
+     MODULE PROCEDURE FINDLOCI_R8Ki
      MODULE PROCEDURE FINDLOCI_IntKi
+     MODULE PROCEDURE FINDLOCI_SiKi
   END INTERFACE
 
 
@@ -836,8 +837,8 @@ END SUBROUTINE InsertDOFrows
 !------------------------------------------------------------------------------------------------------
 !> Returns index of val in Array (val is an integer!)
 ! NOTE: in the future use intrinsinc function findloc
-FUNCTION FINDLOCI_ReKi(Array, Val) result(i)
-   real(ReKi)    , dimension(:), intent(in) :: Array !< Array to search in
+FUNCTION FINDLOCI_R8Ki(Array, Val) result(i)
+   real(R8Ki)    , dimension(:), intent(in) :: Array !< Array to search in
    integer(IntKi), intent(in)               :: val   !< Val
    integer(IntKi)                           :: i     !< Index of joint in joint table
    i = 1
@@ -866,6 +867,21 @@ FUNCTION FINDLOCI_IntKi(Array, Val) result(i)
    enddo
    i=-1
 END FUNCTION
+
+FUNCTION FINDLOCI_SiKi(Array, Val) result(i)
+   real(SiKi), dimension(:), intent(in) :: Array !< Array to search in
+   integer(IntKi), intent(in)               :: val   !< Val
+   integer(IntKi)                           :: i     !< Index of joint in joint table
+   i = 1
+   do while ( i <= size(Array) )
+      if ( Val == Array(i) ) THEN
+         return ! Exit when found
+      else
+         i = i + 1
+      endif
+   enddo
+   i=-1
+END FUNCTION
 !------------------------------------------------------------------------------------------------------
 SUBROUTINE RigidTransformationLine(dx,dy,dz,iLine,Line)
    real(ReKi),               INTENT(IN)  :: dx,dy,dz
@@ -1022,17 +1038,17 @@ END SUBROUTINE GetOrthVectors
 !> Element stiffness matrix for classical beam elements
 !! shear is true  -- non-tapered Timoshenko beam 
 !! shear is false -- non-tapered Euler-Bernoulli beam 
-SUBROUTINE ElemK_Beam(A, L, Ixx, Iyy, Jzz, Shear, kappa, E, G, DirCos, K)
-   REAL(ReKi), INTENT( IN) :: A, L, Ixx, Iyy, Jzz, E, G, kappa
+SUBROUTINE ElemK_Beam(A, L, Ixx, Iyy, Jzz, Shear, kappa_x, kappa_y, E, G, DirCos, K)
+   REAL(ReKi), INTENT( IN) :: A, L, Ixx, Iyy, Jzz, E, G, kappa_x, kappa_y
    REAL(FEKi), INTENT( IN) :: DirCos(3,3) !< From element to global: xg = DC.xe,  Kg = DC.Ke.DC^t
    LOGICAL   , INTENT( IN) :: Shear
    REAL(FEKi), INTENT(OUT) :: K(12, 12) 
    ! Local variables
    REAL(FEKi)                            :: Ax, Ay, Kx, Ky
    REAL(FEKi)                            :: DC(12, 12)
 
-   Ax = kappa*A
-   Ay = kappa*A
+   Ax = kappa_x*A
+   Ay = kappa_y*A
 
    K(1:12,1:12) = 0.0_FEKi