anchor.diff:  !**********************************************************************************************************************************
   ! LICENSING
  -! Copyright (C) 2013-2015  National Renewable Energy Laboratory
  +! Copyright (C) 2013-2016  National Renewable Energy Laboratory
   !
   !    This file is part of the NWTC Subroutine Library.
   !
  
positive.code: !----------------------------------------------------------------------------------------------------------------------------------
  !> This routine writes mesh information in VTK format.
  !! see VTK file information format for XML, here: http://www.vtk.org/wp-content/uploads/2015/04/file-formats.pdf
  SUBROUTINE MeshWrVTK ( M, FileRootName, VTKcount, ErrStat, ErrMsg )
        
     TYPE(MeshType),  INTENT(IN)           :: M             !< mesh to be written
     CHARACTER(*),    INTENT(IN)           :: FileRootName  !< Name of the file to write the output in (excluding extension)
     INTEGER(IntKi),  INTENT(IN)           :: VTKcount      !< Indicates number for VTK output file (when 0, the routine will also write reference information)
  
     INTEGER(IntKi),  INTENT(OUT)          :: ErrStat       !< Indicates whether an error occurred (see NWTC_Library)
     CHARACTER(*),    INTENT(OUT)          :: ErrMsg        !< Error message associated with the ErrStat
  
     ! local variables
     INTEGER(IntKi)                        :: Un            ! fortran unit number
     INTEGER(IntKi)                        :: i,j           ! loop counters
     CHARACTER(1024)                       :: FileName
        
     INTEGER(IntKi)                        :: ErrStat2 
     CHARACTER(ErrMsgLen)                  :: ErrMsg2
  
