anchor.diff:  
      INTEGER(IntKi)              :: idx_qp      !< index of current quadrature point
      
  -   DO idx_qp=1,p%nqp
  -      !m%qp%betaC(:,:,idx_qp,nelem) = MATMUL( diag(p%beta(i)), temp_b,m%qp%Stif(:,:,idx_qp,nelem))
  -      DO j=1,6
  -         DO i=1,6
  -            m%qp%betaC(i,j,idx_qp,nelem) = p%beta(i)*m%qp%Stif(i,j,idx_qp,nelem)
  -         END DO
  -      END DO
  -   END DO
  -
      
      IF (.NOT. fact) then ! skip all but Fc and Fd terms
      
  
positive.code: 
           ! Setup the temporary matrix for modifying the stiffness matrix. RR0 is changing with time.
        tempR6 = 0.0_BDKi
        tempR6(1:3,1:3) = m%qp%RR0(:,:,idx_qp,nelem)       ! upper left   -- translation
        tempR6(4:6,4:6) = m%qp%RR0(:,:,idx_qp,nelem)       ! lower right  -- rotation
  !NOTE: Bauchau has the lower right corner multiplied by H
  
  
           !> Modify the Mass matrix as
           !! \f$ \begin{bmatrix}
           !!        \left(\underline{\underline{R}} \underline{\underline{R}}_0\right)      &  0             \\
           !!                      0  &  \left(\underline{\underline{R}} \underline{\underline{R}}_0\right)
  
negative.code: 
     INTEGER(IntKi)                :: idx_qp         !< index counter for quadrature point
     INTEGER(IntKi)                :: temp_id2       !< Index to last node of previous element
     INTEGER(IntKi)                :: i,j            !< generic counters
     REAL(BDKi)                    :: tempR6(6,6)
     REAL(BDKi)                    :: tempBeta6(6,6)
  
  
     ! see Bauchau 2011 Flexible Multibody Dynamics p 692-693, section 17.7.2
  
  