JGonzalez/fpakc
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Merge branch 'feature/vtk' into 'development'

Browse source 
Implementation of vtu file format for mesh reading and output

See merge request JorgeGonz/fpakc!38
This commit is contained in:
Jorge Gonzalez 2023-02-07 15:11:16 +00:00
commit de2aad02d1
14 changed files with 1148 additions and 60 deletions
Download patch file Download diff file Expand all files Collapse all files

BIN
doc/user-manual/fpakc_UserManual.pdf
View file

Binary file not shown.

19
doc/user-manual/fpakc_UserManual.tex
View file

@ -38,6 +38,7 @@
\newacronym{io}{I/O}{input/output}
\newacronym{mcc}{MCC}{Monte-Carlo Collisions}
\newacronym{cs}{CS}{Coulomb Scattering}
\newacronym{vtu}{VTU}{Visualization Toolkit for unstructured grids}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newglossaryentry{openmp}{name={OpenMP},description={Shared-memory parallelization}}
@ -437,25 +438,27 @@ make
\begin{itemize}
\item \textbf{Cart}: Cartesian coordinates.
Available for \textbf{geometry.dimension} $1$, $2$ and $3$.
For \Gls{gmsh} mesh format, the coordinates $x$, $y$ and $z$ correspond to $x$, $y$ and $z$ respectively.
The coordinates $x$, $y$ and $z$ correspond to $x$, $y$ and $z$ respectively.
\item \textbf{Cyl}: Cylindrical coordinates ($z \hyphen r$) with symmetry axis at $r = 0$.
Only available for \textbf{geometry.dimension} $2$.
For \Gls{gmsh} mesh format, the coordinates $x$ and $y$ correspond to $z$ and $r$ respectively.
The coordinates $x$ and $y$ correspond to $z$ and $r$ respectively.
\item \textbf{Rad}: One-dimensional radial space ($r$).
Only available for \textbf{geometry.dimension} $1$.
For \Gls{gmsh} mesh format, the coordinates $x$ corresponds to $r$.
The coordinates $x$ corresponds to $r$.
\end{itemize}
\item \textbf{meshType}: Character.
Format of mesh file.
The output will be written in the same format as specified here.
Accepted formats are:
\begin{itemize}
\item \textbf{gmsh2}: \Gls{gmsh} file format in version 2.0.
\item \textbf{gmsh2}: \Gls{gmsh} file format in version 2.0. This has to be in ASCII format.
\item \textbf{vtu}: \Gls{vtu} file format. This has to be in ASCII format.
\end{itemize}
\item \textbf{meshFile}: Character.
Mesh filename.
This file is searched in the path \textbf{output.path} and must contain the file extension.
\item \textbf{volume}: Real.
Units of $\unit{m^-3}$.
Units of $\unit{m^{-3}}$.
Used to set a fictitious volume for the $0$ dimension.
Ignored in the other cases.
\end{itemize}
@ -536,7 +539,7 @@ make
Density of neutral background.
Required parameter.
\item \textbf{velocity}: Real.
Units in $\unit{m^{-3}}$.
Units in $\unit{m s^{-1}}$.
Array of dimension $3$.
Mean velocity of neutral background.
Required parameter.
@ -604,7 +607,8 @@ make
\item \textbf{part/s}: Particles (real) per second.
\end{itemize}
\item \textbf{v}: Real.
Module of velocity vector, in $\unitfrac{m}{s}$.
Units of $\unit{m s^{-1}}$.
Module of velocity vector.
\item \textbf{n}: Real.
Array dimension $3$.
Direction of injection.
@ -618,6 +622,7 @@ make
\item \textbf{Delta}: Dirac's delta distribution function. All particles are injected with velocity \textbf{v} times the value of \textbf{n} in the specified direction.
\end{itemize}
\item \textbf{T}: Real.
Units of $\unit{K}$
Array dimension $3$.
Temperature in each direction.
\item \textbf{physicalSurface}: Integer.

2
src/makefile
View file

