anchor.diff:     p%elem_total = InputFileData%member_total      
      ! Number of nodes per elelemt
      p%node_elem  = InputFileData%order_elem + 1   
  -   ! Number of Gauss points
  -   p%ngp        = p%node_elem
  +   ! Quadrature method: 1 Gauss 2 Trapezoidal
  +   p%quadrature = 1
  +   IF(p%quadrature .EQ. 1) THEN
  +       ! Number of Gauss points
  +       p%ngp        = p%node_elem
  +   ELSEIF(p%quadrature .EQ. 2) THEN 
  +       CALL AllocAry(p%ntp,InputFileData%member_total,'Number of trapezoidal point array',ErrStat2,ErrMsg2)
  +          CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  +       p%ntp(:) = 0
  +       p%ntp(:) = InputFileData%kp_member(:)
  +   ENDIF
      ! Degree-of-freedom (DoF) per node
      p%dof_node   = 6
      ! Total number of (finite element) nodes
  
positive.code:    INTEGER(IntKi)          :: j                ! do-loop counter
     INTEGER(IntKi)          :: k                ! do-loop counter
     INTEGER(IntKi)          :: m                ! do-loop counter
     INTEGER(IntKi)          :: id0
     INTEGER(IntKi)          :: id1
     INTEGER(IntKi)          :: temp_int
  
negative.code:          return
        end if
     DEALLOCATE(temp_w)
     IF(p%quadrature .EQ. 1) THEN
         ! temp_L2: the DistrLoad mesh node location
         ! temp_L2: physical coordinates of Gauss points and two end points
         CALL AllocAry(temp_L2,6,p%ngp*p%elem_total+2,'temp_L2',ErrStat2,ErrMsg2)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         CALL AllocAry(p%Gauss,6,p%ngp*p%elem_total+2,'p%Gauss',ErrStat2,ErrMsg2)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         ! Temporary Gauss point intrinsic coordinates array
         CALL AllocAry(temp_GL,p%ngp,'temp_GL',ErrStat2,ErrMsg2)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         ! Temporary Gauss weight function array
         CALL AllocAry(temp_w,p%ngp,'GL weight array',ErrStat2,ErrMsg2)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
            if (ErrStat >= AbortErrLev) then
               call cleanup()
               return
            end if
         temp_L2(:,:) = 0.0D0
         temp_GL(:) = 0.0D0
         temp_w(:) = 0.0D0
         CALL BD_GaussPointWeight(p%ngp,temp_GL,temp_w,ErrStat2,ErrMsg2)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
            if (ErrStat >= AbortErrLev) then
               call cleanup()
               return
            end if
         DEALLOCATE(temp_w)
     ELSEIF(p%quadrature .EQ. 2) THEN
         CALL AllocAry(temp_L2,6,InputFileData%kp_total,'temp_L2',ErrStat2,ErrMsg2)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         CALL AllocAry(p%Gauss,6,InputFileData%kp_total,'p%Trapezoidal',ErrStat2,ErrMsg2)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         ! Temporary Gauss point intrinsic coordinates array
         temp_L2(:,:) = 0.0D0
     ENDIF
  
     DO i=1,p%elem_total
         IF(i .EQ. 1) THEN
  