[Paraview] fortran routines to Write UNFORMATTED files in Ensight
format for rectilinear grid
Stéphane Montésino
Stephane.Montesino at hmg.inpg.fr
Tue Feb 14 12:06:18 EST 2006
Now it works .... with unformatted Format
--
Stephane MONTESINO
PhD Student tel:+33 4 76 82 52 91
LEGI fax:+33 4 76 82 70 22
BP 53
38041 Grenoble Cedex email: stephane.montesino at hmg.inpg.fr
France http://www.legi.hmg.inpg.fr
-------------- next part --------------
!
! ******************************************************************************
! subrout WriteRectEnsightGeo
! ******************************************************************************
!
! writes mesh data in Ensight's ascii format for rectilinear geometry
!
! n1,n2,n3....: number of nodes in the x1,x2,x3 direction
! x1,x2,x3....: coordinates
!
subroutine WriteRectEnsightGeo(imin,imax,jmin,jmax,kmin,kmax,
& x1,x2,x3,FileName,WriteBinary)
implicit none
INTEGER,INTENT(IN)::imin,imax,jmin,jmax,kmin,kmax
REAL*8,DIMENSION(imax),INTENT(IN)::x1
REAL*8,DIMENSION(jmax),INTENT(IN)::x2
REAL*8,DIMENSION(kmax),INTENT(IN)::x3
LOGICAL ,INTENT(IN)::WriteBinary
CHARACTER(LEN=80) ,INTENT(IN)::FileName
! character(LEN=80)::buffer
character(LEN=80)::binary_form
character(LEN=80)::file_description1,file_description2
character(LEN=80)::node_id,element_id
character(LEN=80)::part,description_part,block
integer::FileUnit,i,j,k,npart,isize,jsize,ksize
integer::reclength
FileUnit = 40
binary_form ='C Binary'
file_description1='Ensight Model Geometry File Created by '
file_description2='WriteRectEnsightGeo Routine'
node_id ='node id off'
element_id ='element id off'
part ='part'
npart =1
description_part ='3D periodic channel'
block ='block rectilinear'
isize=imax-imin+1
jsize=jmax-jmin+1
ksize=kmax-kmin+1
reclength=80*8+4*(4+isize+jsize+ksize)
if (WriteBinary) then
open (unit=FileUnit,file=trim(FileName)//'.geo',
& form='UNFORMATTED',access="direct",recl=reclength)
write(unit=FileUnit,rec=1) binary_form
& ,file_description1
& ,file_description2
& ,node_id
& ,element_id
& ,part,npart
& ,description_part
& ,block
& ,isize,jsize,ksize
& ,(sngl(x1(i)),i=imin,imax)
& ,(sngl(x2(j)),j=jmin,jmax)
& ,(sngl(x3(k)),k=kmin,kmax)
else
open (unit=FileUnit,file=trim(FileName)//'.geo')
write(FileUnit,'(A)') 'Ensight Model Geometry File Created by '
write(FileUnit,'(A)') 'WriteRectEnsightGeo Routine'
write(FileUnit,'(A)') 'node id off'
write(FileUnit,'(A)') 'element id off'
write(FileUnit,'(A)') 'part'
write(FileUnit,'(i10)')npart
write(FileUnit,'(A)') '3D periodic channel'
write(FileUnit,'(A)') 'block rectilinear'
write(FileUnit,'(3i10)') isize,jsize,ksize
write(FileUnit,'(E12.5)') (x1(i),i=imin,imax)
write(FileUnit,'(E12.5)') (x2(j),j=jmin,jmax)
write(FileUnit,'(E12.5)') (x3(k),k=kmin,kmax)
endif
end subroutine
!
! ******************************************************************************
!
! WriteEnsightSca writes result data in Ensight's format
!
! m1,m2,m3....: size of the variable in the x1,x2,x3 direction
! ndv.........: number of dimension of the variable (1=>scalar 3=>vector)
! var.........: data to be written
! Varname.....: word used to build filenames
! imin,imax...: range of writting data in the x1 direction
! jmin,jmax...: range of writting data in the x2 direction
! kmin,kmax...: range of writting data in the x3 direction
!
! ******************************************************************************
subroutine WriteEnsightVar(ndv,m1,m2,m3,var,VarName,WriteBinary,
& imin,imax,jmin,jmax,kmin,kmax)
implicit none
INTEGER ,INTENT(IN)::m1,m2,m3,ndv
INTEGER ,INTENT(IN)::imin,imax,jmin,jmax,kmin,kmax
REAL*8,DIMENSION(ndv,m1,m2,m3),INTENT(IN)::var
CHARACTER*80,INTENT(IN)::Varname
LOGICAL ,INTENT(IN)::WriteBinary
character(len=80):: VarFileName,buffer
character(len=80):: part,block
integer::FileUnit,i,j,k,npart,m,reclength
FileUnit = 40
part ='part'
npart=1
block='block rectilinear'
reclength=80*3+4*(1+(imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)*ndv)
if (ndv.eq.1)VarFileName = trim(Varname)//'.scl'
if (ndv.eq.3)VarFileName = trim(Varname)//'.vec'
! write(*,'(5x,A)') VarFileName
if(WriteBinary) then
open (unit=FileUnit,file=VarFileName,
& form='UNFORMATTED',access="direct",recl=reclength)
write(unit=FileUnit,rec=1) VarFileName
& ,part,npart,block
& ,((((SNGl(var(m,i,j,k))
& ,i=imin,imax)
& ,j=jmin,jmax)
& ,k=kmin,kmax)
& ,m=1,ndv)
else
open (unit=FileUnit,file=VarFileName)
write(buffer,'(a,a)') Varname
write(FileUnit,'(A)') buffer
write(FileUnit,'(A)') 'part'
write(FileUnit,'(I10)')npart
write(FileUnit,'(A)') 'block rectilinear'
do m=1,ndv
write(FileUnit,'(e12.5)')
& (((SNGl(var(m,i,j,k)),i=imin,imax),j=jmin,jmax),k=kmin,kmax)
enddo
endif
close(FileUnit)
end subroutine
!
! ******************************************************************************
! EnsightCase helps to write a Ensight's case file
!
! VarName.....: Name of the variable
! GeoName.....: Name of the geometrie
! VarType.....: 1 => Scalar 3 => Vector
! ntini.......: filename start number
! nstop.......: filename end number
! nprint......: filename increment
!
! nfile.......: number of result files (time steps)
!
subroutine EnsightCase(VarName,GeoName,VarType,ntini,nstop,nprint)
implicit none
INTEGER,INTENT(IN)::VarType,nstop,ntini,nprint
CHARACTER(LEN=80),INTENT(IN)::Varname
CHARACTER(LEN=80),INTENT(IN)::GeoName
integer::FileUnit,i,nfile
write(*,'(/2A)') ' Creating case file for Ensight and Paraview: '
& ,Varname
nfile=(nstop-ntini+1)/nprint
FileUnit = 40
open(FileUnit,file=trim(Varname)//'.case')
write(FileUnit,10) trim(GeoName)//'.geo'
10 format(
&'FORMAT' ,/ ,
&'type: ensight gold',//,
&'GEOMETRY' ,/ ,
&'model: ',A ,//,
&'VARIABLE')
if (nfile.eq.1) then
if(VarType.eq.1)
& write(FileUnit,15)trim(Varname),trim(Varname)//'.scl'
if(VarType.eq.3)
& write(FileUnit,25)trim(Varname),trim(Varname)//'.vec'
else
if(VarType.eq.1)
& write(FileUnit,15)trim(Varname),trim(Varname)//'**********.scl'
if(VarType.eq.3)
& write(FileUnit,25)trim(Varname),trim(Varname)//'**********.vec'
write(FileUnit,45) nfile,ntini,nprint
write(FileUnit,'(f15.3)') (ntini+float(i)*nprint,i=1,nfile)
endif
close(FileUnit)
15 format('scalar per node: ',A,' ', A)
25 format('vector per node: ',A,' ', A)
45 format(
&/,'TIME ' ,
&/,'time set: 1 ' ,
&/,'number of steps:' ,i4 ,
&/,'filename start number:',i10
&/,'filename increment:' ,i4
&/,'time values: '
&)
end subroutine
!
! ******************************************************************************
! **** subrout VTR_XML_3D_scalar ****
! ******************************************************************************
! Export 3D scalar data upon a Non uniform rectilinear Grid
!
! n1m,n2m,n3m..: size of the data and the grid
! y1,y2,y3.....: coordinates
! scal_data....: data to write
! name.........: name of the data
!
subroutine VTR_XML_3D_scalar(n1m,n2m,n3m,y1,y2,y3,scal_data,name)
!
implicit none
!
INTEGER ,INTENT(IN)::n1m,n2m,n3m
REAL*8,DIMENSION(n1m) ,INTENT(IN)::y1
REAL*8,DIMENSION(n2m) ,INTENT(IN)::y2
REAL*8,DIMENSION(n3m) ,INTENT(IN)::y3
REAL*8,DIMENSION(n1m,n2m,n3m),INTENT(IN)::scal_data
CHARACTER(LEN=30) ,INTENT(IN)::name
integer::i,j,k,FileUnit
!
FileUnit=40
open(FileUnit,file=trim(name)//'.vtr',FORM='UNFORMATTED')
write(FileUnit,*)'<VTKFile type="RectilinearGrid" version="0.1"',
& ' byte_order="LittleEndian">'
write(FileUnit,*)' <RectilinearGrid WholeExtent=',
& '"1 ',n1m,' 1 ',n2m,' 1 ',n3m,'">'
write(FileUnit,*)' <Piece Extent=',
& '"1 ',n1m,' 1 ',n2m,' 1 ',n3m,'">'
write(FileUnit,*)' <Coordinates>'
write(FileUnit,*)' <DataArray type="Float32"',
& ' Name="X_COORDINATES"',
& ' NumberOfComponents="1">'
write(FileUnit) (y1(i),i=1,n1m)
write(FileUnit,*)' </DataArray>'
write(FileUnit,*)' <DataArray type="Float32"',
& ' Name="Y_COORDINATES"',
& ' NumberOfComponents="1">'
write(FileUnit) (y2(j),j=1,n2m)
write(FileUnit,*)' </DataArray>'
write(FileUnit,*)' <DataArray type="Float32"',
& ' Name="Z_COORDINATES"',
& ' NumberOfComponents="1">'
write(FileUnit) (y3(k),k=1,n3m)
write(FileUnit,*)' </DataArray>'
write(FileUnit,*)' </Coordinates>'
write(FileUnit,*)' <PointData Scalars="scalar">'
write(FileUnit,*)' <DataArray Name="'//trim(name)//'"',
& ' type="Float32"',
& ' NumberOfComponents="1"',
& ' format="ascii">'
write(FileUnit) (((scal_data(i,j,k),i=1,n1m),j=1,n2m),k=1,n3m)
write(FileUnit,*)' </DataArray>'
write(FileUnit,*)' </PointData>'
write(FileUnit,*)' </Piece>'
write(FileUnit,*)' </RectilinearGrid>'
write(FileUnit,*)'</VTKFile>'
close(FileUnit)
!
return
end
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Stephane.Montesino.vcf
Type: text/x-vcard
Size: 308 bytes
Desc: not available
Url : http://public.kitware.com/pipermail/paraview/attachments/20060214/8a832142/Stephane.Montesino.vcf
More information about the ParaView
mailing list