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