From aca84d631228b0837fe848fe59ff92ef59432eb1 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Fri, 3 Feb 2023 20:14:53 +0100 Subject: [PATCH 01/15] First output in VTU format Testing new VTU format. For now, species information is ALWAYS output in .vtu (to test, this will be an independent format in the future). A .pvd file is produced to do time-series. Still to implement other outputs (electromagnetic, average, collisions...) Still to implement reading a mesh from .vtu file --- src/makefile | 1 + src/modules/makefile | 2 +- src/modules/mesh/inout/makefile | 5 +- src/modules/mesh/inout/vtk/makefile | 7 + .../mesh/inout/vtk/moduleMeshInputVTK.f90 | 3 + .../mesh/inout/vtk/moduleMeshOutputVTK.f90 | 224 ++++++++++++++++++ src/modules/solver/moduleSolver.f90 | 2 + 7 files changed, 242 insertions(+), 2 deletions(-) create mode 100644 src/modules/mesh/inout/vtk/makefile create mode 100644 src/modules/mesh/inout/vtk/moduleMeshInputVTK.f90 create mode 100644 src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 diff --git a/src/makefile b/src/makefile index 2a421f8..32c6569 100644 --- a/src/makefile +++ b/src/makefile @@ -5,6 +5,7 @@ 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)/moduleMeshInputVTK.o $(OBJDIR)/moduleMeshOutputVTK.o \ $(OBJDIR)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \ $(OBJDIR)/moduleMeshInput0D.o $(OBJDIR)/moduleMeshOutput0D.o \ $(OBJDIR)/moduleMesh3DCart.o \ diff --git a/src/modules/makefile b/src/modules/makefile index 254d015..743144f 100644 --- a/src/modules/makefile +++ b/src/modules/makefile @@ -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 diff --git a/src/modules/mesh/inout/makefile b/src/modules/mesh/inout/makefile index 2f73e73..81e15bc 100644 --- a/src/modules/mesh/inout/makefile +++ b/src/modules/mesh/inout/makefile @@ -1,4 +1,7 @@ -all: gmsh2.o 0D.o +all: vtk.o gmsh2.o 0D.o + +vtk.o: + $(MAKE) -C vtk all gmsh2.o: $(MAKE) -C gmsh2 all diff --git a/src/modules/mesh/inout/vtk/makefile b/src/modules/mesh/inout/vtk/makefile new file mode 100644 index 0000000..18a13a0 --- /dev/null +++ b/src/modules/mesh/inout/vtk/makefile @@ -0,0 +1,7 @@ +all: moduleMeshInputVTK.o moduleMeshOutputVTK.o + +moduleMeshInputVTK.o: moduleMeshOutputVTK.o moduleMeshInputVTK.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ + +%.o: %.f90 + $(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@ diff --git a/src/modules/mesh/inout/vtk/moduleMeshInputVTK.f90 b/src/modules/mesh/inout/vtk/moduleMeshInputVTK.f90 new file mode 100644 index 0000000..a044b4f --- /dev/null +++ b/src/modules/mesh/inout/vtk/moduleMeshInputVTK.f90 @@ -0,0 +1,3 @@ +MODULE moduleMeshInputVTK + +END MODULE moduleMeshInputVTK diff --git a/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 b/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 new file mode 100644 index 0000000..efaf6e3 --- /dev/null +++ b/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 @@ -0,0 +1,224 @@ +MODULE moduleMeshOutputVTK + + CONTAINS + + SUBROUTINE writeFileHeader(self, fileID) + USE moduleMesh + IMPLICIT NONE + + CLASS(meshParticles), INTENT(in):: self + INTEGER, INTENT(in):: fileID + + WRITE(fileID,"(A)") '' + WRITE(fileID,"(2X, A)") '' + WRITE(fileID,"(4X, A,ES20.6E3,A)") '' + WRITE(fileID,"(6X, A, I10, A, I10, A)") '' + + END SUBROUTINE writeFileHeader + + SUBROUTINE writeFileFooter(fileID) + IMPLICIT NONE + + INTEGER, INTENT(in):: fileID + + WRITE(fileID,"(6X, A)") '' + WRITE(fileID,"(4X, A)") '' + WRITE(fileID,"(2X, A)") '' + + END SUBROUTINE writeFileFooter + + 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 VTK output', 'getCellType') + + END SELECT + + END FUNCTION getCellType + + SUBROUTINE writeFileMesh(self, fileID) + USE moduleMesh + USE moduleRefParam + IMPLICIT NONE + + CLASS(meshParticles), INTENT(in):: self + INTEGER, INTENT(in):: fileID + CHARACTER(LEN=25):: nodeFormat + CHARACTER(LEN=25):: cellFormat + INTEGER:: e + INTEGER:: offset + + !Write nodes coordinates + WRITE(fileID, "(8X, A)") '' + WRITE(fileID, "(10X,A)") '' + WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))" + WRITE(fileID, nodeFormat) (self%nodes(e)%obj%getCoordinates()*L_ref, e = 1, self%numNodes) + WRITE(fileID, "(10X, A)") '' + WRITE(fileID, "(8X, A)") '' + + WRITE(fileID, "(8X, A)") '' + !Write nodes connectivity of each cell + WRITE(fileID, "(10X,A)") '' + DO e = 1, self%numCells + WRITE(cellFormat, "(A,I1,A)") "(",self%cells(e)%obj%nNodes,"(I10))" + WRITE(fileID, cellFormat, advance="no") self%cells(e)%obj%getNodes(self%cells(e)%obj%nNodes) - 1 !Array starts on 0 + + END DO + WRITE(fileID, "(10X, A)") '' + !Write offset of each cell + offset = 0 + WRITE(fileID, "(10X,A)") '' + DO e = 1, self%numCells + WRITE(cellFormat, "(A,I1,A)") "(I10)" + offset = offset + self%cells(e)%obj%nNodes + WRITE(fileID, cellFormat, advance="no") offset + + END DO + WRITE(fileID, "(10X, A)") '' + !Write type of each cell + WRITE(fileID, "(10X,A)") '' + DO e = 1, self%numCells + WRITE(cellFormat, "(A,I1,A)") "(I10)" + WRITE(fileID, cellFormat, advance="no") getCellType(self%cells(e)%obj) + + END DO + WRITE(fileID, "(10X, A)") '' + WRITE(fileID, "(8X, A)") '' + + END SUBROUTINE writeFileMesh + + SUBROUTINE writeNodeData(self, fileID, output) + USE moduleMesh + USE moduleOutput + IMPLICIT NONE + + CLASS(meshParticles), INTENT(in):: self + INTEGER, INTENT(in):: fileID + TYPE(outputFormat):: output(1:self%numNodes) + CHARACTER(LEN=25):: nodeFormat + INTEGER:: n + + WRITE(fileID,"(A)") '' + WRITE(fileID,"(A)") '' + WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" + WRITE(fileID,nodeFormat) (output(n)%density, n = 1, self%numNodes) + WRITE(fileID,"(A)") '' + WRITE(fileID,"(A)") '' + WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))" + WRITE(fileID,nodeFormat) (output(n)%velocity(1:3), n = 1, self%numNodes) + WRITE(fileID,"(A)") '' + WRITE(fileID,"(A)") '' + WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" + WRITE(fileID,nodeFormat) (output(n)%pressure, n = 1, self%numNodes) + WRITE(fileID,"(A)") '' + WRITE(fileID,"(A)") '' + WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" + WRITE(fileID,nodeFormat) (output(n)%temperature, n = 1, self%numNodes) + WRITE(fileID,"(A)") '' + WRITE(fileID,"(A)") '' + + END SUBROUTINE writeNodeData + + SUBROUTINE printOutputVTK(self,t) + USE moduleMesh + USE moduleRefParam + USE moduleSpecies + USE moduleOutput + USE moduleCaseParam + IMPLICIT NONE + + CLASS(meshParticles), 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 = 60 + + DO i = 1, nSpecies + WRITE(tstring, iterationFormat) t + fileName= 'OUTPUT_' // tstring// '_' // species(i)%obj%name // '.vtu' + + WRITE(*, "(6X,A15,A)") "Creating file: ", fileName + OPEN (fileID, file = path // folder // '/' // fileName) + + fileNameCollection = 'OUTPUT_Collection_' // species(i)%obj%name // '.vtu' + 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)") '' + WRITE (fileID + 1, "(2X, A)") '' + CLOSE(fileID + 1) + + END IF + + OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection, ACCESS='APPEND') + WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)"), '' + + IF (t == tFinal) THEN + WRITE (fileID + 1, "(2X, A)") '' + WRITE (fileID + 1, "(A)") '' + + END IF + CLOSE(fileID + 1) + + CALL writeFileHeader(self, fileID) + + CALL writeFileMesh(self, fileID) + + DO n = 1, self%numNodes + CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj) + + END DO + + CALL writeNodeData(self, fileID, output) + + CALL writeFileFooter(fileID) + + CLOSE(fileID) + + END DO + + END SUBROUTINE printOutputVTK + +END MODULE moduleMeshOutputVTK diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index e557495..813aeb3 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -514,6 +514,7 @@ MODULE moduleSolver USE moduleSpecies USE moduleCompTime USE moduleProbe + USE moduleMeshOutputVTK !TEMPORARY TO TEST VTK OUTPUT IMPLICIT NONE INTEGER, INTENT(in):: t @@ -527,6 +528,7 @@ MODULE moduleSolver CALL outputProbes(t) CALL mesh%printOutput(t) + CALL printOutputVTK(mesh, t) IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t) CALL mesh%printEM(t) WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", t, "/", tFinal From ceab516a5b9659ba538ae0f1dd430e9d67e03b28 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sat, 4 Feb 2023 08:54:36 +0100 Subject: [PATCH 02/15] Correction on collection file The extension of the collection file has to be .pvd so that ParaView properly reads it. --- src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 b/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 index efaf6e3..9b1a442 100644 --- a/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 +++ b/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 @@ -181,7 +181,7 @@ MODULE moduleMeshOutputVTK WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (fileID, file = path // folder // '/' // fileName) - fileNameCollection = 'OUTPUT_Collection_' // species(i)%obj%name // '.vtu' + fileNameCollection = 'OUTPUT_Collection_' // species(i)%obj%name // '.pvd' IF (t == tInitial) THEN !Create collection file WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection @@ -192,9 +192,11 @@ MODULE moduleMeshOutputVTK 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)"), '' + !Close collection file IF (t == tFinal) THEN WRITE (fileID + 1, "(2X, A)") '' WRITE (fileID + 1, "(A)") '' From f1c0c5755f6922687025b8ebdd540176661a57fc Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sat, 4 Feb 2023 12:31:33 +0100 Subject: [PATCH 03/15] Collisions and EM field in .vtu The collisions and EM field information is now available in .vtu files. A collection file .pvd is provided per dataset for time-dependent plotting. Still to do: Write average quantities in .vtu Read mesh from .vtu --- .../mesh/inout/vtk/moduleMeshOutputVTK.f90 | 290 ++++++++++++++---- src/modules/solver/moduleSolver.f90 | 4 +- 2 files changed, 233 insertions(+), 61 deletions(-) diff --git a/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 b/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 index 9b1a442..203de0c 100644 --- a/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 +++ b/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 @@ -1,22 +1,44 @@ MODULE moduleMeshOutputVTK + CHARACTER(LEN=6):: prefix = 'OUTPUT' + CONTAINS - SUBROUTINE writeFileHeader(self, fileID) + 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 + + SUBROUTINE writeHeader(nNodes, nCells, fileID) USE moduleMesh IMPLICIT NONE - CLASS(meshParticles), INTENT(in):: self + INTEGER, INTENT(in):: nNodes, nCells INTEGER, INTENT(in):: fileID WRITE(fileID,"(A)") '' WRITE(fileID,"(2X, A)") '' WRITE(fileID,"(4X, A,ES20.6E3,A)") '' - WRITE(fileID,"(6X, A, I10, A, I10, A)") '' + WRITE(fileID,"(6X, A, I10, A, I10, A)") '' - END SUBROUTINE writeFileHeader + END SUBROUTINE writeHeader - SUBROUTINE writeFileFooter(fileID) + SUBROUTINE writeFooter(fileID) IMPLICIT NONE INTEGER, INTENT(in):: fileID @@ -25,7 +47,7 @@ MODULE moduleMeshOutputVTK WRITE(fileID,"(4X, A)") '' WRITE(fileID,"(2X, A)") '' - END SUBROUTINE writeFileFooter + END SUBROUTINE writeFooter FUNCTION getCellType(cell) RESULT(indexType) USE moduleMesh3DCart @@ -74,12 +96,12 @@ MODULE moduleMeshOutputVTK END FUNCTION getCellType - SUBROUTINE writeFileMesh(self, fileID) + SUBROUTINE writeMesh(self, fileID) USE moduleMesh USE moduleRefParam IMPLICIT NONE - CLASS(meshParticles), INTENT(in):: self + CLASS(meshGeneric), INTENT(in):: self INTEGER, INTENT(in):: fileID CHARACTER(LEN=25):: nodeFormat CHARACTER(LEN=25):: cellFormat @@ -123,104 +145,252 @@ MODULE moduleMeshOutputVTK WRITE(fileID, "(10X, A)") '' WRITE(fileID, "(8X, A)") '' - END SUBROUTINE writeFileMesh + END SUBROUTINE writeMesh - SUBROUTINE writeNodeData(self, fileID, output) + 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) CHARACTER(LEN=25):: nodeFormat INTEGER:: n - WRITE(fileID,"(A)") '' - WRITE(fileID,"(A)") '' + 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(fileID,"(8X,A)") '' + WRITE(fileID,"(10X,A)") '' WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" WRITE(fileID,nodeFormat) (output(n)%density, n = 1, self%numNodes) - WRITE(fileID,"(A)") '' - WRITE(fileID,"(A)") '' + WRITE(fileID,"(10X,A)") '' + WRITE(fileID,"(10X,A)") '' WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))" WRITE(fileID,nodeFormat) (output(n)%velocity(1:3), n = 1, self%numNodes) - WRITE(fileID,"(A)") '' - WRITE(fileID,"(A)") '' + WRITE(fileID,"(10X,A)") '' + WRITE(fileID,"(10X,A)") '' WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" WRITE(fileID,nodeFormat) (output(n)%pressure, n = 1, self%numNodes) - WRITE(fileID,"(A)") '' - WRITE(fileID,"(A)") '' + WRITE(fileID,"(10X,A)") '' + WRITE(fileID,"(10X,A)") '' WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" WRITE(fileID,nodeFormat) (output(n)%temperature, n = 1, self%numNodes) - WRITE(fileID,"(A)") '' - WRITE(fileID,"(A)") '' + WRITE(fileID,"(10X,A)") '' + WRITE(fileID,"(8X,A)") '' - END SUBROUTINE writeNodeData + 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 + CHARACTER(LEN=25):: cellFormat + + WRITE(fileID,"(8X,A)") '' + 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)") '' + WRITE(cellFormat, "(A,I10, A)") "(", self%numCells, "(I10))" + WRITE(fileID, cellFormat) (self%cells(n)%obj%tallyColl(k)%tally(c), n = 1, self%numCells) + WRITE(fileID, "(10X, A)") '' + + END DO + END DO + WRITE(fileID,"(8X,A)") '' + + END SUBROUTINE writeCollOutput + + SUBROUTINE writeEM(self, fileID) + USE moduleMesh + USE moduleRefParam + IMPLICIT NONE + + CLASS(meshParticles), INTENT(in):: self + INTEGER, INTENT(in):: fileID + CHARACTER(LEN=25):: nodeFormat + CHARACTER(LEN=25):: cellFormat + INTEGER:: n + REAL(8):: Xi(1:3) + + Xi = (/ 0.D0, 0.D0, 0.D0 /) + + !Points in nodes + WRITE(fileID,"(8X,A)") '' + !Electric potential + WRITE(fileID,"(10X,A)") '' + WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" + WRITE(fileID,nodeFormat) (self%nodes(n)%obj%emData%phi*Volt_ref, n = 1, self%numNodes) + WRITE(fileID,"(10X,A)") '' + !Magnetic Field + WRITE(fileID,"(10X,A)") '' + WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))" + WRITE(fileID,nodeFormat) (self%nodes(n)%obj%emData%B*B_ref, n = 1, self%numNodes) + WRITE(fileID,"(10X,A)") '' + WRITE(fileID,"(8X,A)") '' + + !Cell Data + WRITE(fileID,"(8X,A)") '' + !Electric field + WRITE(fileID,"(10X,A, A, A)") '' + WRITE(cellFormat, "(A,I10, A)") "(", self%numCells, "(3(ES20.6E3)))" + WRITE(fileID, cellFormat) (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells) + WRITE(fileID,"(10X,A)") '' + WRITE(fileID,"(8X,A)") '' + + 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)") '' + WRITE (fileID + 1, "(2X, A)") '' + 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)") '' + + !Close collection file + IF (t == tFinal) THEN + WRITE (fileID + 1, "(2X, A)") '' + WRITE (fileID + 1, "(A)") '' + + END IF + CLOSE(fileID + 1) + + END SUBROUTINE writeCollection SUBROUTINE printOutputVTK(self,t) USE moduleMesh - USE moduleRefParam USE moduleSpecies - USE moduleOutput - USE moduleCaseParam IMPLICIT NONE CLASS(meshParticles), 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 = 60 DO i = 1, nSpecies - WRITE(tstring, iterationFormat) t - fileName= 'OUTPUT_' // tstring// '_' // species(i)%obj%name // '.vtu' - + fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', t) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (fileID, file = path // folder // '/' // fileName) - fileNameCollection = 'OUTPUT_Collection_' // species(i)%obj%name // '.pvd' - 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)") '' - WRITE (fileID + 1, "(2X, A)") '' - CLOSE(fileID + 1) + CALL writeHeader(self%numNodes, self%numCells, fileID) - END IF + CALL writeMesh(self, fileID) - !Write iteration file in collection - OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection, ACCESS='APPEND') - WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)"), '' + CALL writeSpeciesOutput(self, fileID, i) - !Close collection file - IF (t == tFinal) THEN - WRITE (fileID + 1, "(2X, A)") '' - WRITE (fileID + 1, "(A)") '' - - END IF - CLOSE(fileID + 1) - - CALL writeFileHeader(self, fileID) - - CALL writeFileMesh(self, fileID) - - DO n = 1, self%numNodes - CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj) - - END DO - - CALL writeNodeData(self, fileID, output) - - CALL writeFileFooter(fileID) + CALL writeFooter(fileID) CLOSE(fileID) + !Write collection file for time plotting + fileNameCollection = formatFileName(prefix, 'Collection_' // species(i)%obj%name, 'pvd') + CALL writeCollection(fileID, t, fileName, filenameCollection) + END DO END SUBROUTINE printOutputVTK + SUBROUTINE printCollVTK(self,t) + USE moduleMesh + USE moduleOutput + 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(prefix, 'Collection_Collisions', 'pvd') + CALL writeCollection(fileID, t, fileName, filenameCollection) + + END IF + + END SUBROUTINE printCollVTK + + SUBROUTINE printEMVTK(self, t) + USE moduleMesh + 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(prefix, 'Collection_EMField', 'pvd') + CALL writeCollection(fileID, t, fileName, filenameCollection) + + END IF + + END SUBROUTINE printEMVTK + END MODULE moduleMeshOutputVTK diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index 813aeb3..def5ca4 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -528,9 +528,11 @@ MODULE moduleSolver CALL outputProbes(t) CALL mesh%printOutput(t) - CALL printOutputVTK(mesh, t) + CALL printOutputVTK(mesh, t) !TEMPORARY TO TEST VTK OUTPUT IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t) + IF (ASSOCIATED(meshForMCC)) CALL printCollVTK(meshForMCC,t) !TEMPORARY TO TEST VTK OUTPUT CALL mesh%printEM(t) + CALL printEMVTK(mesh, t) !TEMPORARY TO TEST VTK OUTPUT WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", t, "/", tFinal WRITE(*, "(5X,A21,I10)") "Particles: ", nPartOld IF (t == 0) THEN From 6706c5dd1ca26f9dd69648e4d7fe8a5ca903a6c2 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sat, 4 Feb 2023 15:20:36 +0100 Subject: [PATCH 04/15] Average is written in .vtu The average of the species properties can be written now in .vtu format. No .pvd file is provided as no time series is generated. Still to do: Read a .vtu mesh. Improve gmsh format to use more common functions. --- .../mesh/inout/vtk/moduleMeshOutputVTK.f90 | 111 +++++++++++++++++- src/modules/solver/moduleSolver.f90 | 1 + 2 files changed, 111 insertions(+), 1 deletion(-) diff --git a/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 b/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 index 203de0c..d962ae9 100644 --- a/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 +++ b/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 @@ -164,23 +164,29 @@ MODULE moduleMeshOutputVTK END DO + !Write node data WRITE(fileID,"(8X,A)") '' + !Write density WRITE(fileID,"(10X,A)") '' WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" WRITE(fileID,nodeFormat) (output(n)%density, n = 1, self%numNodes) WRITE(fileID,"(10X,A)") '' + !Write velocity WRITE(fileID,"(10X,A)") '' WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))" - WRITE(fileID,nodeFormat) (output(n)%velocity(1:3), n = 1, self%numNodes) + WRITE(fileID,nodeFormat) (output(n)%velocity, n = 1, self%numNodes) WRITE(fileID,"(10X,A)") '' + !Write pressure WRITE(fileID,"(10X,A)") '' WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" WRITE(fileID,nodeFormat) (output(n)%pressure, n = 1, self%numNodes) WRITE(fileID,"(10X,A)") '' + !Write temperature WRITE(fileID,"(10X,A)") '' WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" WRITE(fileID,nodeFormat) (output(n)%temperature, n = 1, self%numNodes) WRITE(fileID,"(10X,A)") '' + !End node data WRITE(fileID,"(8X,A)") '' END SUBROUTINE writeSpeciesOutput @@ -286,6 +292,69 @@ MODULE moduleMeshOutputVTK 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) + CHARACTER(LEN=25):: nodeFormat + 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)") '' + WRITE(fileIDDeviation,"(8X,A)") '' + !Write density + WRITE(fileIDMean, "(10X,A)") '' + WRITE(fileIDDeviation,"(10X,A)") '' + WRITE(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" + WRITE(fileIDMean, nodeFormat) (outputMean(n)%density, n = 1, self%numNodes) + WRITE(fileIDDeviation,nodeFormat) (outputDeviation(n)%density, n = 1, self%numNodes) + WRITE(fileIDMean,"(10X,A)") '' + WRITE(fileIDDeviation,"(10X,A)") '' + !Write velocity + WRITE(fileIDMean, "(10X,A)") '' + WRITE(fileIDDeviation,"(10X,A)") '' + WRITE(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))" + WRITE(fileIDMean, nodeFormat) (outputMean(n)%velocity, n = 1, self%numNodes) + WRITE(fileIDDeviation,nodeFormat) (outputDeviation(n)%velocity, n = 1, self%numNodes) + WRITE(fileIDMean, "(10X,A)") '' + WRITE(fileIDDeviation,"(10X,A)") '' + !Write pressure + WRITE(fileIDMean, "(10X,A)") '' + WRITE(fileIDDeviation,"(10X,A)") '' + WRITE(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" + WRITE(fileIDMean, nodeFormat) (outputMean(n)%pressure, n = 1, self%numNodes) + WRITE(fileIDDeviation,nodeFormat) (outputDeviation(n)%pressure, n = 1, self%numNodes) + WRITE(fileIDMean, "(10X,A)") '' + WRITE(fileIDDeviation,"(10X,A)") '' + !Write temperature + WRITE(fileIDMean, "(10X,A)") '' + WRITE(fileIDDeviation,"(10X,A)") '' + WRITE(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" + WRITE(fileIDMean, nodeFormat) (outputMean(n)%temperature, n = 1, self%numNodes) + WRITE(fileIDDeviation,nodeFormat) (outputDeviation(n)%temperature, n = 1, self%numNodes) + WRITE(fileIDMean, "(10X,A)") '' + WRITE(fileIDDeviation,"(10X,A)") '' + !End node data + WRITE(fileIDMean, "(8X,A)") '' + WRITE(fileIDDeviation,"(8X,A)") '' + + END SUBROUTINE writeAverage + SUBROUTINE printOutputVTK(self,t) USE moduleMesh USE moduleSpecies @@ -393,4 +462,44 @@ MODULE moduleMeshOutputVTK END SUBROUTINE printEMVTK + SUBROUTINE printAverageVTK(self) + USE moduleMesh + USE moduleSpecies + 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 printAverageVTK + END MODULE moduleMeshOutputVTK diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index def5ca4..68c1270 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -565,6 +565,7 @@ MODULE moduleSolver !Output average values IF (useAverage .AND. t == tFinal) THEN CALL mesh%printAverage() + CALL printAverageVTK(mesh) !TEMPORARY TO TEST VTK OUTPUT END IF From f5be04587a82a49b5a5705f8deacf7f53d6fe9d9 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sat, 4 Feb 2023 15:41:13 +0100 Subject: [PATCH 05/15] First step towards reading .vtu mesh Just setting up the required functions. --- src/makefile | 2 +- .../mesh/inout/0D/moduleMeshOutput0D.f90 | 4 +- .../mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 | 4 +- .../inout/gmsh2/moduleMeshOutputGmsh2.f90 | 2 +- src/modules/mesh/inout/makefile | 6 +-- .../mesh/inout/vtk/moduleMeshInputVTK.f90 | 3 -- src/modules/mesh/inout/{vtk => vtu}/makefile | 4 +- .../mesh/inout/vtu/moduleMeshInputVTU.f90 | 50 +++++++++++++++++++ .../moduleMeshOutputVTU.f90} | 22 ++++---- src/modules/mesh/moduleMesh.f90 | 2 +- src/modules/solver/moduleSolver.f90 | 10 ++-- 11 files changed, 78 insertions(+), 31 deletions(-) delete mode 100644 src/modules/mesh/inout/vtk/moduleMeshInputVTK.f90 rename src/modules/mesh/inout/{vtk => vtu}/makefile (50%) create mode 100644 src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 rename src/modules/mesh/inout/{vtk/moduleMeshOutputVTK.f90 => vtu/moduleMeshOutputVTU.f90} (97%) diff --git a/src/makefile b/src/makefile index 32c6569..3223e41 100644 --- a/src/makefile +++ b/src/makefile @@ -5,7 +5,7 @@ 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)/moduleMeshInputVTK.o $(OBJDIR)/moduleMeshOutputVTK.o \ + $(OBJDIR)/moduleMeshInputVTU.o $(OBJDIR)/moduleMeshOutputVTU.o \ $(OBJDIR)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \ $(OBJDIR)/moduleMeshInput0D.o $(OBJDIR)/moduleMeshOutput0D.o \ $(OBJDIR)/moduleMesh3DCart.o \ diff --git a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 index dfe9605..97ec729 100644 --- a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 +++ b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 @@ -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 diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 index ae1fe05..6766c49 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 @@ -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 diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 index 53485f4..52e3922 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 @@ -144,7 +144,7 @@ MODULE moduleMeshOutputGmsh2 USE moduleOutput IMPLICIT NONE - CLASS(meshGeneric), INTENT(inout):: self + CLASS(meshGeneric), INTENT(in):: self INTEGER, INTENT(in):: t INTEGER:: numEdges INTEGER:: k, c diff --git a/src/modules/mesh/inout/makefile b/src/modules/mesh/inout/makefile index 81e15bc..e686b87 100644 --- a/src/modules/mesh/inout/makefile +++ b/src/modules/mesh/inout/makefile @@ -1,7 +1,7 @@ -all: vtk.o gmsh2.o 0D.o +all: vtu.o gmsh2.o 0D.o -vtk.o: - $(MAKE) -C vtk all +vtu.o: + $(MAKE) -C vtu all gmsh2.o: $(MAKE) -C gmsh2 all diff --git a/src/modules/mesh/inout/vtk/moduleMeshInputVTK.f90 b/src/modules/mesh/inout/vtk/moduleMeshInputVTK.f90 deleted file mode 100644 index a044b4f..0000000 --- a/src/modules/mesh/inout/vtk/moduleMeshInputVTK.f90 +++ /dev/null @@ -1,3 +0,0 @@ -MODULE moduleMeshInputVTK - -END MODULE moduleMeshInputVTK diff --git a/src/modules/mesh/inout/vtk/makefile b/src/modules/mesh/inout/vtu/makefile similarity index 50% rename from src/modules/mesh/inout/vtk/makefile rename to src/modules/mesh/inout/vtu/makefile index 18a13a0..07b471d 100644 --- a/src/modules/mesh/inout/vtk/makefile +++ b/src/modules/mesh/inout/vtu/makefile @@ -1,6 +1,6 @@ -all: moduleMeshInputVTK.o moduleMeshOutputVTK.o +all: moduleMeshInputVTU.o moduleMeshOutputVTU.o -moduleMeshInputVTK.o: moduleMeshOutputVTK.o moduleMeshInputVTK.f90 +moduleMeshInputVTU.o: moduleMeshOutputVTU.o moduleMeshInputVTU.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ %.o: %.f90 diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 new file mode 100644 index 0000000..da6c792 --- /dev/null +++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 @@ -0,0 +1,50 @@ +MODULE moduleMeshInputVTK + !Reads mesh in the VTU format + + CONTAINS + SUBROUTINE initVTK(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 initVTK + + 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 + + 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 + + END SUBROUTINE readInitialVTU + +END MODULE moduleMeshInputVTK diff --git a/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 similarity index 97% rename from src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 rename to src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 index d962ae9..bf4b1b4 100644 --- a/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 @@ -1,4 +1,4 @@ -MODULE moduleMeshOutputVTK +MODULE moduleMeshOutputVTU CHARACTER(LEN=6):: prefix = 'OUTPUT' @@ -90,7 +90,7 @@ MODULE moduleMeshOutputVTK indexType = 1 CLASS DEFAULT - CALL criticalError('Cell not valid for VTK output', 'getCellType') + CALL criticalError('Cell not valid for VTU output', 'getCellType') END SELECT @@ -355,7 +355,7 @@ MODULE moduleMeshOutputVTK END SUBROUTINE writeAverage - SUBROUTINE printOutputVTK(self,t) + SUBROUTINE printOutputVTU(self,t) USE moduleMesh USE moduleSpecies IMPLICIT NONE @@ -389,9 +389,9 @@ MODULE moduleMeshOutputVTK END DO - END SUBROUTINE printOutputVTK + END SUBROUTINE printOutputVTU - SUBROUTINE printCollVTK(self,t) + SUBROUTINE printCollVTU(self,t) USE moduleMesh USE moduleOutput IMPLICIT NONE @@ -426,9 +426,9 @@ MODULE moduleMeshOutputVTK END IF - END SUBROUTINE printCollVTK + END SUBROUTINE printCollVTU - SUBROUTINE printEMVTK(self, t) + SUBROUTINE printEMVTU(self, t) USE moduleMesh IMPLICIT NONE @@ -460,9 +460,9 @@ MODULE moduleMeshOutputVTK END IF - END SUBROUTINE printEMVTK + END SUBROUTINE printEMVTU - SUBROUTINE printAverageVTK(self) + SUBROUTINE printAverageVTU(self) USE moduleMesh USE moduleSpecies IMPLICIT NONE @@ -500,6 +500,6 @@ MODULE moduleMeshOutputVTK END DO - END SUBROUTINE printAverageVTK + END SUBROUTINE printAverageVTU -END MODULE moduleMeshOutputVTK +END MODULE moduleMeshOutputVTU diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 7482842..cb252e8 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -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 diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index 68c1270..b1a91fd 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -514,7 +514,7 @@ MODULE moduleSolver USE moduleSpecies USE moduleCompTime USE moduleProbe - USE moduleMeshOutputVTK !TEMPORARY TO TEST VTK OUTPUT + USE moduleMeshOutputVTU !TEMPORARY TO TEST VTU OUTPUT IMPLICIT NONE INTEGER, INTENT(in):: t @@ -528,11 +528,11 @@ MODULE moduleSolver CALL outputProbes(t) CALL mesh%printOutput(t) - CALL printOutputVTK(mesh, t) !TEMPORARY TO TEST VTK OUTPUT + CALL printOutputVTU(mesh, t) !TEMPORARY TO TEST VTU OUTPUT IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t) - IF (ASSOCIATED(meshForMCC)) CALL printCollVTK(meshForMCC,t) !TEMPORARY TO TEST VTK OUTPUT + IF (ASSOCIATED(meshForMCC)) CALL printCollVTU(meshForMCC,t) !TEMPORARY TO TEST VTU OUTPUT CALL mesh%printEM(t) - CALL printEMVTK(mesh, t) !TEMPORARY TO TEST VTK OUTPUT + CALL printEMVTU(mesh, t) !TEMPORARY TO TEST VTU OUTPUT WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", t, "/", tFinal WRITE(*, "(5X,A21,I10)") "Particles: ", nPartOld IF (t == 0) THEN @@ -565,7 +565,7 @@ MODULE moduleSolver !Output average values IF (useAverage .AND. t == tFinal) THEN CALL mesh%printAverage() - CALL printAverageVTK(mesh) !TEMPORARY TO TEST VTK OUTPUT + CALL printAverageVTU(mesh) !TEMPORARY TO TEST VTU OUTPUT END IF From 63fd8fdb91557bbf95ce729956c1986b5379cc53 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sun, 5 Feb 2023 16:23:37 +0100 Subject: [PATCH 06/15] New functions to read VTU Added a few functions to read .vtu meshes. Output in .vtu changed to output data in 6 columns (seems to be the standard.) --- src/modules/init/moduleInput.f90 | 4 + .../mesh/inout/vtu/moduleMeshInputVTU.f90 | 166 +++++++++++++++++- .../mesh/inout/vtu/moduleMeshOutputVTU.f90 | 78 +++----- 3 files changed, 192 insertions(+), 56 deletions(-) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 0a33903..82d458e 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -860,6 +860,7 @@ MODULE moduleInput SUBROUTINE readGeometry(config) USE moduleMesh USE moduleMeshInputGmsh2, ONLY: initGmsh2 + USE moduleMeshInputVTU, ONLY: initVTU, readVTU !TEMPORARY TO TEST VTU OUTPUT USE moduleMeshInput0D, ONLY: init0D USE moduleErrors USE moduleOutput @@ -873,6 +874,7 @@ MODULE moduleInput LOGICAL:: found CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile REAL(8):: volume + CHARACTER(:), ALLOCATABLE:: meshFileVTU !Temporary to test VTU OUTPUT object = 'geometry' @@ -973,6 +975,8 @@ MODULE moduleInput CALL config%get(object // '.meshFile', meshFile, found) pathMeshParticle = path // meshFile CALL mesh%readMesh(pathMeshParticle) + meshFileVTU = '/home/jorge/fpakc/runs/cylFlow/mesh.vtu' !TEMPORARY TO TEST VTU OUTPUT + CALL readVTU(mesh, meshFileVTU) !TEMPORARY TO TEST VTU OUTPUT DEALLOCATE(meshFile) IF (doubleMesh) THEN !Reads the mesh file for collisions diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 index da6c792..505f10d 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 @@ -1,8 +1,138 @@ -MODULE moduleMeshInputVTK +MODULE moduleMeshInputVTU !Reads mesh in the VTU format + INTERFACE getValueFromLine + MODULE PROCEDURE getIntegerFromLine, getRealFromLine + + END INTERFACE + CONTAINS - SUBROUTINE initVTK(self) + 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 + PRINT *, iStart, iEnd + 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 + PRINT *, iStart, iEnd + 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 @@ -20,7 +150,7 @@ MODULE moduleMeshInputVTK END SELECT self%readMesh => readVTU - END SUBROUTINE initVTK + END SUBROUTINE initVTU SUBROUTINE readVTU(self, filename) USE moduleMesh3DCart @@ -34,6 +164,34 @@ MODULE moduleMeshInputVTK 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 + INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets + + fileID = 10 + + OPEN(fileID, FILE=TRIM(filename)) + + !Find the number of nodes and elements (edges+cells) in the mesh. + line = findLine(fileID, '' WRITE(fileID, "(10X,A)") '' - WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))" - WRITE(fileID, nodeFormat) (self%nodes(e)%obj%getCoordinates()*L_ref, e = 1, self%numNodes) + WRITE(fileID, "(6(ES20.6E3))") (self%nodes(e)%obj%getCoordinates()*L_ref, e = 1, self%numNodes) WRITE(fileID, "(10X, A)") '' WRITE(fileID, "(8X, A)") '' WRITE(fileID, "(8X, A)") '' !Write nodes connectivity of each cell WRITE(fileID, "(10X,A)") '' - DO e = 1, self%numCells - WRITE(cellFormat, "(A,I1,A)") "(",self%cells(e)%obj%nNodes,"(I10))" - WRITE(fileID, cellFormat, advance="no") self%cells(e)%obj%getNodes(self%cells(e)%obj%nNodes) - 1 !Array starts on 0 - - END DO + 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)") '' !Write offset of each cell - offset = 0 + offset(1) = self%cells(1)%obj%nNodes WRITE(fileID, "(10X,A)") '' - DO e = 1, self%numCells - WRITE(cellFormat, "(A,I1,A)") "(I10)" - offset = offset + self%cells(e)%obj%nNodes - WRITE(fileID, cellFormat, advance="no") offset + 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)") '' !Write type of each cell WRITE(fileID, "(10X,A)") '' DO e = 1, self%numCells - WRITE(cellFormat, "(A,I1,A)") "(I10)" - WRITE(fileID, cellFormat, advance="no") getCellType(self%cells(e)%obj) + types(e) = getCellType(self%cells(e)%obj) END DO + WRITE(fileID, "(6(I10))") types WRITE(fileID, "(10X, A)") '' WRITE(fileID, "(8X, A)") '' @@ -156,7 +148,6 @@ MODULE moduleMeshOutputVTU INTEGER, INTENT(in):: fileID INTEGER, INTENT(in):: speciesIndex TYPE(outputFormat):: output(1:self%numNodes) - CHARACTER(LEN=25):: nodeFormat INTEGER:: n DO n = 1, self%numNodes @@ -168,23 +159,19 @@ MODULE moduleMeshOutputVTU WRITE(fileID,"(8X,A)") '' !Write density WRITE(fileID,"(10X,A)") '' - WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" - WRITE(fileID,nodeFormat) (output(n)%density, n = 1, self%numNodes) + WRITE(fileID,"(6(ES20.6E3))") (output(n)%density, n = 1, self%numNodes) WRITE(fileID,"(10X,A)") '' !Write velocity WRITE(fileID,"(10X,A)") '' - WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))" - WRITE(fileID,nodeFormat) (output(n)%velocity, n = 1, self%numNodes) + WRITE(fileID,"(6(ES20.6E3))") (output(n)%velocity, n = 1, self%numNodes) WRITE(fileID,"(10X,A)") '' !Write pressure WRITE(fileID,"(10X,A)") '' - WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" - WRITE(fileID,nodeFormat) (output(n)%pressure, n = 1, self%numNodes) + WRITE(fileID,"(6(ES20.6E3))") (output(n)%pressure, n = 1, self%numNodes) WRITE(fileID,"(10X,A)") '' !Write temperature WRITE(fileID,"(10X,A)") '' - WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" - WRITE(fileID,nodeFormat) (output(n)%temperature, n = 1, self%numNodes) + WRITE(fileID,"(6(ES20.6E3))") (output(n)%temperature, n = 1, self%numNodes) WRITE(fileID,"(10X,A)") '' !End node data WRITE(fileID,"(8X,A)") '' @@ -201,7 +188,6 @@ MODULE moduleMeshOutputVTU INTEGER:: k, c, n CHARACTER(:), ALLOCATABLE:: title CHARACTER (LEN=2):: cString - CHARACTER(LEN=25):: cellFormat WRITE(fileID,"(8X,A)") '' DO k = 1, nCollPairs @@ -209,8 +195,7 @@ MODULE moduleMeshOutputVTU WRITE(cString, "(I2)") c title = 'Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString WRITE(fileID,"(10X,A, A, A)") '' - WRITE(cellFormat, "(A,I10, A)") "(", self%numCells, "(I10))" - WRITE(fileID, cellFormat) (self%cells(n)%obj%tallyColl(k)%tally(c), n = 1, self%numCells) + WRITE(fileID, "(6(I10))") (self%cells(n)%obj%tallyColl(k)%tally(c), n = 1, self%numCells) WRITE(fileID, "(10X, A)") '' END DO @@ -226,33 +211,27 @@ MODULE moduleMeshOutputVTU CLASS(meshParticles), INTENT(in):: self INTEGER, INTENT(in):: fileID - CHARACTER(LEN=25):: nodeFormat - CHARACTER(LEN=25):: cellFormat INTEGER:: n REAL(8):: Xi(1:3) - Xi = (/ 0.D0, 0.D0, 0.D0 /) - !Points in nodes WRITE(fileID,"(8X,A)") '' !Electric potential WRITE(fileID,"(10X,A)") '' - WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" - WRITE(fileID,nodeFormat) (self%nodes(n)%obj%emData%phi*Volt_ref, n = 1, self%numNodes) + WRITE(fileID,"(6(ES20.6E3))") (self%nodes(n)%obj%emData%phi*Volt_ref, n = 1, self%numNodes) WRITE(fileID,"(10X,A)") '' !Magnetic Field WRITE(fileID,"(10X,A)") '' - WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))" - WRITE(fileID,nodeFormat) (self%nodes(n)%obj%emData%B*B_ref, n = 1, self%numNodes) + WRITE(fileID,"(6(ES20.6E3))") (self%nodes(n)%obj%emData%B*B_ref, n = 1, self%numNodes) WRITE(fileID,"(10X,A)") '' WRITE(fileID,"(8X,A)") '' !Cell Data + Xi = (/ 0.D0, 0.D0, 0.D0 /) WRITE(fileID,"(8X,A)") '' !Electric field WRITE(fileID,"(10X,A, A, A)") '' - WRITE(cellFormat, "(A,I10, A)") "(", self%numCells, "(3(ES20.6E3)))" - WRITE(fileID, cellFormat) (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells) + WRITE(fileID, "(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells) WRITE(fileID,"(10X,A)") '' WRITE(fileID,"(8X,A)") '' @@ -303,7 +282,6 @@ MODULE moduleMeshOutputVTU INTEGER, INTENT(in):: speciesIndex TYPE(outputFormat):: outputMean(1:self%numNodes) TYPE(outputFormat):: outputDeviation(1:self%numNodes) - CHARACTER(LEN=25):: nodeFormat INTEGER:: n DO n = 1, self%numNodes @@ -320,33 +298,29 @@ MODULE moduleMeshOutputVTU !Write density WRITE(fileIDMean, "(10X,A)") '' WRITE(fileIDDeviation,"(10X,A)") '' - WRITE(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" - WRITE(fileIDMean, nodeFormat) (outputMean(n)%density, n = 1, self%numNodes) - WRITE(fileIDDeviation,nodeFormat) (outputDeviation(n)%density, n = 1, self%numNodes) + 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)") '' WRITE(fileIDDeviation,"(10X,A)") '' !Write velocity WRITE(fileIDMean, "(10X,A)") '' WRITE(fileIDDeviation,"(10X,A)") '' - WRITE(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))" - WRITE(fileIDMean, nodeFormat) (outputMean(n)%velocity, n = 1, self%numNodes) - WRITE(fileIDDeviation,nodeFormat) (outputDeviation(n)%velocity, n = 1, self%numNodes) + 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)") '' WRITE(fileIDDeviation,"(10X,A)") '' !Write pressure WRITE(fileIDMean, "(10X,A)") '' WRITE(fileIDDeviation,"(10X,A)") '' - WRITE(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" - WRITE(fileIDMean, nodeFormat) (outputMean(n)%pressure, n = 1, self%numNodes) - WRITE(fileIDDeviation,nodeFormat) (outputDeviation(n)%pressure, n = 1, self%numNodes) + 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)") '' WRITE(fileIDDeviation,"(10X,A)") '' !Write temperature WRITE(fileIDMean, "(10X,A)") '' WRITE(fileIDDeviation,"(10X,A)") '' - WRITE(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" - WRITE(fileIDMean, nodeFormat) (outputMean(n)%temperature, n = 1, self%numNodes) - WRITE(fileIDDeviation,nodeFormat) (outputDeviation(n)%temperature, n = 1, self%numNodes) + 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)") '' WRITE(fileIDDeviation,"(10X,A)") '' !End node data From 7b470b7f584e98c3f7928c2bd9b9b47806e814c1 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sun, 5 Feb 2023 18:32:38 +0100 Subject: [PATCH 07/15] Subroutines to read .vtu information All subroutines to read .vtu information is ready. Now it is time to create the input and generate a mesh for fpakc. --- .../mesh/inout/vtu/moduleMeshInputVTU.f90 | 33 ++++++++++++++++--- 1 file changed, 29 insertions(+), 4 deletions(-) diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 index 505f10d..b283d67 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 @@ -6,6 +6,11 @@ MODULE moduleMeshInputVTU END INTERFACE + INTERFACE readDataBlock + MODULE PROCEDURE readIntegerBlock, readRealBlock + + END INTERFACE + CONTAINS FUNCTION findLine(fileID, text) RESULT(line) USE moduleErrors @@ -93,7 +98,6 @@ MODULE moduleMeshInputVTU DO WHILE (iStart < nData) iStart = iStart + 1 iEnd = iStart - 1 + block - PRINT *, iStart, iEnd IF (iEnd > nData) THEN iEnd = nData @@ -167,7 +171,8 @@ MODULE moduleMeshInputVTU INTEGER:: fileID, error, found CHARACTER(LEN=256):: line INTEGER:: numNodes, numElements - INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets + INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types + REAL(8), ALLOCATABLE, DIMENSION(:):: coordinates fileID = 10 @@ -182,17 +187,37 @@ MODULE moduleMeshInputVTU !Get the IDs of the cells to identify physical surfaces line = findLine(fileID, 'Name="CellEntityIds"') ALLOCATE(entitiesID(1:numElements)) - CALL readIntegerBlock(fileID, numElements, entitiesID) + CALL readDataBlock(fileID, numElements, entitiesID) REWIND(fileID) !Get the offsets to read connectivity line = findLine(fileID, 'Name="offsets"') ALLOCATE(offsets(1:numElements)) - CALL readIntegerBlock(fileID, numElements, offsets) + 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) + 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. + END SUBROUTINE readVTU SUBROUTINE readInitialVTU(filename, density, velocity, temperature) From 43a74217956507c81535d35d313150c911aa2285 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sun, 5 Feb 2023 19:35:49 +0100 Subject: [PATCH 08/15] Creating of nodes and edges in .vtu format Moving forward making vtu an independent format. Now fpakc can generate nodes and edges from vtu input. Next step is cells. Some minor corrections in gmsh2 format to unify statements. The reading of meshes needs a good overhaul. Testing all geometries with vtu is gonna be fun... --- src/modules/init/moduleInput.f90 | 48 +++++- .../mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 | 47 ++--- .../mesh/inout/vtu/moduleMeshInputVTU.f90 | 163 +++++++++++++++++- src/modules/solver/moduleSolver.f90 | 5 - 4 files changed, 220 insertions(+), 43 deletions(-) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 82d458e..33bd894 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -860,8 +860,13 @@ MODULE moduleInput SUBROUTINE readGeometry(config) USE moduleMesh USE moduleMeshInputGmsh2, ONLY: initGmsh2 - USE moduleMeshInputVTU, ONLY: initVTU, readVTU !TEMPORARY TO TEST VTU OUTPUT + USE moduleMeshInputVTU, ONLY: initVTU USE moduleMeshInput0D, ONLY: init0D + USE moduleMesh3DCart + USE moduleMesh2DCyl + USE moduleMesh2DCart + USE moduleMesh1DRad + USE moduleMesh1DCart USE moduleErrors USE moduleOutput USE moduleRefParam @@ -959,7 +964,34 @@ 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 + + !Get the format of mesh CALL config%get(object // '.meshType', meshFormat, found) SELECT CASE(meshFormat) CASE ("gmsh2") @@ -969,14 +1001,22 @@ 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 CALL config%get(object // '.meshFile', meshFile, found) pathMeshParticle = path // meshFile CALL mesh%readMesh(pathMeshParticle) - meshFileVTU = '/home/jorge/fpakc/runs/cylFlow/mesh.vtu' !TEMPORARY TO TEST VTU OUTPUT - CALL readVTU(mesh, meshFileVTU) !TEMPORARY TO TEST VTU OUTPUT DEALLOCATE(meshFile) IF (doubleMesh) THEN !Reads the mesh file for collisions diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 index 6766c49..0df1289 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 @@ -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) + 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(p(1:2)) - - READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2) - !Associate boundary condition procedure. - bt = getBoundaryId(boundaryType) - 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) + 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(p(1:1)) - - READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1) - !Associate boundary condition - bt = getBoundaryId(boundaryType) - 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 diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 index b283d67..3a72d12 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 @@ -124,7 +124,6 @@ MODULE moduleMeshInputVTU DO WHILE (iStart < nData) iStart = iStart + 1 iEnd = iStart - 1 + block - PRINT *, iStart, iEnd IF (iEnd > nData) THEN iEnd = nData @@ -173,6 +172,9 @@ MODULE moduleMeshInputVTU INTEGER:: numNodes, numElements INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types REAL(8), ALLOCATABLE, DIMENSION(:):: coordinates + INTEGER:: n, e, c + INTEGER, ALLOCATABLE:: p(:) + INTEGER:: bt fileID = 10 @@ -200,6 +202,8 @@ MODULE moduleMeshInputVTU 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 @@ -217,7 +221,164 @@ MODULE moduleMeshInputVTU 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) + + 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 + CALL self%nodes(n)%obj%init(n, r) + + END DO + + !Count the number of edges + SELECT TYPE(self) + TYPE IS(meshParticles) + SELECT CASE(self%dimen) + CASE(3) + !Edges are triangles, type 5 in VTK + self%numEdges = COUNT(types==5) + + CASE(2) + !Edges are segments, type 3 in VTK + self%numEdges = COUNT(types==3) + + CASE(1) + !Edges are nodes, type 1 in VTK + self%numEdges = COUNT(types==1) + + END SELECT + + self%numCells = numElements - self%numEdges + !Allocate array of edges + ALLOCATE(self%edges(1:self%numEdges)) + + TYPE IS(meshCollisions) + self%numCells = numElements + + END SELECT + + !Allocates array of cells + ALLOCATE(self%cells(1:self%numCells)) + + !Read edges + e = 0 + c = 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(meshEdge3DCartTria:: self%edges(e)%obj) + 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 + 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)) + + SELECT CASE(self%geometry) + CASE("Cyl") + ALLOCATE(meshEdge2DCyl:: self%edges(e)%obj) + + CASE("Cart") + ALLOCATE(meshEdge2DCart:: self%edges(e)%obj) + + END SELECT + + !Allocate edge + CALL self%edges(e)%obj%init(n, p, bt, entitiesID(n)) + DEALLOCATE(p) + + END IF + + CASE(1) + IF (types(n) == 3) THEN + e = e + 1 + ALLOCATE(p(1:1)) + p(1) = connectivity(offsets(n)) + + !Associate boundary condition procedure. + bt = getBoundaryId(entitiesID(n)) + SELECT CASE(self%geometry) + CASE("Rad") + ALLOCATE(meshEdge1DRad:: self%edges(e)%obj) + + CASE("Cart") + ALLOCATE(meshEdge1DCart:: self%edges(e)%obj) + + END SELECT + + CALL self%edges(e)%obj%init(n, p, bt, entitiesID(n)) + DEALLOCATE(p) + + END IF + + END SELECT + + END DO + + END SELECT + END SUBROUTINE readVTU SUBROUTINE readInitialVTU(filename, density, velocity, temperature) diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index b1a91fd..e557495 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -514,7 +514,6 @@ MODULE moduleSolver USE moduleSpecies USE moduleCompTime USE moduleProbe - USE moduleMeshOutputVTU !TEMPORARY TO TEST VTU OUTPUT IMPLICIT NONE INTEGER, INTENT(in):: t @@ -528,11 +527,8 @@ MODULE moduleSolver CALL outputProbes(t) CALL mesh%printOutput(t) - CALL printOutputVTU(mesh, t) !TEMPORARY TO TEST VTU OUTPUT IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t) - IF (ASSOCIATED(meshForMCC)) CALL printCollVTU(meshForMCC,t) !TEMPORARY TO TEST VTU OUTPUT CALL mesh%printEM(t) - CALL printEMVTU(mesh, t) !TEMPORARY TO TEST VTU OUTPUT WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", t, "/", tFinal WRITE(*, "(5X,A21,I10)") "Particles: ", nPartOld IF (t == 0) THEN @@ -565,7 +561,6 @@ MODULE moduleSolver !Output average values IF (useAverage .AND. t == tFinal) THEN CALL mesh%printAverage() - CALL printAverageVTU(mesh) !TEMPORARY TO TEST VTU OUTPUT END IF From e3eeb97f4863117a8d138c0ca9b2843e5f7cfdda Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sun, 5 Feb 2023 21:33:03 +0100 Subject: [PATCH 09/15] First working version! First complete implementation of .vtu format. Still a lot of things to improve but right now fpakc can read a vtu mesh and write the output in vtu. Still to test: Multiple geometries. Double mesh. --- .../mesh/inout/vtu/moduleMeshInputVTU.f90 | 125 +++++++++++++++++- 1 file changed, 124 insertions(+), 1 deletion(-) diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 index 3a72d12..f9588a7 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 @@ -302,7 +302,6 @@ MODULE moduleMeshInputVTU !Read edges e = 0 - c = 0 SELECT TYPE(self) TYPE IS(meshParticles) DO n = 1, numElements @@ -378,7 +377,131 @@ MODULE moduleMeshInputVTU 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(meshCell3DCartTetra:: self%cells(c)%obj) + + CALL self%cells(c)%obj%init(c, p, self%nodes) + DEALLOCATE(p) + + END SELECT + + CASE(2) + SELECT CASE(self%geometry) + CASE("Cyl") + 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(meshCell2DCylTria:: self%cells(c)%obj) + + 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(meshCell2DCylQuad:: self%cells(c)%obj) + + CALL self%cells(c)%obj%init(c, p, self%nodes) + DEALLOCATE(p) + + END SELECT + + CASE("Cart") + 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(meshCell2DCartTria:: self%cells(c)%obj) + + 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(meshCell2DCartQuad:: self%cells(c)%obj) + + CALL self%cells(c)%obj%init(c, p, self%nodes) + DEALLOCATE(p) + + END SELECT + + END SELECT + + CASE(1) + SELECT CASE(self%geometry) + CASE("Rad") + 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(meshCell1DRadSegm:: self%cells(c)%obj) + + CALL self%cells(c)%obj%init(c, p, self%nodes) + DEALLOCATE(p) + + END SELECT + + CASE("Cart") + 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(meshCell1DCartSegm:: self%cells(c)%obj) + + CALL self%cells(c)%obj%init(c, p, self%nodes) + DEALLOCATE(p) + + END SELECT + + END SELECT + + END SELECT + + END DO + !Call mesh connectivity + CALL self%connectMesh + END SUBROUTINE readVTU SUBROUTINE readInitialVTU(filename, density, velocity, temperature) From 402dac9068ffa6642734fbe7386a5da687499f02 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sun, 5 Feb 2023 21:44:09 +0100 Subject: [PATCH 10/15] Update to the manual for .vtu Manual has been updated to account for new vtu option in the mesh format. --- doc/user-manual/fpakc_UserManual.pdf | Bin 181644 -> 181522 bytes doc/user-manual/fpakc_UserManual.tex | 16 ++++++++++------ 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/doc/user-manual/fpakc_UserManual.pdf b/doc/user-manual/fpakc_UserManual.pdf index 49536bf078a47bc18c9474931ed1ab9d0ad34bfb..154eb5b6131721f0d957daba28b534627dd05e2b 100644 GIT binary patch delta 27246 zcmZUaQ*b7LvaR#Swr$&XGO=yj_8&VFW8zF~+fF97ZQHtgSDo{4>%4V8^+Q*!>Rzk9 z*ciB^Xt+8cq(-#;A1%*;W5PTBKgDWCBhnfhg_1(Sxutt35SbIox8?i``(Nb?jk zJOlZY5bL(b57>5w?YX&wnX8+Nxv@R$e=0{48(0o5VrJt17A)-Ctc{D12(7@4B>Yk5 zr3W1%mUuQq`Nj)_VrZcZbZZfkb@Jv0vaWluNQj^UNI2v}mC#SHcRTatV!BG+PD_jS ze5(l3p?*0}{RfwOGmR^*ENA2_l%MpUkIN68bhG|mff&}BP5YzVJmOIBT3m$`!_)sV^Q%9T{+}gLDYq= z*G=J80RASo4_!;|RzH(TiqA*wLg|c$FCX^B%} zCba482GwiYJmc$C?@Qdk1#sv_vJjoK7F}+%;N<+rKJ~M8v-b{GC_fG73Ux=}clU$@cK(;FU%5Cte((KqjU>g-<*#hh%Oi5ri&7ibPtec)+acR_eyV~>q zf^SsC<mFSe)uRCC^iPApISeR0EdfvtRzLKMvCgV4}P zfi`ndp#LW>$+(481<5Amkwed43pXLr-k-tU))=)}7Y6wC7oKKS?n)x@Ptyz{q_$3T z;aY(7@@@_1bMQN?215!$$fAojla@0u>?(6t-E{UxDJ-u`5`Hl25d>Rf*aj;4g*wtB z10@2e!=vO`AY?m(%GG4C`~A0#hsSCInrkd;4)yrD^(ROl+`3P)w;V1z8;o9K%|69z zUYA2>v>aGHQc0w<;II^UKp|#y^dwKn(Dj5NHy=W=Bwo2=B5TB(N@lpdqP)a8u6BJ^ zmG9k$`p=_efLn`%Wg$^oD+dP|VQbkDNu8r25>3_pwL#E#r<$%(!l;9zN-49&9k4Bv z2@s>37ljACYisO|UKhS8|mxOX;CwNT5xu zx`CDaTi44hMN#98MI0e;%MaPi?Og=5E)v8j`uo1i`WzV%m;AQK!x^zk?2(!M%{w>f z{6HaC3HcqCHA(vew>{N^AM^$haO+wpU?AU6^$U%6nM3{;cH-(*f5(v3?pJ-VV(h$QIK`>v3a^ruDx|`@|$+$2H z5qZt-old3Lz4mX?_ttlAQiKMW^e2QZ=rwLyPoB+ki_JvSs}QPoO~4k<>UGYiKN7P4 zemegGHh*UVq-NpIa-fuwnDvmUh-OANh%Xt>>3bd@v~%2|yq3FPoGG@b6(}FGC9uv~ zvI4sO7s7z!v9ZNd{}w+P0Ro!)3{ zfEanN|Hf)<|9Xj{67LD@PPiS;-UVy!Vln-@NlXoKAPI!ZFI)gR^9B{$+r|ZfYd;pDDVp~T9b(}B%*kFPSS{)j(M!DgH&xk~HI6dldO!;wA}aYrUaj7#Yg*KQ=qjH8ew7E zaZ4f(qlQ&lXI~Sc7zeWg+_xvIKIl9wYkddKbD< z1X8g%UOdW{b}#fX2CxS|z}E^m!Hi8*Rw3}RM-eS5sVuTd!Dsm`#f`)l>z4B{ zclc;{5W8zj;^YBTY&VFaocyR^}>kOv?{nE;${j zN>3#FlzyqJ-AF8(1;56)SfWyAgddc@e_=vk6t%#)DLDN#%a(VGYlze_bkD(l60@Wx z8#JY^?jHJsj!}CW1w2B~%`^1MI1lYT%yt{>>C2VIPSN391Q=z;<*ZY4qkw82!*8ra zImK)+sg+I;qM`aYyeL(kCs(PPh^YTvoI^CmKAg&5lG!Bub+cY2p?qj|mY^jfw)>PX zHmyjm>O}9F6rR8@ezGGekuE}(KtK~r!54lvJY>Qka=pP00%TLsn>Q)wvVgG#UyV6^ zGz-bXOE>I3oe&M*;O{Q08{PhX;oq}UtXnGJ?Y1RU$GTzs6vd6<_Md0oE6^JC*1hZ~ zRoheTz|+%Pxbcg~Jc(jg>vT*OlE|Js_8HUdogN|lWQJ^0CGW(-yWxavrVFp%=)J1x zHHh^w7BX#V015&_&~m*lc*eKZT08$4v%&log*5x#y?0XrIg3Cf!Z;QWUWkGsgTk71 z@mPkeYghgqobjiwX>L4-8T3!>%DFR$>Vx>C`tvP)yvLl61gdTq(gsg;udNM;)UE<6 zFF$#Q?DEWqRdhpLelZ&|5_kJo0MuTxT8Z-ffq>;=1<*4E4p%;x0S0^ZqaJ4RFB=Gp zO$hf9YDSfa9^wi}nwE(ujMbkcOx1g225iPZ&0wcX&%MqQ{nvcG#M`y1aP#bsUYHc7 zxauV60@B-40(_?fLHUuT(h%rVWe-Dr@JvI+?UsjVT=KBD7A;Zdc#OG;91W~L;Zh5K zj!#y)Ie^FF9KYi;Rz0odq$ILM@=7uOswoi&W}$JXJsVFyZLCIC({tcN369xpAl9NuR+s@nSbYb;EI7ToTl4}mIvc*Vw}t!qM?0sECjPx zO2b6gHzy6~4K$+|cjp*~*F+YJ7fa5LD7A}&Py?-HDUhrlo>f0FkCE}isKsS}7c+v_ zQBGASY*C-3 zc8J08ye{7CMcXJD40W)m_x2uPBV=JV`g1m zTMiVBIS1csuLdt(Sm+3ah!^JVly?bL^9#GIHqIp(tb*Y#9Fy$aogc#NV^v9);%ir^ zX9s^?m);&bl3*2)tigGlwpn<|ujf`QR%;%~VRs3d`6v+;3eH62m$Zl2fxKN6owQyl zN2fy^R8k2a+ZBBr53Dj(XT_Az)LH=PeE<$+lEIiA{^k7*0$?h?CuR1EqY=&wW`zJ{ zV`WV7(S&T1#(opLs2j zCf?XP7?I8?@lnCF@B?BN69e}kevKs6qleanzE6+MNS|Z>eu>jp3}>V(XKC?7Nd#i? z%Lk!Bn7~mW9g^M5!Fra#K9^sg%U!CXx)w27!pyNrWXXqzk8-JM4@(Pvk-l=i;cn+} zLNi>~U*xNJq-7+J7G+tFX_}zowr0qfkxTC`4jG9JF=~NIHdq>K`s9nLw*{*v;1Ff_ zc4i~~8xX1*eI9@pHALP}oY^^?>;%R|t;vbbi%S%9@ZbTyo$cV}Nnwy!`W1Ce{kp+i z!SWZ)Wtivd!cZ5Zqp@Xgdh_Zy>s95nvN(Nmck{xUmsZH>=Mzk3EMBv2wO+62hgaw1 zLPu#mlHNys39nEv1p1{I%*R3JPke>9AJEr+y*3Ef-J$Qt=1Zztxg6*fPXpH)tT<-p z;o!lPcgl$G=STnH0HmSraratu2?^;4_0TR8DGKlrV-#1>-F01HeBjce%=YKgaIB^H zppu>`aeico@{Hcz);-X3@?igMqLztL52{zN2xgZ{TB29e$$+tK^Sm=>>EQOG%PG76 z+6;8Vg;uP)r6FeU%;c-3c>{9u0c;@Yy7df|P&Moi35)-+JH$XUx-08T$AFt%twm1? zF6=B6C(n!{jyAU~Ba%_VKd^oUngwzN>nl5f-&#TT*Tv4UxL@@PE?GL5WD9y#gH`Sc zS}>Ff*eo_2fDVA$bRt5OfY4QZ*{0x{qF-$LP1QYUa|z;2UofYOIsjqWRN7WU?-+iu z8VW2qgf79kTc6m4=G%TG#GRg~PR`5=T9Hb56IIPoJ@0LrOH!VYODDDC>||^mj3={_ zg0yiOw}ZC$)JK=waF2k}2#ujpN0al#Luu?p znL*a=x$US&XrUZXtg?X8(SU2O6P?#6_RnFrHJ_wo%9=fC+wM9!xOaa;@IJ3LPqoL@ zJST@{41p0okSh(SC?m(_%=DsHq5g=_o{(w*fM;1?2vtyKdqCJ%=k9^7Lsb|AJ!!Ww zLZVq8H%?_3g#|H$q8?L($yHDbW|8-VE^EIGAuZ2i9A3-0O2q;CtS;rX9`k?IW|Pu# zcRhCubRFC53tCT3dno#I5;e8@4j7Ya(j#9=^Eue%q$arFzOHs`On^)DLg^+ zP77mrjWES!ZkNAkPflM5P$iJl+7FWsxoOjGwW4yJd4XQ=`}emv;?h$S5Q;-H0z!HE zzd=tP$o4JQD@LP(E~SbdvnTds-Iw<*`K}P7@??Mit4O_VrJ;#|dA#Nxq1C$*Wzjid zW}?bkYRkiV!(A&9(;(Xsw02>&xm2-XBWt@P$n*gZ-B8fmrDi4atAFn=DIrcOE(C}o zZ09QeS^=IoV%906Dv+3rxQI%x@G;cQZwiHWaPFTt3ye_l#p|NOBUQJVN`)bC!VRC> zY8$feFxig$q*wnbEGMsv=bXjUA`El<-HMg*5FPiSy@2bB?#`gdGOH?YcJlTVX&}s$ zFk4)dHAY(hTMhI=T1vybj|wGOM4Wz!vLDs4G6LwhaDQayl%B>^8bj0+ki0m1s5P~~ zfR)dPKQPvof{##?NZNbrExnA0WGP?ce*bb23CA|9k4l@Cy!$pp+F5jamaBetuvhEv zE;pJuUv$+JU=`XmkVPnWC)>laJSvn7(qk9&9K(fvC`SEh&DRnTdx_y4BfvdOmG_@N zNC`+zAxS9UhDgOGVlYnm<#%{mie_)H#s*Bwh~6RKr9aTr*>M z{6{dTyxduQNNV3xC9(cG>vqWOA`4z-Sust618_-ZH&Hqpa##2ejWl;}3_*rTFpk5; zm1<%8+cqdGvREq63-C5&Mn%#3H|<01cRw&4sqUMUU~s@DCJm;%*=B)*dd_7(MKM=p zPZ827J+$sg3{g#7fitepNoOKij-0RpDAvEg%q?WH2WGOotnd#PC4vQX+fhQtPeCw zD>-o-%Rvhxrw9|y6T438>RnR?LARTeoIkZDIgJAJ(ju)!~d z5bzNxZPo*+xl1w8%nwzU{zu>Ss1Ll0s?TN6WEN0h5uj(mzKu{fWiUp^CQ37B(h&~f zTsK|^xq(Cp)3b>#4+T20sH3nh>7Bm@`=?Qv^KB3!6(vqAQRixv(17nMT&$_$$(BUm z&rqv}ed2s-FPUm4H8|B^3l`ELp;AI#jwC0tC*+ZVUN#!c*#XyJiJNhl{e}3m-;kXF- zu?Vq?pZFC8q?YcC85bmkHo!R*jtHG;nty{2oE#oPRJL=5&NeIy(lzUiF(ajV%6o5@Bq3LgPDlmx$c_6^i_mW+A%DFU}(t`A4^sMQo@Jx{iQO;u@&w~mgx2M zB)1Art2q2^o6xPhk~_qQc|zNcQrCLkL=z}GX5Jq2nO>9F#5CXmkTz-S#UPi@3f zO&23k15#QQ%B6b#Z^FKr-|j(b@=NDPFh^CcAu*OY^BwsoJ#*%c49bcQ6s%q*Y!*3? zth~S=@rfKsONndb{(1PUG%J@p??;GmAvKV&OIbBfB=ilU1+YDPw1IDFzM4)rS2BJv z&Aq$)&VmWO?QCW20lc)1et-5pT5lA@5Rgl~I!q$x8)ly8*IO-xT1JYK)w1*yG4}TV zE-AMxTU%%>MDEM~n!5o)3mML+GwA(NVl4C@D;7&%tiJo#M5q-VxRu=te0uw#f&vQw z3#O{_IJf)C=OU`b&l|_{n${=R)?CvSBoYHDiHRqC$_7Mi1FeQ$^U(mZZy8TwfW;(a zWcl&9ZxsT52BkO!bMaNL9ZB-RF`FfS*0Aq?yVQ|glP2vH)QN6#GUVVyY`I^{bhuVzE;xfNZWd1;#hAy+*~QCsbd~?*Peo zb;31f)**@P=J&VLw|K2D^o@on+HG-nt`puxH`NESigE;RCGm+DS%(+@S1Z@YQL6sV zUUenKF^Vu~+p|Fxk|%D=+$b1RGy6gj*+iCqIOF>)=ue5g$bU2%l4k5^{3(aw(*lMI z4V|Hw;8ebXwz)>dFi-}QdTVTL^n-Uac^M&1 zw^ehoe!ztC7Gze04QHwVEb@7kbXBa}GFFbXvCiPnwyp*11Nj$ZfFUl?KsVJEMa^5n zSbIzJ)p$ZeffxAHk1VM&vI@gxN0Nr40TX| zS6j7y#xMiDvh8<%Ge7ZEeN$+{*K#vIVh`&by=W>-P(wS^+EJB%~pxHispV$K2c?V)k}Er_s=0B=H^v1>m?BusKaS~vQ97|MYrN3#4py$U!g`;=SHTm zhzUN=VVsupV)}-`sqyPL!A5;X zfO}+1B&Apmbv=A2xUh0Q7~epkbP{E0b*ixuW^vmA_E_6bC} z7;2rF*O<64O!Lzd>NdO^s6;z@?Plz(3)lMf6ulduf*|7+m$e-%hEIYR5ltUT2Riu= z#|aEn)5=8yns_(DrBDVWMWXdo2xt@l^3$#0-8WBV+=dm<L4xD&&;6fj=Llbjvqp_W{fz3c|Rtf5qp7<48Ri@m*R#I-C zyd+1;ViM}ho?Sgax=8`?!gd7S0MUQ2ba(0VXZt;K2uEPD_RmLvO`aKCF7VdlSiH;S z8Yh4!roL-L8<1Os4pwzHRu9oD#sXsE*EjxX{>U6a1s?8F<}P7_@UW#}{jOAM&%SLa z)BMR00Mn!3koPge9-QT%+3J8C%lL(F=Y(gzDaMf}vIoA|Zf+5a(irKm zFsW8zhL>mV1r0XO0Nu`^0-Sh}5kr0~Fm0@n_$Bp6LY+(S`ZG3BQS*vOy4V>nb+ZHOmv>au zlyox_Fy64@j~0xdP&u@~r8lvR;)xylhUrRO);f!`5baC{Xdc6gE-*0uNuY|C-sPTl z7w61ARd&K?W2;J;tWb|UT;SuoYn!>i4sgdUGH|ouD+iumZ zU#i;=HBTdP%IC&4PtSwa?GbIueEMiezF~kqEN(Di`&@>47nV8sLGvBDAU+6#nf6P* zRIRnEi5^$y>96hUOoI-8ZSb9y0J3>Jl^(6t=A&@-wtCayTwpW&8ISFwhLgviUckt{cW<+&riU_+QO`<1c|M|1h zsgr*}u{U+>bOFPMM2l2TNs{^-_(*8nd>3>Zik(fhU3>f|(?i%RJ1qWHxN3sLzs3A> zyxG9!TuqhiA61{-em&6hf30|Cr~VhzWIq|Aqt}q;NZk09pr(IIeenms6XS3`dC&0v z#@=+&CSth!(4fKO@BZ%CIz@fjKRVZMooBGQV0FR`sIOk z2gcM7#uP1G8;KLwWz=ie^U=7aYAgA{wL|b_lP-1nkThz z9FxP$Xbkw}7-}fsL-t#}(IoW7_3~_5bD44f({M4hxywH?pF(nT;4m4grktZ}HL<#M zqLl^L+KcjYGr44hxDp4O$a(rLO~-oY2XQs`7U_`Zuhb3pmE}xCeWYy;;j;r9EbYo{ z>p$;IGa3p=UPzh8cz$0~t7|S&dWg3dToghqX(olcb!G1?A?O9d*NYvbwdupp{qKM= z9Ov1P&+ZZr^v}~BsJGAax!!)@50v2APs#tu?^(J2H^ArMW?@ZExs(7(r{fQ!JZm5S znvAm%@=W{~K-hR)gZ!*x&xRywA}%kOgDyZZ`T_jd{W`?`le#o9nLm0UTuDi{qSfZV-3%;Hf++hVgOJch_zkLD}ud3VX zrl-I6??0e4(RS>fF8P5AHmrN44+h^K_e8h6Sm8ZPVcB6JjEs5xzQAp}H9H*Fa7Exo zpCLDz&3n=3yGtWSBd3O9deKqoT{-i@;1nj`qmPO+WcG#iwq3P3GZfbs!A{aYAFi*) z8+3yo*gO^n3J)d<`q9<1T0B8py1kjhs~Ij`M4>RM!qF#f-+my^q%fS`HfOaXT^$?c zJ}ZXD5RNINC6sFNGA~=yktsknjNnkOO#9UuE0=d~>d}VGVilUjAMT-1)SAm5tT&@k zMqbg`=XRp)?8>6jqA(lFIbiFqUEAh{kx>Zf^cGRsP$}AeMiF}y_JYW_O{<088sMX+ z_>91)$77WVbOwI%xs4M{?xjy@Ol~_*AN2j63rfujdpz3pTn*V>x7NbGZcx?esN+1~ zi2aa1L%=LGPVYjk>EI& z?zbSDLeJy1EIqkRu;!LVf};lGuW?Kii&S+(1wXbX@Kihwv^AcM2?=;WSA`m(X1#n_ zot>@&EHv_7UccaR@KKP@zzr@X;3Sk1m-$%8T@Xi4?!p``)G^;=tuvhEE5y(z(wFGV zUc;$vVFLj_$ZdtOP0i=o(QGdKloI#-! z=$~aJHP!PfU9K@1PX>TtFd6#uGdp|wC~ zig@zMW5#jK<-NmgW0g7QN2~dd>8kg2wuAeYi~wD|@qsSf87uDkOdlNciK6y z7;X*2AHSXwLh@lOr`m%ld0;!mqoT6L&ISyY9Ui7>lePs3#=;`3d3DxxF^w05VaiI>T|{`^fW!q4#iaVc44IBEQB9|t(;;QP`IM9 zvZa)%HIT#qD_kq>A`YBnff!wAkDYsAunw#b^`;8uofmmE!l}=wU!YVYF>*dYAr0y` zowx-uA&64PNS&CIFq2U67=p)Z2e3zjd{9M-1#>fbwh`T-C6LrG1&m6Am{2bTnqTqS zvgA$6WUQ|!B1@8Z9-4{GffrWj8qI(x5g$7tkBq1I3q!6N7SUaC$SM{3CUPp+{;n6| zHq&g#=_nW%p7(d5tJESU|xLEM4*!VNj04Q@yA$QF9 z8#43Ix@mGm2DD;CBWmc6_8*8nuWkFjHfa4qq}BATp6`HpHEe-84&1e4&6JHO#70Yk$^O^OLAa|TE}SV$=16qJI$A>lOS zI7zNcAPW;Wldan$&CBHHOo2ih?8I{a*;#|YmK{uf*`dsry2w7P-tCo30|QV^J>LdI zndL9(j3i?UJ!?d44T+6|4V0n{lvQ&G!68`8G8-zy5CU1tIqWh!ux-y<)@d;3Ulq~X z?TQ%lH)On_@vgzt_Aje(w6`XQ@5Qmyfc8DiMeYAc2jYy`uV**KLm&(Wgk3bC_U%bc zhsnGasiSv$7lt^N{we`p#-vss+;)it(8F`RQ9L{M-`DD_y=NSXZfKhjUY1ya}D z%Msk)61X5AW zQg-tvj-B0pXeJD2e)l=O{XLcmy2nE%SbLTA>PPR$YIl#in@#Qrtr$y&7+z(a+Wlf; zu}r>!FVXl;!HS)zFFLk_=9e0Y!))2WKUtF(bH-!@;Ms570rL(OcuhNSm~POVQ2ce!U`tBtD#OiYhE%+7$FG6kAgo0XqL}zvB}s< z)WlCw^oxGK|H6l97OqV)3Azg{4SmA95R>p<{9Xa*@ahGA3@*;s8%NbrJ_IO#=-hf> zfzh+B{hxquBG=V$`H%?hn6Azn4AG#0=$PlwgaN;?)a8OC{CJU-Ec$JRPib)wm-1 z3Ie6~IMQ#)DUrnwkhdw-_sHS4V9ZI;$8C~F`-jTK=ZeaFw`ou_rFA_mIS`K%l`{bW?^nJyBCm{8wsg_ah;Vf_>y&y!Q+>P z7>%DJl~+XU=TGDV<^AQ&e+!H&Y%r2e9P8W!NlC?g4m`|9?NevLs$kw>;uR|D6C_#o zU#8J~1g3C|5m#;r;E0kB$j6=NYeU>)YQIU3TfziHQk<6^`yDI2&64Z6$V6cgZ<-mT z)ZHc(*nAf~YiFVE)RVO@7I(9?jK9L8HHw{}Y=t%yV3M79Dwu)(LQx4VU!!kvvI3D3 z(`f<}X}gbtgISa*qiDW`F*)@hXX~OI{5k{&0;li(JuXQ8EmjQ}mb$F8Cf4{1`FFIg zEDS&=D_%z}@?rK1Sd;upW0hc8B691Sxj|<&_m%mM8nHW?1SOtPA&+ATkNct42;%+8 zg0unpCk&fga0UCp1uTL#VvQ#SOBA?#2ZNGBgI+nhb~Z|u!A!q3w`)whhc;k#RE#4( z32Zfnqj5cb8_1u9aC~!*jm)%!&0S`=@6v8@XMmQsFl%C7wbx00thg^dD5inRC|6@{ z&T2{B?e3p#N_fPM?D9=SEq`buBegHKjMFh=(Znv3U^D8VCB(`4`_Ijsj}x5Iqo1WD z#VnSfUnP~pW;hBofg%99lMWx7y5wwI2{^Gen=(M>cGp{ify#L+<$1q;luuiu=x%uK0wyBO3ZPNL;__F!fwhK zx=is;UW;SYhLZ-_l9~<^m|D*kK3+m2-ZrO{rAeOg;j{}refLwpsK;g7r?~0?rOdnvBIqO?@LMFTfIh1l z8$^#=&J#4r*$W~khXal?&3>ns;i*z-KZ?t!{P zUd4z`_m@kMdOo2@h5C>aVBy*C!e|zUVU#LFm|`|r4tm`Ek_78)Hc7SJ#j5s;GJ{hP zclNXx9_ut@Y|C5y!1!3GjqK8Omw$vN7<2HwWxPnXZ{Z%uj+jOe-V;n%=5Qa^8ma==$v&vGK2cK|^N|vLP{`a7prmgHA1S@L_X2l8-se{pd{otF5oxl6tFju!Rl@ zm&cSN7YKeTVDx38U>;)4EYVTflgG~YZ_w%eusw#KLUWb(#x{sirgz41sX~ocy|FQP z1c!*#(vk%BOay-@5N7f}#{eRF86i4!NxyO~3G#n#3WoXGgs4d2&+)gmB1Mj(<@PYJ zo0G9OeXXdW`?!YiA(dpa;~H`u4y`YfLJ&r0_>pWum$bln;VA{7iS=*V2Yr;rqr`DQ zX>dRZ?pFLM&RH-mNGcEB!KEXwWMQ%yVB+Pu#U$}O9n)SJsF)X@k?U+eib{J(r%Chf zQ!DpLefR;-cPmggI9M8;o;+|kneZxrvDhr%^C?xMgkT*1idNPBSp=UC2v_2LdN zw^905B&KZXWqIUhJ+>t>od4{1iT*o>_IeE`ThShv>bS=&z4+2Wb8^Vin0fVD7=F!V8xTITG>^7|4CKf;LJ zFdoAUTFd_zTv{B3I!k9wG{}N)E+%-8Nl*SSm!f?;AU10i3IZtGx92{N-Vsi~X6U7p zTKl|5FHN4AS1X|KqCMZxM5-tRQ;a_<&k$O>@xUxtIS8&z;7`5?e74rOA3;p4&DM#7 z_YyxEJom8Ql(jgskyMq0m?d(7-Mx9@hZ2=p<+Tq)81?dc9Dc~YYhVZ)Y79f!_u|V7 z3*GP=5SK=O#IuwcKS<80=eRj|i`7oUtH4}iOjoz;G-54 zkf#g$ixRbJRx=$GCdNKtXHB!^fYLw1Q~JT@B+yMb3EQk)Uo&V?6HZqaNBP)0Brt$z z0>$XxGCGK0ivq!IMQ+{OqujwMC=^{8GrTQH)E<3NhgS2dwP3kjQt`y^lo)~-&n&oM z+)fLN(aOv;FLb!4D1s0600+34Z++uSU{++K*}Ml&kcE9=@CURZgkdj8MsL1osDl-( zp{A}Z>!ocP{ifrrJiu)xktu&VH}HzgGdM(a5M8{58V*3XqSPe)`xb?8<=pCG0dO;wj4CJ|2^~ER3!*W@f$gQS>SPu z0+~uW49b@Y2}WXd3xDc~wLmNYSe_oWYN;BeJO|Fxpux}FHTxVeBqoA?^#-&5A6pqv z8|!~pg|V^1Q|0fNRqF2B$C0TGQ}bbL~v2x+g$&=-yi*t9ql9ZI>? zY%l${4~d2x33Q+Eu|?K-ue97Ik!O$}aTyQ^NlerLudv^37SLf$%+wU7W*SU5hb$Q0 zZ25&qW!S(W25C5Lfo$BBSE&V{+F+|S*>YvyTZC^p$}NP6AhK!SfXbI*-@~zah`uV^ zU_Leo&u9qsp8%~E{FncvGl62y<6=y{BS&{6kp05ta^MlUu|PZjjx2@WOQ#X-su0aX zZg5^XH0c}1Nji6Xb`%2quf2vabHnj}?KN)B|KVfUxw(0o(E;21Kn?k1d=B{TncCCB zER(LZqZu+KrD04GK&GAw{(xN)u7{;1SSao8%me$E@l#2%k}GyA!^FI=V#AKdw93gL z*B6mO=Y-|2t*h&+-z7*~3?7Dvos~qc*M4>FX6NB6e3_SxH!9A{r4HOz1{q8~!O~-P zh}h1LVN*9>+iNoAKn>UH9E6l{OV^zb8!{+P;ES-b00c!Ev&F%pjR@H^8WhXPYg zWea6qa5mw1+O8j$Un`!VuW`@$i`lq@WA2lkJ|aa%o-ep8%c3)|!bj*ga_)}MN!`$v<`}_3zvVY6tLeo!fU!BCD>MV2we?GIDE2U> zfMqBN!B8@_?&pC;3kk};oBqF4k7eI;sCH5NF6Z8NGVz$G-`W6P=e%iGqZq72!XIV7 zmHBqn1*!sIK^pQ@ZDck-X*{=6y?P+646Zo>PFU;yG0J2kSP+XiOI}~%n!MBTH1St(V|q0wT|-?ktP8(7`i;p5!ZiWkx)MI{6q1T>Gis z&EhLKR`H0`4P1^GgvM<6D7)Ag4r%msVVPw=&Eh(+V07@E`Z#!Yf50UuZvH1u91%C_ z1R2!{rRk5o=-QA(%3rl<4O(WVCIoUJ-(Two>T`UC1+XG^bkz`uogUaM*S?D^3n zpd$FJR+m3?MJTiD(-Tj6^s&Oz^3E^S zq{k&_ji^M`!_rEs<5~MVn5Durhbbu2R3_IzVT4M^o>HHfG!O~G1t474u7Qqt58nPzE!ZL(1lYOHcw}1`3IR}sa zKElQ)j2*C!>@Mw1L>mck&~DdwWJEIr^%vmJx>2AcUT^t4KD7d+c+%3~nL5dlk8W zc~ENxZB*24O7OP3)mf8_@`cSxrG*b!mxuv&k9QR{i@k}Uv6dff&9FfCCy ze3?9$%w~I;_5}9sd~L;(e||wHk(&ee=IZ-NdP>v7PNa;J@@~7po#vJ9XKiS}rHK1} zx~0jke~ZPLqWqBiSi+N#NDH=->;2Ks|0gA@38f&Ni;KxXE{ANGn1ECG?r7KR?Luri=F|@D^Ng+ zO+ERij_WP6i%L<87%QzZ1{Q(EW%^23qh@weepJzCWge|Jj0S$kUaQ6sHMSs5Y0pTV zQffYhn8`xNO`bxWhDH~uK14T!zE5Dnz>)(^HW0;a0#+ZYB_oceh*lU$6;(gPZNh(t zeis`Fx~Hgyypv2fjuue)i6{dMCWrh%`Q|;-I* zY35xq><6pU4_BH$?8Uze!PHy@{o9E;{Qd3z4;&Wm@^D?-=}$l>p_&oUrtlXWZ?`P8 zbAmlsnY}P6`CmCm7FFO7J0f&u?C4(UHxI%1Y-In36@CnxGiWj_Oq`$NGs*%FH?+P} z&zICi!C)N0ECLi9Kg|^$39JeU!C%+?R6I_bBjRt9$pcp4`)j%_>LcVlONfAHPg7>U z+EiMaBSVN?rY9B`FMcPG-+rw+>-d(f*}Muy+U#}hN7=pdjlC%dswU94Ey<5|Qfp^C zju02gz)g>TLz?}JSLZ4T`b!LGOz79*Kh8z*^SI5*Hg)|(E)eK2!<`cvCbP^Qo_qvx9;VQPLVaoYo22&7tITstiwz-97T zp_98yD!qnv)JunVrv80a)1qhAt@hg)hS~_O3v$hZMPNr^tAbYMAWc~-JdbATPkcI& z*w!hrj`!g>VP~u zR`iY@-TSEHV?O{m-C}N7f?-9s$nNqt`w?-kJ$L6QR1ke2V`ncW|CpRUz%_U8R2tJ? zes5gR2r8d6bZ~(>GQ6o}@;Dnzf=C$cGVGzo*ly$*pBF0_-ghO;A;p*_$U3zuUSojgcjQ_}2VgX<0i<`~WWLEtOEw3%#{{~IWL!2MIGbVwoO z6{hah%iE2#SG491kP~9}>DCYoZmYS#- zZMp<_Y5M6nb6);&!SF9&lM?tpkWs44x4{3!dfa=l29bHaYRX9gZ6xB z>}XQILIITb>`0?2vQahFVa3S9aZNh6H+F!@^CcSS zJ-abES4CI`^-h%j-usy1VZS!v{(=t?1Ks6?^K=9$fWFYt@mI2PYQN!3XJMtpZmq-= zP{ zHc}%AWQ$mkzM-ZG6czgt`x3`2(xV)Zn%`|>S2#JM(RbrKTX-ipeF>tjojGga9=03O z3Q2PEv2Giw`uMsAb9$)`Oc3H7GHJk;4s&GYP_oeC7ExM)k3=XfpX*09X!X3s(<~4<<=z=-wmH5n(3%4RJjkOSn@=vx+ zy&aF%W!FB^)2e!-GS&#R1N~o%%PY4!h&dlaYAWBw>N41`toT7`tEfgy&SJHZsLW)s zDpnFE28mQ^$X|psubU=SEZTt%_lw;t74~WvhbcvqXD`n5{^o-6%{gN;Op?bga5CPE04=L`~Kz; z@hC}3WJ%UkeD)=UqzFYDrjTfpQWTSl4@0GfvP6p#5fj4LnhM`RF@x-D*Nv)ytytvuEsuSU7)&lx4N$>~cpB3BN%t$xFBX~>jc zGsAl%lw?17cWY=-pPbUCi;mZ6*gn+{zirh=fArLaCl1(V(7|VeEdyu2TNga7*Z$Q- zkMvsodqhqR%s6SHS7EJY8WmMUy1nS0W&hEBR=Th>;zUf5S#-@|wch%3 z^m&(xBYF3#ojuf2_uMxM+T&g|jvJbjIV7?8iFcoOhA(syUZ18f`ReJq^^3Rl*lSXC zCu02bff28z9hVHc_EQ%7)_CbgW$pXmo|{w{r;~8*8kJE0rD%&!x7CU(X)Dw-U5o3S z?oYg`?d)@Kgof+oEe|5xk1g?z9Ut6o{_9a6miN8-C~dk$h{OHNyWYJV{^~jJz-Xtl zb$QGBp1ir%_CemwNwtrxs$a!*>3^wP>IM5LW2()bbXn?AYkx;_zwwsVezwq{B5zb# z-pTd{uRJ^1Glt^F4c5x*T`~Wa`!9W#7gQa2X8)qt(syLc|E}MuNKdc5kePhw-`eXB zE>}DVj0iEPxez_=<}GP*fOELVqvP2najstmIc!LF+z^;KZhqRBPYZVKKVt4JWN)j9 zuQK>2{!K=8hp;n$^f%moknA|So7I<21@+@obFFn1fy#+ZwWo=-y}cbt47Uwh5ISZ| zP?)0DIGRomrA#=SKgvvzaaBI)t@!s=%{Tja})lX3j zPBJo7SR{EFDm*2UZdGZn+ZZ^>d0MY?ua7OBwYe;GUS>+|ZBujmE*>Mthn}0a#+jK^ z^jWdy>gXM2$D?YU+?3^u+$U83@3mQOyH4db*YyrP_Uho^oo}B``$hhIxUkUc#;3#H z^|R+(X76b2s|(DiE&Dg^zt3+huDrHTKVKL^^g3TDm3PuUy; zY3xeh87x*iOWmFX7#!_pP4r3`yY#dbwby>A?~2iCZfmFRyfUj>@p8u*Zq8YQN)JqZ z?`qv6d9VJFpS>SwIb`POSfnrDyRSd@s%Ut8TkYz@#zO`_xHk4-=D^cqjqcsBu>86W_k* z+jn2@k;ijSjP>1EvU|R5!hZ=@Z3hQhrgmHuQ;|7&{g=SRn{;S>N{5D(%h#(#*GU-< z+Kti}o#9&PcXL%)x-MTrzin7K?3Puz_FmV$8J736+da8qWR)pZYuALI=|eU6FCjF? zzSa+U<8@=3zm=_*9aAs(Sv#VZ1|e{3Fe@8GR)ZzuAf!)j(e>rhbiu*+NfFLPoK zcCzR)Al&Utap0C(AKT}P?m8?N?B9`e3x|HY@c34zAO+xDLNTvP7M@W#azy6<0h zi7}c{#l>&2+&JUpmY%)({`9ff3aVO>-5^{Z$DNb`aBPc z+qd$ry}RZ<+Y>U%7&TM z{Ssm{$Jd{a*4cSgEB(@T&gaD+yH{%z2y?#@wKds<_om0fel$%#ORIDI>kK^js zK2~%h9E|E0Pj>B_5&twHY(RMW`4Iy)jJ~R6dYbYqd8R0Et@E+1%IwXlA04*p_WD(e z-=zGOZGX3YJN|N6&F>Uz6xjkTH)H|nKCugc%zwrW_g!iv@? z{#07G2C)Nu0e;~dY6k1<;#zy?s}$?y6KSb-Fq?r+wnQb^|7S;wcbR#t}bqX zxBgXs&GtD<{#rjlQ)3sElkL~pQr|{Ph+L$XAdalZb;8J7(_j z8+gG=XU#?>zglnv<@UaBlE&D%XJM9cAkgOLlBF8TQlEc`1b)G69a%QEd<>XE=b zxi`EIm1QIrXVlI%znS>Vr}$>tDKCxw#N7$&O%FMRuMgkQRqe*|+ygPklZHRKa`B^9 zdoQ~-XAgOu++mwUR2DiXhnd@1tdgwlv^$-z*|?)?%!GTtO`p(yPW;aNP?wzwA^K-` z^Y=4CMg(!JGBln{8fI#CtNfacTcGi>TNi44`p?~|x%TScZaw;S)~g%j%N!XqB0VI} zu}kVjrLtG&aUB(C%x>))KF3^>x zZ{vbBO@rFD`TKHrJ>9n3)p|Tu`&;dwx5G!bKd`E!X*&a*2lFJ-j*LCEfBV6alT6)C z+w}K$-r0fdx?<(LpUhb&HJey%{X_4j1Uoso#9f%V#ILd;*k#p|v43v2Bs3d#x}^Tu z+4@PS-@TX}r-Ot#ZRhNarFG)o!xKuUY1wVEklNfS^I!jk9}@KIfq*zxTc=@2!!@I{ zs3B)9Z+&!I$`xKXQugYY|n1pu1>k#K=^sqd+eob(q?%0}Dw|WFrwpl0zzwfqeTF}Xytbv=qzS5lg z2Y<%zOxuS`1-r}l8ua2LrbuPKuoG-jJqqgGH@iD+%rKwfR;Tf3LutvH)e}Dzl}xP2 zuD{~D?CayTbFWjQ%}j&By7&9+pL1;W!iV$z@JpS3Ii$2W!h$-KoA5p~!aU77c*drR zH?ei3mqDKz?XFKZ$8;+_KO|~WrSopJz9ySXCnsjN`J%pm*+Bc;nu&^FS0(F0k>(6- zPLfFhI}~K=?|XJA3I{MGK{5PX(ag(h4h6{Y5zGkqn6Vr^$4m6~Bq?tyHz8B8NukbC zBtyeF##su3<}bMuvlK&rk_3V_l(!Lv3k;R=za&qS@_))Df0}88ATk1p(pZvZRB0q9 zk}RQ0gS{e;(}0V-PHM$OBv=wTkPr!if&r(BFd{*6sDL6t(;~W#5*U0PCsM7}0e=FR z2t1+%LL@=Gtrj+}BPkl5W8}|w&4?re-%uJ)VvI?S!Y)UOl7QWgqHx9(jWebgfkhcp zoQMg-(D->Y|4my`tDDdShi#@w5w*KW(-gjr;jr7$9Bw*VR5gjLe~}^a+i^7JNFu|C z*o`?7b4-!pD9lU@ugbWQuOLd;?O1}wZpV^5b~~2FZ^r^>)FdoN;yiff+eM@nM`8s5 zW5N*(=G8pQV+V39jX5yK1SL_cgV zL5dWja4174LKV0YK~fBh-i`n!n6(6!VUPo19f<~FnZ}`31z{LOg9HJ%;3Y(@6fL3w z5h4(e`>_D=2HQ)c0+K98p=JUBX|`EF>09~$BD_?q9tbp!YF^l~j%V?80v>L_?XbO& zW{|xSA&QiUFJy3(ijW~u0Yy$k4utar`E5=9w156;h%i9o7@XGl`iATkp&b5pp7bqygD@--GDvVuG)Dog@#RqUK$79{y&Q=I zUD(T`!I=dA$Kn@cQ21$bL!)>nD3-_f@+8EHM*TP22{sETj)4Xi zMTMbB3?&2vc1&oR!D0YG117$oCp8sMU}6Hrw!jJ~bO7vFgaaJ2E%t2E7PtrC?mLqv z9-u4^)dDL4LNRsG0(u}Y7TFf5n-63$0U5~3C}V;lP-v3`Lp9blP^h)?euhR%30TO8 zXiO&wY1pk|b{#G+b<`c1%41O)_!Ko+2Y!k2;%s$yAI zd^0IvT4JFxLrg}nG{lH+Q2ho3iWFM%ZZmut60nLagw2?}36??2NdP{m3){<~Pz6l` zI*@e|!SX0$NRM!)>O!E=Ke#b3AU*|57?{+m-k4woG#Uc~K`vQU;0sK!Iss~a7LNz4 zh}K6W6tPHZhLFe+XvY9FvXzO)enSG;lF^#`AxNaS`l0N^oRv%JcTk}taag?#yl4*b zPoOcV!hy!(v616YT!3tDn($GQ*vMGn=KVS?+AYv_&6ENjca1(SL z>}}mrHCj?mI9PK?3qN8kHkZQ9LPIu%P>2WvG+byNicD;sgf|2Nl#>`D`NqFol^wL1 zv8)*s7zE8g(8HKu%_ay#GtM0C%XFz~2}&qr+a| zJ0hAJF`39A2ZFlf3AzdssHc7;jEF)J&V$xYb+23?psk7t5J6$sLTE@;X>gMy&wNOj zBtB09nb6z=`9yVe*Jv52Mle>;L!t01XmY^(9vFp>U;~Z8Dnt^-HHgw+GvAD~RzgE4#*_#F^zRw zBxF8pFAX46JrDRw6dEPuLZMMLD*jGmGbX6;DKzTHG%Pv51Ode~ti#J7xCx6!30TNr z(TNAGVrqp_2J>T@MYRe%R7c;k{(%XL4!nQ{#ki`vWE#o@wADCQfMdQ#2Z*u)s#VyF z0TP4}lAfGUbTx+PgqTY{qqU#N=R+mcgc&e6YDxI350nNhp&r1l0P`=D#=vw*l?JRx z)VPpVsm7$NiWz7{5hn1IL!#xP+^V7Nk$fPUM1jVj&p@yj(BC3p6J`XF!rhbsWVS4n zbq1Sh8u1L737VFC}13gCEG@IfBD`Ad<(E z6HMKj2J3MDmY5FH5{wB;ptgg21u+=C2}@zk61>%6?K99+imWEgz$^vB2Ru0-ss_dY zIgByh%|RT(e1@j688O5IQbJQIU<-z!2onGa){D?TsAn0-@2CbL=Ad;235x#Ke6?yF zSZL2twZe1`gCYxNM;%k16Co)LsxI&^RAnwAKkrw~&ikZrhBb6bco{K`m7J2(9p&o}GP&sg*&175U{{e$z)J*^Y delta 27041 zcmV)0K+eCCiVKX33$T9%HwJ$)=6fMZ&*Q!f#XY%zQwmSiTKV=}jr}_DbPH!h+QH0&>BsMn5G|m6r z?sRbDB-HRd{DgS~X!8~H2`x!4hG*NI9D^4^3{SUa6k@Vw*keah7;os`YXhnOc6+oL z6?U7~_d|GDsnJTGmfD73l^MR-PGgYJgdd)i?HXQVf7)vi=+mb-b8C&T!w*|I8WU`P z94|gp?N_kyTH@lF;g8sZ5)x(%(}XP&)>-;lZq_P4?4Q-&FgvATqqQAwYZFiE=DvX4 z@vuz_8-q<6y_^5;y^G&VKU!~xoBBQJ;Meu(?cxO-*bPd89hH>BzP`cjwQ$gT5#?hW z>u{yh#tIl@-dH!e`R`7v6B`S{S-Pr!)6I)-t)@SpY=s|%H+Fcl^&qweW;kwnakbMp z>CH|UP@>=$;+m%dOEA2IpRkQ|qC`S}fPGO4gpEF&Ke>&HfoSeuRx!QJQl#~uN*80+Mj?~siGz)8fz`0kSyV7Pd0q6t7?i@)7W7s}Gf&*< zY^Kbln}N|uJgq|~s#JR;QzISMzkgH5h+3!o81$^yAiD#ZcwP_f=GiD27Aqt5s3u82 z!_|kS0T~(UXQ^8Ylc?QF=oe*wNWo@mFJ_rZ`|MzlbE|eT!&hVtZfj%DsjpzwLfJ4% zEt6hwW}FrG7uYUnk!lS3AFJ?_&xm1om0keDs#xtn@ty?H4m08KPb&xNF7#*ZS>lG) z4BTH7Hz!AJ8{r0dGg`vj{qP^yMUVz;^sBnS$vJpcb{i$0^j3}5p|sV1c=dg2+N~OK zykrQhyuQk#g|EJbNqVrPnrQWFYza)Ka^&G%9aU|UOy!Qw3YrHjp=@U=jDr;t3mqxK z>dc5}=%zm1SJkK@PE6JQtNQL;x2~)ff$7$;LTTF81^6my1`KVM7Ge;H$6f4I)E4qC z?vT;rWF-Ybr6=qL`iEtI2Mv*=r-Uh9Rv9EuXBCG1*=INHL>8g+MqUh~(b=G}h6aAD zinPCJ8k;4aJpBzeYXH)Df>l`KRSoH$YQ&4JM1=?RKz85-eg(M7)9r2@Lzj_Vw8_XozGkXJo4?#WtzDF(JhmAdGgeo9hrF#elMCi zDCK*ZIXDxD-{{O$nK|p85xc?{J$TLTS=&IvkJbxDGzgY=0bP#(`kf5mI-O|^hy*p1 zMm!XBh1mvz?t6h z=N*qD5T0OvDDu@25RrT#1$?O~kdD7nY1c*QtC=)j@wTAk#Yi!1yP#h)#JOXDKGHG?bk-bL97VyiDcWe2_Qp>mXDjHST_Vt6g; zet><|fG?f2Sp&;c2un~VRC$_^SQ^w$fH0zyj$@;LxaOu2JT6vf)_jI(J%(yvoEH4w z$(#vDn~{mt&0v-<yTj^aA z-Bpi&(AXS8T`DCaS;!9iRr#p?Iol~ADOjIL6Ps**9BA7Wk*;j{okPnXsgNd0Jf;#( zNWU5oE+D*f8?Xq$7#Stif!7uxbkwv)$r00pg!8z?NY^9r1Hiy z*i*2fIzqnk^aX4Q$bD2J{JYa=Qb4+4zd{s$DhW0Y&=SZGES`toECB*ireetPhpX#< zeq3y})V{Y%$h}YVZpCZ_>^t~*1)zKQHo7)rES3D8-NxZE4Xd!99!ApTjo zGYJ$-!|_!8evN_(QriVOCpb|u&L+A>SwyX`x=xJxzB81O;W3ofUgJz=6|d-tdJiVx} zhu;48JmF$bu+_0Ty=xv?&)!}5l{6M{`lMmPL$QwYX142|QiV_GZg6o;zw;e*UMA8| z?Y>X`6u~l^>mgr#n2kO7VGM9fB)bzE?>0!3_!!0xf8SnzyBhVvenyv2-+zpM6E3{< zlH~5}+im5rnFP^H24w!=@^nxD>4y#<6y`*N(}{Z{2N%$CcGd;}x>1gTn1!&kAR|v=^G3{7_V6r^cGY@C5?j0r+gaA%H8lVva3>iZ9unOq>*swmMOj$T(iw;c zoE(1KLJTF9X_WME#~HQi9p#FDMw55v1Ptu`)I?{|8GsnIk(_b~Ox8#)$7}#gEHR-3 zw3MVpNgkpF!zHs2WPl6&Z`^XQ`_8Y3_si6RGz~qTMah^(RR~2-l6Onb4!W=>DB?*{ z0Yx+FJUg(!3vp@g82Q9alLEY%Z9@Boi)Unwtjxb5(sZWz^UK~EbD|o5$45-;n=~fm zBhH}D(+W{Os1Bda2^{A3Vb1mG%e1D_TjGYFCE`I&qbW9#vOn3iYr6Vm%7?>#v8{c& zhbl-mLdsB?|EI0S5I_t1kWDC6P=mYTJequBVqbPzkY-hsU^cc^QkY8!HzKSz3(U|J zvCL7a5wtN>E)f7q+PRW{2ef(k44=Xrf};?HXvtx&3jwT2b1P(biNq!`7Q(h-QU*L% z%2?m74Thsr@C`)&wm@{ge1=!AInBO8agcDjad45Osip#P%d0x<1KE8R*yMo3FzGM+ zb4$WSj%yMt^$XN&H1e6W0h;8LVOE+&xRNv}&JPJ>X{oQNLW}Bd1mU@@ETz+BUBZe- zAf4SLIYKJFLiW}$FOXr@1e{Sz;jMsWV<`w-ZS^7<`8i1)!{M&6>&eh% zB{jzOE+X=RrzvovKy0Dpp8VuJ+RT3)su}xKcXcNAnGl$|B!jZj4G)9x;eB3^5M?^=TZasfE=eFqiHGbp) zWzW0?M%~rzPy&7KDB+z+f=FSI2{L+*YLdOA zYlT-$byt5v{sI$}d<7S@=d+x{6r(TgTC)(%hfw$zq!vc6bzm)$rw!8M8huFyIa<<> z9+Q?8hyN{#xB>Cv`_q3NVDD%zK7Y@f94_hIj5`Im4e`*jI?KhADX4w8&5(3`Pi=;w zJNG^}gLm*%t@BbF$yL8}b=sLxYUq?h0v_70R;OSgpItB7tzNgnz#2ia*w2sh1n?vI z*~5Op*Y$$BFt;oRza{R?#Lv1^;tFR`QjS3{J3}vo zwD(E%#+UiCfM*{219VBbCfTOC@PH5i#E6^Vp=nhf>(3p{_f+w1>ke{F$}xz9DuJ)! zc^YvFC$n|JRpCm)?_^T7_0&HDJ_?d(+)NIr05a~ zJePkAkgV0n6+mJQeKjKzGl_fG%YKP#BlQ%G@Jt2V9JeG&G^k{z02ZddPwFM+56ja? zisIZXd8B^`_<{>I2hmi{&VQ$>z5L3{7MU|05apO-eei^*fQ)veWb8u3F-rgjx9E7m zNcRh!$%i*MAgmN~$;ZlRQY`6ojJ?Pysms}mK#=%(+(n-t$qrP3D}4*@g0||hWFxyp zDLv;v`C6$j;C6Z$R)vLxP?!dwS2p9dh~NGv-Q!g`FB)fx$6TSFouW#?Pt zuBKx?ntrb(EP_Kbt2N2HUcJ0>f+sVFc^6jo176P|u0 z5iKdJIBq8J>|<*t4Emn;mK!wzL=^<`Gw6T8_{n+2heY|rff1`Zea7`mbF&D9CM|D$I+rCsnU}NL5_bF z$q_0qmOW3NR5IW9tWq)$dU{Ie4qloSZ&* zi=1{Tvz^bd#QpRwtcfSabvj;ws9G9oL>}t2u5%3vd&hTfg8Pw#m_5t=Ma85^rtYY2 zvahFd2sbB?3YsfC2Gm|GL&SfT0zQtIOQS|cu^g<@o#D*qRrbs4HdUFnoTxq8SIb3C z*<3K!eLRk$rD~3&oKsGRCuM5jbWRGszns)!!jD{%XIG8Ud<7f<7%x)wDPZNAI3NgX zPBf+?yJwy9k0sHQuiz^Je+Skk7v`(U0xmq{)x}FqUy;W+%qfc61xbH^eQWRVJ|gUn zx01J*o&}{RCn-}9mHXxIbR!7wwvmSU1#}GglA`AF*EpE=tR|K51=UeSq?JL^)Sg_3 zv?`=H86Wr{<>bw*5P77!_Z4zS1{2S)FQTF4eG+l4bznF_Jmr*+lEGwZ_45%~Px64A z>l6#qSErxuXp(#y{BVCrm-%^cOih^r5$5OE1~^8S?_pdNP+xkunzFT>aloww5Q;gf zTMj)e4lL2ayJ>mRP6ZOBCY{`IWJ};8&HaVFym0)ux5V<=3rfQ4`;VKneoefc9eMyL z?GJ7LEc|Fa5w&uWz2CU zXi91H_u*ibf3-T52Vpr@8=(v%u2uQTEc)Aiqb%}bi)DwPN75ZHinvJA`&2?6{OZlp zPccg+m2J}o{i%Nl%tSMn5S^%#&=6daTZ}%O&)~|-6#aR9+Cj1ytBVIxlw82Keol10 z=H*@A$HJEHL4+nriBjA&-!QO8HqCff?*3mniK$HEQ0XB#UR(Pkxj9x^JCa*Y0j=>! z`DQdWyxGr2SK6VuS`)#XY@L1D3~?rnmxvwW&nW4M2Sp{?&XiKDuggiR(YN|;4gJAL1bYJplIgQ zQp?LuB-Zzf2048Xv&X8BY)F|)sqkF0^%kq2!do^$#T%*a+FOUZZfagcHR<_#1+{#? zndG82kOT)O>+C3(T_P_wXWW!V z9rKt@fI^`9oh!ZMgU9bW&cbZ+FM$>~E$Z=R>Hgx+EJXh&Bz^f$k>L;){|ZcZM;VFr zXQe{0%CGeM@H$apM$(CN9bX`)nt1owR1TGY@ZpQAe*mtY=)q-fWOHYZMwqj46a?t41ElOQjEHr-xP28(7O zVFo}9nn`4i9q^^Ehf>8Dne=7hx{L!G>6cD2&rxom_NB=ojC;~p;*bwXi!Og9NV6l} z@k=&#%l1crpZvOO#mfG4PM;ic=|?$NM_kS`Kp5GOdF_cvFr}*TjLPSSN7>qiS=_<^ zwdN06!sBI#r6ok_w=Usv@Rf<@7^r*+%31aaqg0oVQKm3f*+tz>>IuDX!aldm?mM^)ei>v8aFB(peUD^UL6+c0e@xZU6QhhwqW!jT@VIOOIdsO?IS; zONyRbA?=xUcj;T!QnN>r&%^?moz_Op`s@;_-ujxV*RK$OTP1{BQxN?W)8Z%$3ni|0 zSWYZyyVa!tn<@!zcLWW1Uxx8#frfMy$XS@7+*Pu#x1o4Py=3pIJM4e<&^9raE)sF( zqmHF$@VOx8`}#URAaC@jNP0tDdSDZ2>q8*LL@3a9JTFSIKrc-os&F7m%+>ov>bYYdD+Tn_H-%^L1Nn$M+a0$K9DXCmAYMcO2>yiy_&9R0PI}P z(0Y1p3*n9T%Y&r(blrb=^+ZW*@qmCoR7x6<-5#J}%Jikni~ym9B($4J6dZvXJkGCn zNI=D~MMsm#WMGIsxprmc_GGbOB$g0cg0aE|Q*xH!bM zI*-1TdQyLQACv-kG_aGOIx&aN0V0_Df%K^+=>Hed7Z`zZ=+1@_t@QvS*7J`ELqzIJI88^KJhd4qU&JtM&R|J@qGvJ03V0O+T z4~>8y%0(&xP!@mPqW^&TA!(sg`mgj!i0y@+FWvvdmlkiF(Z;G&uSXxStN~ewOIH(Q z2ECidJ8fPq&EGPs+YJRV+r&^EZpvD;QUazjtrXsFvhfikD{7SPjA?4yWHpmWj2Ph_ z4QDH?2zSZ6{Y*K5Y(el;bv(OFrnkN(aYhgGYr?)C)+c{Sf+xJ4IT5d%dvssUVyvoh zR)8ldH2`Zoa%{oxPy>J%A8;~#89jsF34x&sYF&t0!+f_@kv^1X8;#eNQZ0XrUxKfb z+sFMgjo_x+m8%;>px$o))MK|{Hs%or8{9_q_E5xUP&fg1ZkZWl#4MuQG&7Yc*~VuX z9>Rch3qpTxjnwjW5f1{@wvK6JlKv9lW3eWab+?jh}mjr`JH>D7+Vj-_(rGA^!voglS)YQ^SdQjoxZw^`_uJPJ} z<=ubX4l5c|IiYxmo4f7q;co86+m~R4Ji3)Hpedk#vvHeuJ^^>?Ht$sK(UkYR+q`jO zlzN!ky!~)U?jdgTPPfCuUSR<61U`fV?y|v^e3#3-N zS?w;k6>z$dUeH7bwz4HE;wj7q-v<((6ip$tjY=-;_76ly7si&-j11r?n#Ji8)>9A> z{$#h~%Ry;YD2iaK6?eg$7;KFjwY@1oJOhmw*$1r9Qv>)QtMfKbET$>h-i4XuB! zHLoGD5UB!sh-9$ZPDl#f-t?AiMci1;oG%^o<^0U5oSCxf+|9pZ{-W~R`!{=d${f3v{aoU3?(csV7x0o7(4$Lvxvd>g8wGXcl+@+{*!s-z8`z~q zFMj@&#=ZStZ(#p<-oP%ks66_4JKEZ}TToMZ)a@^|RZH5c*1~?;{k+Oo!5Ck{*u-Vl z^NfaEEgtwCeV8q4*W>)0Zp}USCZ+ZDeo+a%eN`>3xF)Pk#|^{%{YGN4MGb!$7;^0! z;a87dgae3EF~avb_H^f)Y8P?-FY2aR-abUo5+iP5^2r0&)EZRm$kx;4OBjtbJQIHI z0aIg;)eMY}IOkYPEwJx@{%cS82l5UWdIOswCOjA1n ze)P|Re17%l-{8Oh1Kvdawv(sq9sxCz#_UCZ$&wtm5%rnkNAyiq4tELr!U!K^TMjRh zLl#G;kq_M2a*IPz98$D@z@P0*pa3LLs9t79422GA)}XNF^73V(e%!E$Zy5eJ|NqGk z_rgCneDm*nKi=c>KiK@o%^xoAefEv<8$Wp?jM`j0zsKkHR&1=ZliytYaPM;9<5o|9 z&NHnB6Mj1H<-ace2Cs8r==)O0NjT`@;wilM$8pOh%a}8K_{#`?;Y_IE^Ze`Uacd@R ztR5~%%P-*@>)Q(KjET5Rg&lFmn8De>)TFh} z3}X3Mz*1ocITvaYe%OVlPIxoC8pQ+~A%^c^epV>$hG+1LIjhw01g8hTj2qsLTQwQw zyomE_j?Ya}d$+mx^4`U7FY~Sn<%X9qmKWYT?0BnSQQZ3BQ90Q+xNzm5Y5K(&@$b+j zcQA(!<5o=Gan9qwj?saYn}qcFlc4Ste=PR~i_@@6UI|zm9D-C*_Le5$6u9?kgrDF< zagKI)4PD8JYY}>d0gUD*Zdv@fyI(CmJ^l%ILBgP7(%OUEEq6#THSAvDW61@y{T>7g zjzOx_kMuH)WTcCI8Rp=vt0J(l14^*TNHfHz&e`EPJY*V9(!&Ds0FTOB!G+=6zkdB^{J$ktLGF)N4t7fvOuL;2)wg5b>Z443b6JHRps9B)phJl9xDu z<-#owOrk3G4$KeaX+AL%Pb^O0f53^we#vNivW;1Oai&|a43jjQEuSsQtA)fWt*G!m zMT2HQTO1MKD|`qx!`tC;R0?at7C|%UPS2A?<4yFZ{DSudY8U9QeG@LOW)iq{^}<#^ZI}N0@{Z9TA7$ zlc`1gRGX9~7?+_yw}hduJ*w0pF3xAMvZ{!crEL$1mDY5C9eJy$3sF-PD_EY;^zXA+ zaq_6P26OjOtHSPNr7ANbe_n4s*)TGI0{Cr6!$UHrwMhyEhLtD`QH^PwQIgY$yE4Gq z@{O|_&_h!q8PJkC3 zy^vX%^8ZV5tPnrG*Wu|X&E&_WScIhq9vSbje z!KmE5z_!prwdVU&e{m##e|oYsY6wr%8Pb^BUk=~BA(rK8lptgpgyKC;gEQR?I{+;4 z#+sn71A+Yp4j3%B=hP|^QQA+l<$>c0kr{*s0F<}E#%lkQUgB`vquVu-w6jJmB z_ayWzg_|s($QSv*i4}+k7bODEllzHj-j^6MFzMmC5{P+(f0qe29fcov39pYx4h+c8 zk(}KA5~fxJhq@Kn&fG;%*ILrs`^%m#Ftz!jbL1aF3|_< z%WhwOKW?yYf52dReob0OUfUDlZ;W_QUy8@btjJs+jJ=KCa;+xjIOu?nBtzJ)|H5hj zaWNR}b0lT7*ImxVLj+NXi3ibO!HUTy#NE_jG-mf?YA#|KYE$&s&UPLUrVX|k*`8*$ z1W*)mrTi)&aRn~;dO69KKP)OxG=QiWQIV%T7fKeHe?mrA5jEGPgU6)n<37gU7Eha4%->GqmJ0`Lu5SRL_-uJ#q`` zti?e;hi9@Qd%7p`go*+-#emoRZpf{dveyRP-;H1jL;mP{Ir=U6HF85G*nm)?B>iUDSNN`i|v9TantP$32}j(8DrGOHB5haI(s%G|CvA2Jv*)`WE{ z(B;+4zJSe=`C;c70EK3M(!d#0`sxkD+_JX|%K&yl3_;W_>hMGL0sj zf1_=YVq$35%_$&6GBH^=$6YqA6lPT%JRTFw`(5EXcTXvI6$w2y_qMevg{@6 z!2eZ_%7P_Dw(epCyecj10ZBB?C|1RUr^2MZX#9`@`-Fh@X`@`xQh!K z0o9*(p$^DfWaTGim?iijHuFOzf3Y)#@`fs_lJZ78qunaKqJ)K)3{~^TrO(oo_>LV( zm@e&|u~k~9n||*9(IJwct*M)gkv;!Rn;CeHH=J(5Q&3B*IrX zUx>i1h(RI}e?O<7WOzSHv~Ru`4a7>6X!95rwJ>0zf(RBRK?~W@f1k@F-Ukix z%uvPL?`Z-rqeb^i>PsAgAIlIBVJbGIpBAX>Q=$6P_uO;Mq}Oitt=Wk-!eY@tPR#vwikzyxB@Y!Z`G_Eq{ex@0XR0dze%7xpx3 zcZQuqc12seYUir>)~HFTXiKcOX(qt{nyM6Qf4#IT+|{JuTyWA`XQO9FaDXrMLZ$vm z-^$3f`Sts~1XnR~f3{N2evu)O;e{Tq7`cAf7NPy=VOwc`(09izr*v`3`J9vRK8x8y za*1nxH_Tor=dyif?_zWk%(Sb3(gH{sYwiL_C1gNr@%s4Y$tc7GL0Qh5%a1BY3?T3s zVw9R+R=J^gE>_LeSpPCM!w)swZQ zIJqLEtTdumch|M0IBDb7YfHO#6*%}(qrzPUCsUgEyt@i^%y$*upvYJO7;e6pfMvqZ z!(|0+$fHF&f4snuZ#R&e4Ee7s7a3GQ;0G@<1Zd~5I5yLgSkw6!eqs`V`5D1zT~%dP z5PB=|G70xpl@B$tPr7!zTM3<}!`-?=b%3#*@5GbX!Smj9`*A=YY2pZ0hZNA66%-Z? zaaCW`PvxQ?cWTnyHRLRGfI1qhxSTXg(26tocL02he+S!~rdqSW9-<3T4ZkP`4{Lcw zJBS7USlC@N!dcnEzPr`r%t|>*@qaKD$JYHImc|q0F=5X3O~@RXIWU77U^fRjS-;_c zbLUL=If&uMLTi#pdc`3Y3mZV{ykm91dyrxyk(w{Kpc3M)KOoadnKG?cC_M7ogM(gK zg8cC_e@OY#9+PXuZIUpT7OwSO^V<=io-*M0J{+njkX7FI2IY|5wkH~h1_#!~=_xs~ z=51R8fbg+fUO=5(_F{1_gssa1U`-GV!NTH}{3bXF<9U}kFoCNQn{94Sf$BklU>Y3d z83K{pU((w70-1;i#Gjn^EA3Q8L!F!I|2Rj^e{aJh%)DGL&qpjZNk8s#YE*;O&o&C+ z636h#=_6+eC=Jx?^s_AOFB3x*F9yimk#aKdxbQ5ZhzCJ2<-V6EujvaWR)v+rx!oU%pQJ-AeU#8;? zfBQVAU!L~Ml&`x^pljt>#jF47Onunv{sGZTq3NVUtz~#;<$P(#-9W{l_Pp;Dj-4v* z*%>TNX@8k6$1q9M_o}M00rP(m1sYBqc*X5AEg zArfILmp5|7i{vNOe$7lREn4cw7yZb7!Ep#iGwfLNg@DyJq+s4yQpd7u^8T2R^R$Eq z6!Qf=k|Q-@-H>68Xb_mAzZK8m3eWJ<81lJM<#OBon7@?pZUm$2xu73Nm`MjMf4X=Q zzlsp@T*t`-qXHLjb234L7&x(#t2!1L6+U;jyaApF2ApELebgfuFeKKZqBDq4RvE6y zeW(-iOSiolLKJK5qKrZDIl#dgiRBSgf6Dv&m`CD#jg-#gyHgmBKqfd}GEu>988#`g zev>C|JYSL(&BbW4YpCG{osFh9f2)Z`Ryngc0afzhZQWkjL3Nw26M&gSwWeJxw@cZ6 z+N+K61H%z%hOn2phKF3gj-HC)O{QGyV5!9Km#S7gO>=#LWXyVD%8T}@N}HmpONWB5 z1S#+5?DvELPR>SolG1sy@Q~zgIPTVKs)_wPYRb;_d=U}uU3t1fUuDY0fA)JIO%F3J z9AxZ_>le8>4arJyb$u92ans(c(h4sOtN^42OhsPJopguK#qhYfhhhArB~@M0>mqlD z+qhQlh&H=TeBQj=5kD59pyN+VVuUE^5BShy`~fjpJ$^P1nKt#O#{h#LS0<^_P1iaZn)|l5V~h(fPJZ zzMk7+g=IMcoY#lZAxL+};(Cth?lZhaYJLfe?cXKO%nd0=3w0W>0~LiclYLIoFKKjR zNAOS$gL4Ds_Kb_-f8ES;{n4w$1C(4EMe`&cO)6I1yWABQK))~4f99b={RMY$X?Hh1ghFmDxJ zgy>d~V$}!=95!z-!6KAQ+{65Rn))Mb;ZBtrSeDT4Wn8>Hl5Y`#xCgy$Qe61Os8SZr zFpqM3E@>Iad z_RPyucN|SPHGyrp<7mR+_fu~7T_zGuc!t`%V*BVe^tqQVrp#gAvOueiGHAtdL|cwr z@kn>ca(NHu?H$^<9pdJ+inKli631_M<>H)^Vh)Md7Vru&^HHOXNHlzLaqpk--~Rzg zgV|q|AORf$IX1OF0p2`+RZWi+HxRw|ujr#2(O#F!WtWcxfy62VNG2!5fh-UW%UZxtn(UMBqdMzjzdLp>&4Ynp5)zv$tXRCV`D+-R%@lKgHaCN_ffD&-t_Pp3 zaA5QIFnSxX8O*UPDQh;L%8E12U+pdtU0t7Tjy{GGUC3&1zTf@a?e4ua4~LN~Iq~NC zz}BUp!TdfTS(1>?_hWAFBebe3HSFGWr~A$#zsO?V7GE(6#gd{}=TkJF2X94kq*6?@ z@K4#i9!7zioY@?IhLK@QDDinQkSI=C+Z&XF_?iU3ECHHWsF#avYF7Jez?FT@s^(AZ zx?~**Dr4)(0j}nJe@ZbrRBnecNJ8;y=$xV^!>3gy+`oTL*26(P*X4W*?{NlfB|L4*+CyB%yfDI#WFntO-oqR<4^C zsRUe9(cUZ>f2sGY<&Qiqql4z!i&Of({O8Qq^ldPa9PO}{-A+FrkHVL59Z9r&!cfM!E1Ge)w7MKbSE0u(E& zLEZSB9VFB;cSt3mj-A)*JyAv?H&t=hwIk-CBdHM|t;I=Dlcr7aZ;P`)S;)~doKK6n z12x=d4De;ONe15MnT)4MR@5csc67HBq!~stQsasC=y>(iR$L~AB+zN{GbXA_vNCRh z9S-~Ek*tU}EmO}6qP~nxop8D{Sah{HJ>nJHw8EBDA&nkQ0u{zHApRRXJlb_H@cS2| z7Wf305kUbP12{G~w=h8g3UUZDE;cndF*GulLbCxg2sAD57Ge;2nI`0S4BipmnGT(ssc4Ix6j%E z4*>}ZHZd=AY+-X~x8BFd#HFGPhyh0h%iYNKHgfQB{`#3j(K?tqTHxe=|5S zK0b4Fa%Ev{3V57#RcTZdR~oINp{t@uTvCZ9rm8)Oi6~KlB!-FO0!9o1DoRuou*GJR zu5OxU)ebaBL-WzHBM7}BO#=o5RE%QWg3(CSn1zgEOhz4}6LidS@`^luOf@ExIdf)C zEvM?7`fj~nv?`NFN=_C zlNCz2R(wg5tyd+aD73kJW267z^)B7}$`wh$|cRYU|aiFlHjLQEy55l<0Mixd!^L@Dt- zxsE*LvEJi|#7|=IoZzYRf12X8-RqRpA^nD$O@&ZrXm5Jg6BC6}AE6Ii!b~fPh#{0P zBoWr35(e%;o3@kh-CC6bb3yS&Q8r!nqZH39@ZzZ}HG5(>+(N|_Kp4wYt)?^31}-=Q zt;UnyFSw}H`EZOwk<`yu&t1IOxG9pw`TzIZ&7gs@cG^eK>Fq**f0b+>M(x8AAs(eN zi1Gaq)7aZ@7M_Zq`sX|(MH5CIps6fQ#c@m`#5M5)>SJZ+;Ry7=g^puwNBf()pcjth zI(Kx$LO3KsGDM_rRc(pQQGoc1$&a`B)5pJDm!LP1jv48@`&Eqz$|~5Rl3m51rg=v} zlhxYbY*)8!d!1j9e^bnsA&C?}`#*NwW;zS);Dl@77~e@Y6>O~U=6pt{feW1!&khpJ zyn!~@CBl`_KQm)a-YWfUuRf6*Yo!uP^&xCFFKsd15iz>pHn_$yl+>%necc=~-KOr> zr7q&^s3}vrNT2!&8;I(p2Tj%eESe%kFM{ixMm&qQOMOObf3sZZU6(M&MV`4Qym>=1 z%Apvhx1#bA>!Q+)D?tf6;Z;zZM0mb~y1h3)iCax!DSoCXeKB5h0`;OJXkGuOw|)x8 zZ?&?f0w~=FbR+Dwa!5~^p;BL5X4Dt6F^l=`$g_Sh1HX*}u|JmM1QduDqOfZ>9G4Dp z=A8#qYk;nSf2t-0w^>Z4EeKhfBkgqiMkApWdM~kT6@4OvZOI? z1U00B5h5R#oP)ZOeQD?g{Oj~eoyo8ds-enM-7E?SvaA<=6!Kl<2M^J_>k?tI@Qa)4 z%i|~^z>6qhm&A8=Qwi0 z+z0gg=N0~30)=Bsmo9=O^y2R{7d|=Gcf!I}Z0$%o3bbvo@vq#u%=J`MMv+bdr}wjA z1B`EQ>b98e-6zm=v}9kmiGAF*1y+>OCr$`5e-DczG8)EU?%&!U`-h#XsA)S{QSC@P zs<;F+q8^PP22FeA0hVmos*$H~rf3GtV5)Cw>F#$80Nws`KHjiOFV9HgveOdMHbXQd zcjY|A?a`UqEyy8&e>2rq%K5iA&8wz$|hLv;*9ui@S@G?_d4<@Lf>-Rt&i!VzH zaG;l6TJ?kF5F$kiXesU;d)_P81a=l(f0s6)aIZ{hA(rNC!4rYTlcXY`-j3836f$mZ zMHeN)G_;$so51oY`aYh7rGYQv$$0t?bCCp1JcbBcm7$`zNYty5OJB4~D+gNsQX5*% zS*UUgQ-l0(o7>n7=2lYVzfBo{6 zubf3aJ3G1FzT;3S5-(mxeThL21+=a%zp2JtX*cWIk~i+yxm(979mx%%S^BRET@uxm zSXap9@6P362dwMT)EHvljZK1$^ve5D=t(s3;ps0A>^Z-tgWc4+8D4`(*aNv1eO;ji zY|vcQT-QiH)H8wb$YF-V%}&nZe`>F@0eu@@(&obUJT;fYEA?4GJ27JB1rK~!W9sc` zr3cMeow98K>YDoIK?k?Ecd900wEi4S&JxgfMds zA`3o(VUdKd;ftp8&6hvue>f`U+o0C+T0;~phEP}qOA8`%ut#`KvOCpuVygAI`o;KG zU6r82qYUGK@@zIky*cqkpn1w6a&oqn9y3!(rTKwu8ZUL4M&O`$d=JiXuF|N$QZaVi zP}k6}_mP*29Jr5CZ%Tx4w3x}1mgW_0(C^w)7$v3&jU%K^$AwF+e=W!l${`$Ir~duqjk*JGAPt>$nX&hp9-NXpHd z^Ueb_y9G_Ud%1RT=3zg|@LE!vAVdMz=GCD~p=Qz~1E zkWz%nPGrfNUA9P3GRpX!C*pnY_4)kxV?K_1?>YA@-}AlaX*itNHcejtsJy7*WnsnN zC9`J*_Y~^I#wA6)7@@wRBJt$!JslbXO&ac*)7|M6XFYC&R}=}d+LNjGi^Pn|e*?Et zmKknS92>mv#IN}>%-v(gIxb4xAbFaU|=ia;*yX$;xN1Mhrr<|%Iq#s&T z{qaSP%;UPqCTJRdBY(zr&dyFGVd@wh6k%%{95u;=ru0oHmeA)!MzQC)n#)E5*$|B; zLa4AXka@{DPKms&S|1$HX{O%UcfN^sFLay@?ejB9qkiEs)5+CZUR{>h>>zY*pV2G1 znf2DUFgGuIN?z66V{<*?zC5hDYcn+V-phhjLFcxn$o_M%r`o=5)2UO(B`y4K#K&9k z30)IvE_)P3y>E0%epvtE%lyABgdI~CJr?9edk+*8_Gum(-Mw`wW=oEf&((32zKb91P1{hS_IB{>b*`fd56$-EN+TLpT(NpsR_2~f-@BM~ z_@?ufjdH)jm4|F5rR7xywx56SW4p{HH;Z=DmAUoR7v4E<7#AMfKP%5KGpi(Ki^ITC zl|5&k<{UrzDnr*Zx@GZPufqKk_qNYGdvlu3NBxf-J`!;qj9tsU{ND^J(ok$X zy3^l%e#*GrnK!atmyfZ1GNSLA$WhxV*BF;m8M;Lt+t!Ydvr}g*nB8)Fg4#~rY)#7h z(XA7T$pP=jhU|ED-naV!TdUfewicU{>^j$wMuFYq~Uinc@As@XcA;+I#8Im?WRy3U;qp*)RD@f>!UO zr1P<&wWFjmtUM5bKz6WvE3TWjIL}L-`e)Ej=Ry@>0vJ`XZvXnJ^s*b%9r8A z^`CsJJ7pVu9PKY8u3DEmOgn3J+y1NDd9JQ7d{vX+c-sDGNUge{`RC_F(>_MtOFArdsac9*So`K%re*2Ci`!@TW^)p)8gr!v<34nE_*e5Sfp-7i29TQd8+O4lq+mb znjt^fPRG&UMb78KfN7rVJMLWSV|9o0F890bT%UOO-^kxgtSu%O&JKuVJ3byw*>?`U zlF)gfK$TQklpf1zxWu%OCt^e=v5>fHU@w(`9cG=R7GdJ94^Id#k-(Wn!<&pEZ&mHJxwCf>^+aUnNJ-s(kqQgTSy2}`eH{`o!^9a{eM zX}ZrS&nefd+NV9Z*5X)oP=3Ry>BEXg`xD}(rG}@xW7GfUj#T=OIy^alPNmsIGnb9S zW1f$(wI4`^{ysD_+tc{ zPxW41PFirgd~{XypYzXq59r>;b*QCE_wA!Oa9DKO$lA*jm_ye(X|>R5 zS$v>vkGgiD>;JfL()_~9ibl(69;DN_*2Ydd&yQgYj<2)Z68v=6aJSSmDYG8T-G1lB z%!7x0%vWbzOMIr^t-Z!#8~sHm7a!Z#?E$rUz*5(}H7$2W>Bh-FYnvROq!IMg>fyN5 zJgYLDH)&l~+z-lGdAH4yt5+AT?QZ&KxcuFPC!WI+cBj8>)4*KN>bUTcm%6aJz9zl# zmQ`7tR@)ON9r7cDU zr}ZD*Ey;ly__yMgFQt3XYII4Q!fJlwim}#PR|oFf{h=y%U&d{&-StONHlBXBx9_{B z6@T6}U(_{^tkv3cXzY3;4TDvep3FbsGM--~Cx`u6Sm(X=wEK09x4 zf7wob^4l{X+QAt9XEPx5O-NGE7erp4ax|95mB3Bw$XssUl#Q%u_Bh%CG(`*RLabM0PE@`09Si zY18nMX#s_WW0rQ=HqkizuFGMK;W_Tyy}d@l${1&_Rl4gAx*BY*>ESeWTh7;+?T+cx-jq>-D1`h9*(cmy~Lh8YWtq$;qN938MLm2G?Wwn7>@;EZ z*-ZI}3-9Zb@-q04SI$dnyQqczlkMN1?zO(Vsgw|o2t<^HgWwFTh@&|G2i!D$~U**AfIGA$jY7Fbv)|SB#X9IrR#1Tk6ZMz&5~8~M+P?#wdzAO)jjT{-|w?**z)Om zySdR9r)ryzc$xh=v@WvT)p^d6GS&(be+L(fc#kZuoPz ze=8l@Y>(EshT{wU#)rm~N3L(x$1O)aDCeb+`=z1${hk*QqG4X6{BHV`$u8Ato-ay+ zEhk^tW|%%PV(OAE=Mp=_pRwPc^IkUbkD>bByyc6s?FXH zHy!?;SNDD`GR)M!{`+J%)Ah*ts&=vQho=TOR%ht$3SY72^XQxl)@ECJ_P*D1+kJ-) zvD@aYEq}08)AX|C`0+4BJAAl!IjU9`jq}hY==hsG$GmWajx;W;)}9@D;jV$v+Wj@ux( zP1=p5%?RF%U={CEWTf(+O;Tm2`gzj?MG@YNiMN+G59l#VFd^WViNFyY%TY8*6KsD| z#mZD!zq}#aWwT|iBoPnC?Ms*SX#PXq)pVJWnw%pkeQ7tLAEd7-{g>oviWgO7gdkEP zj2xw3&Fct(5jk2_Msgy>Bj>?h{uk%5Ow)OCr2|ERrjP?+vq1BzN3!ByfCbLSLzsV>9^_EDo6sDMDX!LR-$;$!v zrDJ}&Rip@9F-78vDN1Cpy^M$%!;tuSH2d9K?UfY^oT6r@tdkTI#3q`guyY|zunW-) z7A=wH@hH$Dejkoxzv<|sWD`RZ*dKJ83K!42$bxc>F$$Kx6k|Al1-&vm|~W zmg2u9hFr-`mZ31Gu`G+*&WX4ZfyGjnrv>ahj%HM{RLWhEW3j>%IYGcumxu2d=vSrXc)cuZZszN#dt5@}|y3h$4lTE8x%1 zh0W_&0y9SB7+l8Vg(Qk3;#E?*KN?w5st{};bU4HqSO?$nRNV?P_7{733NaAY!N*G2 zEBU|K{UiZmw6LJ1ptpl{lz`w4$QTO2EReC^fmvCH3XV0!9~MG3pf@4Guc(Ytj@>uX z0gQ%0fDsna0=^fD2s=;4af++iQs=@v1WgjVBn=b-Jp-r*Qkx_JFoWPI$N*7P_Ywlj ze+Px19YI>iA%;qk#;S?{#)cXQot;3P!IB)F5;Uip5)cPa%yX#K;9Rr^l8i{Ph=H(I zdWWhY6i^=U115t;fyrnBYX$&O7hFam8i%-$^fb4e- zcsdJ|;;;R(J|?h1z$y;NGZt$G*oQ(bN7kvch(-fSn^m~NU&*n;r$|ovA}{^=QGN^| zqOF7^AnsTuvlQ-EmJ%^6V43gb_cObo>-_KIs+f=$7gOLVs8d|1+iyYdC zK`~HJIz&a>#T!ub1VI8oorgg)AzzAH8$3TkUd9?VX+o9_SpwWaPCl`tetxM1)GN{vl zL*RwM1DQf?7YH=ck}8xiK*JEe=l0M003?%mN&ut83JCBTO*;Z+0E(K!a38GRU|NAb zo*;866wE+9fE0Egy!8^mqp>;{nt6>fNC3IU?1YhlN>NL-1&@3z2_Wk+=K@c~5=ca1 z+XNCoXt0}z&F}di2?UBcJncme8`+Q?gL_8gNkj)Yh7i4x)Bz>%pip(I6n<>_0rwEG zk`Vy}RrkVqBBfZE50S>IL?X(lgMe~kyhh@yi0#k>DVVXV_Gyq|PbMHh^tUT2MH+;d z?a+v@ssniYg<=lNdH^_#|1dD`SHUQWQ8VyhL>(@H59&0jsL>RY5*`T*cmW$}HSq() zGZ1Npl4F}md;|mrqAdatf`@yjQN-%aI zz%&Cr5Xuon49P$gsU;SG=zwQEY{sqTfIFgz1dR>NZTKlpAZ@C9CA2}|6V~x0 zCr1{7g=qJLj{l36PxA237b2NXVlz}wf(-{iCjOur^6%4Nuu#Bo2}lHbW62n_nh74> zR9HwTY@`MercineV{M?vn5!uo>q?Lwlp^scKnCt{3fP#cGbEd^l;9asz#>R-*h_$< zAjgpNVAS^g$)pS!lz_bppnOoQ0@47UL=NNy8l7uOrH4Zc4MQvCy;5_K=KRQwK}O;S z(j@wRAYk%>Lk3NuOp^cuF)#2y&<0?T0riYv6)0guMQ9O33v06^43F_6uo?vYK#i0y z4%sfiZxv@W`2wp!P!HX%X|Hm$Xb$ZT(t%i-2!Q&y5)Py6g76DGhAEUvVCU!tFh|1)CnKnqq7++J8p|nL^G%4?F+39Rbuhcd_TqCw=rZWH9Z>w^W0q1A zWdY-m1Hq|S>}dGJppvFwGz~GvWE`LXavm*k$QuosD17Kqedoe~4B)6roFo}MeNhNz zL>A3$@G3UshmUOpz7EQgMVcrKYXCD5t4SEZW9|oDg8&!wOMyji3>b|W4|GQ#}rlqB3c*OLuh^VlcvozsF1H*F7zI`1=Icok7 DaJaz+ diff --git a/doc/user-manual/fpakc_UserManual.tex b/doc/user-manual/fpakc_UserManual.tex index 8897869..62e8f50 100644 --- a/doc/user-manual/fpakc_UserManual.tex +++ b/doc/user-manual/fpakc_UserManual.tex @@ -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,26 @@ 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. Accepted formats are: \begin{itemize} \item \textbf{gmsh2}: \Gls{gmsh} file format in version 2.0. + \item \textbf{vtu}: \Gls{vtu} file 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 +538,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 +606,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 +621,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. From a22099ee87092b08ab518bcc536ff80b8a86bdb8 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Mon, 6 Feb 2023 19:18:04 +0100 Subject: [PATCH 11/15] Improvement to vtu read subroutine Some improvements to reduce code repetition when reading vtu mesh. --- .../mesh/inout/vtu/moduleMeshInputVTU.f90 | 191 +++++++++--------- 1 file changed, 94 insertions(+), 97 deletions(-) diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 index f9588a7..3001eac 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 @@ -169,7 +169,7 @@ MODULE moduleMeshInputVTU REAL(8):: r(1:3) !3 generic coordinates INTEGER:: fileID, error, found CHARACTER(LEN=256):: line - INTEGER:: numNodes, numElements + INTEGER:: numNodes, numElements, numEdges INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types REAL(8), ALLOCATABLE, DIMENSION(:):: coordinates INTEGER:: n, e, c @@ -223,6 +223,7 @@ MODULE moduleMeshInputVTU !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)) @@ -238,6 +239,7 @@ MODULE moduleMeshInputVTU 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) @@ -266,41 +268,44 @@ MODULE moduleMeshInputVTU 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) - SELECT CASE(self%dimen) - CASE(3) - !Edges are triangles, type 5 in VTK - self%numEdges = COUNT(types==5) - - CASE(2) - !Edges are segments, type 3 in VTK - self%numEdges = COUNT(types==3) - CASE(1) - !Edges are nodes, type 1 in VTK - self%numEdges = COUNT(types==1) - - END SELECT + self%numEdges = numEdges - self%numCells = numElements - self%numEdges !Allocate array of edges ALLOCATE(self%edges(1:self%numEdges)) - TYPE IS(meshCollisions) - self%numCells = numElements - END SELECT !Allocates array of cells ALLOCATE(self%cells(1:self%numCells)) - !Read edges + !Read edges only if mesh is for tracking particles, not for collisions e = 0 SELECT TYPE(self) TYPE IS(meshParticles) @@ -309,9 +314,8 @@ MODULE moduleMeshInputVTU CASE(3) IF (types(n) == 5) THEN e = e + 1 - ALLOCATE(meshEdge3DCartTria:: self%edges(e)%obj) - ALLOCATE(p(1:3)) + ALLOCATE(p(1:3)) p(1) = connectivity(offsets(n) - 2) p(2) = connectivity(offsets(n) - 1) p(3) = connectivity(offsets(n)) @@ -320,6 +324,9 @@ MODULE moduleMeshInputVTU 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) @@ -335,6 +342,7 @@ MODULE moduleMeshInputVTU !Associate boundary condition procedure. bt = getBoundaryId(entitiesID(n)) + !Allocate edge SELECT CASE(self%geometry) CASE("Cyl") ALLOCATE(meshEdge2DCyl:: self%edges(e)%obj) @@ -344,7 +352,7 @@ MODULE moduleMeshInputVTU END SELECT - !Allocate edge + !Init edge CALL self%edges(e)%obj%init(n, p, bt, entitiesID(n)) DEALLOCATE(p) @@ -358,6 +366,8 @@ MODULE moduleMeshInputVTU !Associate boundary condition procedure. bt = getBoundaryId(entitiesID(n)) + + !Allocate edge SELECT CASE(self%geometry) CASE("Rad") ALLOCATE(meshEdge1DRad:: self%edges(e)%obj) @@ -367,6 +377,7 @@ MODULE moduleMeshInputVTU END SELECT + !Init edge CALL self%edges(e)%obj%init(n, p, bt, entitiesID(n)) DEALLOCATE(p) @@ -387,112 +398,98 @@ MODULE moduleMeshInputVTU 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(self%geometry) - CASE("Cyl") - 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)) + 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) - 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(meshCell2DCylQuad:: self%cells(c)%obj) - - CALL self%cells(c)%obj%init(c, p, self%nodes) - DEALLOCATE(p) - - END SELECT - - CASE("Cart") - 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)) + CASE("Cart") ALLOCATE(meshCell2DCartTria:: self%cells(c)%obj) - CALL self%cells(c)%obj%init(c, p, self%nodes) - DEALLOCATE(p) + END SELECT - 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)) + !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) - CALL self%cells(c)%obj%init(c, p, self%nodes) - DEALLOCATE(p) - END SELECT + !Init cell + CALL self%cells(c)%obj%init(c, p, self%nodes) + DEALLOCATE(p) + + END SELECT CASE(1) - SELECT CASE(self%geometry) - CASE("Rad") - 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)) + 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) - CALL self%cells(c)%obj%init(c, p, self%nodes) - DEALLOCATE(p) - - END SELECT - - CASE("Cart") - 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)) + CASE("Cart") ALLOCATE(meshCell1DCartSegm:: self%cells(c)%obj) - CALL self%cells(c)%obj%init(c, p, self%nodes) - DEALLOCATE(p) - END SELECT + !Init cell + CALL self%cells(c)%obj%init(c, p, self%nodes) + DEALLOCATE(p) + END SELECT END SELECT From 515e5c7744a3b6d2f73113085bf35109bc490451 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Mon, 6 Feb 2023 19:54:54 +0100 Subject: [PATCH 12/15] Fix an issue and starting to read information from .vtu initial files For some reason the connectivity for collision meshes was not being properly assigned. Also, the first subroutine to read information from .vtu files as initial states has been added. It is currently giving wrong results. --- src/modules/init/moduleInput.f90 | 5 +++ .../mesh/inout/vtu/moduleMeshInputVTU.f90 | 38 +++++++++++++++++++ 2 files changed, 43 insertions(+) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 33bd894..da43e98 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -991,6 +991,11 @@ MODULE moduleInput 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) diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 index 3001eac..c62dec1 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 @@ -508,6 +508,44 @@ MODULE moduleMeshInputVTU 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, ' Date: Tue, 7 Feb 2023 15:19:13 +0100 Subject: [PATCH 14/15] Finishing implementation of vtu mesh format I have tested all geometries and cases and all seem to work perfectly with .vtu meshes. Input and output works great, starting from a previous case also works. Everything I was able to test is okay. Still, what I want to do is to change a few things in the output (e.g., change OUTPUT prefix to Step) and try to improve reading gmsh meshes to make them more compact as vtu is now. --- src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 index f58564e..763517f 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 @@ -334,7 +334,7 @@ MODULE moduleMeshInputVTU CASE(2) IF (types(n) == 3) THEN - e = e+1 + e = e + 1 ALLOCATE(p(1:2)) p(1) = connectivity(offsets(n) - 1) p(2) = connectivity(offsets(n)) @@ -359,7 +359,7 @@ MODULE moduleMeshInputVTU END IF CASE(1) - IF (types(n) == 3) THEN + IF (types(n) == 1) THEN e = e + 1 ALLOCATE(p(1:1)) p(1) = connectivity(offsets(n)) From 2c559135013a6a6c06418f57ea2fc089a2eba9e9 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Tue, 7 Feb 2023 16:02:36 +0100 Subject: [PATCH 15/15] First version of vtu file format After some testing and making things a bit better and more general, I am quite happy with the implementation of vtu and it seems that it is working (at least as good as Gmsh2). There are some procedures that might be useful for other XML-like formats that might be moved in the future to the common module (I am thinking right now in the implementation of a general format like XDMF3). --- doc/user-manual/fpakc_UserManual.pdf | Bin 181522 -> 181687 bytes doc/user-manual/fpakc_UserManual.tex | 5 +-- src/makefile | 1 + .../inout/gmsh2/moduleMeshOutputGmsh2.f90 | 20 +++++------ src/modules/mesh/inout/makefile | 5 ++- .../mesh/inout/moduleMeshInoutCommon.f90 | 27 +++++++++++++++ .../mesh/inout/vtu/moduleMeshOutputVTU.f90 | 32 ++++-------------- 7 files changed, 51 insertions(+), 39 deletions(-) create mode 100644 src/modules/mesh/inout/moduleMeshInoutCommon.f90 diff --git a/doc/user-manual/fpakc_UserManual.pdf b/doc/user-manual/fpakc_UserManual.pdf index 154eb5b6131721f0d957daba28b534627dd05e2b..da9471c2e3131c41080d8a832dc6ef1febb6fe46 100644 GIT binary patch delta 16188 zcmaiZcRbba`@fl0$jZpbF6-=JW@fKsWlP8^D+OOJj1BKR4A>@{zLakjO^f>S-$pSp(C#Ux9Md>0Ymk-h4?n5$AcELwO4_ zVm{H7DV=`h;OF(T9oIX4J^3Yj*Nge6;dF+axj#-|SgvblVQFEtIPOK|$aQu>k-l7Z~G_LEQ4E+3>HN zjtk+qt?%}$!EfKDt}dkFjdNLV(xWc7C$j-H)3zZ=YO~r3Y3t?FQs_3; z!(FNuPpUHxcP$T8K77MT(GLX|OTb~mouOF>x+-rH~OZrr|0%>7a6PTV$bn=fn5 zqAi#5`CAQV=j2-tel<_8_MUpIJI(!`uirKp_3q-M!NeV#1HU?a`(haVw7nFZ5^7o{ zV|7WN)oppBCTJ4(TxRMd^>eM&11BHj!d7)|{FHXmmFoxYXxgMu!NXn0x6jtEuIbzr z6Wm{iXExNROs>aCVOimCAr^=F8Z~aMZBCK*?Q1P<>C{!zpqh-es>b%8oeE8g4|g}4 z0>a;(KgChHzWWUIwc7QI$RsD-bY3}hGOzXd;q&j}O%gw?nrG+6zpGK;ER)Jxm~j?L zzAm<$&pc<0i9<9wR#ou0+ooBQ#`JH7uVlHJ`Z*zg?s*Cw@-$UU-r7H;QjR$)v8dC_ zQ$~YDh`rCg%!<63mzc2ezBC%$6&c*+yUfQq@y>D_nmc}8PE&Ff+4Ch;Q@Zq?w=YT^r(Knj za7XuXUfGA0_Tj*2=T)bQt29YZXVsU_*r~Q&zadfC-{J2K>6t0we7(z2dE$0W#hFo~ z@*6p82vxDi7JGXAJHgD~jV7O^H{Ce7KXGNR*Qo>26iF9oeWk=Z+b1sA;U3$qy@C`{ zuQU5qyP~%golbm_h{#9^sg1cA(!_=nv;Gor;5OOGe`<0W?G^OcW3y-tr7r!sl)co} zvH0uCz47X+$oD?I^v1saNh&mM{xnHDR%*zhmAXlWH0_T)`0ri@HCnKgyOJUfT;-va zS%nZQK+{(4T43iB}n`A>%y6%MrvJy>vFH-UJh1P z(}!G?{^DE~TpMqH>s~aMipY7&3o_Q|K{-(=hX{TdWzDaRuXui-iXm$Q=NcS^5-B@{ zhHY45-kgX*A~r8y$!m8RF;4Z_T9>n_I>D1AWE_UH7>nn6s?s#oy2_YKQAt|1DfaTp zL)Y7@tn(dWm_R&f~OIjX+Ac__?; z^zU}ys5vzfM#1029yAz{-_%MMc(a-q_k%XCKPZw$_ru6)u|(O6cfy>51NCngD178! z@}pcsZ}Y!tondM2zkdz4{m~|R zg!nnr0%iHb%~PK6rwpISE`|i3@O~R;D7)41EpKhU!#V_^z2~oEMk2ZOy2CEBo~TvK zq@q94g9lG(o!QpX z6E3`8*9%$B@QqHt^Rb*QYAjUO!0lJPm<2>zbB|?Mpl>}!XHLfN>gi6d5EjRLsBua& zWQ2!)j>_isO5#||etP|*(J)lu{rB!G8n~Hr{1n!a4veyjX=LwyrPDPLCu`iqWw*`O zef?n9Q(8l__vp(_3ZY=BN3B#R26vvxx!h$9+9Xp>s57)u+6iac>de2Wbfh^n`!$W1 zWq}CaJ=mJuQ20K=+F!o)+sthnw`;EH`mD@)p|-O%lDz^#My$M17q3#D=~ixX3N~79 zx!=XuZ^-S}MY}&3@{#8at5%3sI}z^c=uaY!Vq$Aa$0YVtw{sBrbls2B5uNkd)U1!Q zmueHR;<^s0QnuAx7m=*mG)K2#=D~w7(TYa4{zGdsXz(rVFebM7X#ca^ylUSm9WKzV z@{&q=MJKR(&_faxjd=RIU&M+fm z@a)=7mf~Y>pLjtfv&&RMd6v;9?Du9HPGRbLS=so-jqAJzHvFNQ9 zKVy6xuzb=k`}WC?)p*LJl6hvCiahJpZ8+8CtW9)glIznmlW?jK{l4%!Bx%GDO zymn=f-bZHnyns;ddCU@4O_7x=7q{&elG~@o7teV{q#T_ej zO{*f1f|K{(@33yp?EPF(x4bSnlaD=7DzmGdZBK2_5mHJ!wr##WDEZvy!x;vXj_ssy zD#p$`1JlkkG6FP@hQ&SvX9OHbElcY>$@%h=KI=lT9co%}Zr`hN?|n41M{EeUX5rj4 z%^GJ+Vf%Ld*V?m|tI>Yemla=B%Hl#8#rwZyE`5ofF}JHoXVMcSW4L$ISHXdVIgpl{ zd6g_`B$A??0-f_?O_P?$^YK|_7~I-XF_|=noWZtiY#}>;DKfrv4Qk51o=Zi8sP*AK zQC-sFR79k@qwOv!E#7!vke+`WH-1R@K#Ej8C~5kUH@%{o9lL!-^po574(+O>Z3ZF=IC z)|0x9GjeeqA!jKrJM-G%5^)or-+Ghdy4#(v9)!|eiN6+_oCEuHTPo&OAj^&h>WM+Z zsV`C64qK6t{G#RQpjT@qLG@x=w z9Em8wR=pG`x}3(X-#e^-@3O!dFKQC6tA*+XeT|wlu0fvO93Lsweo)U{U6l^Gr^b0x zrBgsO&~8y8;jVj@&1d>qjk%sD-=Ik+L*dI=WBD8Mw!ROedT+#JI4qJBe)kUxAVYHc zL{GhH%6e%_1N~&70X1X|&*{gVJlzsTpG<7hU1TX5^{hOfIN*~SqMNDbqNT1NLyhNFYA%^CJORuMaXwU47r){4F2X{0Jm z_7z5c*&&AY!2-i4>H3?^aq1@pPO*-Y`fXm;x1fKI|Kz4mJr~a8lFR4#DwiQMGbhly zt5#fW@GCv5-=zM+D(B?ci}~USAwlu47qblwI*fCF{p_es$hEPVt{>Vq_$69fN-JcO z5YBcNMfxGMzs=3Qw=33(_nKjC(x}?pwnN<0Wh&Df=Pgm^5?)Y<(B8*=vj0)3v*{m{ zy>PyRBZ*yWsXxXe<&r&3-q^)e2;{6%%Hzr&8Xl1wWTKosPNB=LYQCbln_7dH48z>9 zqT1HmV=i9>o{F;yC8pkgem(Ep`dL-l%)DPBc)#GSPsKBOcVwCJ8U*6Qbng`>+xlSi zV&4dx+A>jGZ}Ql&Mh%1t;N%1vzfhiy*Cstp%_6U*I*=te>=Yo|uN|1*sxEPcV{z>Y zM=M1kEanq$utu|i(K}vZ>jPTa*^1P(kJ|S))MC49cqVB1F1VQ%KT2TRR4MPHrX}Xm zieMdwVV*vdnJn4w*NT{+DErn?)IjIkU#S=4>Ye%0fZlZ?=VhK%SUnC$m1`SqGn-k4 zy)8@E{i)kmc`5OgR-wvCj*Ocx1T>e53#MYq@8syOXU5Qbx^P#X{^?hF%5;iCi|qV* z2y(weHuHIZ&vGWe2E(NVj_JG7zBg4gi7H=Lo9*kFn{*l?QMSY>+>2BQ<>5!&q#~Qu z5f9FPS5rIjk|j5FRgYT`w`O}#p4_D_;3shNi=BL$_OO;grF^X^s@8LPl~(=c{FhcEZlF;-nr$_XY|ahKvrA zcLYBe%XG7v<+*S2GJTbqU651X0iv;}vbz^J{B@&WNV%Si4r-a+J8c|2i-9;(8ek zjdsHowbD7h5_1GZy8@T7Nb0#J<9p`7aD`9Gy*DMoTPQt|&kc5}H_oZWxUI0OJfly@ zs-=xN;m*l$k?y>@)u{rdh29TNC)ux$SL%o>)5not4waBY+i@sdkv(|FFwd8NIIaKY z5I*+dhsIPjtgl+f`$Vr*gxsRAPGf*eZwk?w{pAwcG=>C*Ca)IUtMg<5_d-v1To-TT z$Hqnp$Qmj*WYTg)*Gq?=(9_>0X=9Z4)EO7TKYsb+Fqr2-?Mj6o50>G)%qRD~#;h8i zsnpqX@}AnpCN&JCXyxnYGAv{^F6gi4LQFRx5+WJQO^?p(u)obyHVm1Q^>v3?t?4_# z59-!fA2gCFO=+9G$MG=;_*E>m26XKu*ySn4F8L`mN(WRCMdcb$pWaF(=Q7`@ewpj> zi~h@0#L=hb`E>2(pKMYFyt&X0FzNhxZT^pyCK8g8QYr4G9|MM5@$Mt}t_6MK zLbbQZc~E2HbQpt}%QyQ^=NKD^Ja;E6$wR(wzm`JFOomW9QpY{HWGbH3CwXNx>N1Cr zeD$W|D;5)*nW^>_W5I%g%Q>$EzBR-kgI|C9`hfQBgk=9|eWSMl?`cdJewB5}KIt7z z^y;juIzUkp336H*n{c?+Lq3h=eP@@CvkBelPV;WO5nuU?6VfT^og?U7w zf05926Ax-*tCv!zkHa*~%$HhD}? zOlzFF4so-zZhO@e6H0|jL`l{&9r0%m1sW{lwT7~iRo^Aq-xhE5}B4* zH`yNMXIlF@CEV9hEwOWyxCV1_z(KK(e=Sl&R&qk!uQ^2^K87+rY@@l7#Axks$Euw* z-mC%XPH##pbjQ!u?9$e~!c^eN*;=q|r5ku6-_o)(oUsqBGcC=7?G&`UivQol7*D6ZjaK~%j+_~ewLJ%!NyRv|s9Qqt%nQi+c~ ziZ7al_%A+TqkOU%E=}(MVH=Y}@49EPmgiv>u6^+<6~!-N6`n%&k9ejO=`>XOb3dUe zN=NtgCY9gg)&`QRtlxk6yv>ql>OF>P?SD<4e_k#?(Y3UfD$EusNoLZoHJM1PFSu1e z<`&xK>N_V8;8PX^t-AjzKa_s!!u%aJPhn3F=>C|`mWlzI^`7u-eNb<^lBU(8({EyL zJSLe*6bPd{HC0$XUF2I(tgSm3RIV=`&o3$c+)eW3XPo?i#o%YLM|?v*^XZpr`s@zk3R_dxGG| zjLMl0&dTAW-{%t5h@2rmET=h>@gr>^SkMcS@chcPV9Yd?1g+Nwz=^KMeV#XbDa#Ef!J zvQoNh8%5TTmsRa=VPD?^5V^3g2rEXb`jT&$K~|k_F_X$iY;a9y_YRDqv5}BN0{PT zPpMOCVlExa?$*nFawdkhqOZZ*=fC~P$Uj@eVwuUANJPJ0matcB7H&DcP&Z4NL!{U%fwa-vz#>D$bI;Dz_}NB2ut9eBz-;gWm9%q3m#0 zYuq^^fpfU+r`;)-9hLSjSI8$$QPZB8EZM6{F;uNX$!p6JD{YFG#=aPYhsOMLA%aKO z8-5e2S|0f*&*EFrS%O(k@YJ-Ulvg_Q`Oc$2^$<(Vi!$vuXbPP(7S`w-t(6%(q^2To zU}cbFdZE}x)gjAKw;$5q7ldEEt*P2Lh!sZRPC-p^Vox+OmN@Qs+l{ng3m-f+w)*i@ zztXp@Wmxz{-~*$2K;*X~$<6eWf|B?1wQd{YDj0oERVUCis&2J+d`>M2;G6Ul z(9_i?x`_LAM5V?N($mIzaj3uWNtXz6jaD!1oON-Sv z+rYl-Wj=gkv)${u6I!qFPxJfN@ zatx4#pKGe+7b{aWeJE5LcZgigGryCVArqUF1`R>t@(KkNuip^;`8LktZV4R)&FZ=n z@mD7kiEk1Sau(**^Jcs3_uh3^)~JSwmAdNqH0xxN9JQ^i9gJ&O$!|)t>}0xl_Ryfi zHns5bS+`~`9mX8Kt%LKK57}>+798+|rZ#o}z zGjwgp!kJ&v4h~7`ju6f;=aoR=H7!iPr)dW&WbbE*G%to_eCM>ce)z*Jhy^3I^5NCN zk;**3(k~%?|04>9Ta;BTT4O2G3C40FVOIij-gQqzQp0n#J46(O+_{o(KhiSM`=N-i zwD>8(GqaFx=X%Y<)7{S6g&O$H&E^I*429u=Q1Kl9ih@HhRh45AvlJL{Bo>8&;^nJ| zVbm}zLL3H1!jZ57|0-gx6N&~ZI!4C$kRf6R3J4E`2Z6!jZ;cW&;seKssq5&5iPtdy zJwr@GQowF6;Vj`Q;U?kcYA12i&fU@NnuNQAhn=^bE54P6#2JS`!DMCsBa~03*A+(t zCL{LJ+??WERWp4m!Tcf_wP<0zApX$E$ETxN$qbU8_V6lHQH`BKIIUd_K_Iwa(*_TU z3?r;AUhEL=5@U@ai&5iWU`^8NUMAzL%51XXlFMfoO2T?v*A^>S$-HCPyfJq$tHJO_ zn@lZnI1m>Pk#OrC;r)QsH-67AEkS%%iVfW`Ms6nXe&w(=n61=9ukUP!T;yDqma1yv zfGD-D@HQ-wG0C8^^A;COLQc+7C;fXAb3i51aVvVU-q?bLp;$e!ptGn0Z}65y*{sqm zMd)f1OznyOk@b(Q%*;%?3di(5Cbr08X=yubZxN2)Bu)DA!vwK7eCg%Q3nmWE56n=_ z4R#~niIg9}Gz13X!-VHQvPsL+>ug6wdj*-im9Kcdz4bA&@WsPELm%G*`VJKv#)>eX zbYsmX4QIYoh9f1S)ADQlD#i>l%>`B&n+=*M-)WZD=|^p=SJ@Z&1gP!pPEMq ze@%A4W!Ynj&6n!eys-i`>Mc+F&LhVZ1_~BFOFWEVd^HnFTuwuy$4j&>$i8FmEq;f$oWTb_+5`Z8ZN+qa++cAPVvw1FYt zN~YjApIdywe@ZR?sLIoI3v^ zlf34dnm_V}`8%qXQXRQi{T<%PcP?EtILG^=E|s#=v2%4a_tzKG@3v*2`IF8@ zS?w1WNmO`WFt4K>wxVmIR&O?bxRL8>#KP~+AuQTCz*JDy!J(5KfQUfR`3J0ivO4sZ zDwYVX{CN?r_#{%1{T0#JgGBeN{LmQV-tXZLzlHwT<6H}Az9|~3DCSX;6N7LEuw%jwkNI9pYc@&8DRjU+$Sc9h#_(A_;aL340VuBK2cPv3jRbDNXEf z<>ejcMs=c1o0NAIFBV?jxe&e;JHzy_xzLqAW_{F~JsauRWMSRv-iFJ)7o|3>wN1%&JT}n$`E-_e3 z*mFu-vl8vMZ@Wu~$F}L&75iH=H*Xz8(>&tjz*t`^=6%i(N8-Zecm%&K?nJ|!YEqKg z{+ZNy5VG*$QKg8@?G38MSi|*=v8p)7tF8*TUeftcp7weSZ<74UrK8fHg)JseiyC8- zuGEZ&KHk!}$hR^?W;57i#AVex-V;4}O`H2oAn{I&m`k(bu64E)?)=NE6HYd8#)iz% z`cY1LS~Zy;wclUslG9$y*c4B)HWRK+>T+UeC7rc#w(S;)-joph9;lvY*Jo==ySj2E z9{0(r(qNn6nRhpXXL%mmjcl8MyL4&!-?)gW7hwooSMZfu;9_acaCoa!`*k@X>!C2(u%_#(DG4ZCpO*eP-`ukai9g1-==~znDZya`sd*3~v z^*NZmw!Qx#MD-_bO!ehq^dkXXF~G zOM8ULirnF?tQRe1SyhOM_D>hzDT!~cX$ii@SR zno(P!8y6ZMp;K$JceoM&of)4?ri@Zcte9t@9w0HlT$u1CF@Ih);f=cD&(8d=#s^M= zI8pZp{N@odDPQkzas;AZON7?*g$%yyc$i1zTYYYSk5$Z-`d-G+RmMEkVlo#IJMrt^ zySE!mKC9arehTWET4ArPtR%S*I38m9{G7k7CXy*&R@$Ae1ya$TH?p#Syhl&>}U5M!#KaxTL&;Y&6J?`QZgnHF_cpDm}0>^c8Tbud9Kyh1^l z7MD>tld#D+D!}heJl7t(Tf}h%9ys0W3Ks1e@C`{H#XSI*jHyV=r!godr z>&N{hLYnMxie_rHXp2wrQnSCZwT=1^M(IYCMg8bpxOBqW%y}Q-slKSP&cg!IGcC3L zlowvFQ#4As$3*dRCgUr*`s`)TlJ7si?1aO&a8(F8liquLt5zfbqCV}6^KGtd{M2B{ zb6c0ikIx+rWKP?JrXOw%y9?(c$8x{ognvHUSbwsY{xUBjD`0%U=h;=K+x1m4(|J9b zMt#j2dr~L1c4sfWEBj&NC~|l8eS^dYogJ^11ObUMs&8?t7i6tmr3T+-z0JG=YkvHq zPv(2)P%d4}mv22{(`_5`!#0V1$semTp1uquRUOi`Po1Asc_nM>DJa$cc4Bcqv~X_B ztGd&RMgEZS^gK6x9WM(>R}?kc0E(1EVkD6W{8lyz^v~~NIV2r|=xc~;HrL@cV%M%i z5n@p2b&MDU_=rJ~P+L1YJ2=YL79k{okIo?xtdq$laU!8WLy<@fa9vp*3CDlG)D7j4 z2opPEArJ(N=a?VDBw#uIGYSo2K|oL-8WIik*MEtC1=?)zIY0)&b5JN0Pzm=pGZc;h zA4H%*@g90fBgTek)4)@!@5G)Jw@nzwp7LkQmS-2s9j| zF$w_+JuCu)0=Wr+1wV|0{;Aw=|3k1y7z*?-U_p>}NF-Q0BpL#W86dep7DQq(U`23H z80Zlc3R(Qqh8Ks53nUd*572fh#p17t503u4AV5fD78 zhy;U!LkO>iA|ar#g~EV~33P|TQCJwk%}^*7jv#x3>;)=3=)L_85IL3L`gs@1ccBb5GWS$n-~A8 zn1lz1fD)!46a~YAu>^x4%r+oSkc3h|kthTTEEy6@C(a1lfK2{iD#$Z4+CvXwLp6{)YyIqw!PafQ^2T6eZo!Wnb8DL-tVKIavpfC(rd=$vB5F`q0 zU!XpM(jAHUGsAw{2rwLyup9KxZX zgasc8K>^k2Z$^LyCv1U_X&~dmpr8r>g=3L_wmtv25@@52eg2mrA^*!WC{Q;EJd6U; z`1qdx%c)};7+FwgFtVV4LE&HC{~8jUA`qbcfZ+G5fbtDY8~<|r79WZRAen#$g8-F| zfCd4&R01CW4GB0u|K`7AU=!4=6Ac8XG zcq$+i=D43Bcm(i2FjwJu8%TKZ{q-awIN*O>1S$Zv85EX49s3`*BM8G15El$-0PRnR zJ`M)p8~|y@?^^xgTp*_)g!p=#X9OV(Far|=NNB)8f71X}Q9zz)Hv*uJ#2+;PpbmwC zfEn00fbs#Lfhr*syoLa_6$rwn2zV3&s^*Fk#^;2soVMFPeyU~2nY zAt)SZxIo5$f(9YrJ~SMk&_u$B1KSe{0%Hgg=z0hs0fK^HK%MaMCQbmzZC)O4|J%Y(hP-T2sc@PCIC2tObB}9_+i490$9fY%U>`r7LNOy z2%G{C^=}#gW5B@uHx0ND_<(3|LMs5ge~+`r@DAwvf0&Blr{pL%>2K2ul+11~7?&-arC1 z6OU^pp~oT71i*QGGhy?8ys9UZ08j}9CrF`V8sRz#fP4hFk^@`v-#+}ag#uava8)>7 zcoB*Z05RBXz&e=l91;eGPN2v5llONE1>icUkpv1OILU#1E=U+8(6)je2CNRE0nHmB z4F$!4`vqXWBRmL<$)G?1)*c|J16J`M0RhK>Tn7gl(|-i~qX--?{0Z|9$PAE65Ge2l z2Uwkhq6`U441^I1c;xRi0PGL!C4q85U>YE_{%UoB+=W7c;sY{22G~{-@InDlBUA*a z@&9aA|LE#~8V?2ZKgNZBOCbyZ90GKK!k|D&B0LAIGYEx-A;1U*tTRBJ9WW-s{&OGR zuY&|_se)Aj3WGA>WJM%gNeQE%fR+a)L||)!LMvd=in7xG-xKE4va-}3p4RT3KJIqc Ssew)hg@sTH2q^2QQ2#$a8L1-x delta 15937 zcmaiacRbba`+tt?A~VU}Bb?XS$IJ{7S;@%WNn~YovR9IkBZNq0Qz2QIQT7a3B{N%+ zhTlu?&*%HD^ZGma=G#2C7=o0Mc%lY%J`hbL2Am!4HoBB}_c zJU~2>Ugz2Jr95iIR2F9?zA?dyrvkN`8B_*pfc+m&55UR&P~5qX!RXxlm8joSF$%G(LtVH%y?l>z6PO80Wb zl?v+ZaYx)W*SNCrH6Ks7@H^b4ANmzq?u|{Z(Xoz|gjd|jci^_E!j}`uy8A~wpBaB! znYsS+$TI^DjOqA_h1&<`iT4-EWVtIft!>LRLw!H_u2}?H?lt9|fYm+fwd|nQtyoVi*EBUedT&BF7H^bn9_si=W#Y5X42kx(juD4y8T9Jz?4n9FI zd&ztScgz3T@!o9O8^R^m#mf$lkE}CRCY(#wSnw^A)(S4`lYT4Vdp@DoC+2>7@X@s} zwpMDd6*?y)jvw2#F?@{O0x_KJwc&Xs#-W6$*`|WecgsS1$J9mH85!iKj7Qy){rb(&rmyph z+X{J12d)g={3xVi_34Waw`A=XR_*>!p&`}6DVwuaGIs?mCr8c?VLA*~7mKp5zhZjh z^QZ}D!Uqf_XqwzVaqQ^Cuseehx*&?#z9pJD-xx*NxY?b( zOq)IR_E5>dz2Tc0dMeH>wn`UL>j;Me1n>)cS+dj~af3Q0-hpwp?8;6sS-4}Fvw8t_ z^_inh?sL`-t9L^kMK6%`J(>LFEQ@bK0W)TX#JDP~g6!_D8wG2?EwhI-wy0B9?M=ML zW*XO=Mn9|$2^F6FLcVdjn$FGhc;E6M@z-+|L;9w~Ln46@8Z`C$bUXp%QbfYEN_eaL zK=a!tX0FFo4!vyV%8B9*iV$Vj9tha4UGX^F#LTNUX{yB3!KXh+OoK6gqoeS9ZA))@ zi`*$QhYH!Hqh@(zSyue<{3>&Rz%_}YrYBMSoXaIcno{4YqXlz#pR0K|o3}JNgoR55 z_WQxcmqL5qDe?H%a<-lzcY%iy*llB9@>toSZ)%IU889f`*@?`{j$Ufc4WN%3q%7{9$blz^wLGJ$Q60~Bq4}JS!k_E z;JOga&iLMKJ*5EHPh)?gs_&2hYA!;Ls1>cj_O$s1vd)muEv%(MJ51+nlQbGn?pDL! ze>blFn$xxbYrkr=tN2XVb}ZS!3Fq~>{M;CoA3iF4=#tC{+}6M`$lmh0 zztuOpZTa!mv@73LA2JGO{|^NGXR%aWF21Sl6^VCDcfLPS->n&N(m(yuxV~R+n8#-cFr`@@mkj8rElwWoy&(tT$IzM*0d**yhiE z{fxOH@_1hTmEIX6!!@;GyO7;jPe&%v)QFGydq3GISRYo>d7ZnOW0#@hMTmPEWoqet zpZ$x9qo_c$t!Qyw%RLs>(uL?xN6oS<{Zg008fG)Rn(JR>YO*JrVg|{q(&FxqiQ*Vx z7v_$Pk3?UM9hWwI{Fa-U$^Xh%y~3rWA|o3mqnyD$dTaLCuIwSL_uOaJ#){O=&)O~C zF;P;P{*tZItt)77_0+Y*h&yrw)%huPHKLNh9WJpDA=!wH))pH{rPXoNA*8-}zs-YM z2O8JI?_w9W?LyN}5-sZ%-YQ3o%hf$8dp-Hw^0DEJ!ON58<^yYik!kbMD8m}}Bt^B1&N=^fvxe?A?jLZ<7X$;rn&T(e zFq9P%5#{3zDzez0#2-k3$Aqs+~qV>!;#<`ux?H%xW3N z@x=|UE{8Ox&mm;w1Ye#V?)gq|ig{X=orgTefUEv_gK8$p@X`7HkMfR#1-{9296Fhj zH0>sYOcW^dpB-Ol*{CobDyE8EhVtWS6>ELB24w)8O&8aaFc~ zce2!gt{=7S$0lUUVZ_HGuQ}NzFT2nboE`O8lY*t|$yUp&BhoE-)sC~w1gsA2x`gVNkJ6~sQ4YdrVtjDZ5hI(8I>CCv49qaBd<0H;G zc~vb(6?HR0IjNANs#-tnMoJc!lEir8BS$ND7P$*}N2$eB?0;qXaHt%TTp(Un&JW>u z&pO{X%gp@!hPc=`yw@t5F^^mCZSs@0etsR3O zY;(sr!$nok-p-RcTqfLAd}l(4DZRyndOW5BHLWjZ8&rc~zkRzAd+g^jP0rRQWUc>Q z)OT2`Qqy8vqgvvlG};zP_ocFUyRR3mlE|5wA2zJdrx0Sdheg?zjz1xcI?W=T$x$6* zD(BNKzDPugbu-Y|Ip>$=E5|M!o~bTAoHsD=GK?brv`meTXP9AN_~&BlOuptK8eDF> zdd+ly3)?G%K2bxq`dm5c@Xna#Wcw#Igpzv6iS2Hefv?%+INdzMinjYG&)}Z#kCYWd z(&**X-6*Cg*WW33I*;TNSS0A1_4Fg!rM^&npouKCQV2hPXiTiw9P^P_oQ_vlk6bZPk1-`U%xp)rb)HQF;i*A%`L%g&ipw+=m95v ztnND&Z!=FCf4Y=tyAS;-C!958sS)WLbUtHs-@}Ka2Lwd!{GKqRmDHfI>#T zx+6=mgcdn+62AveCzx_HvS^?8KJ)&j7YDU&*}537q;HyRsn|E1ESU?XM4h>%9+R=# zEO*`;)RkCjDevfMle#PX4GYSrpkL6tkpMHcj9d@aqr+7p-BvJ6wHJQVfWqC5N zt8D8cR^ep#GMuPdDs|#xM-qbn&NFrYhM@$Xz}A?zN94{@Pu-oku3sL<*BL)q*szb;;|)Lm~u#0)M%1H^I=h ze&R$6{bYFs%A^{3u0N_ZKbaTlae!Sc;hTmBi=y zO~?p~>EuU?sO)BFcRlKwYG(~-j+$$FwR7|I`g>u{)0`VnIr^~P_Gsrjm;+pC?@q0A zV5*{t@7xI)$G$?{rc>#)=SS=9-xS+*YHDxPZ%9hKoJ0+ne7QQqXwj!uQfeVcnfue$ zx>2l3Jbx%4tl}(zvQ5kc`dZC5$uPF)ALS=H+xA|E+NPd5)-?V1JttPs8?b&lc>l5=YNh$@Klm^N7 zrxXpiG8-Q1NIcX2Zs&MspS6rWS8j2yLdf@!`9f~2j1`-z7GZt(r}+HBS1u;CV=2|M ziL*HIly#HnOb;2qBOf-`PpE2ky`|7?>5*5=vbYVK-xB!XI9BjFCit_)!|jZ>GiSYr zKRC*+&_`#Tz5S-(?s}D&SS-!TG|HEN~@pIOHNuZ#=vtRtFIu9eq+Q$giN*JF^99S|F}){E zxWnCVR~=8?ii!JvaeyO$#amJr=rz!>oB#GD>e6ZHJH2@i(+TbT<+lxC%bFU-_OC=3 zlay3X4GA|gxfe#d4{UzAe3h6Y{WzBA`HiH(-dU5b0mr%E5Q}A*RT2H+EL(Zw)@Ro1vGTZ9kt}cy(i0R?=a%=P7kS3WI-4iN z!=?IpaCGFX2!91EyC;8gM^5i~O;hF<)qI$AiL_`)Nsr~+J%!->d=J$YjSs&13FWKl z>n(Q6=|{ANbyJMdWS<2tzkS>rx{BO(DR>)o}ugCERszz3m#D9+jtY(~s<1 zmz{`ERZhxx!`qEShO%G$T+K}`u_-k96)pp4(LT=)$-0^PGrE*liD6}5ROk&=&tc+T zT{V51q%Fu@Kv1c#6M0hKc2rQ^PlY?mqZnnJ&~~XrAm}VpvB;K{>lg2^7q?gkaz)g7 z8RX}Z-zR2TelI8Jm#TCz1kWh;?#61~y)Q)+ZoN|IAZC2p!Ybv_ynFk7=5TH~<5(5B z`S&*@FY?&y&&<5f%Jil>o~em(@afTO=PoM}C}B4li5`nP${RIgA|qXJUsomhSoN#V z*8IYrno$>d3)Ih_>Nw=o84>+hziZ1P^Ssi>CavmYiwnD-ZxCEKA5+NtD?P3-KUA`j zY;z$y%t-UisoHIejp*`B$qR7#bAbFQfPFpE$X3#tg`FqGUfT3@C}zxh?9dt|W4O6X zOnw;QDqzHj7&2el4GBsSv6mg^W_XzJc1SeS2u4Y1}D!xisbAOH46W6;>X1y z+Y5C$nF6;(1FpiSydEO90X)OB-}cuW9*w z*p=eOrp0)_)l7X6RIi5Hcwc=k5(O_viWp0JEwi<3PToGb)pcK-Y0Kum>XMxGph}v> zGW`!>gMcN!krkTWkYO4I=8BPl-dX&&&htgigEy}%ZJe68NO%*t6du8GH2c=kTZ2D1 z3m))0dP}7)^KrZM>7ye}yT#r#Z(-0riXyIkmL4v^fM&_qv)H!B=7ML++6oPuY0Iu> z1fC=_6KYPQ|CQ-|;q<8@?&yoLz)f=Iq_^LPoXGgOx0XH(J#lXz!F2*LV`GV41-?#k zoRcnL6B~~<=(h;{+^#Pd7$1yD6Of!_0lsmlWU5|<*}~gt5_TW0#+b7V+Z+weSQCsX zH4Ugdi~Do#e6S0g{&+X(YfW1Sx+v3<_qaWLL-5Bbc-oW}6wCA`-%sizx27}B=zFX*{iR0j`X>cQa<(rAIKup?wVMAd7_vnvf^j7((`3b&spAM_MAiIyM;NuQzg&G z=N|P=J~WP(;A#y^()3moW|MuJ=Yme9lXJOVtI0cbU(m2AS<8m2;Pw$=^2ly}`>mCy zSu2Mxy>poj5l*~J5ZS0QP%B})>~Qaih$d~nqVeg(Y~C_Br5^H{%LD7RTMbrukLN{x zdhC&@#~X2v!VxWM85Miq7QV%s{A3w_8ZADl>W!V38T2yP5-7+&dR^dO=D~4maB9&0eEww;{yR;`I69Z?4I${#GcDy zLBg#^Qf{>7T6$1-as8~@G2vmL7MJT!c3~;J$*djG)K>IlVz*!uy~bMe|#k6-eoI8 zM(Yv1op}7pq;GQg#ZOt}0#{83KYkkgIq~3P)Gu9)y&t9bW{x=!G#*KC3s~PcPh0py zxgoxknbkq?z6N%W!ZRh-H%o`2ZyPJ1(6A@VDfZ3T{)*i-FV7S9)^4=mukJQ3v`vLgL^!G!g!s3_+wVBIBqEd``B`MQL!(NLf>d=-E%^M2^DY ziGAZ_d{w{t$v)%$`_Vcb?C}W)XCHfMN9oJbSESvg-LKe7U$ysicDIxEB$gk6T_&KB zSOo=8iUT41GQ5og44EC-viJm%M$Z*wiYt96tE71hDGG**3t5Xy4jT>@MDZ&7X$l`x zrgowwqY)&Y<_+PWjSjoXb~wtj$b=J zc4i>IcKZ?DqPRv@bJ$$uaZ?o4$?63y!tNl7n zfz{h8akh;6g8q(2`1mvxsszpNS!BP_ID3KWoMt8myPW&v`|MYDE@F7w?@c$L`Lao~(q2?e`xT=1Hb( zCO@|CqE7=PaDvIaDa;a4MlBL8EFZXRB!TQD z%f_uJ*o>~>v8g194(1z$Kj^gxA%c|qjK5?SO%@#%Y3`ZrppCRFvg`9cZT}K^AtDnV z)2Gbf7p_Uq#26%6t;=7e`sv8~aN{qR>SgPN$sYy(_|hLZww-%KS3G6Yu)#ZkGcr9t z%qGh2dO>R-t!X3O#n&Y^7k5bXeiLIREqt)l*%DML8?n#WxC*EH+ogWc|y7&45GwCm z|2(bHuq!3ST{84?n(uL}3}+2N&TX_feQq6DVP8a}X}3GNf8M+B8|ry5RmshO7m0!5 zohDNqaoijXm%Pm7#x+rkGNvnu)IU`S*4)?qB0oP94-)_SBwVYwH}cuH+M>}wv(7%U zs#mfZDoS+@??*fMwj`G>Sj?ulMFyls z7ankNNc)dSXm2DwUZr(6Cm!(-y}ejcX>M0%c-w<-AEm zmRQ%lcp{I=!~(DB_1%GVMYWkKQh z7{{?teH6^NrMz*9<@%?W-#^?SbXVG2K4Nu}7({vIWNg#VRA)|o3R7TdN}ZaC&7KoH zb>mv*%z`Fckj_kpf)!OZe^~Q?XjH`W%FekAay6dQ%%8)*)?V_pJ${Hwo3G^DQXapS zYXDwQFE`x6|M7YT@i~o4^s%rjNLh28NGap~OEK6bUOH6SQ4Bp%n+s*iLeM)Rq#pWU zLC9Q@?0YV~N+``MM0pGJ-yC>XD>hR&r~z; zRTw1*lb3eU5p*A^PESuuV|dtz=?+>?Ev8!d+S3tQZAa@~WP&@+*1OeUbbo(VqX&>R z~=djcD+2+%!I4UY#d{&v@UZzJoAmMqpdrR^=juvnga{yFS|z#vEky)Z^dr zwPScI+w}xPs|l|AeZEo!V`KicFNg2@YusRTQkCaaJ;JfS6m;q$SE<&U=}uSnvf%3a z^6IU!s^bflS+VO&d3z7vqofoi@UHSj!P~c$duq--q2$S~Sr?9UJ1x4CxKqevE&j5@ zPl~Jfd9%Cr6|xGAw*{}CeO-3-0GblYL2Bt$tblenfUe|_%s+H1)i ze}Pv*+Fp0K8+G|_v0k$F84Nv9u(-FP*Z$$`QyWfk?+%w6Wou_1Qj$79On8tKGubdRI)U|#mzI>vT zXBl|;GDp3z>%-GU(}9skmq#DjvaW7hRiPcSFnvcTR=SfjS+}aiG z(2DKEfbs=R)MjHxVL)Y;*<^Q#>qUlVe3X?c!4{SmY^X#~L#QEiFT=JrK~&CbeEM`! z4A;ckyvTAtxz$R@Sr4xgC=XTj7L-E{d30wD28My+*9%W*as{~tH&#Zc` z!#cgDnk=7UN;Y=?Y+|?K=NA{A_hWGM#rr^*Y9NYRbwatEzQv;0_YQ-+)3P^>zt(Ot z{^2q_dY5<65zA2!mx9P8Fupna!l?SBxaYK6l&@3qc;PK8*h`i-m7cpeqeJ$xl6UXQ zF25TlKa(Ectmzc}j17}#*cP3Zq*8?7la^GwD@u8WTWt5eO%ZVM=9SICx{(4mL)O{k zhngW%x7f8 z)a6}4f!345V@G6~pJsfDK3TmOb(AelBdj%F! zEvrs~ZBFIzsqpYt!1VL+Je{vPm+s@P2h>KkN=9rxJ

!!;UqadsaQVY+!jf?eTkl z_s3=nXE5T$Iiks3AMClect118uQeh&9obk-`EILQ7Dz0N1ZXS?1*{Psbe?45qwR9s zv||m=EBsIQGAPgDn46xaFL2%~0OG^04Tcj6E(tbMyAX$&TCWSk^X;0Pof z0sqg5DxW-W9_ad_~dzdwfK;b?&53P^(^z()Mw=f5I=G=Qk|7{*LMBJrdjA#f;C>F{9xjvze; z$6-k{gW~}ZP||Y%0*QgpP@piNK!iVK0w z4&)Dq;}IDAZ)yG%1_j#kkAuI@A+b=jkZ=s-E)u{)41+`>pyyBo3`9F525K!N9%>R4 z0M!nIL_?UMQ20MJ{;dcG_>Z#@{%{Et1_7}X3j048Qg==JRkA#E-FPEaV9c8VT`P3>Ff3cr2V~4L+-ggdWFY zNo@uOB^(MJ_8`W_3t%_`0*fWl4~$_j_57_f0N}BtmI6(TMUxl>JO`jifCZ399EyYo z{NTBRcoG6J{=dM-p%HjUkm2B9f~+cp;Ry&hNn`;C6q;00IGEEQ(BjZ26a;P@K$`r( zy8wpN1aJV4`i%zBy$HrjK;lWm49o&ps59ZfpeMl)40%#105k?ciYq7u5=|-v9*6p! zR}ON*|KWwlVGz)RSTqSt;Di4oDnPe@+y(PDp~>>y#`F{CK79EBn6kic_T49Ph-4lFi*7?%4VYolgnu;+r~hyu@GA_P+Q2mq4V!GJ9_vMpIA5?0vi_OOzH(#00FTx_P=fCzx5Rd zHaHL;$AZ=VkN1ef6)-7r$fN?+E+jO4f;}R!vH~VR06)e+QV|@BB?TWmhb0|kKpKwJ zO+h2!p&1v1mqdIx5(Q3(yu{v0n1CS`PdZrua2$aAQ@h{x$HJk`2qsqusyMJ6Bc7=O zeF6zJ1DLNM5s&*HX%+`I*MGnZ3IGa)fjBQ%QX$m@Y}Ec~AOG?R9QJ=O#zFHt4i9Av z_DB$+!5#^c3-L(!pP~LwO$Fl`;*5B379rIR?3*CYh{r&!iN`?eEI4IDg~75wJrAr5 zjLqQ8fPs1w9uMiEaaaV7#0>zjnZqNAw$(7a1RRn#zG4~}52ZQ5@=U>8DU^o;p;3Z5%p44DqJtH*)$OZ$cxBvi7+QbVrFd+hDX#s#A zl9nU@EWD(%CfJGooq~SL4~|`=aRY#p9EtmaLSg>_^UoLoz^RACEx{}dNk(8{14|pJ zCBOmmPv8G73h-2p88uV(poi~wL5lg0?BGAV4J)L_&>7~@HU^FTn-#103m9Hb2Z z$5e4B)Ch3yC1rvDX9!|pEts-#IPh`E$B4g_I&dWn8bhi)=!E~n)e+EX5FBW+5LZWl zUmWuuETjJPf3RzX_&Neoiouo@Ixhj>4voZ$2Lh4~U~nwh-2N*J;h*&&fBalg1|;OTZza&>$h>`ays}QYT1*q)sr3NFx@!3#r+l zBIw_3*8%?jrXM5<(x+oVR3Wf}Wf|hPNc11D{zesK0$K8b!hkEKf06lZe~<`vcmEO( zeByvc$^`5X{+{v<4j%m11`+|myO00`g#n?4@F0K-jlYQp4?>&<2WC1V+asvP+%!1g%C_PkdX>R8Zttl!9kDI{^Q#TeSyfeqBZh2cMS0sL|Yk~F|gHL1^l8?L{%vr-^52ynea zsxbnjl0Jw4u*CCaFlGs8f5>K|mL_(^