From 43a74217956507c81535d35d313150c911aa2285 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sun, 5 Feb 2023 19:35:49 +0100 Subject: [PATCH] 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