diff --git a/modules-local/beamdyn/src/BeamDyn_Subs.f90 b/modules-local/beamdyn/src/BeamDyn_Subs.f90 index 89dd81e459..6cc143f969 100644 --- a/modules-local/beamdyn/src/BeamDyn_Subs.f90 +++ b/modules-local/beamdyn/src/BeamDyn_Subs.f90 @@ -339,15 +339,10 @@ SUBROUTINE BD_CrvExtractCrv(R, cc, ErrStat, ErrMsg) CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None !Local variables - REAL(BDKi) :: pivot(4) ! Trace of the rotation matrix and diagonal elements - REAL(BDKi) :: sm0 - REAL(BDKi) :: sm1 - REAL(BDKi) :: sm2 - REAL(BDKi) :: sm3 - REAL(BDKi) :: em - REAL(BDKi) :: temp + integer :: ipivot(1) + REAL(BDKi) :: sm(0:3) + REAL(BDKi) :: em REAL(BDKi) :: Rr(3,3) ! (possibly) corrected rotation matrix - INTEGER :: i ! case indicator INTEGER(IntKi) :: ErrStat2 ! Temporary Error status CHARACTER(ErrMsgLen) :: ErrMsg2 ! Temporary Error message @@ -376,54 +371,51 @@ SUBROUTINE BD_CrvExtractCrv(R, cc, ErrStat, ErrMsg) ! mjs--determine whether R is a valid rotation matrix and correct it if not call BD_CheckRotMat(R, Rr, ErrStat2, ErrMsg2) CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) - if (ErrStat >= AbortErrLev) return + if (ErrStat >= AbortErrLev) return - ! mjs--find max value of T := Tr(Rr) and diagonal elements of Rr - pivot = (/ Rr(1, 1), Rr(2, 2), Rr(3, 3), Rr(1,1) + Rr(2,2) + Rr(3,3) /) - i = maxloc(pivot, 1) - - !> - Condition 1: \f$ \underline{\underline{R(3,3)}} > T \f$ - !! This implies that \f$ c_3^2 > c_0^2 \f$ - IF (i == 3) THEN - sm0 = Rr(2,1) - Rr(1,2) ! 4 c_0 c_3 t_{r0} - sm1 = Rr(1,3) + Rr(3,1) ! 4 c_1 c_3 t_{r0} - sm2 = Rr(2,3) + Rr(3,2) ! 4 c_2 c_3 t_{r0} - sm3 = 1.0_BDKi - Rr(1,1) - Rr(2,2) + Rr(3,3) ! 4 c_3^2 t_{r0} - temp = SIGN( 2.0_BDKi*SQRT(ABS(sm3)), sm0 ) - - !> - Condition 2: \f$ \underline{\underline{R(2,2)}} > T \f$ - !! This implies that \f$ c_2^2 > c_0^2 \f$ - ELSEIF (i == 2) THEN - sm0 = Rr(1,3) - Rr(3,1) ! 4 c_0 c_2 t_{r0} - sm1 = Rr(1,2) + Rr(2,1) ! 4 c_1 c_2 t_{r0} - sm2 = 1.0_BDKi - Rr(1,1) + Rr(2,2) - Rr(3,3) ! 4 c_2^2 t_{r0} - sm3 = Rr(2,3) + Rr(3,2) ! 4 c_2 c_3 t_{r0} - temp = SIGN( 2.0_BDKi*SQRT(ABS(sm2)), sm0 ) - - !> - Condition 3: \f$ \underline{\underline{R(1,1)}} > T \f$ - !! This implies that \f$ c_1^2 > c_0^2 \f$ - ELSEIF (i == 1) THEN - sm0 = Rr(3,2) - Rr(2,3) ! 4 c_0 c_1 t_{r0} - sm1 = 1.0_BDKi + Rr(1,1) - Rr(2,2) - Rr(3,3) ! 4 c_1^2 t_{r0} - sm2 = Rr(1,2) + Rr(2,1) ! 4 c_1 c_2 t_{r0} - sm3 = Rr(1,3) + Rr(3,1) ! 4 c_1 c_3 t_{r0} - temp = SIGN( 2.0_BDKi*SQRT(ABS(sm1)), sm0 ) - - !> - Condition 4: all diagonal terms are less than the trace - ELSE + ! Note that 0 <= c_0 <= 2 since 0 <= dot_product(cc,cc) <= 16 + + ! Let's first find out which parameter is largest so we can divide by it later + sm(0) = 1.0_BDKi + Rr(1,1) + Rr(2,2) + Rr(3,3) ! 4 c_0 c_0 t_{r0} + sm(1) = 1.0_BDKi + Rr(1,1) - Rr(2,2) - Rr(3,3) ! 4 c_1 c_1 t_{r0} + sm(2) = 1.0_BDKi - Rr(1,1) + Rr(2,2) - Rr(3,3) ! 4 c_2 c_2 t_{r0} + sm(3) = 1.0_BDKi - Rr(1,1) - Rr(2,2) + Rr(3,3) ! 4 c_3 c_3 t_{r0} + ipivot = maxloc(sm) - 1 ! our sm array starts at 0, so we need to subtract 1 here to get the correct index + + select case (ipivot(1)) + case (3) + !! We need the condition that c_3 is not zero. + sm(0) = Rr(2,1) - Rr(1,2) ! 4 c_0 c_3 t_{r0} + sm(1) = Rr(1,3) + Rr(3,1) ! 4 c_1 c_3 t_{r0} + sm(2) = Rr(2,3) + Rr(3,2) ! 4 c_2 c_3 t_{r0} + !sm(3) = 1.0_BDKi - Rr(1,1) - Rr(2,2) + Rr(3,3) ! 4 c_3 c_3 t_{r0} + + case (2) + !! We need the condition that c_2 is not zero. + sm(0) = Rr(1,3) - Rr(3,1) ! 4 c_0 c_2 t_{r0} + sm(1) = Rr(1,2) + Rr(2,1) ! 4 c_1 c_2 t_{r0} + !sm(2) = 1.0_BDKi - Rr(1,1) + Rr(2,2) - Rr(3,3) ! 4 c_2 c_2 t_{r0} + sm(3) = Rr(2,3) + Rr(3,2) ! 4 c_3 c_2 t_{r0} + + case (1) + !! We need the condition that c_1 is not zero. + sm(0) = Rr(3,2) - Rr(2,3) ! 4 c_0 c_1 t_{r0} + !sm(1) = 1.0_BDKi + Rr(1,1) - Rr(2,2) - Rr(3,3) ! 4 c_1 c_1 t_{r0} + sm(2) = Rr(1,2) + Rr(2,1) ! 4 c_2 c_1 t_{r0} + sm(3) = Rr(1,3) + Rr(3,1) ! 4 c_3 c_1 t_{r0} + + case (0) + !! We need the condition that c_0 is not zero. !ipivot = 0 - sm0 = 1.0_BDKi + Rr(1,1) + Rr(2,2) + Rr(3,3) ! 4 c_0^2 t_{r0} - sm1 = Rr(3,2) - Rr(2,3) ! -4 c_0 c_1 t_{r0} - sm2 = Rr(1,3) - Rr(3,1) ! -4 c_0 c_2 t_{r0} - sm3 = Rr(2,1) - Rr(1,2) ! -4 c_0 c_3 t_{r0} - temp = SIGN( 2.0_BDKi*SQRT(ABS(sm0)), sm0 ) - ENDIF - - em = sm0 + temp - em = 4.0_BDKi/em - cc(1) = em*sm1 - cc(2) = em*sm2 - cc(3) = em*sm3 + !sm(0) = 1.0_BDKi + Rr(1,1) + Rr(2,2) + Rr(3,3) ! 4 c_0 c_0 t_{r0} + sm(1) = Rr(3,2) - Rr(2,3) ! 4 c_1 c_0 t_{r0} + sm(2) = Rr(1,3) - Rr(3,1) ! 4 c_2 c_0 t_{r0} + sm(3) = Rr(2,1) - Rr(1,2) ! 4 c_3 c_0 t_{r0} + end select + + em = sm(0) + SIGN( 2.0_BDKi*SQRT(sm(ipivot(1))), sm(0) ) + em = 4.0_BDKi/em ! 1 / ( 4 t_{r0} c_{ipivot} ), assuming 0 <= c_0 < 4 and c_{ipivot} > 0 + cc = em*sm(1:3) END SUBROUTINE BD_CrvExtractCrv !------------------------------------------------------------------------------ diff --git a/modules-local/beamdyn/src/Driver_Beam_Subs.f90 b/modules-local/beamdyn/src/Driver_Beam_Subs.f90 index 9363b6326e..19927d800f 100644 --- a/modules-local/beamdyn/src/Driver_Beam_Subs.f90 +++ b/modules-local/beamdyn/src/Driver_Beam_Subs.f90 @@ -27,26 +27,22 @@ module BeamDyn_driver_subs TYPE , PUBLIC :: BD_DriverInternalType REAL(ReKi) , DIMENSION(1:6) :: DistrLoad !< Constant distributed load along beam axis, 3 forces and 3 moments [-] REAL(ReKi) , DIMENSION(1:6) :: TipLoad !< Constant point load applied at tip, 3 forces and 3 moments [-] - INTEGER(IntKi) :: NumPointLoads !< Number of constant point loads applied along beam axis, 3 forces and 3 moments [-] - REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: MultiPointLoad !< Constant point loads applied along beam axis (index 1= Relative position along blade span; indices 2-7 = Fx, Fy, Fz, Mx, My, Mz) [-] -! INTEGER(IntKi) :: NumPointLoads !< Number of multi-point loads in the driver input file -! REAL(DbKi) :: MultiPointLoads{:}{:} !< The array of multipoint loads Index 1: [1, NumPointLoads]; - !< Index 2: [1,7] (index of Loads 1 = Relative position along blade span; - !! 2-7 = Fx, Fy, Fz, Mx, My, Mz ) - - TYPE(MeshType) :: mplMotion ! Mesh for blade motion at multipoint loads locations - TYPE(MeshType) :: mplLoads ! Mesh for multipoint loads - TYPE(MeshMapType) :: Map_BldMotion_to_mplMotion - TYPE(MeshMapType) :: Map_mplLoads_to_PointLoad - TYPE(MeshType) :: y_BldMotion_at_u_point ! Intermediate mesh to transfer motion from output mesh to input mesh - TYPE(MeshMapType) :: Map_y_BldMotion_to_u_point - - TYPE(MeshType) :: RotationCenter - TYPE(MeshMapType) :: Map_RotationCenter_to_RootMotion + INTEGER(IntKi) :: NumPointLoads !< Number of constant point loads applied along beam axis, 3 forces and 3 moments, from the driver input file [-] + REAL(BDKi) , DIMENSION(:,:), ALLOCATABLE :: MultiPointLoad !< Constant point loads applied along beam axis, size (NumPointLoads,7); (dimension 2: index 1=Relative position along blade span; indices 2-7 = Fx, Fy, Fz, Mx, My, Mz) [-] + + TYPE(MeshType) :: mplMotion ! Mesh for blade motion at multipoint loads locations + TYPE(MeshType) :: mplLoads ! Mesh for multipoint loads + TYPE(MeshMapType) :: Map_BldMotion_to_mplMotion + TYPE(MeshMapType) :: Map_mplLoads_to_PointLoad + TYPE(MeshType) :: y_BldMotion_at_u_point ! Intermediate mesh to transfer motion from output mesh to input mesh + TYPE(MeshMapType) :: Map_y_BldMotion_to_u_point - REAL(DbKi) :: t_initial - REAL(DbKi) :: t_final - REAL(ReKi) :: w ! magnitude of rotational velocity vector + TYPE(MeshType) :: RotationCenter + TYPE(MeshMapType) :: Map_RotationCenter_to_RootMotion + + REAL(DbKi) :: t_initial + REAL(DbKi) :: t_final + REAL(R8Ki) :: w ! magnitude of rotational velocity vector END TYPE @@ -83,6 +79,7 @@ SUBROUTINE BD_ReadDvrFile(DvrInputFile,dt,InitInputData,DvrData,& INTEGER(IntKi) :: UnEc CHARACTER(1024) :: FTitle ! "File Title": the 2nd line of the input file, which contains a description of its contents + CHARACTER(1024) :: PriPath ! Path name of the primary file INTEGER(IntKi) :: i !------------------------------------------------------------------------------------ @@ -98,6 +95,8 @@ SUBROUTINE BD_ReadDvrFile(DvrInputFile,dt,InitInputData,DvrData,& CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) IF (ErrStat >= AbortErrLev) RETURN + CALL GetPath( DvrInputFile, PriPath ) ! Input files will be relative to the path where the primary input file is located. + !-------------------------- HEADER --------------------------------------------- CALL ReadCom(UnIn,DvrInputFile,'File Header: Module Version (line 1)',ErrStat2,ErrMsg2,UnEc) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) @@ -270,6 +269,8 @@ SUBROUTINE BD_ReadDvrFile(DvrInputFile,dt,InitInputData,DvrData,& !---------------------- BEAM SECTIONAL PARAMETER ---------------------------------------- CALL ReadVar ( UnIn, DvrInputFile, InitInputData%InputFile, 'InputFile', 'Name of the primary input file', ErrStat2,ErrMsg2, UnEc ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + IF ( PathIsRelative( InitInputData%InputFile ) ) InitInputData%InputFile = TRIM(PriPath)//TRIM(InitInputData%InputFile) + if (ErrStat >= AbortErrLev) then call cleanup() return @@ -586,20 +587,20 @@ SUBROUTINE Init_RotationCenterMesh(DvrData, InitInputData, RootMotionMesh, ErrSt DvrData%w = TwoNorm( InitInputData%RootVel(4:6) ) - if (EqualRealNos(DvrData%w,0.0_ReKi)) then - DvrData%w = 0.0_ReKi + if (EqualRealNos(DvrData%w,0.0_R8Ki)) then + DvrData%w = 0.0_R8Ki ! the beam is not rotating, so pick an orientation call eye(orientation, ErrStat2, ErrMsg2) else z_hat = InitInputData%RootVel(4:6) / DvrData%w - if ( EqualRealNos( z_hat(3), 1.0_ReKi ) ) then + if ( EqualRealNos( z_hat(3), 1.0_R8Ki ) ) then call eye(orientation, ErrStat2, ErrMsg2) - elseif ( EqualRealNos( z_hat(3), -1.0_ReKi ) ) then + elseif ( EqualRealNos( z_hat(3), -1.0_R8Ki ) ) then orientation = 0.0_ReKi - orientation(1,1) = -1.0_ReKi - orientation(2,2) = 1.0_ReKi - orientation(3,3) = -1.0_ReKi + orientation(1,1) = -1.0_R8Ki + orientation(2,2) = 1.0_R8Ki + orientation(3,3) = -1.0_R8Ki else Z_unit = (/0, 0, 1/) diff --git a/modules-local/orcaflex-interface/src/OrcaDriver_Subs.f90 b/modules-local/orcaflex-interface/src/OrcaDriver_Subs.f90 index e1306b5b48..e07b6105c1 100644 --- a/modules-local/orcaflex-interface/src/OrcaDriver_Subs.f90 +++ b/modules-local/orcaflex-interface/src/OrcaDriver_Subs.f90 @@ -1907,12 +1907,13 @@ SUBROUTINE PointsForce_OutputWrite(ProgInfo, OutUnit, OutFileName, InputFileName ! Temporary local variables INTEGER(IntKi) :: ErrStatTmp !< Temporary variable for the status of error message - CHARACTER(*), PARAMETER :: RoutineName = 'PointsForce_OutputWrite' + CHARACTER(*), PARAMETER :: RoutineName = 'PointsForce_OutputWrite' CHARACTER(2048) :: ErrMsgTmp !< Temporary variable for the error message INTEGER(IntKi) :: LenErrMsgTmp !< Length of ErrMsgTmp (for getting WindGrid info) REAL(ReKi) :: rotdisp(3) !< Rotational displacement (euler angles) INTEGER(IntKi) :: I !< Generic counter - + REAL(ReKi) :: outputArray(13) + CHARACTER(47) :: PointsOutputFmt !< Format specifier for the output file for wave elevation series CHARACTER(3) :: AngleUnit !< Units for the angle @@ -1972,32 +1973,30 @@ SUBROUTINE PointsForce_OutputWrite(ProgInfo, OutUnit, OutFileName, InputFileName WRITE (OutUnit,'(A)', IOSTAT=ErrStatTmp ) '' ENDIF - rotdisp = GetSmllRotAngs ( u%PtfmMesh%Orientation(:,:,1), ErrStatTmp, ErrMsgTmp ) CALL SetErrStat( ErrStatTmp, ErrMsgTmp, ErrStat, ErrMsg, RoutineName ) - IF ( AnglesInDegrees ) THEN - CALL WrNumAryFileNR( OutUnit, (/ REAL(Time, ReKi), & - u%PtfmMesh%TranslationDisp(1,1), u%PtfmMesh%TranslationDisp(2,1), u%PtfmMesh%TranslationDisp(3,1), & - rotdisp(1)*R2D, rotdisp(2)*R2D, rotdisp(3)*R2D, & - u%PtfmMesh%TranslationVel(1,1), u%PtfmMesh%TranslationVel(2,1), u%PtfmMesh%TranslationVel(3,1), & - u%PtfmMesh%RotationVel(1,1)*R2D, u%PtfmMesh%RotationVel(2,1)*R2D, u%PtfmMesh%RotationVel(3,1)*R2D /), & - '3x,ES10.3E2', ErrStatTmp, ErrMsgTmp ) - CALL WrNumAryFileNR( OutUnit, y%WriteOutput, '3x,ES10.3E2', ErrStatTmp, ErrMsgTmp ) - WRITE (OutUnit,'(A)', IOSTAT=ErrStatTmp ) '' + outputArray = (/ REAL(Time, ReKi), & + REAL(u%PtfmMesh%TranslationDisp(1,1), ReKi), & + REAL(u%PtfmMesh%TranslationDisp(2,1), ReKi), & + REAL(u%PtfmMesh%TranslationDisp(3,1), ReKi), & + rotdisp(1)*R2D, rotdisp(2)*R2D, rotdisp(3)*R2D, & + u%PtfmMesh%TranslationVel(1,1), u%PtfmMesh%TranslationVel(2,1), u%PtfmMesh%TranslationVel(3,1), & + u%PtfmMesh%RotationVel(1,1)*R2D, u%PtfmMesh%RotationVel(2,1)*R2D, u%PtfmMesh%RotationVel(3,1)*R2D /) ELSE - CALL WrNumAryFileNR( OutUnit, (/ REAL(Time, ReKi), & - u%PtfmMesh%TranslationDisp(1,1), u%PtfmMesh%TranslationDisp(2,1), u%PtfmMesh%TranslationDisp(3,1), & - rotdisp(1), rotdisp(2), rotdisp(3), & - u%PtfmMesh%TranslationVel(1,1), u%PtfmMesh%TranslationVel(2,1), u%PtfmMesh%TranslationVel(3,1), & - u%PtfmMesh%RotationVel(1,1), u%PtfmMesh%RotationVel(2,1), u%PtfmMesh%RotationVel(3,1) /), & - '3x,ES10.3E2', ErrStatTmp, ErrMsgTmp ) - CALL WrNumAryFileNR( OutUnit, y%WriteOutput, '3x,ES10.3E2', ErrStatTmp, ErrMsgTmp ) - WRITE (OutUnit,'(A)', IOSTAT=ErrStatTmp ) '' + outputArray = (/ REAL(Time, ReKi), & + REAL(u%PtfmMesh%TranslationDisp(1,1), ReKi), & + REAL(u%PtfmMesh%TranslationDisp(2,1), ReKi), & + REAL(u%PtfmMesh%TranslationDisp(3,1), ReKi), & + rotdisp(1), rotdisp(2), rotdisp(3), & + u%PtfmMesh%TranslationVel(1,1), u%PtfmMesh%TranslationVel(2,1), u%PtfmMesh%TranslationVel(3,1), & + u%PtfmMesh%RotationVel(1,1), u%PtfmMesh%RotationVel(2,1), u%PtfmMesh%RotationVel(3,1) /) ENDIF - - + + CALL WrNumAryFileNR( OutUnit, outputArray, '3x,ES10.3E2', ErrStatTmp, ErrMsgTmp ) + CALL WrNumAryFileNR( OutUnit, y%WriteOutput, '3x,ES10.3E2', ErrStatTmp, ErrMsgTmp ) + WRITE (OutUnit,'(A)', IOSTAT=ErrStatTmp ) '' END SUBROUTINE PointsForce_OutputWrite