WriteDataToVTK

public subroutine WriteDataToVTK(dim1, vecdim, nVal, NPlot, nElems, VarNames, Coord, Values, FileString)

Uses

- MODgvec_Globals
Help

Graph Key

Nodes of different colours represent the following:

Solid arrows point from a submodule to the (sub)module which it is descended from. Dashed arrows point from a module or program unit to modules which it uses.

Subroutine to write 3D point data to VTK format

Arguments

Type	Intent		Name
integer,	intent(in)	::	dim1	dimension of the data (either 1:lines,2=quads or 3=hexas)
integer,	intent(in)	::	vecdim	dimension of coordinates
integer,	intent(in)	::	nVal	Number of nodal output variables
integer,	intent(in)	::	NPlot(dim1)	Number of output points per element : (nPlot+1)**dim1
integer,	intent(in)	::	nElems	Number of output elements
character(len=*),	intent(in)	::	VarNames(nVal)	Names of all variables that will be written out
real(kind=wp),	intent(in)	::	Coord(vecdim,1:PRODUCT(Nplot+1),nElems)
real(kind=wp),	intent(in)	::	Values(nVal,1:PRODUCT(Nplot+1),nElems)	Statevector
character(len=*),	intent(in)	::	FileString	Output file name

Calls

Help

Called by

Help

Source Code

SUBROUTINE WriteDataToVTK(dim1,vecDim,nVal,NPlot,nElems,VarNames,Coord,Values,FileString)
! MODULES
USE MODgvec_Globals
! IMPLICIT VARIABLE HANDLING
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
INTEGER,INTENT(IN)            :: dim1                    !! dimension of the data (either 1:lines,2=quads or 3=hexas)
INTEGER,INTENT(IN)            :: vecdim                  !! dimension of coordinates
INTEGER,INTENT(IN)            :: nVal                    !! Number of nodal output variables
INTEGER,INTENT(IN)            :: NPlot(dim1)             !! Number of output points per element : (nPlot+1)**dim1
INTEGER,INTENT(IN)            :: nElems                  !! Number of output elements
REAL(wp),INTENT(IN)           :: Coord(vecdim,1:PRODUCT(Nplot+1),nElems)      ! CoordinatesVector
CHARACTER(LEN=*),INTENT(IN)   :: VarNames(nVal)          !! Names of all variables that will be written out
REAL(wp),INTENT(IN)           :: Values(nVal,1:PRODUCT(Nplot+1),nElems)   !! Statevector
CHARACTER(LEN=*),INTENT(IN)   :: FileString              !! Output file name
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
INTEGER,PARAMETER     :: kindFloat=8  !set floating point accuracy single (4) double (8), should be equal or lower than input data!
REAL(KIND=kindFloat)  :: FLOATdummy
CHARACTER(LEN=7)      :: strfloat
INTEGER               :: INTdummy
INTEGER               :: sizefloat,sizeInt
INTEGER            :: i,j,k,iVal,iElem,Offset,nBytes,nVTKPoints,nVTKCells,ivtk
INTEGER            :: Vertex(2**dim1,PRODUCT(Nplot)*nElems)  ! ?
INTEGER            :: ProdNplot,ProdNplot_p1,NPlot_p1(dim1),CellID,PointID,ElemType  ! ?
CHARACTER(LEN=35)  :: StrOffset,TempStr1,TempStr2  ! ?
CHARACTER(LEN=300) :: Buffer
CHARACTER(LEN=255) :: tmpVarName,tmpVarNameY,tmpVarNameZ
INTEGER            :: StrLen,iValVec,nValVec,VecOffset(0:nVal)
LOGICAL            :: isVector,maybeVector
CHARACTER(LEN=1)   :: strvecdim
CHARACTER(LEN=1)   :: lf
!===================================================================================================================================
WRITE(UNIT_stdOut,'(A)',ADVANCE='NO')'   WRITE DATA TO VTX XML BINARY (VTU) FILE "'//TRIM(FileString)//'" ...'
ivtk=GETFREEUNIT()
NPlot_p1  =(Nplot(:)+1)
ProdNPlot  =PRODUCT(Nplot(:))
ProdNPlot_p1  =PRODUCT(Nplot_p1(:))