negative.code:                
        
     END SUBROUTINE MeshWrVTK
  !----------------------------------------------------------------------------------------------------------------------------------
  !> This routine writes line2 mesh surface information in VTK format.
  !! see VTK file information format for XML, here: http://www.vtk.org/wp-content/uploads/2015/04/file-formats.pdf
  SUBROUTINE MeshWrVTK_Ln2Surface ( M, FileRootName, VTKcount, NumSegments, Radius, ErrStat, ErrMsg )
        
     TYPE(MeshType),  INTENT(IN)           :: M             !< mesh to be written
     CHARACTER(*),    INTENT(IN)           :: FileRootName  !< Name of the file to write the output in (excluding extension)
     INTEGER(IntKi),  INTENT(IN)           :: VTKcount      !< Indicates number for VTK output file (when 0, the routine will also write reference information)
     INTEGER(IntKi),  INTENT(IN)           :: NumSegments   !< Number of segments to split the circle into
     REAL(SiKi),      INTENT(IN)           :: Radius(:)     !< Radius of each node
     
     INTEGER(IntKi),  INTENT(OUT)          :: ErrStat       !< Indicates whether an error occurred (see NWTC_Library)
     CHARACTER(*),    INTENT(OUT)          :: ErrMsg        !< Error message associated with the ErrStat
  
  
     ! local variables
     INTEGER(IntKi)                        :: Un            ! fortran unit number
     INTEGER(IntKi)                        :: i,j           ! loop counters
     CHARACTER(1024)                       :: FileName
     REAL(SiKi)                            :: angle
     REAL(SiKi)                            :: delX
     REAL(SiKi)                            :: delY
  
     INTEGER(IntKi)                        :: firstPntEnd, firstPntStart, secondPntStart, secondPntEnd  ! node indices for forming rectangle 
     
     
     
     INTEGER(IntKi)                        :: ErrStat2 
     CHARACTER(ErrMsgLen)                  :: ErrMsg2
     CHARACTER(*),PARAMETER                :: RoutineName = 'MeshWrVTK_Ln2Surface'
  
     
     ErrStat = ErrID_None
     ErrMsg  = ""
  
     IF (.NOT. M%Initialized) RETURN
     IF (.NOT. ALLOCATED(M%TranslationDisp) ) RETURN
     IF (.NOT. ALLOCATED(M%Orientation) ) RETURN
        
     !.................................................................
     ! write the data that potentially changes each time step:
     !.................................................................
        
     ! PolyData (.vtp)  Serial vtkPolyData (unstructured) file
     FileName = TRIM(FileRootName)//'.t'//TRIM(Num2LStr(VTKcount))//'.vtp'
        
     CALL GetNewUnit( Un, ErrStat2, ErrMsg2 )      
     CALL OpenFOutFile ( Un, TRIM(FileName), ErrStat2, ErrMsg2 )
        call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
        if (ErrStat >= AbortErrLev) return
           
        ! Write a VTP mesh file (Polygonal VTK file) with positions and polygons (surfaces)
        ! (note alignment of WRITE statements to make sure spaces are lined up in XML file)
        WRITE(Un,'(A)')         '<?xml version="1.0"?>'
        WRITE(Un,'(A)')         '<VTKFile type="PolyData" version="0.1" byte_order="LittleEndian">' ! bjj note: we don't have binary data in this file, so byte_order shouldn't matter, right?
        WRITE(Un,'(A)')         '  <PolyData>'
        WRITE(Un,'(2(A,i7),A)') '    <Piece NumberOfPoints="', M%Nnodes*NumSegments, '" NumberOfVerts="  0" NumberOfLines="', M%ElemTable(ELEMENT_LINE2)%nelem, '"'
        WRITE(Un,'(A,i7,A)')    '           NumberOfStrips="  0" NumberOfPolys="',  M%ElemTable(ELEMENT_LINE2)%nelem*NumSegments, '">'
           
  ! points (nodes, augmented with NumSegments):   
        WRITE(Un,'(A)')         '      <Points>'
        WRITE(Un,'(A)')         '        <DataArray type="Float32" NumberOfComponents="3" format="ascii">'
        
        DO i=1,M%Nnodes
           DO j=1,NumSegments
              angle = TwoPi*(j-1.0_ReKi)/NumSegments
              delX = radius(i)*COS(angle)
              delY = radius(i)*SIN(angle)
              WRITE(Un,VTK_AryFmt) M%Position(:,i) + M%TranslationDisp(:,i) + delX*M%Orientation(1,:,i) + delY*M%Orientation(2,:,i)
           END DO
        END DO
     
        WRITE(Un,'(A)')         '        </DataArray>'
        WRITE(Un,'(A)')         '      </Points>'
    
  
     if ( M%ElemTable(ELEMENT_LINE2)%nelem > 0) then   
           ! Using rectangle to render surfaces (for line2 meshes only):
        WRITE(Un,'(A)')         '      <Polys>'
        
        WRITE(Un,'(A)')         '        <DataArray type="Int32" Name="connectivity" format="ascii">'      
        DO i=1,M%ElemTable(ELEMENT_LINE2)%nelem
           firstPntStart  = (M%ElemTable(ELEMENT_LINE2)%Elements(i)%ElemNodes(1)-1)*NumSegments
           firstPntEnd    = firstPntStart + NumSegments - 1
           secondPntStart = (M%ElemTable(ELEMENT_LINE2)%Elements(i)%ElemNodes(2)-1)*NumSegments
           secondPntEnd   = secondPntStart + NumSegments - 1
           DO j=1,NumSegments-1
              WRITE(Un,'(4(i7))') firstPntStart + (j-1), firstPntStart + j, &
                                  secondPntStart + j,    secondPntStart + (j-1)
           END DO
           WRITE(Un,'(4(i7))')  firstPntEnd, firstPntStart, secondPntStart, secondPntEnd
        END DO      
        WRITE(Un,'(A)')         '        </DataArray>'
  
        WRITE(Un,'(A)')         '        <DataArray type="Int32" Name="offsets" format="ascii">'
        DO i=1,M%ElemTable(ELEMENT_LINE2)%nelem
           DO j=1,NumSegments
              WRITE(Un,'(i7)') 4*NumSegments*(i - 1) + 4*j
           END DO
        END DO
        WRITE(Un,'(A)')         '        </DataArray>'
  
        WRITE(Un,'(A)')         '      </Polys>'      
              
     end if ! do this only for line2 elements    
  
        WRITE(Un,'(A)')         '    </Piece>'
        WRITE(Un,'(A)')         '  </PolyData>'
        WRITE(Un,'(A)')         '</VTKFile>'
        CLOSE(Un)         
                       
     END SUBROUTINE MeshWrVTK_Ln2Surface
  !----------------------------------------------------------------------------------------------------------------------------------
  !> This routine writes point mesh surfaces information in VTK format.
  !! see VTK file information format for XML, here: http://www.vtk.org/wp-content/uploads/2015/04/file-formats.pdf
  SUBROUTINE MeshWrVTK_PointSurface ( M, FileRootName, VTKcount, ErrStat, ErrMsg, NumSegments, Radius, verts )
        
     TYPE(MeshType),  INTENT(IN)           :: M             !< mesh to be written
     CHARACTER(*),    INTENT(IN)           :: FileRootName  !< Name of the file to write the output in (excluding extension)
     INTEGER(IntKi),  INTENT(IN)           :: VTKcount      !< Indicates number for VTK output file (when 0, the routine will also write reference information)
     INTEGER(IntKi),  INTENT(IN), OPTIONAL :: NumSegments   !< Number of segments to split the circle into
     REAL(SiKi),      INTENT(IN), OPTIONAL :: Radius        !< Radius of each node
     REAL(SiKi),      INTENT(IN), OPTIONAL :: verts(:,:)    !< X-Y-Z verticies (3xn) of points that define a volume around each node
     !bjj: we don't need this to be limited to 8, I guess...
     
     INTEGER(IntKi),  INTENT(OUT)          :: ErrStat       !< Indicates whether an error occurred (see NWTC_Library)
     CHARACTER(*),    INTENT(OUT)          :: ErrMsg        !< Error message associated with the ErrStat
  
  
     ! local variables
     INTEGER(IntKi)                        :: Un            ! fortran unit number
     INTEGER(IntKi)                        :: i,j,k         ! loop counters
     INTEGER(IntKi)                        :: NumberOfPoints
     INTEGER(IntKi)                        :: NumberOfPolys
     INTEGER(IntKi)                        :: NumSegments1   
     INTEGER(IntKi)                        :: NumSegments2   
     CHARACTER(1024)                       :: FileName
     REAL(SiKi)                            :: angle, r, ratio
     REAL(SiKi)                            :: xyz(3)
  
     INTEGER(IntKi)                        :: firstPntEnd, firstPntStart, secondPntStart, secondPntEnd  ! node indices for forming rectangle 
     
     
     
     INTEGER(IntKi)                        :: ErrStat2 
     CHARACTER(ErrMsgLen)                  :: ErrMsg2
     CHARACTER(*),PARAMETER                :: RoutineName = 'MeshWrVTK_PointSurface'
  
     
     ErrStat = ErrID_None
     ErrMsg  = ""
  
     IF (.NOT. M%Initialized) RETURN
     IF (.NOT. ALLOCATED(M%TranslationDisp) ) RETURN
     IF (.NOT. ALLOCATED(M%Orientation) ) RETURN
        
     if (present(verts)) then
        NumberOfPoints = size(verts,2)
        if (size(verts,2)==8) then
           NumberOfPolys = 6 ! per node
        elseif (size(verts,2)==4) then
           NumberOfPolys = 1 ! per node
        else
           ! it would be nice if we could add this sometime, but ...
           call SetErrStat(ErrID_Fatal,'When verticies are specified, there must be exactly 8.',ErrStat,ErrMsg,RoutineName)
           RETURN
        end if
        
        NumberOfPoints = M%Nnodes * NumberOfPoints
        NumberOfPolys  = M%Nnodes * NumberOfPolys
        
     elseif (present(Radius) ) then
        if (present(NumSegments)) then ! a volume
           NumSegments1 = NumSegments
           NumSegments2 = max(4,NumSegments1)
        else                           ! a plane
           NumSegments1 = 20
           NumSegments2 = 1
        end if
        
        NumberOfPoints = M%Nnodes *  NumSegments2   * NumSegments1
        NumberOfPolys  = M%Nnodes * (NumSegments2-1)* NumSegments1
     else
        call SetErrStat(ErrID_Fatal,'Incorrect number of arguments.',ErrStat,ErrMsg,RoutineName)
        RETURN
     end if
     
     !.................................................................
     ! write the data that potentially changes each time step:
     !.................................................................
        
     ! PolyData (.vtp)  Serial vtkPolyData (unstructured) file
     FileName = TRIM(FileRootName)//'.t'//TRIM(Num2LStr(VTKcount))//'.vtp'
        
     CALL GetNewUnit( Un, ErrStat2, ErrMsg2 )      
     CALL OpenFOutFile ( Un, TRIM(FileName), ErrStat2, ErrMsg2 )
        call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
        if (ErrStat >= AbortErrLev) return
           
        ! Write a VTP mesh file (Polygonal VTK file) with positions and polygons (surfaces)
        ! (note alignment of WRITE statements to make sure spaces are lined up in XML file)
        WRITE(Un,'(A)')         '<?xml version="1.0"?>'
        WRITE(Un,'(A)')         '<VTKFile type="PolyData" version="0.1" byte_order="LittleEndian">' ! bjj note: we don't have binary data in this file, so byte_order shouldn't matter, right?
        WRITE(Un,'(A)')         '  <PolyData>'
        WRITE(Un,'(2(A,i7),A)') '    <Piece NumberOfPoints="', NumberOfPoints, '" NumberOfVerts="  0" NumberOfLines="', 0, '"'
        WRITE(Un,'(A,i7,A)')    '           NumberOfStrips="  0" NumberOfPolys="',  NumberOfPolys, '">'
           
  ! points (nodes, augmented with NumSegments):   
        WRITE(Un,'(A)')         '      <Points>'
        WRITE(Un,'(A)')         '        <DataArray type="Float32" NumberOfComponents="3" format="ascii">'
        
        if ( present(verts) ) then
           do i=1,M%Nnodes
              do j=1,NumberOfPoints
                 WRITE(Un,VTK_AryFmt) M%Position(:,i) + M%TranslationDisp(:,i) + MATMUL(M%Orientation(:,:,i),verts(:,j))
              end do    
           end do
           
        else                           
           DO i=1,M%Nnodes
              DO j=1,NumSegments2
                 ratio  = 2.0_SiKi*REAL(j-1,SiKi)/REAL(NumSegments2-1,SiKi) - 1.0_SiKi  ! where we are in [-1, 1]
                 xyz(1) = radius*ratio
                 
                 ! now calculate the radius of the y-z circle we're going to create at this x:              
                 !r = acos( radius / xyz(1) ) or r = sqrt( radius**2 - xyz(1)**2 ) = radius*sqrt(1 - ratio**2)
                 r = radius*sqrt(abs(1.0_SiKi - ratio**2))  ! note the abs in case ratio**2 gets slightly larger than 1 
                 
                 DO k=1,NumSegments1
                    angle = TwoPi*(k-1.0_ReKi)/NumSegments1
                    xyz(2) = r*COS(angle)
                    xyz(3) = r*SIN(angle)
                    !WRITE(Un,VTK_AryFmt) M%Position(:,i) + M%TranslationDisp(:,i) + MATMUL(M%Orientation(:,:,i),xyz)
                    WRITE(Un,VTK_AryFmt) M%Position(:,i) + M%TranslationDisp(:,i) + xyz
                 END DO            
              END DO
           END DO         
        end if
        
        WRITE(Un,'(A)')         '        </DataArray>'
        WRITE(Un,'(A)')         '      </Points>'
    
                    
        WRITE(Un,'(A)')         '      <Polys>'
        
        WRITE(Un,'(A)')         '        <DataArray type="Int32" Name="connectivity" format="ascii">'   
        
        if ( present(verts) ) then
           if (size(verts,2)==8) then
              
                 ! Write points for the 6 corners of the box
              
              do i=1,M%Nnodes
                 firstPntStart  = (i-1)*size(verts,2)
  
                 WRITE(Un,'(4(i7))') firstPntStart  ,firstPntStart+1,firstPntStart+2,firstPntStart+3    ! bottom
                 WRITE(Un,'(4(i7))') firstPntStart+4,firstPntStart+5,firstPntStart+6,firstPntStart+7    ! top
                 
                 WRITE(Un,'(4(i7))') firstPntStart+7,firstPntStart+6,firstPntStart+1,firstPntStart      ! sides
                 WRITE(Un,'(4(i7))') firstPntStart+6,firstPntStart+5,firstPntStart+2,firstPntStart+1
                 WRITE(Un,'(4(i7))') firstPntStart+5,firstPntStart+4,firstPntStart+3,firstPntStart+2
                 WRITE(Un,'(4(i7))') firstPntStart+4,firstPntStart+7,firstPntStart  ,firstPntStart+3                  
              end do
              
           elseif (size(verts,2)==4) then
              
                 ! Write points for the 4 corners of the polygon
              
              do i=1,M%Nnodes
                 firstPntStart  = (i-1)*size(verts,2)
                 
                 WRITE(Un,'(4(i7))') firstPntStart  ,firstPntStart+1,firstPntStart+2,firstPntStart+3    ! bottom
              end do
              
           end if
        else
              
           DO i=1,M%Nnodes
              DO j=1,NumSegments2-1
                 firstPntStart  = (j-1)*NumSegments1 + (i-1)*NumSegments1*NumSegments2
                 firstPntEnd    = firstPntStart + NumSegments1 - 1
                 secondPntStart = j*NumSegments1 + (i-1)*NumSegments1*NumSegments2
                 secondPntEnd   = secondPntStart + NumSegments1 - 1
                 DO k=1,NumSegments1-1
                    WRITE(Un,'(4(i7))') firstPntStart + (k-1), firstPntStart + k, &
                                        secondPntStart + k,    secondPntStart + (k-1)
                 END DO
                 WRITE(Un,'(4(i7))')  firstPntEnd, firstPntStart, secondPntStart, secondPntEnd
              END DO         
           END DO      
           
        end if
        
        WRITE(Un,'(A)')         '        </DataArray>'
        
        WRITE(Un,'(A)')         '        <DataArray type="Int32" Name="offsets" format="ascii">'      
        do i=1,NumberOfPolys
           WRITE(Un,'(i7)') 4*i
        end do      
        WRITE(Un,'(A)')         '        </DataArray>'
        WRITE(Un,'(A)')         '      </Polys>'      
              
  
        WRITE(Un,'(A)')         '    </Piece>'
        WRITE(Un,'(A)')         '  </PolyData>'
        WRITE(Un,'(A)')         '</VTKFile>'
        CLOSE(Un)         
                       
     END SUBROUTINE MeshWrVTK_PointSurface
     
  !-------------------------------------------------------------------------------------------------------------------------------
  !> This routine writes mesh information in text form. It is used for debugging.
  