@ -5,6 +5,8 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module
$(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \
$(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o $(OBJDIR)/moduleMath.o \
$(OBJDIR)/moduleProbe.o $(OBJDIR)/moduleAverage.o \
$(OBJDIR)/moduleMeshInoutCommon.o \
$(OBJDIR)/moduleMeshInputVTU.o $(OBJDIR)/moduleMeshOutputVTU.o \
$(OBJDIR)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \
$(OBJDIR)/moduleMeshInput0D.o $(OBJDIR)/moduleMeshOutput0D.o \
$(OBJDIR)/moduleMesh3DCart.o \

51
src/modules/init/moduleInput.f90
View file

@ -860,7 +860,13 @@ MODULE moduleInput
SUBROUTINE readGeometry(config)
USE moduleMesh
USE moduleMeshInputGmsh2, ONLY: initGmsh2
USE moduleMeshInputVTU, ONLY: initVTU
USE moduleMeshInput0D, ONLY: init0D
USE moduleMesh3DCart
USE moduleMesh2DCyl
USE moduleMesh2DCart
USE moduleMesh1DRad
USE moduleMesh1DCart
USE moduleErrors
USE moduleOutput
USE moduleRefParam
@ -873,6 +879,7 @@ MODULE moduleInput
LOGICAL:: found
CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile
REAL(8):: volume
CHARACTER(:), ALLOCATABLE:: meshFileVTU !Temporary to test VTU OUTPUT
object = 'geometry'
@ -957,7 +964,39 @@ MODULE moduleInput
END SELECT
!Get the type of mesh
!Link the procedure to connect meshes
SELECT CASE(mesh%dimen)
CASE(3)
mesh%connectMesh => connectMesh3DCart
CASE(2)
SELECT CASE(mesh%geometry)
CASE("Cyl")
mesh%connectMesh => connectMesh2DCyl
CASE("Cart")
mesh%connectMesh => connectMesh2DCart
END SELECT
CASE(1)
SELECT CASE(mesh%geometry)
CASE("Rad")
mesh%connectMesh => connectMesh1DRad
CASE("Cart")
mesh%connectMesh => connectMesh1DCart
END SELECT
END SELECT
IF (doubleMesh) THEN
meshColl%connectMesh => mesh%connectMesh
END IF
!Get the format of mesh
CALL config%get(object // '.meshType', meshFormat, found)
SELECT CASE(meshFormat)
CASE ("gmsh2")
@ -967,6 +1006,16 @@ MODULE moduleInput
END IF
CASE ("vtu")
CALL initVTU(mesh)
IF (doubleMesh) THEN
CALL initVTU(meshColl)
END IF
CASE DEFAULT
CALL criticalError('Mesh format ' // meshFormat // ' not defined.', 'readGeometry')
END SELECT
!Reads the mesh file

2
src/modules/makefile
View file

@ -11,7 +11,7 @@ common.o:
output.o: moduleSpecies.o common.o
$(MAKE) -C output all
mesh.o: moduleCollisions.o moduleBoundary.o output.o
mesh.o: moduleCollisions.o moduleBoundary.o output.o common.o
$(MAKE) -C mesh all
solver.o: moduleSpecies.o moduleProbe.o common.o output.o mesh.o

4
src/modules/mesh/inout/0D/moduleMeshOutput0D.f90
View file

@ -42,13 +42,13 @@ MODULE moduleMeshOutput0D
USE moduleOutput
IMPLICIT NONE
CLASS(meshGeneric), INTENT(inout):: self
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
CHARACTER(:), ALLOCATABLE:: fileName
INTEGER:: k
fileName='OUTPUT_Collisions.dat'
IF (t == 0) THEN
IF (t == tInitial) THEN
OPEN(20, file = path // folder // '/' // fileName, action = 'write')
WRITE(20, "(A1, 14X, A5, A20)") "#","t (s)","collisions"
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName

41
src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90
View file

@ -2,7 +2,7 @@ MODULE moduleMeshInputGmsh2
!Reads a mesh in the Gmsh v2.0 format
CONTAINS
!Inits a mesh to use Gmsh2 format
!Init a mesh to use Gmsh2 format
SUBROUTINE initGmsh2(self)
USE moduleMesh
USE moduleMeshOutputGmsh2
@ -23,7 +23,7 @@ MODULE moduleMeshInputGmsh2
END SUBROUTINE initGmsh2
!Reads a Gmsh 2 format
!Read a Gmsh 2 format
SUBROUTINE readGmsh2(self, filename)
USE moduleMesh3DCart
USE moduleMesh2DCyl
@ -71,30 +71,26 @@ MODULE moduleMeshInputGmsh2
SELECT CASE(self%dimen)
CASE(3)
ALLOCATE(meshNode3Dcart::self%nodes(n)%obj)
self%connectMesh => connectMesh3DCart
CASE(2)
SELECT CASE(self%geometry)
CASE("Cyl")
ALLOCATE(meshNode2DCyl:: self%nodes(n)%obj)
self%connectMesh => connectMesh2DCyl
CASE("Cart")
ALLOCATE(meshNode2DCart:: self%nodes(n)%obj)
self%connectMesh => connectMesh2DCart
END SELECT
r(3) = 0.D0
CASE(1)
SELECT CASE(self%geometry)
CASE("Rad")
ALLOCATE(meshNode1DRad:: self%nodes(n)%obj)
self%connectMesh => connectMesh1DRad
CASE("Cart")
ALLOCATE(meshNode1DCart:: self%nodes(n)%obj)
self%connectMesh => connectMesh1DCart
END SELECT
r(2:3) = 0.D0
@ -111,7 +107,7 @@ MODULE moduleMeshInputGmsh2
!Reads total number of elements (no nodes)
READ(10, *) totalNumElem
!conts edges and volume elements
!Count edges and volume elements
SELECT TYPE(self)
TYPE IS(meshParticles)
self%numEdges = 0
@ -151,7 +147,7 @@ MODULE moduleMeshInputGmsh2
END SELECT
!Allocates arrays
!Allocates array of cells
ALLOCATE(self%cells(1:self%numCells))
SELECT TYPE(self)
@ -179,45 +175,30 @@ MODULE moduleMeshInputGmsh2
END SELECT
CASE (2)
SELECT CASE(self%geometry)
CASE("Cyl")
ALLOCATE(p(1:2))
READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2)
!Associate boundary condition procedure.
bt = getBoundaryId(boundaryType)
SELECT CASE(self%geometry)
CASE("Cyl")
ALLOCATE(meshEdge2DCyl:: self%edges(e)%obj)
CASE("Cart")
ALLOCATE(p(1:2))
READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2)
!Associate boundary condition procedure.
bt = getBoundaryId(boundaryType)
ALLOCATE(meshEdge2DCart:: self%edges(e)%obj)
END SELECT
CASE(1)
SELECT CASE(self%geometry)
CASE("Rad")
ALLOCATE(p(1:1))
READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1)
!Associate boundary condition
bt = getBoundaryId(boundaryType)
SELECT CASE(self%geometry)
CASE("Rad")
ALLOCATE(meshEdge1DRad:: self%edges(e)%obj)
CASE("Cart")
ALLOCATE(p(1:1))
READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1)
!Associate boundary condition
bt = getBoundaryId(boundaryType)
ALLOCATE(meshEdge1DCart:: self%edges(e)%obj)
END SELECT
@ -231,9 +212,9 @@ MODULE moduleMeshInputGmsh2
END SELECT
!Read and initialize volumes
!Read and initialize cells
DO e = 1, self%numCells
!Reads the volume according to the geometry
!Read the cell according to the geometry
SELECT CASE(self%dimen)
CASE(3)
READ(10, *) n, elemType

22
src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
View file

@ -85,6 +85,7 @@ MODULE moduleMeshOutputGmsh2
USE moduleRefParam
USE moduleSpecies
USE moduleOutput
USE moduleMeshInoutCommon
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
@ -93,13 +94,11 @@ MODULE moduleMeshOutputGmsh2
TYPE(outputFormat):: output(1:self%numNodes)
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
CHARACTER (LEN=iterationDigits):: tstring
time = DBLE(t)*tauMin*ti_ref
DO i = 1, nSpecies
WRITE(tstring, iterationFormat) t
fileName= 'OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh'
fileName = formatFileName(prefix, species(i)%obj%name, 'msh', t)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (60, file = path // folder // '/' // fileName)
@ -142,16 +141,16 @@ MODULE moduleMeshOutputGmsh2
USE moduleCaseParam
USE moduleCollisions
USE moduleOutput
USE moduleMeshInoutCommon
IMPLICIT NONE
CLASS(meshGeneric), INTENT(inout):: self
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: numEdges
INTEGER:: k, c
INTEGER:: n
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
CHARACTER (LEN=iterationDigits):: tString
CHARACTER(:), ALLOCATABLE:: title
CHARACTER (LEN=2):: cString
@ -169,9 +168,8 @@ MODULE moduleMeshOutputGmsh2
IF (collOutput) THEN
time = DBLE(t)*tauMin*ti_ref
WRITE(tString, iterationFormat) t
fileName='OUTPUT_' // tString// '_Collisions.msh'
fileName = formatFileName(prefix, 'Collisions', 'msh', t)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (60, file = path // folder // '/' // fileName)
@ -203,6 +201,7 @@ MODULE moduleMeshOutputGmsh2
USE moduleRefParam
USE moduleCaseParam
USE moduleOutput
USE moduleMeshInoutCommon
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
@ -210,16 +209,14 @@ MODULE moduleMeshOutputGmsh2
INTEGER:: n, e
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
CHARACTER (LEN=iterationDigits):: tstring
REAL(8):: Xi(1:3)
Xi = (/ 0.D0, 0.D0, 0.D0 /)
IF (emOutput) THEN
time = DBLE(t)*tauMin*ti_ref
WRITE(tstring, iterationFormat) t
fileName='OUTPUT_' // tstring// '_EMField.msh'
fileName = formatFileName(prefix, 'Collisions', 'msh', t)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (20, file = path // folder // '/' // fileName)
@ -256,6 +253,7 @@ MODULE moduleMeshOutputGmsh2
USE moduleSpecies
USE moduleOutput
USE moduleAverage
USE moduleMeshInoutCommon
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
@ -265,11 +263,11 @@ MODULE moduleMeshOutputGmsh2
INTEGER:: fileMean=10, fileDeviation=20
DO i = 1, nSpecies
fileName= 'Average_mean_' // species(i)%obj%name // '.msh'
fileName = formatFileName('Average_mean', species(i)%obj%name, 'msh')
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (fileMean, file = path // folder // '/' // fileName)
fileName= 'Average_deviation_' // species(i)%obj%name // '.msh'
fileName = formatFileName('Average_deviation', species(i)%obj%name, 'msh')
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (filedeviation, file = path // folder // '/' // fileName)

8
src/modules/mesh/inout/makefile
View file

@ -1,7 +1,13 @@
all: gmsh2.o 0D.o
all: vtu.o gmsh2.o 0D.o
vtu.o: moduleMeshInoutCommon.o
$(MAKE) -C vtu all
gmsh2.o:
$(MAKE) -C gmsh2 all
0D.o:
$(MAKE) -C 0D all
%.o: %.f90
$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@

27
src/modules/mesh/inout/moduleMeshInoutCommon.f90 Normal file
View file

@ -0,0 +1,27 @@
MODULE moduleMeshInoutCommon
CHARACTER(LEN=4):: prefix = 'Step'
CONTAINS
PURE FUNCTION formatFileName(prefix, suffix, extension, t) RESULT(fileName)
USE moduleOutput
IMPLICIT NONE
CHARACTER(*), INTENT(in):: prefix, suffix, extension
INTEGER, INTENT(in), OPTIONAL:: t
CHARACTER (LEN=iterationDigits):: tString
CHARACTER(:), ALLOCATABLE:: fileName
IF (PRESENT(t)) THEN
WRITE(tString, iterationFormat) t
fileName = prefix // '_' // tString // '_' // suffix // '.' // extension
ELSE
fileName = prefix // '_' // suffix // '.' // extension
END IF
END FUNCTION formatFileName
END MODULE moduleMeshInoutCommon

7
src/modules/mesh/inout/vtu/makefile Normal file
View file

@ -0,0 +1,7 @@
all: moduleMeshInputVTU.o moduleMeshOutputVTU.o
moduleMeshInputVTU.o: moduleMeshOutputVTU.o moduleMeshInputVTU.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
%.o: %.f90
$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@

552
src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 Normal file
View file

@ -0,0 +1,552 @@
MODULE moduleMeshInputVTU
!Reads mesh in the VTU format
INTERFACE getValueFromLine
MODULE PROCEDURE getIntegerFromLine, getRealFromLine
END INTERFACE
INTERFACE readDataBlock
MODULE PROCEDURE readIntegerBlock, readRealBlock
END INTERFACE
CONTAINS
FUNCTION findLine(fileID, text) RESULT(line)
USE moduleErrors
IMPLICIT NONE
INTEGER, INTENT(in):: fileID
CHARACTER(*):: text
CHARACTER(LEN=256):: line
INTEGER:: error, found
error = 0
found = 0
!Reads all the file to find the 'text' string in a line
DO WHILE(error == 0 .AND. found == 0)
READ(fileID, "(A)", IOSTAT=error) line
found = INDEX(line, text)
IF (found > 0) THEN
EXIT
END IF
END DO
!If no line is found, return an error
IF (found == 0) THEN
CALL criticalError('String ' // text // ' not found in file.', 'findLine')
END IF
END FUNCTION findLine
SUBROUTINE getIntegerFromLine(line, label, valueInteger)
USE moduleErrors
IMPLICIT NONE
CHARACTER(LEN=256), INTENT(in):: line
CHARACTER(*), INTENT(in):: label
INTEGER, INTENT(out):: valueInteger
INTEGER:: labelStart, valueStart, valueEnd
labelStart = 0
labelStart = INDEX(line, label)
IF (labelStart == 0) THEN
CALL criticalError('Label '// label // ' not found in line', 'getValueFromLine')
END IF
valueStart = INDEX(line(labelStart:256), '"') + labelStart
valueEnd = INDEX(line(valueStart:256), '"') + valueStart - 2
READ(line(valueStart:valueEnd), '(I10)') valueInteger
END SUBROUTINE getIntegerFromLine
SUBROUTINE getRealFromLine(line, label, valueReal)
USE moduleErrors
IMPLICIT NONE
CHARACTER(LEN=256), INTENT(in):: line
CHARACTER(*), INTENT(in):: label
REAL(8), INTENT(out):: valueReal
INTEGER:: labelStart, valueStart, valueEnd
labelStart = 0
labelStart = INDEX(line, label)
IF (labelStart == 0) THEN
CALL criticalError('Label '// label // ' not found in line', 'getValueFromLine')
END IF
valueStart = INDEX(line(labelStart:256), '"') + labelStart
valueEnd = INDEX(line(valueStart:256), '"') + valueStart - 2
READ(line(valueStart:valueEnd), '(F16.14)') valueReal
END SUBROUTINE getRealFromLine
SUBROUTINE readIntegerBlock(fileID, nData, array)
IMPLICIT NONE
INTEGER, INTENT(in):: fileID
INTEGER, INTENT(in):: nData
INTEGER, INTENT(out):: array(1:nData)
INTEGER:: iStart, iEnd, block
iStart = 0
iEnd = 0
block = 6 !Assumes block of data in 6 columns
DO WHILE (iStart < nData)
iStart = iStart + 1
iEnd = iStart - 1 + block
IF (iEnd > nData) THEN
iEnd = nData
END IF
READ(fileID, *) array(iStart:iEnd)
iStart = iEnd
END DO
END SUBROUTINE readIntegerBlock
SUBROUTINE readRealBlock(fileID, nData, array)
IMPLICIT NONE
INTEGER, INTENT(in):: fileID
INTEGER, INTENT(in):: nData
REAL(8), INTENT(out):: array(1:nData)
INTEGER:: iStart, iEnd, block
iStart = 0
iEnd = 0
block = 6 !Assumes block of data in 6 columns
DO WHILE (iStart < nData)
iStart = iStart + 1
iEnd = iStart - 1 + block
IF (iEnd > nData) THEN
iEnd = nData
END IF
READ(fileID, *) array(iStart:iEnd)
iStart = iEnd
END DO
END SUBROUTINE readRealBlock
SUBROUTINE initVTU(self)
USE moduleMesh
USE moduleMeshOutputVTU
IMPLICIT NONE
CLASS(meshGeneric), INTENT(inout), TARGET:: self
IF (ASSOCIATED(meshForMCC, self)) self%printColl => printCollVTU
SELECT TYPE(self)
TYPE IS(meshParticles)
self%printOutput => printOutputVTU
self%printEM => printEMVTU
self%readInitial => readInitialVTU
self%printAverage => printAverageVTU
END SELECT
self%readMesh => readVTU
END SUBROUTINE initVTU
SUBROUTINE readVTU(self, filename)
USE moduleMesh3DCart
USE moduleMesh2DCyl
USE moduleMesh2DCart
USE moduleMesh1DRad
USE moduleMesh1DCart
USE moduleBoundary
IMPLICIT NONE
CLASS(meshGeneric), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8):: r(1:3) !3 generic coordinates
INTEGER:: fileID, error, found
CHARACTER(LEN=256):: line
INTEGER:: numNodes, numElements, numEdges
INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types
REAL(8), ALLOCATABLE, DIMENSION(:):: coordinates
INTEGER:: n, e, c
INTEGER, ALLOCATABLE:: p(:)
INTEGER:: bt
fileID = 10
OPEN(fileID, FILE=TRIM(filename))
!Find the number of nodes and elements (edges+cells) in the mesh.
line = findLine(fileID, '<Piece')
CALL getValueFromLine(line, 'NumberOfPoints', numNodes)
CALL getValueFromLine(line, 'NumberOfCells', numElements)
REWIND(fileID)
!Get the IDs of the cells to identify physical surfaces
line = findLine(fileID, 'Name="CellEntityIds"')
ALLOCATE(entitiesID(1:numElements))
CALL readDataBlock(fileID, numElements, entitiesID)
REWIND(fileID)
!Get the offsets to read connectivity
line = findLine(fileID, 'Name="offsets"')
ALLOCATE(offsets(1:numElements))
CALL readDataBlock(fileID, numElements, offsets)
REWIND(fileID)
!Get the connectivity of elements to nodes
line = findline(fileID, 'Name="connectivity"')
ALLOCATE(connectivity(1:MAXVAL(offsets)))
CALL readDataBlock(fileID, MAXVAL(offsets), connectivity)
!Shift connectivity to start in 1
connectivity = connectivity + 1
REWIND(fileID)
!Get the type of elements
line = findline(fileID, 'Name="types"')
ALLOCATE(types(1:numElements))
CALL readDataBlock(fileID, numElements, types)
REWIND(fileID)
!Get nodes coordinates
line = findline(fileID, 'Name="Points"')
ALLOCATE(coordinates(1:3*numNodes))
CALL readDataBlock(fileID, 3*numNodes, coordinates)
REWIND(fileID)
CLOSE(fileID)
!All relevant information from the .vtu file has been read. Time to build the mesh.
self%numNodes = numNodes
ALLOCATE(self%nodes(1:self%numNodes))
SELECT TYPE(self)
TYPE IS(meshParticles)
ALLOCATE(self%K(1:self%numNodes, 1:self%numNodes))
ALLOCATE(self%IPIV(1:self%numNodes, 1:self%numNodes))
self%K = 0.D0
self%IPIV = 0
END SELECT
DO n = 1, self%numNodes
!Get the coordinates for each direction
r(1) = coordinates(3*(n-1)+1)
r(2) = coordinates(3*(n-1)+2)
r(3) = coordinates(3*(n-1)+3)
!Allocate node
SELECT CASE(self%dimen)
CASE(3)
ALLOCATE(meshNode3Dcart::self%nodes(n)%obj)
CASE(2)
SELECT CASE(self%geometry)
CASE("Cyl")
ALLOCATE(meshNode2DCyl:: self%nodes(n)%obj)
CASE("Cart")
ALLOCATE(meshNode2DCart:: self%nodes(n)%obj)
END SELECT
r(3) = 0.D0
CASE(1)
SELECT CASE(self%geometry)
CASE("Rad")
ALLOCATE(meshNode1DRad:: self%nodes(n)%obj)
CASE("Cart")
ALLOCATE(meshNode1DCart:: self%nodes(n)%obj)
END SELECT
r(2:3) = 0.D0
END SELECT
!Init node
CALL self%nodes(n)%obj%init(n, r)
END DO
!Count the number of edges
SELECT CASE(self%dimen)
CASE(3)
!Edges are triangles, type 5 in VTK
numEdges = COUNT(types==5)
CASE(2)
!Edges are segments, type 3 in VTK
numEdges = COUNT(types==3)
CASE(1)
!Edges are nodes, type 1 in VTK
numEdges = COUNT(types==1)
END SELECT
self%numCells = numElements - numEdges
SELECT TYPE(self)
TYPE IS(meshParticles)
self%numEdges = numEdges
!Allocate array of edges
ALLOCATE(self%edges(1:self%numEdges))
END SELECT
!Allocates array of cells
ALLOCATE(self%cells(1:self%numCells))
!Read edges only if mesh is for tracking particles, not for collisions
e = 0
SELECT TYPE(self)
TYPE IS(meshParticles)
DO n = 1, numElements
SELECT CASE(self%dimen)
CASE(3)
IF (types(n) == 5) THEN
e = e + 1
ALLOCATE(p(1:3))
p(1) = connectivity(offsets(n) - 2)
p(2) = connectivity(offsets(n) - 1)
p(3) = connectivity(offsets(n))
!Associate boundary condition procedure.
bt = getBoundaryId(entitiesID(n))
!Allocate edge
ALLOCATE(meshEdge3DCartTria:: self%edges(e)%obj)
!Init edge
CALL self%edges(e)%obj%init(n, p, bt, entitiesID(n))
DEALLOCATE(p)
END IF
CASE(2)
IF (types(n) == 3) THEN
e = e + 1
ALLOCATE(p(1:2))
p(1) = connectivity(offsets(n) - 1)
p(2) = connectivity(offsets(n))
!Associate boundary condition procedure.
bt = getBoundaryId(entitiesID(n))
!Allocate edge
SELECT CASE(self%geometry)
CASE("Cyl")
ALLOCATE(meshEdge2DCyl:: self%edges(e)%obj)
CASE("Cart")
ALLOCATE(meshEdge2DCart:: self%edges(e)%obj)
END SELECT
!Init edge
CALL self%edges(e)%obj%init(n, p, bt, entitiesID(n))
DEALLOCATE(p)
END IF
CASE(1)
IF (types(n) == 1) THEN
e = e + 1
ALLOCATE(p(1:1))
p(1) = connectivity(offsets(n))
!Associate boundary condition procedure.
bt = getBoundaryId(entitiesID(n))
!Allocate edge
SELECT CASE(self%geometry)
CASE("Rad")
ALLOCATE(meshEdge1DRad:: self%edges(e)%obj)
CASE("Cart")
ALLOCATE(meshEdge1DCart:: self%edges(e)%obj)
END SELECT
!Init edge
CALL self%edges(e)%obj%init(n, p, bt, entitiesID(n))
DEALLOCATE(p)
END IF
END SELECT
END DO
END SELECT
!Read cells
c = 0
DO n = 1, numElements
SELECT CASE(self%dimen)
CASE(3)
SELECT CASE(types(n))
CASE(10)
!Thetraedron
c = c + 1
ALLOCATE(p(1:4))
p(1) = connectivity(offsets(n) - 3)
p(2) = connectivity(offsets(n) - 2)
p(3) = connectivity(offsets(n) - 1)
p(4) = connectivity(offsets(n))
!Allocate cell
ALLOCATE(meshCell3DCartTetra:: self%cells(c)%obj)
!Init cell
CALL self%cells(c)%obj%init(c, p, self%nodes)
DEALLOCATE(p)
END SELECT
CASE(2)
SELECT CASE(types(n))
CASE(5)
!Triangular element
c = c + 1
ALLOCATE(p(1:3))
p(1) = connectivity(offsets(n) - 2)
p(2) = connectivity(offsets(n) - 1)
p(3) = connectivity(offsets(n))
!Allocate cell
SELECT CASE(self%geometry)
CASE("Cyl")
ALLOCATE(meshCell2DCylTria:: self%cells(c)%obj)
CASE("Cart")
ALLOCATE(meshCell2DCartTria:: self%cells(c)%obj)
END SELECT
!Init cell
CALL self%cells(c)%obj%init(c, p, self%nodes)
DEALLOCATE(p)
CASE(9)
!Quadrilateral element
c = c + 1
ALLOCATE(p(1:4))
p(1) = connectivity(offsets(n) - 3)
p(2) = connectivity(offsets(n) - 2)
p(3) = connectivity(offsets(n) - 1)
p(4) = connectivity(offsets(n))
!Allocate cell
SELECT CASE(self%geometry)
CASE("Cyl")
ALLOCATE(meshCell2DCylQuad:: self%cells(c)%obj)
CASE("Cart")
ALLOCATE(meshCell2DCartQuad:: self%cells(c)%obj)
END SELECT
!Init cell
CALL self%cells(c)%obj%init(c, p, self%nodes)
DEALLOCATE(p)
END SELECT
CASE(1)
SELECT CASE(types(n))
CASE(3)
!Segment element
c = c + 1
ALLOCATE(p(1:2))
p(1) = connectivity(offsets(n) - 1)
p(2) = connectivity(offsets(n))
!Allocate cell
SELECT CASE(self%geometry)
CASE("Rad")
ALLOCATE(meshCell1DRadSegm:: self%cells(c)%obj)
CASE("Cart")
ALLOCATE(meshCell1DCartSegm:: self%cells(c)%obj)
END SELECT
!Init cell
CALL self%cells(c)%obj%init(c, p, self%nodes)
DEALLOCATE(p)
END SELECT
END SELECT
END DO
!Call mesh connectivity
CALL self%connectMesh
END SUBROUTINE readVTU
SUBROUTINE readInitialVTU(filename, density, velocity, temperature)
IMPLICIT NONE
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density
REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity
REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: temperature
INTEGER:: fileID
CHARACTER(:), ALLOCATABLE:: line
INTEGER:: numNodes
REAL(8), ALLOCATABLE:: velocityBlock(:)
INTEGER:: n
fileID = 10
OPEN(fileID, file = TRIM(filename))
line = findLine(fileID, '<Piece')
CALL getValueFromLine(line, 'NumberOfPoints', numNodes)
!Read the species density
line = findLine(fileID, 'Name="Density')
ALLOCATE(density(1:numNodes))
CALL readDataBlock(fileID, numNodes, density)
REWIND(fileID)
!Read the species velocity
line = findLine(fileID, 'Name="Velocity')
ALLOCATE(velocityBlock(1:3*numNodes))
CALL readDataBlock(fileID, 3*numNodes, velocityBlock)
ALLOCATE(velocity(1:numNodes, 1:3))
DO n = 1, numNodes
velocity(n, 1) = velocityBlock(3*(n-1)+1)
velocity(n, 2) = velocityBlock(3*(n-1)+2)
velocity(n, 3) = velocityBlock(3*(n-1)+3)
END DO
DEALLOCATE(velocityBlock)
REWIND(fileID)
!Read the species temperature
line = findLine(fileID, 'Name="Temperature')
ALLOCATE(temperature(1:numNodes))
CALL readDataBlock(fileID, numNodes, temperature)
REWIND(fileID)
END SUBROUTINE readInitialVTU
END MODULE moduleMeshInputVTU

461
src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 Normal file
View file

@ -0,0 +1,461 @@
MODULE moduleMeshOutputVTU
CONTAINS
SUBROUTINE writeHeader(nNodes, nCells, fileID)
USE moduleMesh
IMPLICIT NONE
INTEGER, INTENT(in):: nNodes, nCells
INTEGER, INTENT(in):: fileID
WRITE(fileID,"(A)") '<?xml version="1.0"?>'
WRITE(fileID,"(2X, A)") '<VTKFile type="UnstructuredGrid">'
WRITE(fileID,"(4X, A,ES20.6E3,A)") '<UnstructuredGrid>'
WRITE(fileID,"(6X, A, I10, A, I10, A)") '<Piece NumberOfPoints="', nNodes, '" NumberOfCells="', nCells, '">'
END SUBROUTINE writeHeader
SUBROUTINE writeFooter(fileID)
IMPLICIT NONE
INTEGER, INTENT(in):: fileID
WRITE(fileID,"(6X, A)") '</Piece>'
WRITE(fileID,"(4X, A)") '</UnstructuredGrid>'
WRITE(fileID,"(2X, A)") '</VTKFile>'
END SUBROUTINE writeFooter
FUNCTION getCellType(cell) RESULT(indexType)
USE moduleMesh3DCart
USE moduleMesh2DCyl
USE moduleMesh2DCart
USE moduleMesh1DRad
USE moduleMesh1DCart
USE moduleMesh0D
USE moduleErrors
IMPLICIT NONE
CLASS(meshCell), INTENT(in):: cell
INTEGER:: indexType
indexType = 0
SELECT TYPE(cell)
TYPE IS(meshCell3DCartTetra)
indexType = 10
TYPE IS(meshCell2DCylQuad)
indexType = 9
TYPE IS(meshCell2DCartQuad)
indexType = 9
TYPE IS(meshCell2DCylTria)
indexType = 5
TYPE IS(meshCell2DCartTria)
indexType = 5
TYPE IS(meshCell1DRadSegm)
indexType = 3
TYPE IS(meshCell1DCartSegm)
indexType = 3
TYPE IS(meshCell0D)
indexType = 1
CLASS DEFAULT
CALL criticalError('Cell not valid for VTU output', 'getCellType')
END SELECT
END FUNCTION getCellType
SUBROUTINE writeMesh(self, fileID)
USE moduleMesh
USE moduleRefParam
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: fileID
INTEGER:: e
INTEGER:: offset(1:self%numCells), types(1:self%numCells)
!Write nodes coordinates
WRITE(fileID, "(8X, A)") '<Points>'
WRITE(fileID, "(10X,A)") '<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
WRITE(fileID, "(6(ES20.6E3))") (self%nodes(e)%obj%getCoordinates()*L_ref, e = 1, self%numNodes)
WRITE(fileID, "(10X, A)") '</DataArray>'
WRITE(fileID, "(8X, A)") '</Points>'
WRITE(fileID, "(8X, A)") '<Cells>'
!Write nodes connectivity of each cell
WRITE(fileID, "(10X,A)") '<DataArray type="Int32" Name="connectivity" format="ascii">'
WRITE(fileID, "(6(I10))") (self%cells(e)%obj%getNodes(self%cells(e)%obj%nNodes) - 1, e = 1, self%numCells) !Array starts on 0
WRITE(fileID, "(10X, A)") '</DataArray>'
!Write offset of each cell
offset(1) = self%cells(1)%obj%nNodes
WRITE(fileID, "(10X,A)") '<DataArray type="Int32" Name="offsets" format="ascii">'
DO e = 2, self%numCells
offset(e) = offset(e - 1) + self%cells(e)%obj%nNodes
END DO
WRITE(fileID, "(6(I10))") offset
WRITE(fileID, "(10X, A)") '</DataArray>'
!Write type of each cell
WRITE(fileID, "(10X,A)") '<DataArray type="Int32" Name="types" format="ascii">'
DO e = 1, self%numCells
types(e) = getCellType(self%cells(e)%obj)
END DO
WRITE(fileID, "(6(I10))") types
WRITE(fileID, "(10X, A)") '</DataArray>'
WRITE(fileID, "(8X, A)") '</Cells>'
END SUBROUTINE writeMesh
SUBROUTINE writeSpeciesOutput(self, fileID, speciesIndex)
USE moduleMesh
USE moduleOutput
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: fileID
INTEGER, INTENT(in):: speciesIndex
TYPE(outputFormat):: output(1:self%numNodes)
INTEGER:: n
DO n = 1, self%numNodes
CALL calculateOutput(self%nodes(n)%obj%output(speciesIndex), output(n), self%nodes(n)%obj%v, species(speciesIndex)%obj)
END DO
!Write node data
WRITE(fileID,"(8X,A)") '<PointData>'
!Write density
WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Density (m^-3)" NumberOfComponents="1">'
WRITE(fileID,"(6(ES20.6E3))") (output(n)%density, n = 1, self%numNodes)
WRITE(fileID,"(10X,A)") '</DataArray>'
!Write velocity
WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Velocity (m s^-1)" NumberOfComponents="3">'
WRITE(fileID,"(6(ES20.6E3))") (output(n)%velocity, n = 1, self%numNodes)
WRITE(fileID,"(10X,A)") '</DataArray>'
!Write pressure
WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Pressure (Pa)" NumberOfComponents="1">'
WRITE(fileID,"(6(ES20.6E3))") (output(n)%pressure, n = 1, self%numNodes)
WRITE(fileID,"(10X,A)") '</DataArray>'
!Write temperature
WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Temperature (K)" NumberOfComponents="1">'
WRITE(fileID,"(6(ES20.6E3))") (output(n)%temperature, n = 1, self%numNodes)
WRITE(fileID,"(10X,A)") '</DataArray>'
!End node data
WRITE(fileID,"(8X,A)") '</PointData>'
END SUBROUTINE writeSpeciesOutput
SUBROUTINE writeCollOutput(self,fileID)
USE moduleMesh
USE moduleCollisions
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: fileID
INTEGER:: k, c, n
CHARACTER(:), ALLOCATABLE:: title
CHARACTER (LEN=2):: cString
WRITE(fileID,"(8X,A)") '<CellData>'
DO k = 1, nCollPairs
DO c = 1, interactionMatrix(k)%amount
WRITE(cString, "(I2)") c
title = 'Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString
WRITE(fileID,"(10X,A, A, A)") '<DataArray type="Float64" Name="',title, '" NumberOfComponents="1">'
WRITE(fileID, "(6(I10))") (self%cells(n)%obj%tallyColl(k)%tally(c), n = 1, self%numCells)
WRITE(fileID, "(10X, A)") '</DataArray>'
END DO
END DO
WRITE(fileID,"(8X,A)") '</CellData>'
END SUBROUTINE writeCollOutput
SUBROUTINE writeEM(self, fileID)
USE moduleMesh
USE moduleRefParam
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: fileID
INTEGER:: n
REAL(8):: Xi(1:3)
!Points in nodes
WRITE(fileID,"(8X,A)") '<PointData>'
!Electric potential
WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Potential (V)" NumberOfComponents="1">'
WRITE(fileID,"(6(ES20.6E3))") (self%nodes(n)%obj%emData%phi*Volt_ref, n = 1, self%numNodes)
WRITE(fileID,"(10X,A)") '</DataArray>'
!Magnetic Field
WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Magnetic Field (T)" NumberOfComponents="3">'
WRITE(fileID,"(6(ES20.6E3))") (self%nodes(n)%obj%emData%B*B_ref, n = 1, self%numNodes)
WRITE(fileID,"(10X,A)") '</DataArray>'
WRITE(fileID,"(8X,A)") '</PointData>'
!Cell Data
Xi = (/ 0.D0, 0.D0, 0.D0 /)
WRITE(fileID,"(8X,A)") '<CellData>'
!Electric field
WRITE(fileID,"(10X,A, A, A)") '<DataArray type="Float64" Name="Electric Field (V m^-1)" NumberOfComponents="3">'
WRITE(fileID, "(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells)
WRITE(fileID,"(10X,A)") '</DataArray>'
WRITE(fileID,"(8X,A)") '</CellData>'
END SUBROUTINE writeEM
SUBROUTINE writeCollection(fileID, t, fileNameStep, fileNameCollection)
USE moduleCaseParam
USE moduleOutput
USE moduleRefParam
IMPLICIT NONE
INTEGER:: fileID
INTEGER, INTENT(in):: t
CHARACTER(*):: fileNameStep, fileNameCollection
IF (t == tInitial) THEN
!Create collection file
WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection
OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection)
WRITE (fileID + 1, "(A)") '<VTKFile type="Collection">'
WRITE (fileID + 1, "(2X, A)") '<Collection>'
CLOSE(fileID + 1)
END IF
!Write iteration file in collection
OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection, ACCESS='APPEND')
WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") '<DataSet timestep="', DBLE(t)*tauMin*ti_ref,'" file="', fileNameStep,'"/>'
!Close collection file
IF (t == tFinal) THEN
WRITE (fileID + 1, "(2X, A)") '</Collection>'
WRITE (fileID + 1, "(A)") '</VTKFile>'
END IF
CLOSE(fileID + 1)
END SUBROUTINE writeCollection
SUBROUTINE writeAverage(self, fileIDMean, fileIDDeviation, speciesIndex)
USE moduleMesh
USE moduleOutput
USE moduleAverage
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: fileIDMean, fileIDDeviation
INTEGER, INTENT(in):: speciesIndex
TYPE(outputFormat):: outputMean(1:self%numNodes)
TYPE(outputFormat):: outputDeviation(1:self%numNodes)
INTEGER:: n
DO n = 1, self%numNodes
CALL calculateOutput(averageScheme(n)%mean%output(speciesIndex), outputMean(n), &
self%nodes(n)%obj%v, species(speciesIndex)%obj)
CALL calculateOutput(averageScheme(n)%deviation%output(speciesIndex), outputDeviation(n), &
self%nodes(n)%obj%v, species(speciesIndex)%obj)
END DO
!Write node data
WRITE(fileIDMean, "(8X,A)") '<PointData>'
WRITE(fileIDDeviation,"(8X,A)") '<PointData>'
!Write density
WRITE(fileIDMean, "(10X,A)") '<DataArray type="Float64" Name="Density, mean (m^-3)" NumberOfComponents="1">'
WRITE(fileIDDeviation,"(10X,A)") '<DataArray type="Float64" Name="Density, deviation (m^-3)" NumberOfComponents="1">'
WRITE(fileIDMean, "(6(ES20.6E3))") (outputMean(n)%density, n = 1, self%numNodes)
WRITE(fileIDDeviation,"(6(ES20.6E3))") (outputDeviation(n)%density, n = 1, self%numNodes)
WRITE(fileIDMean,"(10X,A)") '</DataArray>'
WRITE(fileIDDeviation,"(10X,A)") '</DataArray>'
!Write velocity
WRITE(fileIDMean, "(10X,A)") '<DataArray type="Float64" Name="Velocity, mean (m s^-1)" NumberOfComponents="3">'
WRITE(fileIDDeviation,"(10X,A)") '<DataArray type="Float64" Name="Velocity, deviation (m s^-1)" NumberOfComponents="3">'
WRITE(fileIDMean, "(6(ES20.6E3))") (outputMean(n)%velocity, n = 1, self%numNodes)
WRITE(fileIDDeviation,"(6(ES20.6E3))") (outputDeviation(n)%velocity, n = 1, self%numNodes)
WRITE(fileIDMean, "(10X,A)") '</DataArray>'
WRITE(fileIDDeviation,"(10X,A)") '</DataArray>'
!Write pressure
WRITE(fileIDMean, "(10X,A)") '<DataArray type="Float64" Name="Pressure, mean (Pa)" NumberOfComponents="1">'
WRITE(fileIDDeviation,"(10X,A)") '<DataArray type="Float64" Name="Pressure, deviation (Pa)" NumberOfComponents="1">'
WRITE(fileIDMean, "(6(ES20.6E3))") (outputMean(n)%pressure, n = 1, self%numNodes)
WRITE(fileIDDeviation,"(6(ES20.6E3))") (outputDeviation(n)%pressure, n = 1, self%numNodes)
WRITE(fileIDMean, "(10X,A)") '</DataArray>'
WRITE(fileIDDeviation,"(10X,A)") '</DataArray>'
!Write temperature
WRITE(fileIDMean, "(10X,A)") '<DataArray type="Float64" Name="Temperature, mean (K)" NumberOfComponents="1">'
WRITE(fileIDDeviation,"(10X,A)") '<DataArray type="Float64" Name="Temperature, deviation (K)" NumberOfComponents="1">'
WRITE(fileIDMean, "(6(ES20.6E3))") (outputMean(n)%temperature, n = 1, self%numNodes)
WRITE(fileIDDeviation,"(6(ES20.6E3))") (outputDeviation(n)%temperature, n = 1, self%numNodes)
WRITE(fileIDMean, "(10X,A)") '</DataArray>'
WRITE(fileIDDeviation,"(10X,A)") '</DataArray>'
!End node data
WRITE(fileIDMean, "(8X,A)") '</PointData>'
WRITE(fileIDDeviation,"(8X,A)") '</PointData>'
END SUBROUTINE writeAverage
SUBROUTINE printOutputVTU(self,t)
USE moduleMesh
USE moduleSpecies
USE moduleMeshInoutCommon
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n, i, fileID
CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
TYPE(outputFormat):: output(1:self%numNodes)
fileID = 60
DO i = 1, nSpecies
fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', t)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (fileID, file = path // folder // '/' // fileName)
CALL writeHeader(self%numNodes, self%numCells, fileID)
CALL writeMesh(self, fileID)
CALL writeSpeciesOutput(self, fileID, i)
CALL writeFooter(fileID)
CLOSE(fileID)
!Write collection file for time plotting
fileNameCollection = formatFileName('Collection', species(i)%obj%name, 'pvd')
CALL writeCollection(fileID, t, fileName, filenameCollection)
END DO
END SUBROUTINE printOutputVTU
SUBROUTINE printCollVTU(self,t)
USE moduleMesh
USE moduleOutput
USE moduleMeshInoutCommon
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n, i, fileID
CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
CHARACTER (LEN=iterationDigits):: tstring
TYPE(outputFormat):: output(1:self%numNodes)
fileID = 62
IF (collOutput) THEN
fileName = formatFileName(prefix, 'Collisions', 'vtu', t)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (fileID, file = path // folder // '/' // fileName)
CALL writeHeader(self%numNodes, self%numCells, fileID)
CALL writeMesh(self, fileID)
CALL writeCollOutput(self, fileID)
CALL writeFooter(fileID)
CLOSE(fileID)
!Write collection file for time plotting
fileNameCollection = formatFileName('Collection', 'Collisions', 'pvd')
CALL writeCollection(fileID, t, fileName, filenameCollection)
END IF
END SUBROUTINE printCollVTU
SUBROUTINE printEMVTU(self, t)
USE moduleMesh
USE moduleMeshInoutCommon
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: fileID
CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
fileID = 64
IF (emOutput) THEN
fileName = formatFileName(prefix, 'EMField', 'vtu', t)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (fileID, file = path // folder // '/' // fileName)
CALL writeHeader(self%numNodes, self%numCells, fileID)
CALL writeMesh(self, fileID)
CALL writeEM(self, fileID)
CALL writeFooter(fileID)
CLOSE(fileID)
!Write collection file for time plotting
fileNameCollection = formatFileName('Collection', 'EMField', 'pvd')
CALL writeCollection(fileID, t, fileName, filenameCollection)
END IF
END SUBROUTINE printEMVTU
SUBROUTINE printAverageVTU(self)
USE moduleMesh
USE moduleSpecies
USE moduleMeshInoutCommon
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
INTEGER:: n, i, fileIDMean, fileIDDeviation
CHARACTER(:), ALLOCATABLE:: fileNameMean, fileNameDeviation
TYPE(outputFormat):: output(1:self%numNodes)
fileIDMean = 66
fileIDDeviation = 67
DO i = 1, nSpecies
fileNameMean = formatFileName('Average_mean', species(i)%obj%name, 'vtu')
WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameMean
OPEN (fileIDMean, file = path // folder // '/' // fileNameMean)
fileNameDeviation = formatFileName('Average_deviation', species(i)%obj%name, 'vtu')
WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameDeviation
OPEN (fileIDDeviation, file = path // folder // '/' // fileNameDeviation)
CALL writeHeader(self%numNodes, self%numCells, fileIDMean)
CALL writeHeader(self%numNodes, self%numCells, fileIDDeviation)
CALL writeMesh(self, fileIDMean)
CALL writeMesh(self, fileIDDeviation)
CALL writeAverage(self, fileIDMean, fileIDDeviation, i)
CALL writeFooter(fileIDMean)
CALL writeFooter(fileIDDeviation)
CLOSE(fileIDMean)
CLOSE(fileIDDeviation)
END DO
END SUBROUTINE printAverageVTU
END MODULE moduleMeshOutputVTU

2
src/modules/mesh/moduleMesh.f90
View file

@ -374,7 +374,7 @@ MODULE moduleMesh
!Prints number of collisions in each cell
SUBROUTINE printColl_interface(self, t)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(inout):: self
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printColl_interface