IF(kindFloat.EQ.4) THEN
  strfloat='Float32'
ELSEIF(kindFloat.EQ.8)THEN
  strfloat='Float64'
ELSE
  STOP 'kindFloat not implemented in output vtk'
END IF
sizefloat=SIZEOF_F(FLOATdummy)
sizeInt  =SIZEOF_F(INTdummy)


IF(vecdim.LT.dim1) THEN
  WRITE(*,*)'WARNING:data dimension should be <= vecdim! dim1= ',dim1,' vecdim= ',vecdim
  STOP
END IF
! Line feed character
lf = char(10)
WRITE(strvecdim,'(I1)') vecdim

! Write file
OPEN(UNIT=ivtk,FILE=TRIM(FileString),ACCESS='STREAM')
! Write header
Buffer='<?xml version="1.0"?>'//lf;WRITE(ivtk) TRIM(Buffer)
Buffer='<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">'//lf;WRITE(ivtk) TRIM(Buffer)
! Specify file type
nVTKPoints= ProdNPlot_p1*nElems
nVTKCells = ProdNPlot   *nElems
Buffer='  <UnstructuredGrid>'//lf;WRITE(ivtk) TRIM(Buffer)
WRITE(TempStr1,'(I16)')nVTKPoints
WRITE(TempStr2,'(I16)')nVTKCells
Buffer='    <Piece NumberOfPoints="'//TRIM(ADJUSTL(TempStr1))//'" &
       &NumberOfCells="'//TRIM(ADJUSTL(TempStr2))//'">'//lf;WRITE(ivtk) TRIM(Buffer)
! Specify point data
Buffer='      <PointData>'//lf;WRITE(ivtk) TRIM(Buffer)
Offset=0
WRITE(StrOffset,'(I16)')Offset
!accout for vectors:
! if Variable Name ends with an X and the following have the same name with Y and Z
! then it forms a vector variable (X is omitted for the name)

