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:          !!                      0  &  \left(\underline{\underline{R}} \underline{\underline{R}}_0\right)^T
           !!     \end{bmatrix} \f$
        m%qp%Stif(:,:,idx_qp,nelem) = MATMUL(tempR6,MATMUL(p%Stif0_QP(1:6,1:6,temp_id2+idx_qp),TRANSPOSE(tempR6)))
     ENDDO
  
  END SUBROUTINE BD_StifAtDeformedQP
  
negative.code: 
           ! Setup the temporary matrix for modifying the stiffness matrix. RR0 is changing with time.
        tempR6 = 0.0_BDKi
        tempBeta6 = 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
  
           ! Move damping ratio from material frame to the calculation reference frame
           !     This is the following:
           !        tempBEta6=matmul(tempR6,matmul(diag(p%beta),transpose(tempR6)))
        do j=1,6
           do i=1,6
                 ! diagonal of p%beta * TRANSPOSE(tempR6)
              tempBeta6(i,j) = p%beta(i)*tempR6(j,i)
           enddo
        enddo
        tempBeta6 = matmul(tempR6,tempBeta6)
  
  
           !> Modify the Mass matrix so it is in the calculation reference frame
           !! \f$ \begin{bmatrix}
           !!        \left(\underline{\underline{R}} \underline{\underline{R}}_0\right)      &  0             \\
           !!                      0  &  \left(\underline{\underline{R}} \underline{\underline{R}}_0\right)
  