iVal=0 !scalars
iValVec=0 !scalars & vectors
VecOffset(0)=0
DO WHILE(iVal.LT.nVal)
  iVal=iVal+1
  iValVec=iValVec+1
  tmpVarName=TRIM(VarNames(iVal))
  StrLen=LEN(TRIM(tmpVarName))
  maybeVector=(iVal+vecdim-1.LE.nVal)
  isVector=.FALSE.
  IF(maybeVector)THEN
    SELECT CASE(vecdim)
    CASE(2)
      tmpVarNameY=TRIM(VarNames(iVal+1))
      isVector=((iVal+2.LE.nVal).AND.(INDEX(tmpVarName( StrLen:StrLen),"X").NE.0) &
                                .AND.(INDEX(tmpVarNameY(:StrLen),TRIM(tmpVarName(:StrLen-1))//"Y").NE.0))
    CASE(3)
      tmpVarNameY=TRIM(VarNames(iVal+1))
      tmpVarNameZ=TRIM(VarNames(iVal+2))
      isVector=((iVal+2.LE.nVal).AND.(INDEX(tmpVarName( StrLen:StrLen),"X").NE.0) &
                                .AND.(INDEX(tmpVarNameY(:StrLen),TRIM(tmpVarName(:StrLen-1))//"Y").NE.0) &
                                .AND.(INDEX(tmpVarNameZ(:StrLen),TRIM(tmpVarName(:StrLen-1))//"Z").NE.0))

    END SELECT
  END IF !maybevector

  IF(isvector)THEN !variable is a vector!
    tmpVarName=tmpVarName(:StrLen-1)
    Buffer='        <DataArray type="'//strfloat//'" Name="'//TRIM(tmpVarName)//'" NumberOfComponents="'//strvecdim// &
           &'" format="appended" offset="'//TRIM(ADJUSTL(StrOffset))//'"/>'//lf;WRITE(ivtk) TRIM(Buffer)
    Offset=Offset+sizeInt+vecdim*nVTKPoints*sizefloat
    WRITE(StrOffset,'(I16)')Offset
    VecOffset(iValVec)=VecOffset(iValVec-1)+vecdim
    iVal=iVal+vecdim-1 !skip the Y (& Z) components
  ELSE
    Buffer='        <DataArray type="'//strfloat//'" Name="'//TRIM(tmpVarName)// &
           &'" format="appended" offset="'//TRIM(ADJUSTL(StrOffset))//'"/>'//lf;WRITE(ivtk) TRIM(Buffer)
    Offset=Offset+sizeInt+nVTKPoints*sizeFloat
    WRITE(StrOffset,'(I16)')Offset
    VecOffset(iValVec)=VecOffset(iValVec-1)+1
  END IF !isvector
END DO !iVal <=nVal
nValVec=iValVec

Buffer='      </PointData>'//lf;WRITE(ivtk) TRIM(Buffer)
! Specify cell data
Buffer='      <CellData> </CellData>'//lf;WRITE(ivtk) TRIM(Buffer)
! Specify coordinate data
Buffer='      <Points>'//lf;WRITE(ivtk) TRIM(Buffer)
Buffer='        <DataArray type="'//strfloat//'" Name="Coordinates" NumberOfComponents="'//strvecdim// &
'" format="appended" offset="'//TRIM(ADJUSTL(StrOffset))//'"/>'//lf;WRITE(ivtk) TRIM(Buffer)
Offset=Offset+sizeInt+vecdim*nVTKPoints*sizeFloat
WRITE(StrOffset,'(I16)')Offset
Buffer='      </Points>'//lf;WRITE(ivtk) TRIM(Buffer)
! Specify necessary cell data
Buffer='      <Cells>'//lf;WRITE(ivtk) TRIM(Buffer)
! Connectivity
Buffer='        <DataArray type="Int32" Name="connectivity" format="appended" &
         &offset="'//TRIM(ADJUSTL(StrOffset))//'"/>'//lf;WRITE(ivtk) TRIM(Buffer)
Offset=Offset+sizeInt+2**dim1*nVTKCells*sizeInt
WRITE(StrOffset,'(I16)')Offset
! Offset in connectivity data
Buffer='        <DataArray type="Int32" Name="offsets" format="appended" &
         &offset="'//TRIM(ADJUSTL(StrOffset))//'"/>'//lf;WRITE(ivtk) TRIM(Buffer)
Offset=Offset+sizeInt+nVTKCells*sizeInt
WRITE(StrOffset,'(I16)')Offset
! Elem types
Buffer='        <DataArray type="Int32" Name="types" format="appended" &
         &offset="'//TRIM(ADJUSTL(StrOffset))//'"/>'//lf;WRITE(ivtk) TRIM(Buffer)
Buffer='      </Cells>'//lf;WRITE(ivtk) TRIM(Buffer)
Buffer='    </Piece>'//lf;WRITE(ivtk) TRIM(Buffer)
Buffer='  </UnstructuredGrid>'//lf;WRITE(ivtk) TRIM(Buffer)
! Prepare append section
Buffer='  <AppendedData encoding="raw">'//lf;WRITE(ivtk) TRIM(Buffer)
! Write leading data underscore
Buffer='_';WRITE(ivtk) TRIM(Buffer)

! Write binary raw data into append section
! Point data
nBytes = nVTKPoints*sizeFloat
DO iValVec=1,nValVec
  WRITE(ivtk) (vecOffset(iValVec)-vecOffset(iValVec-1))*nBytes
  WRITE(ivtk) REAL(Values(VecOffSet(iValVec-1)+1:VecOffset(iValVec),:,:),kindFloat)
END DO !iValVec
! Point coordinates
nBytes = nVTKPoints * vecdim*sizeFloat
WRITE(ivtk) nBytes
WRITE(ivtk) REAL(Coord(:,:,:),kindFloat)
! Connectivity
SELECT CASE(dim1)
CASE(1)
  CellID = 0
  PointID= 0
  DO iElem=1,nElems
    DO i=1,NPlot(1)
      CellID = CellID+1
      !visuLineElem
      Vertex(:,CellID) = (/ PointID+(i-1), PointID+ i /)
    END DO
    PointID=PointID+NPlot(1)
  END DO
CASE(2)
  CellID = 0
  PointID= 0
  DO iElem=1,nElems
    DO j=1,NPlot(2)
      DO i=1,NPlot(1)
        CellID = CellID+1
        !visuQuadElem
        Vertex(:,CellID) = (/                  &
          PointID+(i-1)+  j   * NPlot_p1(1) ,    & !P4
          PointID+(i-1)+ (j-1)* NPlot_p1(1) ,    & !P1(CGNS=tecplot standard)
          PointID+ i   + (j-1)* NPlot_p1(1) ,    & !P2
          PointID+ i   +  j   * NPlot_p1(1)     /) !P3
      END DO
    END DO
    PointID=PointID+ProdNPlot_p1
  END DO
CASE(3)
  CellID=0
  PointID=0
  DO iElem=1,nElems
    DO k=1,NPlot(3)
      DO j=1,NPlot(2)
        DO i=1,NPlot(1)
          CellID=CellID+1
          !
          Vertex(:,CellID)=(/                                       &
            PointID+(i-1)+( j   +(k-1)*NPlot_p1(2))*NPlot_p1(1),      & !P4(CGNS=tecplot standard)
            PointID+(i-1)+((j-1)+(k-1)*NPlot_p1(2))*NPlot_p1(1),      & !P1
            PointID+ i   +((j-1)+(k-1)*NPlot_p1(2))*NPlot_p1(1),      & !P2
            PointID+ i   +( j   +(k-1)*NPlot_p1(2))*NPlot_p1(1),      & !P3
            PointID+(i-1)+( j   + k   *NPlot_p1(2))*NPlot_p1(1),      & !P8
            PointID+(i-1)+((j-1)+ k   *NPlot_p1(2))*NPlot_p1(1),      & !P5
            PointID+ i   +((j-1)+ k   *NPlot_p1(2))*NPlot_p1(1),      & !P6
            PointID+ i   +( j   + k   *NPlot_p1(2))*NPlot_p1(1)      /) !P7
        END DO
      END DO
    END DO
    !
    PointID=PointID+ProdNPlot_p1
  END DO
END SELECT
nBytes = 2**dim1*nVTKCells*sizeInt
WRITE(ivtk) nBytes
WRITE(ivtk) Vertex(:,:)
! Offset in connectivity
nBytes = nVTKCells*sizeInt
WRITE(ivtk) nBytes
WRITE(ivtk) (Offset,Offset=2**dim1,2**dim1*nVTKCells,2**dim1)
! Elem type
SELECT CASE(dim1)
CASE(1)
  ElemType =3 !VTK_LINE
CASE(2)
  ElemType =9 !VTK_QUAD
CASE(3)
  ElemType =12  !VTK_HEXAHEDRON
END SELECT
WRITE(ivtk) nBytes
WRITE(ivtk) (ElemType,iElem=1,nVTKCells)
! Write footer
Buffer=lf//'  </AppendedData>'//lf;WRITE(ivtk) TRIM(Buffer)
Buffer='</VTKFile>'//lf;WRITE(ivtk) TRIM(Buffer)
CLOSE(ivtk)
WRITE(UNIT_stdOut,'(A)',ADVANCE='YES')"   DONE"
END SUBROUTINE WriteDataToVTK

WriteDataToVTK Subroutine

Contents

Source Code

public subroutine WriteDataToVTK(dim1, vecdim, nVal, NPlot, nElems, VarNames, Coord, Values, FileString)

Uses

Arguments

Calls

Called by

Source Code