diff --git a/runs/1D_Cathode/inputCart.json b/runs/1D_Cathode/inputCart.json index 2d1248e..74eb5f7 100644 --- a/runs/1D_Cathode/inputCart.json +++ b/runs/1D_Cathode/inputCart.json @@ -15,8 +15,8 @@ }, "geometry": { "type": "1DCart", - "meshType": "gmsh", - "meshFile": "mesh.msh" + "meshType": "gmsh2", + "meshFile": "mesh.msh" }, "species": [ {"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e1}, diff --git a/runs/1D_Cathode/inputRadEmission.json b/runs/1D_Cathode/inputRadEmission.json index 012f105..14fb767 100644 --- a/runs/1D_Cathode/inputRadEmission.json +++ b/runs/1D_Cathode/inputRadEmission.json @@ -15,8 +15,8 @@ }, "geometry": { "type": "1DRad", - "meshType": "gmsh", - "meshFile": "mesh.msh" + "meshType": "gmsh2", + "meshFile": "mesh.msh" }, "species": [ {"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e1}, diff --git a/runs/ALPHIE_Grid/inputBaseCase.json b/runs/ALPHIE_Grid/inputBaseCase.json new file mode 100644 index 0000000..7cb1ba0 --- /dev/null +++ b/runs/ALPHIE_Grid/inputBaseCase.json @@ -0,0 +1,76 @@ +{ + "output": { + "path": "./runs/ALPHIE_Grid/", + "triggerOutput": 500, + "cpuTime": false, + "numColl": false, + "EMField": true, + "folder": "base_case" + }, + "geometry": { + "type": "2DCyl", + "meshType": "gmsh2", + "meshFile": "mesh.msh" + }, + "species": [ + {"name": "Argon+", "type": "charged", "mass": 6.633e-26, "charge": 1.0, "weight": 1.0e1}, + {"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e2} + ], + "boundary": [ + {"name": "Ionization Chanber", "physicalSurface": 1, "bTypes": [ + {"type": "transparent"}, + {"type": "transparent"} + ]}, + {"name": "Vacuum Chamber", "physicalSurface": 2, "bTypes": [ + {"type": "transparent"}, + {"type": "transparent"} + ]}, + {"name": "Exterior", "physicalSurface": 3, "bTypes": [ + {"type": "reflection"}, + {"type": "reflection"} + ]}, + {"name": "Grid Extraction", "physicalSurface": 4, "bTypes": [ + {"type": "absorption"}, + {"type": "absorption"} + ]}, + {"name": "Grid Acceleration", "physicalSurface": 5, "bTypes": [ + {"type": "absorption"}, + {"type": "absorption"} + ]}, + {"name": "Axis", "physicalSurface": 6, "bTypes": [ + {"type": "axis"}, + {"type": "axis"} + ]} + ], + "boundaryEM": [ + {"name": "Extraction Grid", "type": "dirichlet", "potential": -150.0, "physicalSurface": 4}, + {"name": "Acceleration Grid", "type": "dirichlet", "potential": -600.0, "physicalSurface": 5}, + {"name": "Ionization Chamber", "type": "dirichlet", "potential": 0.0, "physicalSurface": 1} + ], + "inject": [ + {"name": "Ionization Argon+", "species": "Argon+", "flow": 27.0e-6, "units": "A", "v": 322.0, "T": [ 500.0, 500.0, 500.0], + "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1}, + {"name": "Ionization Electron", "species": "Electron", "flow": 27.0e-6, "units": "A", "v": 87000.0, "T": [ 500.0, 500.0, 500.0], + "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1}, + {"name": "Cathode Electron", "species": "Electron", "flow": 9.0e-5, "units": "A", "v": 87000.0, "T": [2500.0, 2500.0, 2500.0], + "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [-1, 0, 0], "physicalSurface": 2} + ], + "reference": { + "density": 1.0e16, + "mass": 9.109e-31, + "temperature": 2500.0, + "radius": 1.88e-10 + }, + "case": { + "tau": [1.0e-9, 1.0e-11], + "time": 1.0e-6, + "pusher": ["2DCylCharged", "2DCylCharged"], + "WeightingScheme": "Volume", + "EMSolver": "Electrostatic" + }, + "parallel": { + "OpenMP":{ + "nThreads": 24 + } + } +} diff --git a/runs/ALPHIE_Grid/inputIonization_0.10mA.json b/runs/ALPHIE_Grid/inputIonization_0.10mA.json new file mode 100644 index 0000000..e7664bd --- /dev/null +++ b/runs/ALPHIE_Grid/inputIonization_0.10mA.json @@ -0,0 +1,73 @@ +{ + "output": { + "path": "./runs/ALPHIE_Grid/", + "triggerOutput": 500, + "cpuTime": false, + "numColl": false, + "EMField": true, + "folder": "ionization_0.10mA" + }, + "geometry": { + "type": "2DCyl", + "meshType": "gmsh2", + "meshFile": "mesh.msh" + }, + "species": [ + {"name": "Argon+", "type": "charged", "mass": 6.633e-26, "charge": 1.0, "weight": 1.0e1}, + {"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e2} + ], + "boundary": [ + {"name": "Ionization Chanber", "physicalSurface": 1, "bTypes": [ + {"type": "transparent"}, + {"type": "ionization", "neutral": {"ion": "Argon+", "mass": 6.633e-26, "density": 1.0e17, "velocity": [323, 0, 0], "temperature": 300}, + "effectiveTime": 5.0e-6,"energyThreshold": 15.76, "crossSection": "./data/collisions/IO_e-Ar.dat"} + ]}, + {"name": "Vacuum Chamber", "physicalSurface": 2, "bTypes": [ + {"type": "transparent"}, + {"type": "transparent"} + ]}, + {"name": "Exterior", "physicalSurface": 3, "bTypes": [ + {"type": "reflection"}, + {"type": "reflection"} + ]}, + {"name": "Grid Extraction", "physicalSurface": 4, "bTypes": [ + {"type": "absorption"}, + {"type": "absorption"} + ]}, + {"name": "Grid Acceleration", "physicalSurface": 5, "bTypes": [ + {"type": "absorption"}, + {"type": "absorption"} + ]}, + {"name": "Axis", "physicalSurface": 6, "bTypes": [ + {"type": "axis"}, + {"type": "axis"} + ]} + ], + "boundaryEM": [ + {"name": "Extraction Grid", "type": "dirichlet", "potential": -150.0, "physicalSurface": 4}, + {"name": "Acceleration Grid", "type": "dirichlet", "potential": -600.0, "physicalSurface": 5}, + {"name": "Ionization Chamber", "type": "dirichlet", "potential": 0.0, "physicalSurface": 1} + ], + "inject": [ + {"name": "Cathode Electron", "species": "Electron", "flow": 1.0e-4, "units": "A", "v": 87000.0, "T": [2500.0, 2500.0, 2500.0], + "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [-1, 0, 0], "physicalSurface": 2} + ], + "reference": { + "density": 1.0e16, + "mass": 9.109e-31, + "temperature": 2500.0, + "radius": 1.88e-10 + }, + "case": { + "tau": [1.0e-9, 1.0e-11], + "time": 1.0e-6, + "pusher": ["2DCylCharged", "2DCylCharged"], + "WeightingScheme": "Volume", + "EMSolver": "Electrostatic" + }, + "parallel": { + "OpenMP":{ + "nThreads": 24 + } + } +} diff --git a/runs/Argon_Expansion/CX_case.json b/runs/Argon_Expansion/CX_case.json index 21adc6a..463c8af 100644 --- a/runs/Argon_Expansion/CX_case.json +++ b/runs/Argon_Expansion/CX_case.json @@ -8,7 +8,7 @@ }, "geometry": { "type": "2DCyl", - "meshType": "gmsh", + "meshType": "gmsh2", "meshFile": "mesh.msh" }, "species": [ diff --git a/runs/Argon_Expansion/elastic_case.json b/runs/Argon_Expansion/elastic_case.json index 0a10d86..ac9fcc7 100644 --- a/runs/Argon_Expansion/elastic_case.json +++ b/runs/Argon_Expansion/elastic_case.json @@ -8,7 +8,7 @@ }, "geometry": { "type": "2DCyl", - "meshType": "gmsh", + "meshType": "gmsh2", "meshFile": "mesh.msh" }, "species": [ diff --git a/runs/Argon_Expansion/nocoll_case.json b/runs/Argon_Expansion/nocoll_case.json new file mode 100644 index 0000000..55dceb5 --- /dev/null +++ b/runs/Argon_Expansion/nocoll_case.json @@ -0,0 +1,55 @@ +{ + "output": { + "path": "./runs/Argon_Expansion/", + "triggerOutput": 10, + "cpuTime": false, + "numColl": false, + "folder": "Nocoll_case" + }, + "geometry": { + "type": "2DCyl", + "meshType": "gmsh2", + "meshFile": "mesh.msh" + }, + "species": [ + {"name": "Argon", "type": "neutral", "mass": 6.633e-26, "weight": 1.0e8, "ion": "Argon+"}, + {"name": "Argon+", "type": "charged", "mass": 6.633e-26, "weight": 1.0e8, "charge": 1.0, "neutral": "Argon"} + ], + "boundary": [ + {"name": "Injection", "physicalSurface": 1, "bTypes": [ + {"type": "transparent"}, + {"type": "transparent"} + ]}, + {"name": "Exterior", "physicalSurface": 2, "bTypes": [ + {"type": "transparent"}, + {"type": "transparent"} + ]}, + {"name": "Axis", "physicalSurface": 3, "bTypes": [ + {"type": "axis"}, + {"type": "axis"} + ]} + ], + "inject": [ + {"name": "Exhausts Ar", "species": "Argon", "flow": 0.7, "units": "sccm", "v": 300.0, "T": [300.0, 300.0, 300.0], + "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [1, 0, 0], "physicalSurface": 1}, + {"name": "Exhausts Ar+", "species": "Argon+", "flow": 0.3, "units": "sccm", "v": 40000.0, "T": [300.0, 300.0, 300.0], + "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [1, 0, 0], "physicalSurface": 1} + ], + "reference": { + "density": 1.0e19, + "mass": 6.633e-26, + "temperature": 300.0, + "radius": 1.88e-10 + }, + "case": { + "tau": [1.0e-6, 1.0e-6], + "time": 4.0e-3, + "pusher": ["2DCylNeutral", "2DCylNeutral"], + "WeightingScheme": "Volume" + }, + "parallel": { + "OpenMP":{ + "nThreads": 24 + } + } +} diff --git a/runs/cylFlow/input.json b/runs/cylFlow/input.json index 6abd1c7..b61af96 100644 --- a/runs/cylFlow/input.json +++ b/runs/cylFlow/input.json @@ -7,7 +7,7 @@ }, "geometry": { "type": "2DCyl", - "meshType": "gmsh", + "meshType": "gmsh2", "meshFile": "mesh.msh" }, "species": [ diff --git a/src/makefile b/src/makefile index b398134..75795c8 100644 --- a/src/makefile +++ b/src/makefile @@ -3,12 +3,13 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module $(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \ $(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \ $(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \ - $(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o $(OBJDIR)/moduleMath.o \ - $(OBJDIR)/moduleMesh3DCart.o $(OBJDIR)/moduleMesh3DCartRead.o \ - $(OBJDIR)/moduleMesh2DCyl.o $(OBJDIR)/moduleMesh2DCylRead.o \ - $(OBJDIR)/moduleMesh2DCart.o $(OBJDIR)/moduleMesh2DCartRead.o \ - $(OBJDIR)/moduleMesh1DCart.o $(OBJDIR)/moduleMesh1DCartRead.o \ - $(OBJDIR)/moduleMesh1DRad.o $(OBJDIR)/moduleMesh1DRadRead.o + $(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o $(OBJDIR)/moduleMath.o \ + $(OBJDIR)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \ + $(OBJDIR)/moduleMesh3DCart.o \ + $(OBJDIR)/moduleMesh2DCyl.o \ + $(OBJDIR)/moduleMesh2DCart.o \ + $(OBJDIR)/moduleMesh1DRad.o \ + $(OBJDIR)/moduleMesh1DCart.o all: $(OUTPUT) diff --git a/src/modules/mesh/1DCart/makefile b/src/modules/mesh/1DCart/makefile index e078727..f1045e9 100644 --- a/src/modules/mesh/1DCart/makefile +++ b/src/modules/mesh/1DCart/makefile @@ -1,8 +1,5 @@ -all: moduleMesh1DCart.o moduleMesh1DCartRead.o +all: moduleMesh1DCart.o moduleMesh1DCart.o: moduleMesh1DCart.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleMesh1DCartRead.o: moduleMesh1DCart.o moduleMesh1DCartRead.f90 - $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ - diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 34209fa..12898d5 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -308,27 +308,27 @@ MODULE moduleMesh1DCart END SUBROUTINE partialDerSegm !Computes local stiffness matrix - FUNCTION elemKSegm(self) RESULT(ke) + PURE FUNCTION elemKSegm(self) RESULT(localK) IMPLICIT NONE CLASS(meshVol1DCartSegm), INTENT(in):: self - REAL(8):: ke(1:2,1:2) + REAL(8), ALLOCATABLE:: localK(:,:) REAL(8):: Xii(1:3) REAL(8):: dPsi(1:1, 1:2) REAL(8):: invJ(1), detJ INTEGER:: l - ke = 0.D0 + ALLOCATE(localK(1:2,1:2)) + localK = 0.D0 Xii = 0.D0 - DO l = 1, 3 xii(1) = corSeg(l) dPsi = self%dPsi(Xii) detJ = self%detJac(Xii, dPsi) invJ = self%invJac(Xii, dPsi) - ke = ke + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & - RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & - wSeg(l)/detJ + localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & + RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & + wSeg(l)/detJ END DO @@ -522,6 +522,117 @@ MODULE moduleMesh1DCart END FUNCTION invJ1DCart + SUBROUTINE connectMesh1DCart(self) + IMPLICIT NONE + + CLASS(meshParticle), INTENT(inout):: self + INTEGER:: e, et + + DO e = 1, self%numVols + !Connect Vol-Vol + DO et = 1, self%numVols + IF (e /= et) THEN + CALL connectVolVol(self%vols(e)%obj, self%vols(et)%obj) + + END IF + + END DO + + !Connect Vol-Edge + DO et = 1, self%numEdges + CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) + + END DO + + END DO + + END SUBROUTINE connectMesh1DCart + + SUBROUTINE connectVolVol(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol), INTENT(inout):: elemA + CLASS(meshVol), INTENT(inout):: elemB + + SELECT TYPE(elemA) + TYPE IS(meshVol1DCartSegm) + SELECT TYPE(elemB) + TYPE IS(meshVol1DCartSegm) + CALL connectSegmSegm(elemA, elemB) + + END SELECT + + END SELECT + + END SUBROUTINE connectVolVol + + SUBROUTINE connectSegmSegm(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemA + CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemB + + IF (.NOT. ASSOCIATED(elemA%e1) .AND. & + elemA%n2%n == elemB%n1%n) THEN + elemA%e1 => elemB + elemB%e2 => elemA + + END IF + + IF (.NOT. ASSOCIATED(elemA%e2) .AND. & + elemA%n1%n == elemB%n2%n) THEN + + elemA%e2 => elemB + elemB%e1 => elemA + + END IF + + END SUBROUTINE connectSegmSegm + + SUBROUTINE connectVolEdge(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol), INTENT(inout):: elemA + CLASS(meshEdge), INTENT(inout):: elemB + + SELECT TYPE(elemA) + TYPE IS (meshVol1DCartSegm) + SELECT TYPE(elemB) + CLASS IS(meshEdge1DCart) + CALL connectSegmEdge(elemA, elemB) + + END SELECT + + END SELECT + + END SUBROUTINE connectVolEdge + + SUBROUTINE connectSegmEdge(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemA + CLASS(meshEdge1DCart), INTENT(inout), TARGET:: elemB + + IF (.NOT. ASSOCIATED(elemA%e1) .AND. & + elemA%n2%n == elemB%n1%n) THEN + + elemA%e1 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + + END IF + + IF (.NOT. ASSOCIATED(elemA%e2) .AND. & + elemA%n1%n == elemB%n1%n) THEN + + elemA%e2 => elemB + elemB%e1 => elemA + + END IF + + END SUBROUTINE connectSegmEdge END MODULE moduleMesh1DCart diff --git a/src/modules/mesh/1DCart/moduleMesh1DCartRead.f90 b/src/modules/mesh/1DCart/moduleMesh1DCartRead.f90 deleted file mode 100644 index d69cfdc..0000000 --- a/src/modules/mesh/1DCart/moduleMesh1DCartRead.f90 +++ /dev/null @@ -1,264 +0,0 @@ -MODULE moduleMesh1DCartRead - USE moduleMesh - USE moduleMesh1DCart - - !TODO: make this abstract to allow different mesh formats - TYPE, EXTENDS(meshGeneric):: mesh1DCartGeneric - CONTAINS - PROCEDURE, PASS:: init => init1DCartMesh - PROCEDURE, PASS:: readMesh => readMesh1DCart - - END TYPE - - INTERFACE connected - MODULE PROCEDURE connectedVolVol, connectedVolEdge - - END INTERFACE connected - - CONTAINS - !Init 1D mesh - SUBROUTINE init1DCartMesh(self, meshFormat) - USE moduleMesh - USE moduleErrors - IMPLICIT NONE - - CLASS(mesh1DCartGeneric), INTENT(out):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat - - SELECT CASE(meshFormat) - CASE ("gmsh") - self%printOutput => printOutputGmsh - self%printColl => printCollGmsh - self%printEM => printEMGmsh - - CASE DEFAULT - CALL criticalError("Mesh type " // meshFormat // " not supported.", "init1DCart") - - END SELECT - - END SUBROUTINE init1DCartMesh - - !Reads 1D mesh - SUBROUTINE readMesh1DCart(self, filename) - USE moduleBoundary - IMPLICIT NONE - - CLASS(mesh1DCartGeneric), INTENT(inout):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: filename - REAL(8):: x - INTEGER:: p(1:2) - INTEGER:: e, et, n, eTemp, elemType, bt - INTEGER:: totalNumElem - INTEGER:: boundaryType - - !Open file mesh - OPEN(10, FILE=TRIM(filename)) - !Skip header - READ(10, *) - READ(10, *) - READ(10, *) - READ(10, *) - !Read number of nodes - READ(10, *) self%numNodes - !Allocate required matrices and vectors - ALLOCATE(self%nodes(1:self%numNodes)) - 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 - !Read nodes coordinates. Only relevant for x - DO e = 1, self%numNodes - READ(10, *) n, x - ALLOCATE(meshNode1DCart:: self%nodes(n)%obj) - CALL self%nodes(n)%obj%init(n, (/ x, 0.D0, 0.D0 /)) - - END DO - !Skips comments - READ(10, *) - READ(10, *) - !Reads the total number of elements (edges+vol) - READ(10, *) totalNumElem - self%numEdges = 0 - DO e = 1, totalNumElem - READ(10, *) eTemp, elemType - IF (elemType == 15) THEN !15 is physical node in GMSH - self%numEdges = e - - END IF - - END DO - - !Substract the number of edges to the total number of elements - !to obtain the number of volume elements - self%numVols = totalNumelem - self%numEdges - !Allocates arrays - ALLOCATE(self%edges(1:self%numEdges)) - ALLOCATE(self%vols(1:self%numVols)) - - !Go back to the beginning of reading elements - DO e = 1, totalNumelem - BACKSPACE(10) - - END DO - - !Reads edges - DO e = 1, self%numEdges - READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1) - !Associate boundary condition - bt = getBoundaryId(boundaryType) - - ALLOCATE(meshEdge1DCart:: self%edges(e)%obj) - - CALL self%edges(e)%obj%init(n, p(1:1), bt, boundaryType) - - END DO - - !Read and initialize volumes - DO e = 1, self%numVols - READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2) - ALLOCATE(meshVol1DCartSegm:: self%vols(e)%obj) - CALL self%vols(e)%obj%init(n - self%numEdges, p(1:2)) - - END DO - - CLOSE(10) - - !Build connectivity between elements - DO e = 1, self%numVols - !Connectivity between volumes - DO et = 1, self%numVols - IF (e /= et) THEN - CALL connected(self%vols(e)%obj, self%vols(et)%obj) - - END IF - - END DO - - !Connectivity betwen vols and edges - DO et = 1, self%numEdges - CALL connected(self%vols(e)%obj, self%edges(et)%obj) - - END DO - - !Constructs the global K matrix - CALL constructGlobalK(self%K, self%vols(e)%obj) - - END DO - - END SUBROUTINE readMesh1DCart - - SUBROUTINE connectedVolVol(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshVol), INTENT(inout):: elemB - - SELECT TYPE(elemA) - TYPE IS(meshVol1DCartSegm) - SELECT TYPE(elemB) - TYPE IS(meshVol1DCartSegm) - CALL connectedSegmSegm(elemA, elemB) - - END SELECT - - END SELECT - - END SUBROUTINE connectedVolVol - - SUBROUTINE connectedSegmSegm(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemA - CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemB - - IF (.NOT. ASSOCIATED(elemA%e1) .AND. & - elemA%n2%n == elemB%n1%n) THEN - elemA%e1 => elemB - elemB%e2 => elemA - - END IF - - IF (.NOT. ASSOCIATED(elemA%e2) .AND. & - elemA%n1%n == elemB%n2%n) THEN - - elemA%e2 => elemB - elemB%e1 => elemA - - END IF - - END SUBROUTINE connectedSegmSegm - - SUBROUTINE connectedVolEdge(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshEdge), INTENT(inout):: elemB - - SELECT TYPE(elemA) - TYPE IS (meshVol1DCartSegm) - SELECT TYPE(elemB) - CLASS IS(meshEdge1DCart) - CALL connectedSegmEdge(elemA, elemB) - - END SELECT - - END SELECT - - END SUBROUTINE connectedVolEdge - - SUBROUTINE connectedSegmEdge(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemA - CLASS(meshEdge1DCart), INTENT(inout), TARGET:: elemB - - IF (.NOT. ASSOCIATED(elemA%e1) .AND. & - elemA%n2%n == elemB%n1%n) THEN - - elemA%e1 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal - - END IF - - IF (.NOT. ASSOCIATED(elemA%e2) .AND. & - elemA%n1%n == elemB%n1%n) THEN - - elemA%e2 => elemB - elemB%e1 => elemA - - END IF - - END SUBROUTINE connectedSegmEdge - - SUBROUTINE constructGlobalK(K, elem) - IMPLICIT NONE - - REAL(8), INTENT(inout):: K(1:,1:) - CLASS(meshVol), INTENT(in):: elem - REAL(8):: localK(1:2,1:2) - INTEGER:: i, j - INTEGER:: n(1:2) - - SELECT TYPE(elem) - TYPE IS(meshVol1DCartSegm) - localK = elem%elemK() - n = (/ elem%n1%n, elem%n2%n /) - - CLASS DEFAULT - n = 0 - localK = 0.D0 - - END SELECT - - DO i = 1, 2 - DO j = 1, 2 - K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j) - END DO - END DO - - END SUBROUTINE constructGlobalK - -END MODULE moduleMesh1DCartRead diff --git a/src/modules/mesh/1DRad/makefile b/src/modules/mesh/1DRad/makefile index 13e34ab..eeacc76 100644 --- a/src/modules/mesh/1DRad/makefile +++ b/src/modules/mesh/1DRad/makefile @@ -1,8 +1,5 @@ -all: moduleMesh1DRad.o moduleMesh1DRadRead.o +all: moduleMesh1DRad.o moduleMesh1DRad.o: moduleMesh1DRad.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleMesh1DRadRead.o: moduleMesh1DRad.o moduleMesh1DRadRead.f90 - $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ - diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 9678dff..dc47c62 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -312,19 +312,20 @@ MODULE moduleMesh1DRad END SUBROUTINE partialDerRad !Computes local stiffness matrix - PURE FUNCTION elemKRad(self) RESULT(ke) + PURE FUNCTION elemKRad(self) RESULT(localK) USE moduleConstParam, ONLY: PI2 IMPLICIT NONE CLASS(meshVol1DRadSegm), INTENT(in):: self - REAL(8):: ke(1:2,1:2) + REAL(8), ALLOCATABLE:: localK(:,:) REAL(8):: Xii(1:3) REAL(8):: dPsi(1:1, 1:2) REAL(8):: invJ(1), detJ REAL(8):: r, fPsi(1:2) INTEGER:: l - ke = 0.D0 + ALLOCATE(localK(1:2, 1:2)) + localK = 0.D0 Xii = 0.D0 DO l = 1, 3 xii(1) = corSeg(l) @@ -333,13 +334,13 @@ MODULE moduleMesh1DRad invJ = self%invJac(Xii, dPsi) fPsi = self%fPsi(Xii) r = DOT_PRODUCT(fPsi, self%r) - ke = ke + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & - RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & - r*wSeg(l)/detJ + localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & + RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & + r*wSeg(l)/detJ END DO - ke = ke*PI2 + localK = localK*PI2 END FUNCTION elemKRad @@ -532,5 +533,117 @@ MODULE moduleMesh1DRad END FUNCTION invJ1DRad + SUBROUTINE connectMesh1DRad(self) + IMPLICIT NONE + + CLASS(meshParticle), INTENT(inout):: self + INTEGER:: e, et + + DO e = 1, self%numVols + !Connect Vol-Vol + DO et = 1, self%numVols + IF (e /= et) THEN + CALL connectVolVol(self%vols(e)%obj, self%vols(et)%obj) + + END IF + + END DO + + !Connect Vol-Edge + DO et = 1, self%numEdges + CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) + + END DO + + END DO + + END SUBROUTINE connectMesh1DRad + + SUBROUTINE connectVolVol(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol), INTENT(inout):: elemA + CLASS(meshVol), INTENT(inout):: elemB + + SELECT TYPE(elemA) + TYPE IS(meshVol1DRadSegm) + SELECT TYPE(elemB) + TYPE IS(meshVol1DRadSegm) + CALL connectSegmSegm(elemA, elemB) + + END SELECT + + END SELECT + + END SUBROUTINE connectVolVol + + SUBROUTINE connectSegmSegm(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemA + CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemB + + IF (.NOT. ASSOCIATED(elemA%e1) .AND. & + elemA%n2%n == elemB%n1%n) THEN + elemA%e1 => elemB + elemB%e2 => elemA + + END IF + + IF (.NOT. ASSOCIATED(elemA%e2) .AND. & + elemA%n1%n == elemB%n2%n) THEN + + elemA%e2 => elemB + elemB%e1 => elemA + + END IF + + END SUBROUTINE connectSegmSegm + + SUBROUTINE connectVolEdge(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol), INTENT(inout):: elemA + CLASS(meshEdge), INTENT(inout):: elemB + + SELECT TYPE(elemA) + TYPE IS (meshVol1DRadSegm) + SELECT TYPE(elemB) + CLASS IS(meshEdge1DRad) + CALL connectSegmEdge(elemA, elemB) + + END SELECT + + END SELECT + + END SUBROUTINE connectVolEdge + + SUBROUTINE connectSegmEdge(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemA + CLASS(meshEdge1DRad), INTENT(inout), TARGET:: elemB + + IF (.NOT. ASSOCIATED(elemA%e1) .AND. & + elemA%n2%n == elemB%n1%n) THEN + + elemA%e1 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + + END IF + + IF (.NOT. ASSOCIATED(elemA%e2) .AND. & + elemA%n1%n == elemB%n1%n) THEN + + elemA%e2 => elemB + elemB%e1 => elemA + + END IF + + END SUBROUTINE connectSegmEdge + END MODULE moduleMesh1DRad diff --git a/src/modules/mesh/1DRad/moduleMesh1DRadRead.f90 b/src/modules/mesh/1DRad/moduleMesh1DRadRead.f90 deleted file mode 100644 index 4f14b5f..0000000 --- a/src/modules/mesh/1DRad/moduleMesh1DRadRead.f90 +++ /dev/null @@ -1,264 +0,0 @@ -MODULE moduleMesh1DRadRead - USE moduleMesh - USE moduleMesh1DRad - - !TODO: make this abstract to allow different mesh formats - TYPE, EXTENDS(meshGeneric):: mesh1DRadGeneric - CONTAINS - PROCEDURE, PASS:: init => init1DRadMesh - PROCEDURE, PASS:: readMesh => readMesh1DRad - - END TYPE - - INTERFACE connected - MODULE PROCEDURE connectedVolVol, connectedVolEdge - - END INTERFACE connected - - CONTAINS - !Init 1D mesh - SUBROUTINE init1DRadMesh(self, meshFormat) - USE moduleMesh - USE moduleErrors - IMPLICIT NONE - - CLASS(mesh1DRadGeneric), INTENT(out):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat - - SELECT CASE(meshFormat) - CASE ("gmsh") - self%printOutput => printOutputGmsh - self%printColl => printCollGmsh - self%printEM => printEMGmsh - - CASE DEFAULT - CALL criticalError("Mesh type " // meshFormat // " not supported.", "init1DRad") - - END SELECT - - END SUBROUTINE init1DRadMesh - - !Reads 1D mesh - SUBROUTINE readMesh1DRad(self, filename) - USE moduleBoundary - IMPLICIT NONE - - CLASS(mesh1DRadGeneric), INTENT(inout):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: filename - REAL(8):: x - INTEGER:: p(1:2) - INTEGER:: e, et, n, eTemp, elemType, bt - INTEGER:: totalNumElem - INTEGER:: boundaryType - - !Open file mesh - OPEN(10, FILE=TRIM(filename)) - !Skip header - READ(10, *) - READ(10, *) - READ(10, *) - READ(10, *) - !Read number of nodes - READ(10, *) self%numNodes - !Allocate required matrices and vectors - ALLOCATE(self%nodes(1:self%numNodes)) - 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 - !Read nodes coordinates. Only relevant for x - DO e = 1, self%numNodes - READ(10, *) n, x - ALLOCATE(meshNode1DRad:: self%nodes(n)%obj) - CALL self%nodes(n)%obj%init(n, (/ x, 0.D0, 0.D0 /)) - - END DO - !Skips comments - READ(10, *) - READ(10, *) - !Reads the total number of elements (edges+vol) - READ(10, *) totalNumElem - self%numEdges = 0 - DO e = 1, totalNumElem - READ(10, *) eTemp, elemType - IF (elemType == 15) THEN !15 is physical node in GMSH - self%numEdges = e - - END IF - - END DO - - !Substract the number of edges to the total number of elements - !to obtain the number of volume elements - self%numVols = totalNumelem - self%numEdges - !Allocates arrays - ALLOCATE(self%edges(1:self%numEdges)) - ALLOCATE(self%vols(1:self%numVols)) - - !Go back to the beginning of reading elements - DO e = 1, totalNumelem - BACKSPACE(10) - - END DO - - !Reads edges - DO e = 1, self%numEdges - READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1) - !Associate boundary condition - bt = getBoundaryId(boundaryType) - - ALLOCATE(meshEdge1DRad:: self%edges(e)%obj) - - CALL self%edges(e)%obj%init(n, p(1:1), bt, boundaryType) - - END DO - - !Read and initialize volumes - DO e = 1, self%numVols - READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2) - ALLOCATE(meshVol1DRadSegm:: self%vols(e)%obj) - CALL self%vols(e)%obj%init(n - self%numEdges, p(1:2)) - - END DO - - CLOSE(10) - - !Build connectivity between elements - DO e = 1, self%numVols - !Connectivity between volumes - DO et = 1, self%numVols - IF (e /= et) THEN - CALL connected(self%vols(e)%obj, self%vols(et)%obj) - - END IF - - END DO - - !Connectivity betwen vols and edges - DO et = 1, self%numEdges - CALL connected(self%vols(e)%obj, self%edges(et)%obj) - - END DO - - !Constructs the global K matrix - CALL constructGlobalK(self%K, self%vols(e)%obj) - - END DO - - END SUBROUTINE readMesh1DRad - - SUBROUTINE connectedVolVol(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshVol), INTENT(inout):: elemB - - SELECT TYPE(elemA) - TYPE IS(meshVol1DRadSegm) - SELECT TYPE(elemB) - TYPE IS(meshVol1DRadSegm) - CALL connectedSegmSegm(elemA, elemB) - - END SELECT - - END SELECT - - END SUBROUTINE connectedVolVol - - SUBROUTINE connectedSegmSegm(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemA - CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemB - - IF (.NOT. ASSOCIATED(elemA%e1) .AND. & - elemA%n2%n == elemB%n1%n) THEN - elemA%e1 => elemB - elemB%e2 => elemA - - END IF - - IF (.NOT. ASSOCIATED(elemA%e2) .AND. & - elemA%n1%n == elemB%n2%n) THEN - - elemA%e2 => elemB - elemB%e1 => elemA - - END IF - - END SUBROUTINE connectedSegmSegm - - SUBROUTINE connectedVolEdge(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshEdge), INTENT(inout):: elemB - - SELECT TYPE(elemA) - TYPE IS (meshVol1DRadSegm) - SELECT TYPE(elemB) - CLASS IS(meshEdge1DRad) - CALL connectedSegmEdge(elemA, elemB) - - END SELECT - - END SELECT - - END SUBROUTINE connectedVolEdge - - SUBROUTINE connectedSegmEdge(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemA - CLASS(meshEdge1DRad), INTENT(inout), TARGET:: elemB - - IF (.NOT. ASSOCIATED(elemA%e1) .AND. & - elemA%n2%n == elemB%n1%n) THEN - - elemA%e1 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal - - END IF - - IF (.NOT. ASSOCIATED(elemA%e2) .AND. & - elemA%n1%n == elemB%n1%n) THEN - - elemA%e2 => elemB - elemB%e1 => elemA - - END IF - - END SUBROUTINE connectedSegmEdge - - SUBROUTINE constructGlobalK(K, elem) - IMPLICIT NONE - - REAL(8), INTENT(inout):: K(1:,1:) - CLASS(meshVol), INTENT(in):: elem - REAL(8):: localK(1:2,1:2) - INTEGER:: i, j - INTEGER:: n(1:2) - - SELECT TYPE(elem) - TYPE IS(meshVol1DRadSegm) - localK = elem%elemK() - n = (/ elem%n1%n, elem%n2%n /) - - CLASS DEFAULT - n = 0 - localK = 0.D0 - - END SELECT - - DO i = 1, 2 - DO j = 1, 2 - K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j) - END DO - END DO - - END SUBROUTINE constructGlobalK - -END MODULE moduleMesh1DRadRead diff --git a/src/modules/mesh/2DCart/makefile b/src/modules/mesh/2DCart/makefile index 65211cb..06c5288 100644 --- a/src/modules/mesh/2DCart/makefile +++ b/src/modules/mesh/2DCart/makefile @@ -1,8 +1,5 @@ -all : moduleMesh2DCart.o moduleMesh2DCartRead.o +all : moduleMesh2DCart.o moduleMesh2DCart.o: moduleMesh2DCart.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleMesh2DCartRead.o: moduleMesh2DCart.o moduleMesh2DCartRead.f90 - $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ - diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index d78c663..11e6454 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -414,17 +414,18 @@ MODULE moduleMesh2DCart END SUBROUTINE partialDerQuad !Computes element local stiffness matrix - PURE FUNCTION elemKQuad(self) RESULT(ke) + PURE FUNCTION elemKQuad(self) RESULT(localK) IMPLICIT NONE CLASS(meshVol2DCartQuad), INTENT(in):: self + REAL(8), ALLOCATABLE:: localK(:,:) REAL(8):: xi(1:3) REAL(8):: fPsi(1:4), dPsi(1:2,1:4) - REAL(8):: ke(1:4,1:4) REAL(8):: invJ(1:2,1:2), detJ INTEGER:: l, m - ke=0.D0 + ALLOCATE(localK(1:4, 1:4)) + localK=0.D0 xi=0.D0 !Start 2D Gauss Quad Integral DO l=1, 3 @@ -436,7 +437,7 @@ MODULE moduleMesh2DCart fPsi = self%fPsi(xi) detJ = self%detJac(xi,dPsi) invJ = self%invJac(xi,dPsi) - ke = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wQuad(l)*wQuad(m)/detJ + localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wQuad(l)*wQuad(m)/detJ END DO END DO @@ -776,27 +777,28 @@ MODULE moduleMesh2DCart END SUBROUTINE partialDerTria !Computes element local stiffness matrix - PURE FUNCTION elemKTria(self) RESULT(ke) + PURE FUNCTION elemKTria(self) RESULT(localK) IMPLICIT NONE CLASS(meshVol2DCartTria), INTENT(in):: self + REAL(8), ALLOCATABLE:: localK(:,:) REAL(8):: xi(1:3) REAL(8):: fPsi(1:3), dPsi(1:2,1:3) - REAL(8):: ke(1:3,1:3) REAL(8):: invJ(1:2,1:2), detJ INTEGER:: l - ke=0.D0 + ALLOCATE(localK(1:4, 1:4)) + localK=0.D0 xi=0.D0 !Start 2D Gauss Quad Integral DO l=1, 4 - xi(1) = xi1Tria(l) - xi(2) = xi2Tria(l) - dPsi = self%dPsi(xi) - detJ = self%detJac(xi,dPsi) - invJ = self%invJac(xi,dPsi) - fPsi = self%fPsi(xi) - ke = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ + xi(1) = xi1Tria(l) + xi(2) = xi2Tria(l) + dPsi = self%dPsi(xi) + detJ = self%detJac(xi,dPsi) + invJ = self%invJac(xi,dPsi) + fPsi = self%fPsi(xi) + localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ END DO @@ -1017,4 +1019,449 @@ MODULE moduleMesh2DCart END FUNCTION invJ2DCart + SUBROUTINE connectMesh2DCart(self) + IMPLICIT NONE + + CLASS(meshParticle), INTENT(inout):: self + INTEGER:: e, et + + DO e = 1, self%numVols + !Connect Vol-Vol + DO et = 1, self%numVols + IF (e /= et) THEN + CALL connectVolVol(self%vols(e)%obj, self%vols(et)%obj) + + END IF + + END DO + + !Connect Vol-Edge + DO et = 1, self%numEdges + CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) + + END DO + + END DO + + END SUBROUTINE connectMesh2DCart + + !Selects type of elements to build connection + SUBROUTINE connectVolVol(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol), INTENT(inout):: elemA + CLASS(meshVol), INTENT(inout):: elemB + + SELECT TYPE(elemA) + TYPE IS(meshVol2DCartQuad) + !Element A is a quadrilateral + SELECT TYPE(elemB) + TYPE IS(meshVol2DCartQuad) + !Element B is a quadrilateral + CALL connectQuadQuad(elemA, elemB) + + TYPE IS(meshVol2DCartTria) + !Element B is a triangle + CALL connectQuadTria(elemA, elemB) + + END SELECT + + TYPE IS(meshVol2DCartTria) + !Element A is a Triangle + SELECT TYPE(elemB) + TYPE IS(meshVol2DCartQuad) + !Element B is a quadrilateral + CALL connectQuadTria(elemB, elemA) + + TYPE IS(meshVol2DCartTria) + !Element B is a triangle + CALL connectTriaTria(elemA, elemB) + + END SELECT + + END SELECT + + END SUBROUTINE connectVolVol + + SUBROUTINE connectVolEdge(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol), INTENT(inout):: elemA + CLASS(meshEdge), INTENT(inout):: elemB + + SELECT TYPE(elemB) + CLASS IS(meshEdge2DCart) + SELECT TYPE(elemA) + TYPE IS(meshVol2DCartQuad) + !Element A is a quadrilateral + CALL connectQuadEdge(elemA, elemB) + + TYPE IS(meshVol2DCartTria) + !Element A is a triangle + CALL connectTriaEdge(elemA, elemB) + + END SELECT + + END SELECT + + END SUBROUTINE connectVolEdge + + SUBROUTINE connectQuadQuad(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA + CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemB + + !Check direction 1 + IF (.NOT. ASSOCIATED(elemA%e1) .AND. & + elemA%n1%n == elemB%n4%n .AND. & + elemA%n2%n == elemB%n3%n) THEN + elemA%e1 => elemB + elemB%e3 => elemA + + END IF + + !Check direction 2 + IF (.NOT. ASSOCIATED(elemA%e2) .AND. & + elemA%n2%n == elemB%n1%n .AND. & + elemA%n3%n == elemB%n4%n) THEN + elemA%e2 => elemB + elemB%e4 => elemA + + END IF + + !Check direction 3 + IF (.NOT. ASSOCIATED(elemA%e3) .AND. & + elemA%n3%n == elemB%n2%n .AND. & + elemA%n4%n == elemB%n1%n) THEN + elemA%e3 => elemB + elemB%e1 => elemA + + END IF + + !Check direction 4 + IF (.NOT. ASSOCIATED(elemA%e4) .AND. & + elemA%n4%n == elemB%n3%n .AND. & + elemA%n1%n == elemB%n2%n) THEN + elemA%e4 => elemB + elemB%e2 => elemA + + END IF + + END SUBROUTINE connectQuadQuad + + SUBROUTINE connectQuadTria(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA + CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemB + + !Check direction 1 + IF (.NOT. ASSOCIATED(elemA%e1)) THEN + IF (elemA%n1%n == elemB%n1%n .AND. & + elemA%n2%n == elemB%n3%n) THEN + elemA%e1 => elemB + elemB%e3 => elemA + + ELSEIF (elemA%n1%n == elemB%n3%n .AND. & + elemA%n2%n == elemB%n2%n) THEN + elemA%e1 => elemB + elemB%e2 => elemA + + ELSEIF (elemA%n1%n == elemB%n2%n .AND. & + elemA%n2%n == elemB%n1%n) THEN + elemA%e1 => elemB + elemB%e1 => elemA + + END IF + + END IF + + !Check direction 2 + IF (.NOT. ASSOCIATED(elemA%e2)) THEN + IF (elemA%n2%n == elemB%n1%n .AND. & + elemA%n3%n == elemB%n3%n) THEN + elemA%e2 => elemB + elemB%e3 => elemA + + ELSEIF (elemA%n2%n == elemB%n3%n .AND. & + elemA%n3%n == elemB%n2%n) THEN + elemA%e2 => elemB + elemB%e2 => elemA + + ELSEIF (elemA%n2%n == elemB%n2%n .AND. & + elemA%n3%n == elemB%n1%n) THEN + elemA%e2 => elemB + elemB%e1 => elemA + + END IF + + END IF + + !Check direction 3 + IF (.NOT. ASSOCIATED(elemA%e3)) THEN + IF (elemA%n3%n == elemB%n1%n .AND. & + elemA%n4%n == elemB%n3%n) THEN + elemA%e3 => elemB + elemB%e3 => elemA + + ELSEIF (elemA%n3%n == elemB%n3%n .AND. & + elemA%n4%n == elemB%n2%n) THEN + elemA%e3 => elemB + elemB%e2 => elemA + + ELSEIF (elemA%n3%n == elemB%n2%n .AND. & + elemA%n4%n == elemB%n1%n) THEN + elemA%e3 => elemB + elemB%e1 => elemA + + END IF + + END IF + + !Check direction 4 + IF (.NOT. ASSOCIATED(elemA%e4)) THEN + IF (elemA%n4%n == elemB%n1%n .AND. & + elemA%n1%n == elemB%n3%n) THEN + elemA%e4 => elemB + elemB%e3 => elemA + + ELSEIF (elemA%n4%n == elemB%n3%n .AND. & + elemA%n1%n == elemB%n2%n) THEN + elemA%e4 => elemB + elemB%e2 => elemA + + ELSEIF (elemA%n4%n == elemB%n2%n .AND. & + elemA%n1%n == elemB%n1%n) THEN + elemA%e4 => elemB + elemB%e1 => elemA + + END IF + + END IF + + END SUBROUTINE connectQuadTria + + SUBROUTINE connectTriaTria(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemA + CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemB + + !Check direction 1 + IF (.NOT. ASSOCIATED(elemA%e1)) THEN + IF (elemA%n1%n == elemB%n1%n .AND. & + elemA%n2%n == elemB%n3%n) THEN + elemA%e1 => elemB + elemB%e3 => elemA + + ELSEIF (elemA%n1%n == elemB%n2%n .AND. & + elemA%n2%n == elemB%n1%n) THEN + elemA%e1 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n1%n == elemB%n3%n .AND. & + elemA%n2%n == elemB%n2%n) THEN + elemA%e1 => elemB + elemB%e2 => elemA + + END IF + + END IF + + !Check direction 2 + IF (.NOT. ASSOCIATED(elemA%e2)) THEN + IF (elemA%n2%n == elemB%n1%n .AND. & + elemA%n3%n == elemB%n3%n) THEN + elemA%e2 => elemB + elemB%e3 => elemA + + ELSEIF (elemA%n2%n == elemB%n2%n .AND. & + elemA%n3%n == elemB%n1%n) THEN + elemA%e2 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n2%n == elemB%n3%n .AND. & + elemA%n3%n == elemB%n2%n) THEN + elemA%e2 => elemB + elemB%e2 => elemA + + END IF + + + END IF + + !Check direction 3 + IF (.NOT. ASSOCIATED(elemA%e3)) THEN + IF (elemA%n3%n == elemB%n1%n .AND. & + elemA%n1%n == elemB%n3%n) THEN + elemA%e3 => elemB + elemB%e3 => elemA + + ELSEIF (elemA%n3%n == elemB%n2%n .AND. & + elemA%n1%n == elemB%n1%n) THEN + elemA%e3 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n3%n == elemB%n3%n .AND. & + elemA%n1%n == elemB%n2%n) THEN + elemA%e3 => elemB + elemB%e2 => elemA + + END IF + + + END IF + + END SUBROUTINE connectTriaTria + + SUBROUTINE connectQuadEdge(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA + CLASS(meshEdge2DCart), INTENT(inout), TARGET:: elemB + + !Check direction 1 + IF (.NOT. ASSOCIATED(elemA%e1)) THEN + IF (elemA%n1%n == elemB%n1%n .AND. & + elemA%n2%n == elemB%n2%n) THEN + elemA%e1 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n1%n == elemB%n2%n .AND. & + elemA%n2%n == elemB%n1%n) THEN + elemA%e1 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + + END IF + + END IF + + !Check direction 2 + IF (.NOT. ASSOCIATED(elemA%e2)) THEN + IF (elemA%n2%n == elemB%n1%n .AND. & + elemA%n3%n == elemB%n2%n) THEN + elemA%e2 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n2%n == elemB%n2%n .AND. & + elemA%n3%n == elemB%n1%n) THEN + elemA%e2 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + + END IF + + END IF + + !Check direction 3 + IF (.NOT. ASSOCIATED(elemA%e3)) THEN + IF (elemA%n3%n == elemB%n1%n .AND. & + elemA%n4%n == elemB%n2%n) THEN + elemA%e3 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n3%n == elemB%n2%n .AND. & + elemA%n4%n == elemB%n1%n) THEN + elemA%e3 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + + END IF + + END IF + + !Check direction 4 + IF (.NOT. ASSOCIATED(elemA%e4)) THEN + IF (elemA%n4%n == elemB%n1%n .AND. & + elemA%n1%n == elemB%n2%n) THEN + elemA%e4 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n4%n == elemB%n2%n .AND. & + elemA%n1%n == elemB%n1%n) THEN + elemA%e4 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + + END IF + + END IF + + END SUBROUTINE connectQuadEdge + + SUBROUTINE connectTriaEdge(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemA + CLASS(meshEdge2DCart), INTENT(inout), TARGET:: elemB + + !Check direction 1 + IF (.NOT. ASSOCIATED(elemA%e1)) THEN + IF (elemA%n1%n == elemB%n1%n .AND. & + elemA%n2%n == elemB%n2%n) THEN + elemA%e1 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n1%n == elemB%n2%n .AND. & + elemA%n2%n == elemB%n1%n) THEN + elemA%e1 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + + END IF + + END IF + + !Check direction 2 + IF (.NOT. ASSOCIATED(elemA%e2)) THEN + IF (elemA%n2%n == elemB%n1%n .AND. & + elemA%n3%n == elemB%n2%n) THEN + elemA%e2 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n2%n == elemB%n2%n .AND. & + elemA%n3%n == elemB%n1%n) THEN + elemA%e2 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + + END IF + + END IF + + !Check direction 3 + IF (.NOT. ASSOCIATED(elemA%e3)) THEN + IF (elemA%n3%n == elemB%n1%n .AND. & + elemA%n1%n == elemB%n2%n) THEN + elemA%e3 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n3%n == elemB%n2%n .AND. & + elemA%n1%n == elemB%n1%n) THEN + elemA%e3 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + + END IF + + END IF + + END SUBROUTINE connectTriaEdge + END MODULE moduleMesh2DCart diff --git a/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 b/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 deleted file mode 100644 index 96664df..0000000 --- a/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 +++ /dev/null @@ -1,620 +0,0 @@ -MODULE moduleMesh2DCartRead - USE moduleMesh - USE moduleMesh2DCart - - TYPE, EXTENDS(meshGeneric):: mesh2DCartGeneric - CONTAINS - PROCEDURE, PASS:: init => init2DCartMesh - PROCEDURE, PASS:: readMesh => readMesh2DCartGmsh - - END TYPE - - INTERFACE connected - MODULE PROCEDURE connectedVolVol, connectedVolEdge - - END INTERFACE connected - - CONTAINS - !Init mesh - SUBROUTINE init2DCartMesh(self, meshFormat) - USE moduleMesh - USE moduleErrors - IMPLICIT NONE - - CLASS(mesh2DCartGeneric), INTENT(out):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat - - SELECT CASE(meshFormat) - CASE ("gmsh") - self%printOutput => printOutputGmsh - self%printColl => printCollGmsh - self%printEM => printEMGmsh - - CASE DEFAULT - CALL criticalError("Mesh type " // meshFormat // " not supported.", "init2DCartMesh") - - END SELECT - - END SUBROUTINE init2DCartMesh - - !Read mesh from gmsh file - SUBROUTINE readMesh2DCartGmsh(self, filename) - USE moduleBoundary - IMPLICIT NONE - - CLASS(mesh2DCartGeneric), INTENT(inout):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: filename - REAL(8):: x, y - INTEGER:: p(1:4) - INTEGER:: e=0, et=0, n=0, eTemp=0, elemType=0, bt = 0 - INTEGER:: totalNumElem - INTEGER:: boundaryType - - !Read msh - OPEN(10, FILE=TRIM(filename)) - !Skip header - READ(10, *) - READ(10, *) - READ(10, *) - READ(10, *) - !Read number of nodes - READ(10, *) self%numNodes - !Allocate required matrices and vectors - ALLOCATE(self%nodes(1:self%numNodes)) - 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 - !Read node cartesian coordinates (x=x, y=y, z=null) - DO e=1, self%numNodes - READ(10, *) n, x, y - ALLOCATE(meshNode2DCart:: self%nodes(n)%obj) - CALL self%nodes(n)%obj%init(n, (/x, y, 0.D0 /)) - - END DO - !Skips comments - READ(10, *) - READ(10, *) - !Reads Totalnumber of elements - READ(10, *) TotalnumElem - !counts edges and volume elements - self%numEdges = 0 - DO e=1, TotalnumElem - READ(10,*) eTemp, elemType - IF (elemType==1) THEN - self%numEdges=e - END IF - END DO - !Substract the number of edges to the total number of elements - !to obtain the number of volume elements - self%numVols = TotalnumElem - self%numEdges - !Allocates arrays - ALLOCATE(self%edges(1:self%numEdges)) - ALLOCATE(self%vols(1:self%numVols)) - - !Go back to the beggining to read elements - DO e=1, totalNumElem - BACKSPACE(10) - END DO - - !Reads edges - DO e=1, self%numEdges - READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2) - !Associate boundary condition procedure. - bt = getBoundaryId(boundaryType) - - ALLOCATE(meshEdge2DCart:: self%edges(e)%obj) - - CALL self%edges(e)%obj%init(n, p(1:2), bt, boundaryType) - - END DO - - !Read and initialize volumes - DO e=1, self%numVols - READ(10,*) n, elemType - BACKSPACE(10) - - SELECT CASE(elemType) - CASE (2) - !Triangular element - READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3) - ALLOCATE(meshVol2DCartTria:: self%vols(e)%obj) - CALL self%vols(e)%obj%init(n - self%numEdges, p(1:3)) - - CASE (3) - !Quadrilateral element - READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4) - ALLOCATE(meshVol2DCartQuad:: self%vols(e)%obj) - CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4)) - - END SELECT - - END DO - - CLOSE(10) - - !Build connectivity between elements - DO e = 1, self%numVols - !Connectivity between volumes - DO et = 1, self%numVols - IF (e /= et) THEN - CALL connected(self%vols(e)%obj, self%vols(et)%obj) - - END IF - END DO - - !Connectivity between vols and edges - DO et = 1, self%numEdges - CALL connected(self%vols(e)%obj, self%edges(et)%obj) - - END DO - - !Constructs the global K matrix - CALL constructGlobalK(self%K, self%vols(e)%obj) - - END DO - - END SUBROUTINE readMesh2DCartGmsh - - !Selects type of elements to build connection - SUBROUTINE connectedVolVol(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshVol), INTENT(inout):: elemB - - SELECT TYPE(elemA) - TYPE IS(meshVol2DCartQuad) - !Element A is a quadrilateral - SELECT TYPE(elemB) - TYPE IS(meshVol2DCartQuad) - !Element B is a quadrilateral - CALL connectedQuadQuad(elemA, elemB) - - TYPE IS(meshVol2DCartTria) - !Element B is a triangle - CALL connectedQuadTria(elemA, elemB) - - END SELECT - - TYPE IS(meshVol2DCartTria) - !Element A is a Triangle - SELECT TYPE(elemB) - TYPE IS(meshVol2DCartQuad) - !Element B is a quadrilateral - CALL connectedQuadTria(elemB, elemA) - - TYPE IS(meshVol2DCartTria) - !Element B is a triangle - CALL connectedTriaTria(elemA, elemB) - - END SELECT - - END SELECT - - END SUBROUTINE connectedVolVol - - SUBROUTINE connectedVolEdge(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshEdge), INTENT(inout):: elemB - - SELECT TYPE(elemB) - CLASS IS(meshEdge2DCart) - SELECT TYPE(elemA) - TYPE IS(meshVol2DCartQuad) - !Element A is a quadrilateral - CALL connectedQuadEdge(elemA, elemB) - - TYPE IS(meshVol2DCartTria) - !Element A is a triangle - CALL connectedTriaEdge(elemA, elemB) - - END SELECT - - END SELECT - - END SUBROUTINE connectedVolEdge - - SUBROUTINE connectedQuadQuad(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA - CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemB - - !Check direction 1 - IF (.NOT. ASSOCIATED(elemA%e1) .AND. & - elemA%n1%n == elemB%n4%n .AND. & - elemA%n2%n == elemB%n3%n) THEN - elemA%e1 => elemB - elemB%e3 => elemA - - END IF - - !Check direction 2 - IF (.NOT. ASSOCIATED(elemA%e2) .AND. & - elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n4%n) THEN - elemA%e2 => elemB - elemB%e4 => elemA - - END IF - - !Check direction 3 - IF (.NOT. ASSOCIATED(elemA%e3) .AND. & - elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n1%n) THEN - elemA%e3 => elemB - elemB%e1 => elemA - - END IF - - !Check direction 4 - IF (.NOT. ASSOCIATED(elemA%e4) .AND. & - elemA%n4%n == elemB%n3%n .AND. & - elemA%n1%n == elemB%n2%n) THEN - elemA%e4 => elemB - elemB%e2 => elemA - - END IF - - END SUBROUTINE connectedQuadQuad - - SUBROUTINE connectedQuadTria(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA - CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemB - - !Check direction 1 - IF (.NOT. ASSOCIATED(elemA%e1)) THEN - IF (elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n3%n) THEN - elemA%e1 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n1%n == elemB%n3%n .AND. & - elemA%n2%n == elemB%n2%n) THEN - elemA%e1 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n1%n) THEN - elemA%e1 => elemB - elemB%e1 => elemA - - END IF - - END IF - - !Check direction 2 - IF (.NOT. ASSOCIATED(elemA%e2)) THEN - IF (elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n3%n) THEN - elemA%e2 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n2%n) THEN - elemA%e2 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n1%n) THEN - elemA%e2 => elemB - elemB%e1 => elemA - - END IF - - END IF - - !Check direction 3 - IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF (elemA%n3%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n3%n) THEN - elemA%e3 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n3%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n2%n) THEN - elemA%e3 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n1%n) THEN - elemA%e3 => elemB - elemB%e1 => elemA - - END IF - - END IF - - !Check direction 4 - IF (.NOT. ASSOCIATED(elemA%e4)) THEN - IF (elemA%n4%n == elemB%n1%n .AND. & - elemA%n1%n == elemB%n3%n) THEN - elemA%e4 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n4%n == elemB%n3%n .AND. & - elemA%n1%n == elemB%n2%n) THEN - elemA%e4 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n4%n == elemB%n2%n .AND. & - elemA%n1%n == elemB%n1%n) THEN - elemA%e4 => elemB - elemB%e1 => elemA - - END IF - - END IF - - END SUBROUTINE connectedQuadTria - - SUBROUTINE connectedTriaTria(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemA - CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemB - - !Check direction 1 - IF (.NOT. ASSOCIATED(elemA%e1)) THEN - IF (elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n3%n) THEN - elemA%e1 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n1%n) THEN - elemA%e1 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n1%n == elemB%n3%n .AND. & - elemA%n2%n == elemB%n2%n) THEN - elemA%e1 => elemB - elemB%e2 => elemA - - END IF - - END IF - - !Check direction 2 - IF (.NOT. ASSOCIATED(elemA%e2)) THEN - IF (elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n3%n) THEN - elemA%e2 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n1%n) THEN - elemA%e2 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n2%n) THEN - elemA%e2 => elemB - elemB%e2 => elemA - - END IF - - - END IF - - !Check direction 3 - IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF (elemA%n3%n == elemB%n1%n .AND. & - elemA%n1%n == elemB%n3%n) THEN - elemA%e3 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n3%n == elemB%n2%n .AND. & - elemA%n1%n == elemB%n1%n) THEN - elemA%e3 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n3%n == elemB%n3%n .AND. & - elemA%n1%n == elemB%n2%n) THEN - elemA%e3 => elemB - elemB%e2 => elemA - - END IF - - - END IF - - END SUBROUTINE connectedTriaTria - - SUBROUTINE connectedQuadEdge(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA - CLASS(meshEdge2DCart), INTENT(inout), TARGET:: elemB - - !Check direction 1 - IF (.NOT. ASSOCIATED(elemA%e1)) THEN - IF (elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n2%n) THEN - elemA%e1 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n1%n) THEN - elemA%e1 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal - - END IF - - END IF - - !Check direction 2 - IF (.NOT. ASSOCIATED(elemA%e2)) THEN - IF (elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n2%n) THEN - elemA%e2 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n1%n) THEN - elemA%e2 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal - - END IF - - END IF - - !Check direction 3 - IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF (elemA%n3%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n2%n) THEN - elemA%e3 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n1%n) THEN - elemA%e3 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal - - END IF - - END IF - - !Check direction 4 - IF (.NOT. ASSOCIATED(elemA%e4)) THEN - IF (elemA%n4%n == elemB%n1%n .AND. & - elemA%n1%n == elemB%n2%n) THEN - elemA%e4 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n4%n == elemB%n2%n .AND. & - elemA%n1%n == elemB%n1%n) THEN - elemA%e4 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal - - END IF - - END IF - - END SUBROUTINE connectedQuadEdge - - SUBROUTINE connectedTriaEdge(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemA - CLASS(meshEdge2DCart), INTENT(inout), TARGET:: elemB - - !Check direction 1 - IF (.NOT. ASSOCIATED(elemA%e1)) THEN - IF (elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n2%n) THEN - elemA%e1 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n1%n) THEN - elemA%e1 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal - - END IF - - END IF - - !Check direction 2 - IF (.NOT. ASSOCIATED(elemA%e2)) THEN - IF (elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n2%n) THEN - elemA%e2 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n1%n) THEN - elemA%e2 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal - - END IF - - END IF - - !Check direction 3 - IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF (elemA%n3%n == elemB%n1%n .AND. & - elemA%n1%n == elemB%n2%n) THEN - elemA%e3 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n3%n == elemB%n2%n .AND. & - elemA%n1%n == elemB%n1%n) THEN - elemA%e3 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal - - END IF - - END IF - - END SUBROUTINE connectedTriaEdge - - SUBROUTINE constructGlobalK(K, elem) - IMPLICIT NONE - - REAL(8), INTENT(inout):: K(1:,1:) - CLASS(meshVol), INTENT(in):: elem - REAL(8), ALLOCATABLE:: localK(:,:) - INTEGER:: nNodes, i, j - INTEGER, ALLOCATABLE:: n(:) - - SELECT TYPE(elem) - TYPE IS(meshVol2DCartQuad) - nNodes = 4 - ALLOCATE(localK(1:nNodes,1:nNodes)) - localK = elem%elemK() - ALLOCATE(n(1:nNodes)) - n = (/ elem%n1%n, elem%n2%n, & - elem%n3%n, elem%n4%n /) - - TYPE IS(meshVol2DCartTria) - nNodes = 3 - ALLOCATE(localK(1:nNodes,1:nNodes)) - localK = elem%elemK() - ALLOCATE(n(1:nNodes)) - n = (/ elem%n1%n, elem%n2%n, elem%n3%n /) - - CLASS DEFAULT - nNodes = 0 - ALLOCATE(localK(1:1, 1:1)) - localK = 0.D0 - ALLOCATE(n(1:1)) - n = 0 - - END SELECT - - DO i = 1, nNodes - DO j = 1, nNodes - K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j) - END DO - END DO - - END SUBROUTINE constructGlobalK - -END MODULE moduleMesh2DCartRead diff --git a/src/modules/mesh/2DCyl/makefile b/src/modules/mesh/2DCyl/makefile index ca2d7d1..e03a302 100644 --- a/src/modules/mesh/2DCyl/makefile +++ b/src/modules/mesh/2DCyl/makefile @@ -1,8 +1,5 @@ -all : moduleMesh2DCyl.o moduleMesh2DCylRead.o +all : moduleMesh2DCyl.o moduleMesh2DCyl.o: moduleMesh2DCyl.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleMesh2DCylRead.o: moduleMesh2DCyl.o moduleMesh2DCylRead.f90 - $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ - diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index d5a1432..362380d 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -427,18 +427,19 @@ MODULE moduleMesh2DCyl END FUNCTION randposVolQuad !Computes element local stiffness matrix - PURE FUNCTION elemKQuad(self) RESULT(ke) + PURE FUNCTION elemKQuad(self) RESULT(localK) USE moduleConstParam, ONLY: PI2 IMPLICIT NONE CLASS(meshVol2DCylQuad), INTENT(in):: self + REAL(8), ALLOCATABLE:: localK(:,:) REAL(8):: r, xi(1:3) REAL(8):: fPsi(1:4), dPsi(1:2,1:4) - REAL(8):: ke(1:4,1:4) REAL(8):: invJ(1:2,1:2), detJ INTEGER:: l, m - ke=0.D0 + ALLOCATE(localK(1:4, 1:4)) + localK=0.D0 xi=0.D0 !Start 2D Gauss Quad Integral DO l=1, 3 @@ -451,13 +452,13 @@ MODULE moduleMesh2DCyl detJ = self%detJac(xi,dPsi) invJ = self%invJac(xi,dPsi) r = DOT_PRODUCT(fPsi,self%r) - ke = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), & + localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), & MATMUL(invJ,dPsi))* & r*wQuad(l)*wQuad(m)/detJ END DO END DO - ke = ke*PI2 + localK = localK*PI2 END FUNCTION elemKQuad @@ -800,32 +801,33 @@ MODULE moduleMesh2DCyl END SUBROUTINE partialDerTria !Computes element local stiffness matrix - PURE FUNCTION elemKTria(self) RESULT(ke) + PURE FUNCTION elemKTria(self) RESULT(localK) USE moduleConstParam, ONLY: PI2 IMPLICIT NONE CLASS(meshVol2DCylTria), INTENT(in):: self + REAL(8), ALLOCATABLE:: localK(:,:) REAL(8):: r, xi(1:3) REAL(8):: fPsi(1:3), dPsi(1:2,1:3) - REAL(8):: ke(1:3,1:3) REAL(8):: invJ(1:2,1:2), detJ INTEGER:: l - ke=0.D0 + ALLOCATE(localK(1:4, 1:4)) + localK=0.D0 xi=0.D0 !Start 2D Gauss Quad Integral DO l=1, 4 - xi(1) = xi1Tria(l) - xi(2) = xi2Tria(l) - dPsi = self%dPsi(xi) - detJ = self%detJac(xi,dPsi) - invJ = self%invJac(xi,dPsi) - fPsi = self%fPsi(xi) - r = DOT_PRODUCT(fPsi,self%r) - ke = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ + xi(1) = xi1Tria(l) + xi(2) = xi2Tria(l) + dPsi = self%dPsi(xi) + detJ = self%detJac(xi,dPsi) + invJ = self%invJac(xi,dPsi) + fPsi = self%fPsi(xi) + r = DOT_PRODUCT(fPsi,self%r) + localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ END DO - ke = ke*PI2 + localK = localK*PI2 END FUNCTION elemKTria @@ -1047,4 +1049,449 @@ MODULE moduleMesh2DCyl END FUNCTION invJ2DCyl + SUBROUTINE connectMesh2DCyl(self) + IMPLICIT NONE + + CLASS(meshParticle), INTENT(inout):: self + INTEGER:: e, et + + DO e = 1, self%numVols + !Connect Vol-Vol + DO et = 1, self%numVols + IF (e /= et) THEN + CALL connectVolVol(self%vols(e)%obj, self%vols(et)%obj) + + END IF + + END DO + + !Connect Vol-Edge + DO et = 1, self%numEdges + CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) + + END DO + + END DO + + END SUBROUTINE connectMesh2DCyl + + !Selects type of elements to build connection + SUBROUTINE connectVolVol(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol), INTENT(inout):: elemA + CLASS(meshVol), INTENT(inout):: elemB + + SELECT TYPE(elemA) + TYPE IS(meshVol2DCylQuad) + !Element A is a quadrilateral + SELECT TYPE(elemB) + TYPE IS(meshVol2DCylQuad) + !Element B is a quadrilateral + CALL connectQuadQuad(elemA, elemB) + + TYPE IS(meshVol2DCylTria) + !Element B is a triangle + CALL connectQuadTria(elemA, elemB) + + END SELECT + + TYPE IS(meshVol2DCylTria) + !Element A is a Triangle + SELECT TYPE(elemB) + TYPE IS(meshVol2DCylQuad) + !Element B is a quadrilateral + CALL connectQuadTria(elemB, elemA) + + TYPE IS(meshVol2DCylTria) + !Element B is a triangle + CALL connectTriaTria(elemA, elemB) + + END SELECT + + END SELECT + + END SUBROUTINE connectVolVol + + SUBROUTINE connectVolEdge(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol), INTENT(inout):: elemA + CLASS(meshEdge), INTENT(inout):: elemB + + SELECT TYPE(elemB) + CLASS IS(meshEdge2DCyl) + SELECT TYPE(elemA) + TYPE IS(meshVol2DCylQuad) + !Element A is a quadrilateral + CALL connectQuadEdge(elemA, elemB) + + TYPE IS(meshVol2DCylTria) + !Element A is a triangle + CALL connectTriaEdge(elemA, elemB) + + END SELECT + + END SELECT + + END SUBROUTINE connectVolEdge + + SUBROUTINE connectQuadQuad(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA + CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemB + + !Check direction 1 + IF (.NOT. ASSOCIATED(elemA%e1) .AND. & + elemA%n1%n == elemB%n4%n .AND. & + elemA%n2%n == elemB%n3%n) THEN + elemA%e1 => elemB + elemB%e3 => elemA + + END IF + + !Check direction 2 + IF (.NOT. ASSOCIATED(elemA%e2) .AND. & + elemA%n2%n == elemB%n1%n .AND. & + elemA%n3%n == elemB%n4%n) THEN + elemA%e2 => elemB + elemB%e4 => elemA + + END IF + + !Check direction 3 + IF (.NOT. ASSOCIATED(elemA%e3) .AND. & + elemA%n3%n == elemB%n2%n .AND. & + elemA%n4%n == elemB%n1%n) THEN + elemA%e3 => elemB + elemB%e1 => elemA + + END IF + + !Check direction 4 + IF (.NOT. ASSOCIATED(elemA%e4) .AND. & + elemA%n4%n == elemB%n3%n .AND. & + elemA%n1%n == elemB%n2%n) THEN + elemA%e4 => elemB + elemB%e2 => elemA + + END IF + + END SUBROUTINE connectQuadQuad + + SUBROUTINE connectQuadTria(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA + CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemB + + !Check direction 1 + IF (.NOT. ASSOCIATED(elemA%e1)) THEN + IF (elemA%n1%n == elemB%n1%n .AND. & + elemA%n2%n == elemB%n3%n) THEN + elemA%e1 => elemB + elemB%e3 => elemA + + ELSEIF (elemA%n1%n == elemB%n3%n .AND. & + elemA%n2%n == elemB%n2%n) THEN + elemA%e1 => elemB + elemB%e2 => elemA + + ELSEIF (elemA%n1%n == elemB%n2%n .AND. & + elemA%n2%n == elemB%n1%n) THEN + elemA%e1 => elemB + elemB%e1 => elemA + + END IF + + END IF + + !Check direction 2 + IF (.NOT. ASSOCIATED(elemA%e2)) THEN + IF (elemA%n2%n == elemB%n1%n .AND. & + elemA%n3%n == elemB%n3%n) THEN + elemA%e2 => elemB + elemB%e3 => elemA + + ELSEIF (elemA%n2%n == elemB%n3%n .AND. & + elemA%n3%n == elemB%n2%n) THEN + elemA%e2 => elemB + elemB%e2 => elemA + + ELSEIF (elemA%n2%n == elemB%n2%n .AND. & + elemA%n3%n == elemB%n1%n) THEN + elemA%e2 => elemB + elemB%e1 => elemA + + END IF + + END IF + + !Check direction 3 + IF (.NOT. ASSOCIATED(elemA%e3)) THEN + IF (elemA%n3%n == elemB%n1%n .AND. & + elemA%n4%n == elemB%n3%n) THEN + elemA%e3 => elemB + elemB%e3 => elemA + + ELSEIF (elemA%n3%n == elemB%n3%n .AND. & + elemA%n4%n == elemB%n2%n) THEN + elemA%e3 => elemB + elemB%e2 => elemA + + ELSEIF (elemA%n3%n == elemB%n2%n .AND. & + elemA%n4%n == elemB%n1%n) THEN + elemA%e3 => elemB + elemB%e1 => elemA + + END IF + + END IF + + !Check direction 4 + IF (.NOT. ASSOCIATED(elemA%e4)) THEN + IF (elemA%n4%n == elemB%n1%n .AND. & + elemA%n1%n == elemB%n3%n) THEN + elemA%e4 => elemB + elemB%e3 => elemA + + ELSEIF (elemA%n4%n == elemB%n3%n .AND. & + elemA%n1%n == elemB%n2%n) THEN + elemA%e4 => elemB + elemB%e2 => elemA + + ELSEIF (elemA%n4%n == elemB%n2%n .AND. & + elemA%n1%n == elemB%n1%n) THEN + elemA%e4 => elemB + elemB%e1 => elemA + + END IF + + END IF + + END SUBROUTINE connectQuadTria + + SUBROUTINE connectTriaTria(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemA + CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemB + + !Check direction 1 + IF (.NOT. ASSOCIATED(elemA%e1)) THEN + IF (elemA%n1%n == elemB%n1%n .AND. & + elemA%n2%n == elemB%n3%n) THEN + elemA%e1 => elemB + elemB%e3 => elemA + + ELSEIF (elemA%n1%n == elemB%n2%n .AND. & + elemA%n2%n == elemB%n1%n) THEN + elemA%e1 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n1%n == elemB%n3%n .AND. & + elemA%n2%n == elemB%n2%n) THEN + elemA%e1 => elemB + elemB%e2 => elemA + + END IF + + END IF + + !Check direction 2 + IF (.NOT. ASSOCIATED(elemA%e2)) THEN + IF (elemA%n2%n == elemB%n1%n .AND. & + elemA%n3%n == elemB%n3%n) THEN + elemA%e2 => elemB + elemB%e3 => elemA + + ELSEIF (elemA%n2%n == elemB%n2%n .AND. & + elemA%n3%n == elemB%n1%n) THEN + elemA%e2 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n2%n == elemB%n3%n .AND. & + elemA%n3%n == elemB%n2%n) THEN + elemA%e2 => elemB + elemB%e2 => elemA + + END IF + + + END IF + + !Check direction 3 + IF (.NOT. ASSOCIATED(elemA%e3)) THEN + IF (elemA%n3%n == elemB%n1%n .AND. & + elemA%n1%n == elemB%n3%n) THEN + elemA%e3 => elemB + elemB%e3 => elemA + + ELSEIF (elemA%n3%n == elemB%n2%n .AND. & + elemA%n1%n == elemB%n1%n) THEN + elemA%e3 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n3%n == elemB%n3%n .AND. & + elemA%n1%n == elemB%n2%n) THEN + elemA%e3 => elemB + elemB%e2 => elemA + + END IF + + + END IF + + END SUBROUTINE connectTriaTria + + SUBROUTINE connectQuadEdge(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA + CLASS(meshEdge2DCyl), INTENT(inout), TARGET:: elemB + + !Check direction 1 + IF (.NOT. ASSOCIATED(elemA%e1)) THEN + IF (elemA%n1%n == elemB%n1%n .AND. & + elemA%n2%n == elemB%n2%n) THEN + elemA%e1 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n1%n == elemB%n2%n .AND. & + elemA%n2%n == elemB%n1%n) THEN + elemA%e1 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + + END IF + + END IF + + !Check direction 2 + IF (.NOT. ASSOCIATED(elemA%e2)) THEN + IF (elemA%n2%n == elemB%n1%n .AND. & + elemA%n3%n == elemB%n2%n) THEN + elemA%e2 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n2%n == elemB%n2%n .AND. & + elemA%n3%n == elemB%n1%n) THEN + elemA%e2 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + + END IF + + END IF + + !Check direction 3 + IF (.NOT. ASSOCIATED(elemA%e3)) THEN + IF (elemA%n3%n == elemB%n1%n .AND. & + elemA%n4%n == elemB%n2%n) THEN + elemA%e3 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n3%n == elemB%n2%n .AND. & + elemA%n4%n == elemB%n1%n) THEN + elemA%e3 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + + END IF + + END IF + + !Check direction 4 + IF (.NOT. ASSOCIATED(elemA%e4)) THEN + IF (elemA%n4%n == elemB%n1%n .AND. & + elemA%n1%n == elemB%n2%n) THEN + elemA%e4 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n4%n == elemB%n2%n .AND. & + elemA%n1%n == elemB%n1%n) THEN + elemA%e4 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + + END IF + + END IF + + END SUBROUTINE connectQuadEdge + + SUBROUTINE connectTriaEdge(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemA + CLASS(meshEdge2DCyl), INTENT(inout), TARGET:: elemB + + !Check direction 1 + IF (.NOT. ASSOCIATED(elemA%e1)) THEN + IF (elemA%n1%n == elemB%n1%n .AND. & + elemA%n2%n == elemB%n2%n) THEN + elemA%e1 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n1%n == elemB%n2%n .AND. & + elemA%n2%n == elemB%n1%n) THEN + elemA%e1 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + + END IF + + END IF + + !Check direction 2 + IF (.NOT. ASSOCIATED(elemA%e2)) THEN + IF (elemA%n2%n == elemB%n1%n .AND. & + elemA%n3%n == elemB%n2%n) THEN + elemA%e2 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n2%n == elemB%n2%n .AND. & + elemA%n3%n == elemB%n1%n) THEN + elemA%e2 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + + END IF + + END IF + + !Check direction 3 + IF (.NOT. ASSOCIATED(elemA%e3)) THEN + IF (elemA%n3%n == elemB%n1%n .AND. & + elemA%n1%n == elemB%n2%n) THEN + elemA%e3 => elemB + elemB%e1 => elemA + + ELSEIF (elemA%n3%n == elemB%n2%n .AND. & + elemA%n1%n == elemB%n1%n) THEN + elemA%e3 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + + END IF + + END IF + + END SUBROUTINE connectTriaEdge + END MODULE moduleMesh2DCyl diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 deleted file mode 100644 index cbe3faf..0000000 --- a/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 +++ /dev/null @@ -1,643 +0,0 @@ -MODULE moduleMesh2DCylRead - USE moduleMesh - USE moduleMesh2DCyl - - TYPE, EXTENDS(meshGeneric):: mesh2DCylGeneric - CONTAINS - PROCEDURE, PASS:: init => init2DCylMesh - PROCEDURE, PASS:: readMesh => readMesh2DCylGmsh - - END TYPE - - INTERFACE connected - MODULE PROCEDURE connectedVolVol, connectedVolEdge - - END INTERFACE connected - - CONTAINS - !Init mesh - SUBROUTINE init2DCylMesh(self, meshFormat) - USE moduleMesh - USE moduleErrors - IMPLICIT NONE - - CLASS(mesh2DCylGeneric), INTENT(out):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat - - SELECT CASE(meshFormat) - CASE ("gmsh") - self%printOutput => printOutputGmsh - self%printColl => printCollGmsh - self%printEM => printEMGmsh - - CASE DEFAULT - CALL criticalError("Mesh type " // meshFormat // " not supported.", "init2DCylMesh") - - END SELECT - - END SUBROUTINE init2DCylMesh - - !Read mesh from gmsh file - SUBROUTINE readMesh2DCylGmsh(self, filename) - USE moduleBoundary - IMPLICIT NONE - - CLASS(mesh2DCylGeneric), INTENT(inout):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: filename - REAL(8):: r, z - INTEGER:: p(1:4) - INTEGER:: e=0, et=0, n=0, eTemp=0, elemType=0, bt = 0 - INTEGER:: totalNumElem - INTEGER:: boundaryType - - !Read msh - OPEN(10, FILE=TRIM(filename)) - !Skip header - READ(10, *) - READ(10, *) - READ(10, *) - READ(10, *) - !Read number of nodes - READ(10, *) self%numNodes - !Allocate required matrices and vectors - ALLOCATE(self%nodes(1:self%numNodes)) - 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 - !Read nodes cartesian coordinates (x=z, y=r, z=null) - DO e=1, self%numNodes - READ(10, *) n, z, r - ALLOCATE(meshNode2DCyl:: self%nodes(n)%obj) - CALL self%nodes(n)%obj%init(n, (/z, r, 0.D0 /)) - - END DO - !Skips comments - READ(10, *) - READ(10, *) - !Reads Totalnumber of elements - READ(10, *) TotalnumElem - !counts edges and volume elements - self%numEdges = 0 - DO e=1, TotalnumElem - READ(10,*) eTemp, elemType - IF (elemType==1) THEN - self%numEdges=e - END IF - END DO - !Substract the number of edges to the total number of elements - !to obtain the number of volume elements - self%numVols = TotalnumElem - self%numEdges - !Allocates arrays - ALLOCATE(self%edges(1:self%numEdges)) - ALLOCATE(self%vols(1:self%numVols)) - - !Go back to the beggining to read elements - DO e=1, totalNumElem - BACKSPACE(10) - END DO - - !Reads edges - DO e=1, self%numEdges - READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2) - !Associate boundary condition procedure. - bt = getBoundaryId(boundaryType) - - ALLOCATE(meshEdge2DCyl:: self%edges(e)%obj) - - CALL self%edges(e)%obj%init(n, p(1:2), bt, boundaryType) - - END DO - - !Read and initialize volumes - DO e=1, self%numVols - READ(10,*) n, elemType - BACKSPACE(10) - - SELECT CASE(elemType) - CASE (2) - !Triangular element - READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3) - ALLOCATE(meshVol2DCylTria:: self%vols(e)%obj) - CALL self%vols(e)%obj%init(n - self%numEdges, p(1:3)) - - CASE (3) - !Quadrilateral element - READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4) - ALLOCATE(meshVol2DCylQuad:: self%vols(e)%obj) - CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4)) - - END SELECT - - END DO - - CLOSE(10) - - !Build connectivity between elements - DO e = 1, self%numVols - !Connectivity between volumes - DO et = 1, self%numVols - IF (e /= et) THEN - CALL connected(self%vols(e)%obj, self%vols(et)%obj) - - END IF - END DO - - !Connectivity between vols and edges - DO et = 1, self%numEdges - CALL connected(self%vols(e)%obj, self%edges(et)%obj) - - END DO - - !Constructs the global K matrix - CALL constructGlobalK(self%K, self%vols(e)%obj) - - END DO - - END SUBROUTINE readMesh2DCylGmsh - - !Selects type of elements to build connection - SUBROUTINE connectedVolVol(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshVol), INTENT(inout):: elemB - - SELECT TYPE(elemA) - TYPE IS(meshVol2DCylQuad) - !Element A is a quadrilateral - SELECT TYPE(elemB) - TYPE IS(meshVol2DCylQuad) - !Element B is a quadrilateral - CALL connectedQuadQuad(elemA, elemB) - - TYPE IS(meshVol2DCylTria) - !Element B is a triangle - CALL connectedQuadTria(elemA, elemB) - - END SELECT - - TYPE IS(meshVol2DCylTria) - !Element A is a Triangle - SELECT TYPE(elemB) - TYPE IS(meshVol2DCylQuad) - !Element B is a quadrilateral - CALL connectedQuadTria(elemB, elemA) - - TYPE IS(meshVol2DCylTria) - !Element B is a triangle - CALL connectedTriaTria(elemA, elemB) - - END SELECT - - END SELECT - - END SUBROUTINE connectedVolVol - - SUBROUTINE connectedVolEdge(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshEdge), INTENT(inout):: elemB - - SELECT TYPE(elemB) - CLASS IS(meshEdge2DCyl) - SELECT TYPE(elemA) - TYPE IS(meshVol2DCylQuad) - !Element A is a quadrilateral - CALL connectedQuadEdge(elemA, elemB) - - TYPE IS(meshVol2DCylTria) - !Element A is a triangle - CALL connectedTriaEdge(elemA, elemB) - - END SELECT - - END SELECT - - END SUBROUTINE connectedVolEdge - - PURE FUNCTION coincidentNodes(nodesA, nodesB) RESULT(coincident) - IMPLICIT NONE - - INTEGER, DIMENSION(1:2), INTENT(in):: nodesA, nodesB - LOGICAL:: coincident - INTEGER:: i - - coincident = .FALSE. - DO i = 1, 2 - IF (ANY(nodesA(i) == nodesB)) THEN - coincident = .TRUE. - - ELSE - coincident = .FALSE. - EXIT - - END IF - - END DO - - END FUNCTION coincidentNodes - - - SUBROUTINE connectedQuadQuad(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA - CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemB - - !Check direction 1 - IF (.NOT. ASSOCIATED(elemA%e1) .AND. & - elemA%n1%n == elemB%n4%n .AND. & - elemA%n2%n == elemB%n3%n) THEN - elemA%e1 => elemB - elemB%e3 => elemA - - END IF - - !Check direction 2 - IF (.NOT. ASSOCIATED(elemA%e2) .AND. & - elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n4%n) THEN - elemA%e2 => elemB - elemB%e4 => elemA - - END IF - - !Check direction 3 - IF (.NOT. ASSOCIATED(elemA%e3) .AND. & - elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n1%n) THEN - elemA%e3 => elemB - elemB%e1 => elemA - - END IF - - !Check direction 4 - IF (.NOT. ASSOCIATED(elemA%e4) .AND. & - elemA%n4%n == elemB%n3%n .AND. & - elemA%n1%n == elemB%n2%n) THEN - elemA%e4 => elemB - elemB%e2 => elemA - - END IF - - END SUBROUTINE connectedQuadQuad - - SUBROUTINE connectedQuadTria(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA - CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemB - - !Check direction 1 - IF (.NOT. ASSOCIATED(elemA%e1)) THEN - IF (elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n3%n) THEN - elemA%e1 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n1%n == elemB%n3%n .AND. & - elemA%n2%n == elemB%n2%n) THEN - elemA%e1 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n1%n) THEN - elemA%e1 => elemB - elemB%e1 => elemA - - END IF - - END IF - - !Check direction 2 - IF (.NOT. ASSOCIATED(elemA%e2)) THEN - IF (elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n3%n) THEN - elemA%e2 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n2%n) THEN - elemA%e2 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n1%n) THEN - elemA%e2 => elemB - elemB%e1 => elemA - - END IF - - END IF - - !Check direction 3 - IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF (elemA%n3%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n3%n) THEN - elemA%e3 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n3%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n2%n) THEN - elemA%e3 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n1%n) THEN - elemA%e3 => elemB - elemB%e1 => elemA - - END IF - - END IF - - !Check direction 4 - IF (.NOT. ASSOCIATED(elemA%e4)) THEN - IF (elemA%n4%n == elemB%n1%n .AND. & - elemA%n1%n == elemB%n3%n) THEN - elemA%e4 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n4%n == elemB%n3%n .AND. & - elemA%n1%n == elemB%n2%n) THEN - elemA%e4 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n4%n == elemB%n2%n .AND. & - elemA%n1%n == elemB%n1%n) THEN - elemA%e4 => elemB - elemB%e1 => elemA - - END IF - - END IF - - END SUBROUTINE connectedQuadTria - - SUBROUTINE connectedTriaTria(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemA - CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemB - - !Check direction 1 - IF (.NOT. ASSOCIATED(elemA%e1)) THEN - IF (elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n3%n) THEN - elemA%e1 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n1%n) THEN - elemA%e1 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n1%n == elemB%n3%n .AND. & - elemA%n2%n == elemB%n2%n) THEN - elemA%e1 => elemB - elemB%e2 => elemA - - END IF - - END IF - - !Check direction 2 - IF (.NOT. ASSOCIATED(elemA%e2)) THEN - IF (elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n3%n) THEN - elemA%e2 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n1%n) THEN - elemA%e2 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n2%n) THEN - elemA%e2 => elemB - elemB%e2 => elemA - - END IF - - - END IF - - !Check direction 3 - IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF (elemA%n3%n == elemB%n1%n .AND. & - elemA%n1%n == elemB%n3%n) THEN - elemA%e3 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n3%n == elemB%n2%n .AND. & - elemA%n1%n == elemB%n1%n) THEN - elemA%e3 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n3%n == elemB%n3%n .AND. & - elemA%n1%n == elemB%n2%n) THEN - elemA%e3 => elemB - elemB%e2 => elemA - - END IF - - - END IF - - END SUBROUTINE connectedTriaTria - - SUBROUTINE connectedQuadEdge(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA - CLASS(meshEdge2DCyl), INTENT(inout), TARGET:: elemB - - !Check direction 1 - IF (.NOT. ASSOCIATED(elemA%e1)) THEN - IF (elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n2%n) THEN - elemA%e1 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n1%n) THEN - elemA%e1 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal - - END IF - - END IF - - !Check direction 2 - IF (.NOT. ASSOCIATED(elemA%e2)) THEN - IF (elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n2%n) THEN - elemA%e2 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n1%n) THEN - elemA%e2 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal - - END IF - - END IF - - !Check direction 3 - IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF (elemA%n3%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n2%n) THEN - elemA%e3 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n1%n) THEN - elemA%e3 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal - - END IF - - END IF - - !Check direction 4 - IF (.NOT. ASSOCIATED(elemA%e4)) THEN - IF (elemA%n4%n == elemB%n1%n .AND. & - elemA%n1%n == elemB%n2%n) THEN - elemA%e4 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n4%n == elemB%n2%n .AND. & - elemA%n1%n == elemB%n1%n) THEN - elemA%e4 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal - - END IF - - END IF - - END SUBROUTINE connectedQuadEdge - - SUBROUTINE connectedTriaEdge(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemA - CLASS(meshEdge2DCyl), INTENT(inout), TARGET:: elemB - - !Check direction 1 - IF (.NOT. ASSOCIATED(elemA%e1)) THEN - IF (elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n2%n) THEN - elemA%e1 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n1%n) THEN - elemA%e1 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal - - END IF - - END IF - - !Check direction 2 - IF (.NOT. ASSOCIATED(elemA%e2)) THEN - IF (elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n2%n) THEN - elemA%e2 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n1%n) THEN - elemA%e2 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal - - END IF - - END IF - - !Check direction 3 - IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF (elemA%n3%n == elemB%n1%n .AND. & - elemA%n1%n == elemB%n2%n) THEN - elemA%e3 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n3%n == elemB%n2%n .AND. & - elemA%n1%n == elemB%n1%n) THEN - elemA%e3 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal - - END IF - - END IF - - END SUBROUTINE connectedTriaEdge - - SUBROUTINE constructGlobalK(K, elem) - IMPLICIT NONE - - REAL(8), INTENT(inout):: K(1:,1:) - CLASS(meshVol), INTENT(in):: elem - REAL(8), ALLOCATABLE:: localK(:,:) - INTEGER:: nNodes, i, j - INTEGER, ALLOCATABLE:: n(:) - - SELECT TYPE(elem) - TYPE IS(meshVol2DCylQuad) - nNodes = 4 - ALLOCATE(localK(1:nNodes,1:nNodes)) - localK = elem%elemK() - ALLOCATE(n(1:nNodes)) - n = (/ elem%n1%n, elem%n2%n, & - elem%n3%n, elem%n4%n /) - - TYPE IS(meshVol2DCylTria) - nNodes = 3 - ALLOCATE(localK(1:nNodes,1:nNodes)) - localK = elem%elemK() - ALLOCATE(n(1:nNodes)) - n = (/ elem%n1%n, elem%n2%n, elem%n3%n /) - - CLASS DEFAULT - nNodes = 0 - ALLOCATE(localK(1:1, 1:1)) - localK = 0.D0 - ALLOCATE(n(1:1)) - n = 0 - - END SELECT - - DO i = 1, nNodes - DO j = 1, nNodes - K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j) - END DO - END DO - - END SUBROUTINE constructGlobalK - -END MODULE moduleMesh2DCylRead diff --git a/src/modules/mesh/3DCart/makefile b/src/modules/mesh/3DCart/makefile index 0c77a10..0f23c79 100644 --- a/src/modules/mesh/3DCart/makefile +++ b/src/modules/mesh/3DCart/makefile @@ -1,8 +1,5 @@ -all : moduleMesh3DCart.o moduleMesh3DCartRead.o +all : moduleMesh3DCart.o moduleMesh3DCart.o: moduleMesh3DCart.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleMesh3DCartRead.o: moduleMesh3DCart.o moduleMesh3DCartRead.f90 - $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ - diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 494459a..3e3b970 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -417,23 +417,25 @@ MODULE moduleMesh3DCart END SUBROUTINE partialDerTetra - PURE FUNCTION elemKTetra(self) RESULT(ke) + PURE FUNCTION elemKTetra(self) RESULT(localK) IMPLICIT NONE CLASS(meshVol3DCartTetra), INTENT(in):: self + REAL(8), ALLOCATABLE:: localK(:,:) REAL(8):: xii(1:3) REAL(8):: fPsi(1:4), dPsi(1:3, 1:4) - REAL(8):: ke(1:4,1:4) REAL(8):: invJ(1:3,1:3), detJ + ALLOCATE(localK(1:4,1:4)) + localK = 0.D0 + xii = 0.D0 !TODO: One point Gauss integral. Upgrade when possible - ke = 0.D0 - xii = (/ 0.25D0, 0.25D0, 0.25D0 /) - dPsi = self%dPsi(xii) - detJ = self%detJac(xii, dPsi) - invJ = self%invJac(xii, dPsi) - fPsi = self%fPsi(xii) - ke = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*1.D0/detJ + xii = (/ 0.25D0, 0.25D0, 0.25D0 /) + dPsi = self%dPsi(xii) + detJ = self%detJac(xii, dPsi) + invJ = self%invJac(xii, dPsi) + fPsi = self%fPsi(xii) + localK = MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*1.D0/detJ END FUNCTION elemKTetra @@ -456,7 +458,7 @@ MODULE moduleMesh3DCart detJ = self%detJac(xii, dPsi) fPsi = self%fPsi(xii) f = DOT_PRODUCT(fPsi, source) - localF = localF + f*fPsi*1.D0*detJ + localF = f*fPsi*1.D0*detJ END FUNCTION elemFTetra @@ -662,5 +664,369 @@ MODULE moduleMesh3DCart END FUNCTION invJ3DCart + !Selects type of elements to build connection + SUBROUTINE connectVolVol(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol), INTENT(inout):: elemA + CLASS(meshVol), INTENT(inout):: elemB + + SELECT TYPE(elemA) + TYPE IS(meshVol3DCartTetra) + !Element A is a tetrahedron + SELECT TYPE(elemB) + TYPE IS(meshVol3DCartTetra) + !Element B is a tetrahedron + CALL connectTetraTetra(elemA, elemB) + + END SELECT + + END SELECT + + END SUBROUTINE connectVolVol + + SUBROUTINE connectVolEdge(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol), INTENT(inout):: elemA + CLASS(meshEdge), INTENT(inout):: elemB + + SELECT TYPE(elemB) + CLASS IS(meshEdge3DCartTria) + SELECT TYPE(elemA) + TYPE IS(meshVol3DCartTetra) + !Element A is a tetrahedron + CALL connectTetraEdge(elemA, elemB) + + END SELECT + + END SELECT + + END SUBROUTINE connectVolEdge + + SUBROUTINE connectMesh3DCart(self) + IMPLICIT NONE + + CLASS(meshParticle), INTENT(inout):: self + INTEGER:: e, et + + DO e = 1, self%numVols + !Connect Vol-Vol + DO et = 1, self%numVols + IF (e /= et) THEN + CALL connectVolVol(self%vols(e)%obj, self%vols(et)%obj) + + END IF + + END DO + + !Connect Vol-Edge + DO et = 1, self%numEdges + CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) + + END DO + + END DO + + END SUBROUTINE connectMesh3DCart + + !Checks if two sets of nodes are coincidend in any order + PURE FUNCTION coincidentNodes(nodesA, nodesB) RESULT(coincident) + IMPLICIT NONE + + INTEGER, DIMENSION(1:3), INTENT(in):: nodesA, nodesB + LOGICAL:: coincident + INTEGER:: i + + coincident = .FALSE. + DO i = 1, 3 + IF (ANY(nodesA(i) == nodesB)) THEN + coincident = .TRUE. + + ELSE + coincident = .FALSE. + EXIT + + END IF + + END DO + + END FUNCTION coincidentNodes + + SUBROUTINE connectTetraTetra(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA + CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemB + + !Check surface 1 + IF (.NOT. ASSOCIATED(elemA%e1)) THEN + IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN + + elemA%e1 => elemB + elemB%e1 => elemA + + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), & + (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN + + elemA%e1 => elemB + elemB%e2 => elemA + + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN + + elemA%e1 => elemB + elemB%e3 => elemA + + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), & + (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN + + elemA%e1 => elemB + elemB%e4 => elemA + + END IF + + END IF + + !Check surface 2 + IF (.NOT. ASSOCIATED(elemA%e2)) THEN + IF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN + + elemA%e2 => elemB + elemB%e1 => elemA + + ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN + + elemA%e2 => elemB + elemB%e2 => elemA + + ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN + + elemA%e2 => elemB + elemB%e3 => elemA + + ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN + + elemA%e2 => elemB + elemB%e4 => elemA + + END IF + + END IF + + !Check surface 3 + IF (.NOT. ASSOCIATED(elemA%e3)) THEN + IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN + + elemA%e3 => elemB + elemB%e1 => elemA + + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN + + elemA%e3 => elemB + elemB%e2 => elemA + + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN + + elemA%e3 => elemB + elemB%e3 => elemA + + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN + + elemA%e3 => elemB + elemB%e4 => elemA + + END IF + + END IF + + !Check surface 4 + IF (.NOT. ASSOCIATED(elemA%e4)) THEN + IF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN + + elemA%e4 => elemB + elemB%e1 => elemA + + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN + + elemA%e4 => elemB + elemB%e2 => elemA + + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN + + elemA%e4 => elemB + elemB%e3 => elemA + + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN + + elemA%e4 => elemB + elemB%e4 => elemA + + END IF + + END IF + + END SUBROUTINE connectTetraTetra + + SUBROUTINE connectTetraEdge(elemA, elemB) + USE moduleMath + IMPLICIT NONE + + CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA + CLASS(meshEdge3DCartTria), INTENT(inout), TARGET:: elemB + INTEGER:: nodesEdge(1:3) + REAL(8), DIMENSION(1:3):: vec1, vec2 + REAL(8):: normVol(1:3) + + nodesEdge = (/ elemB%n1%n, elemB%n2%n, elemB%n3%n /) + + !Check surface 1 + IF (.NOT. ASSOCIATED(elemA%e1)) THEN + IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elema%n3%n/), & + nodesEdge)) THEN + + vec1 = (/ elemA%x(2) - elemA%x(1), & + elemA%y(2) - elemA%y(1), & + elemA%z(2) - elemA%z(1) /) + vec2 = (/ elemA%x(3) - elemA%x(1), & + elemA%y(3) - elemA%y(1), & + elemA%z(3) - elemA%z(1) /) + normVol = crossProduct(vec1, vec2) + normVol = normalize(normVol) + + IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN + + elemA%e1 => elemB + elemB%e1 => elemA + + ELSE + + elemA%e1 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + + END IF + + END IF + + END IF + + !Check surface 2 + IF (.NOT. ASSOCIATED(elemA%e2)) THEN + IF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + nodesEdge)) THEN + + vec1 = (/ elemA%x(3) - elemA%x(2), & + elemA%y(3) - elemA%y(2), & + elemA%z(3) - elemA%z(2) /) + vec2 = (/ elemA%x(4) - elemA%x(2), & + elemA%y(4) - elemA%y(2), & + elemA%z(4) - elemA%z(2) /) + normVol = crossProduct(vec1, vec2) + normVol = normalize(normVol) + + IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN + + elemA%e2 => elemB + elemB%e1 => elemA + + ELSE + + elemA%e2 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + + END IF + + END IF + + END IF + + !Check surface 3 + IF (.NOT. ASSOCIATED(elemA%e3)) THEN + IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elema%n4%n/), & + nodesEdge)) THEN + + vec1 = (/ elemA%x(2) - elemA%x(1), & + elemA%y(2) - elemA%y(1), & + elemA%z(2) - elemA%z(1) /) + vec2 = (/ elemA%x(4) - elemA%x(1), & + elemA%y(4) - elemA%y(1), & + elemA%z(4) - elemA%z(1) /) + normVol = crossProduct(vec1, vec2) + normVol = normalize(normVol) + + IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN + + elemA%e3 => elemB + elemB%e1 => elemA + + ELSE + + elemA%e3 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + + END IF + + END IF + + END IF + + !Check surface 4 + IF (.NOT. ASSOCIATED(elemA%e4)) THEN + IF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elema%n4%n/), & + nodesEdge)) THEN + + vec1 = (/ elemA%x(3) - elemA%x(1), & + elemA%y(3) - elemA%y(1), & + elemA%z(3) - elemA%z(1) /) + vec2 = (/ elemA%x(4) - elemA%x(1), & + elemA%y(4) - elemA%y(1), & + elemA%z(4) - elemA%z(1) /) + normVol = crossProduct(vec1, vec2) + normVol = normalize(normVol) + + IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN + + elemA%e4 => elemB + elemB%e1 => elemA + + + ELSE + + elemA%e4 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + + END IF + + END IF + + END IF + + END SUBROUTINE connectTetraEdge + END MODULE moduleMesh3DCart diff --git a/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 b/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 deleted file mode 100644 index 01aef8e..0000000 --- a/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 +++ /dev/null @@ -1,544 +0,0 @@ -MODULE moduleMesh3DCartRead - USE moduleMesh - USE moduleMesh3DCart - - TYPE, EXTENDS(meshGeneric):: mesh3DCartGeneric - CONTAINS - PROCEDURE, PASS:: init => init3DCartMesh - PROCEDURE, PASS:: readMesh => readMesh3DCartGmsh - - END TYPE - - INTERFACE connected - MODULE PROCEDURE connectedVolVol, connectedVolEdge - - END INTERFACE connected - - CONTAINS - !Init mesh - SUBROUTINE init3DCartMesh(self, meshFormat) - USE moduleMesh - USE moduleErrors - IMPLICIT NONE - - CLASS(mesh3DCartGeneric), INTENT(out):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat - - SELECT CASE(meshFormat) - CASE ("gmsh") - self%printOutput => printOutputGmsh - self%printColl => printCollGmsh - self%printEM => printEMGmsh - - CASE DEFAULT - CALL criticalError("Mesh type " // meshFormat // " not supported.", "init3DCartMesh") - - END SELECT - - END SUBROUTINE init3DCartMesh - - !Read mesh from gmsh file - SUBROUTINE readMesh3DCartGmsh(self, filename) - USE moduleBoundary - IMPLICIT NONE - - CLASS(mesh3DCartGeneric), INTENT(inout):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: filename - REAL(8):: x, y, z - INTEGER:: p(1:4) - INTEGER:: e = 0, et = 0, n = 0, eTemp = 0, elemType = 0, bt = 0 - INTEGER:: totalNumElem - INTEGER:: boundaryType - - !Read mesh - OPEN(10, FILE=TRIM(filename)) - !Skip header - READ(10, *) - READ(10, *) - READ(10, *) - READ(10, *) - !Read number of nodes - READ(10, *) self%numNodes - !Allocate required matrices and vectors - ALLOCATE(self%nodes(1:self%numNodes)) - 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 - - !Read node cartesian coordinates (x = x, y = y, z = z) - DO e = 1, self%numNodes - READ(10, *) n, x, y, z - ALLOCATE(meshNode3Dcart::self%nodes(n)%obj) - CALL self%nodes(n)%obj%init(n, (/x, y, z /)) - - END DO - - !Skip comments - READ(10, *) - READ(10, *) - - !Reads total number of elements - READ(10, *) totalNumElem - !conts edges and volume elements - self%numEdges = 0 - DO e = 1, totalNumElem - READ(10, *) eTemp, elemType - IF (elemType == 2) THEN - self%numEdges = e - - END IF - - END DO - - !Substract the number of edges to the total number of elements to obtain the number - !of volume elements - self%numVols = totalNumElem - self%numEdges - - !Allocate required arrays - ALLOCATE(self%edges(1:self%numEdges)) - ALLOCATE(self%vols(1:self%numVols)) - - !Go back to the beggining to read each specific element - DO e = 1, totalNumElem - BACKSPACE(10) - - END DO - - !Reads surfaces - DO e = 1, self%numEdges - READ(10, *) n, elemType - BACKSPACE(10) - - SELECT CASE(elemType) - CASE(2) - !Triangular surface - READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1:3) - bt = getBoundaryID(boundaryType) - - ALLOCATE(meshEdge3DCartTria:: self%edges(e)%obj) - - CALL self%edges(e)%obj%init(n, p(1:3), bt, boundaryType) - - END SELECT - - END DO - - !Read and initialize volumes - DO e = 1, self%numVols - READ(10, *) n, elemType - BACKSPACE(10) - - SELECT CASE(elemType) - CASE(4) - !Tetrahedron element - READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:4) - ALLOCATE(meshVol3DCartTetra:: self%vols(e)%obj) - CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4)) - - END SELECT - - END DO - - CLOSE(10) - - !Build connectivy between elements - DO e = 1, self%numVols - !Connectivity between volumes - DO et = 1, self%numVols - IF (e /= et) THEN - CALL connected(self%vols(e)%obj, self%vols(et)%obj) - - END IF - - END DO - - !Connectivity between vols and surfaces - DO et = 1, self%numEdges - CALL connected(self%vols(e)%obj, self%edges(et)%obj) - - END DO - - !Constructs the global K matrix - CALL constructGlobalK(self%K, self%vols(e)%obj) - - END DO - - END SUBROUTINE readMesh3DCartGmsh - - !Selects type of elements to build connection - SUBROUTINE connectedVolVol(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshVol), INTENT(inout):: elemB - - SELECT TYPE(elemA) - TYPE IS(meshVol3DCartTetra) - !Element A is a tetrahedron - SELECT TYPE(elemB) - TYPE IS(meshVol3DCartTetra) - !Element B is a tetrahedron - CALL connectedTetraTetra(elemA, elemB) - - END SELECT - - END SELECT - - END SUBROUTINE connectedVolVol - - - SUBROUTINE connectedVolEdge(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshEdge), INTENT(inout):: elemB - - SELECT TYPE(elemB) - CLASS IS(meshEdge3DCartTria) - SELECT TYPE(elemA) - TYPE IS(meshVol3DCartTetra) - !Element A is a tetrahedron - CALL connectedTetraEdge(elemA, elemB) - - END SELECT - - END SELECT - - END SUBROUTINE connectedVolEdge - - PURE FUNCTION coincidentNodes(nodesA, nodesB) RESULT(coincident) - IMPLICIT NONE - - INTEGER, DIMENSION(1:3), INTENT(in):: nodesA, nodesB - LOGICAL:: coincident - INTEGER:: i - - coincident = .FALSE. - DO i = 1, 3 - IF (ANY(nodesA(i) == nodesB)) THEN - coincident = .TRUE. - - ELSE - coincident = .FALSE. - EXIT - - END IF - - END DO - - END FUNCTION coincidentNodes - - SUBROUTINE connectedTetraTetra(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA - CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemB - - !TODO: Try to find a much clear way to do this - - !Check surface 1 - IF (.NOT. ASSOCIATED(elemA%e1)) THEN - IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), & - (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN - - elemA%e1 => elemB - elemB%e1 => elemA - - ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), & - (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN - - elemA%e1 => elemB - elemB%e2 => elemA - - ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), & - (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN - - elemA%e1 => elemB - elemB%e3 => elemA - - ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), & - (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN - - elemA%e1 => elemB - elemB%e4 => elemA - - END IF - - END IF - - !Check surface 2 - IF (.NOT. ASSOCIATED(elemA%e2)) THEN - IF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & - (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN - - elemA%e2 => elemB - elemB%e1 => elemA - - ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & - (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN - - elemA%e2 => elemB - elemB%e2 => elemA - - ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & - (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN - - elemA%e2 => elemB - elemB%e3 => elemA - - ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & - (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN - - elemA%e2 => elemB - elemB%e4 => elemA - - END IF - - END IF - - !Check surface 3 - IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & - (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN - - elemA%e3 => elemB - elemB%e1 => elemA - - ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & - (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN - - elemA%e3 => elemB - elemB%e2 => elemA - - ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & - (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN - - elemA%e3 => elemB - elemB%e3 => elemA - - ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & - (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN - - elemA%e3 => elemB - elemB%e4 => elemA - - END IF - - END IF - - !Check surface 4 - IF (.NOT. ASSOCIATED(elemA%e4)) THEN - IF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), & - (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN - - elemA%e4 => elemB - elemB%e1 => elemA - - ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), & - (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN - - elemA%e4 => elemB - elemB%e2 => elemA - - ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), & - (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN - - elemA%e4 => elemB - elemB%e3 => elemA - - ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), & - (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN - - elemA%e4 => elemB - elemB%e4 => elemA - - END IF - - END IF - - END SUBROUTINE connectedTetraTetra - - SUBROUTINE connectedTetraEdge(elemA, elemB) - USE moduleMath - IMPLICIT NONE - - CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA - CLASS(meshEdge3DCartTria), INTENT(inout), TARGET:: elemB - INTEGER:: nodesEdge(1:3) - REAL(8), DIMENSION(1:3):: vec1, vec2 - REAL(8):: normVol(1:3) - - nodesEdge = (/ elemB%n1%n, elemB%n2%n, elemB%n3%n /) - - !Check surface 1 - IF (.NOT. ASSOCIATED(elemA%e1)) THEN - IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elema%n3%n/), & - nodesEdge)) THEN - - vec1 = (/ elemA%x(2) - elemA%x(1), & - elemA%y(2) - elemA%y(1), & - elemA%z(2) - elemA%z(1) /) - vec2 = (/ elemA%x(3) - elemA%x(1), & - elemA%y(3) - elemA%y(1), & - elemA%z(3) - elemA%z(1) /) - normVol = crossProduct(vec1, vec2) - normVol = normalize(normVol) - - IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN - - elemA%e1 => elemB - elemB%e1 => elemA - - ELSE - - elemA%e1 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = -elemB%normal - - END IF - - END IF - - END IF - - !Check surface 2 - IF (.NOT. ASSOCIATED(elemA%e2)) THEN - IF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & - nodesEdge)) THEN - - vec1 = (/ elemA%x(3) - elemA%x(2), & - elemA%y(3) - elemA%y(2), & - elemA%z(3) - elemA%z(2) /) - vec2 = (/ elemA%x(4) - elemA%x(2), & - elemA%y(4) - elemA%y(2), & - elemA%z(4) - elemA%z(2) /) - normVol = crossProduct(vec1, vec2) - normVol = normalize(normVol) - - IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN - - elemA%e2 => elemB - elemB%e1 => elemA - - ELSE - - elemA%e2 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = -elemB%normal - - END IF - - END IF - - END IF - - !Check surface 3 - IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elema%n4%n/), & - nodesEdge)) THEN - - vec1 = (/ elemA%x(2) - elemA%x(1), & - elemA%y(2) - elemA%y(1), & - elemA%z(2) - elemA%z(1) /) - vec2 = (/ elemA%x(4) - elemA%x(1), & - elemA%y(4) - elemA%y(1), & - elemA%z(4) - elemA%z(1) /) - normVol = crossProduct(vec1, vec2) - normVol = normalize(normVol) - - IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN - - elemA%e3 => elemB - elemB%e1 => elemA - - ELSE - - elemA%e3 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = -elemB%normal - - END IF - - END IF - - END IF - - !Check surface 4 - IF (.NOT. ASSOCIATED(elemA%e4)) THEN - IF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elema%n4%n/), & - nodesEdge)) THEN - - vec1 = (/ elemA%x(3) - elemA%x(1), & - elemA%y(3) - elemA%y(1), & - elemA%z(3) - elemA%z(1) /) - vec2 = (/ elemA%x(4) - elemA%x(1), & - elemA%y(4) - elemA%y(1), & - elemA%z(4) - elemA%z(1) /) - normVol = crossProduct(vec1, vec2) - normVol = normalize(normVol) - - IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN - - elemA%e4 => elemB - elemB%e1 => elemA - - - ELSE - - elemA%e4 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = -elemB%normal - - END IF - - END IF - - END IF - - END SUBROUTINE connectedTetraEdge - - SUBROUTINE constructGlobalK(K, elem) - IMPLICIT NONE - - REAL(8), INTENT(inout):: K(1:, 1:) - CLASS(meshVol), INTENT(in):: elem - REAL(8), ALLOCATABLE:: localK(:,:) - INTEGER:: nNodes, i, j - INTEGER, ALLOCATABLE:: n(:) - - SELECT TYPE(elem) - TYPE IS(meshVol3DCartTetra) - nNodes = 4 - ALLOCATE(localK(1:nNodes,1:nNodes)) - localK = elem%elemK() - ALLOCATE(n(1:nNodes)) - n = (/ elem%n1%n, elem%n2%n, & - elem%n3%n, elem%n4%n /) - - CLASS DEFAULT - nNodes = 0 - ALLOCATE(localK(1:1, 1:1)) - localK = 0.D0 - ALLOCATE(n(1:1)) - n = 0 - - END SELECT - - DO i = 1, nNodes - DO j = 1, nNodes - K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j) - END DO - END DO - - END SUBROUTINE constructGlobalK - -END MODULE moduleMesh3DCartRead diff --git a/src/modules/mesh/inout/gmsh2/makefile b/src/modules/mesh/inout/gmsh2/makefile new file mode 100644 index 0000000..68ad11b --- /dev/null +++ b/src/modules/mesh/inout/gmsh2/makefile @@ -0,0 +1,7 @@ +all: moduleMeshInputGmsh2.o moduleMeshOutputGmsh2.o + +moduleMeshInputGmsh2.o: moduleMeshOutputGmsh2.o moduleMeshInputGmsh2.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ + +%.o: %.f90 + $(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@ diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 new file mode 100644 index 0000000..7cfdb01 --- /dev/null +++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 @@ -0,0 +1,295 @@ +MODULE moduleMeshInputGmsh2 + + CONTAINS + !Inits a mesh to use Gmsh2 format + SUBROUTINE initGmsh2(self) + USE moduleMesh + USE moduleMeshOutputGmsh2 + IMPLICIT NONE + + TYPE(meshParticle), INTENT(inout):: self + + self%printOutput => printOutputGmsh2 + self%printColl => printCollGmsh2 + self%printEM => printEMGmsh2 + self%readMesh => readGmsh2 + + END SUBROUTINE initGmsh2 + + !Reads a Gmsh 2 format + SUBROUTINE readGmsh2(self, filename) + USE moduleMesh3DCart + USE moduleMesh2DCyl + USE moduleMesh2DCart + USE moduleMesh1DRad + USE moduleMesh1DCart + USE moduleBoundary + IMPLICIT NONE + + CLASS(meshParticle), INTENT(inout):: self + CHARACTER(:), ALLOCATABLE, INTENT(in):: filename + REAL(8):: x1, x2, x3 !3 generic coordinates + INTEGER, ALLOCATABLE:: p(:) !Array for nodes + INTEGER:: e = 0, n = 0, eTemp = 0, elemType = 0, bt = 0 + INTEGER:: totalNumElem + INTEGER:: boundaryType + + !Read mesh + OPEN(10, FILE=TRIM(filename)) + + !Skip header + READ(10, *) + READ(10, *) + READ(10, *) + READ(10, *) + + !Read number of nodes + READ(10, *) self%numNodes + + !Allocate required matrices and vectors + ALLOCATE(self%nodes(1:self%numNodes)) + 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 + + !Read the nodes information + DO e = 1, self%numNodes + READ(10, *) n, x1, x2, x3 + SELECT CASE(self%geometry) + CASE("3DCart") + ALLOCATE(meshNode3Dcart::self%nodes(n)%obj) + CALL self%nodes(n)%obj%init(n, (/x1, x2, x3 /)) + + CASE("2DCyl") + ALLOCATE(meshNode2DCyl:: self%nodes(n)%obj) + CALL self%nodes(n)%obj%init(n, (/x1, x2, 0.D0 /)) + + CASE("2DCart") + ALLOCATE(meshNode2DCart:: self%nodes(n)%obj) + CALL self%nodes(n)%obj%init(n, (/x1, x2, 0.D0 /)) + + CASE("1DRad") + ALLOCATE(meshNode1DRad:: self%nodes(n)%obj) + CALL self%nodes(n)%obj%init(n, (/x1, 0.D0, 0.D0 /)) + + CASE("1DCart") + ALLOCATE(meshNode1DCart:: self%nodes(n)%obj) + CALL self%nodes(n)%obj%init(n, (/x1, 0.D0, 0.D0 /)) + + END SELECT + + END DO + + !Skip comments + READ(10, *) + READ(10, *) + + !Reads total number of elements (no nodes) + READ(10, *) totalNumElem + + !conts edges and volume elements + self%numEdges = 0 + DO e = 1, totalNumElem + READ(10, *) eTemp, elemType + SELECT CASE(self%geometry) + CASE("3DCart") + !Element type 2 is triangle in gmsh + IF (elemType == 2) self%numEdges = e + + CASE("2DCyl","2DCart") + !Element type 1 is segment in Gmsh + IF (elemType == 1) self%numEdges = e + + CASE("1DRad","1DCart") + !Element type 15 is physical point in Gmsh + IF (elemType == 15) self%numEdges = e + + END SELECT + + END DO + !Substract the number of edges to the total number of elements + !to obtain the number of volume elements + self%numVols = TotalnumElem - self%numEdges + !Allocates arrays + ALLOCATE(self%edges(1:self%numEdges)) + ALLOCATE(self%vols(1:self%numVols)) + + !Go back to the beggining to read elements + DO e=1, totalNumElem + BACKSPACE(10) + END DO + + !Reads edges + DO e=1, self%numEdges + !Reads the edge according to the geometry + SELECT CASE(self%geometry) + CASE("3DCart") + READ(10, *) n, elemType, eTemp, boundaryType + BACKSPACE(10) + + !Associate boundary condition procedure. + bt = getBoundaryID(boundaryType) + + SELECT CASE(elemType) + CASE(2) + !Triangular surface + ALLOCATE(p(1:3)) + + READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1:3) + + ALLOCATE(meshEdge3DCartTria:: self%edges(e)%obj) + + CALL self%edges(e)%obj%init(n, p(1:3), bt, boundaryType) + + DEALLOCATE(p) + + END SELECT + + CASE("2DCyl") + 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) + + CALL self%edges(e)%obj%init(n, p(1:2), bt, boundaryType) + + DEALLOCATE(p) + + CASE("2DCart") + 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) + + CALL self%edges(e)%obj%init(n, p(1:2), bt, boundaryType) + + DEALLOCATE(p) + + CASE("1DRad") + 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) + + CALL self%edges(e)%obj%init(n, p(1:1), bt, boundaryType) + + DEALLOCATE(p) + + CASE("1DCart") + 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) + + CALL self%edges(e)%obj%init(n, p(1:1), bt, boundaryType) + + DEALLOCATE(p) + + END SELECT + + END DO + + !Read and initialize volumes + DO e = 1, self%numVols + !Reads the volume according to the geometry + SELECT CASE(self%geometry) + CASE("3DCart") + READ(10, *) n, elemType + BACKSPACE(10) + + SELECT CASE(elemType) + CASE(4) + !Tetrahedron element + ALLOCATE(p(1:4)) + READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:4) + ALLOCATE(meshVol3DCartTetra:: self%vols(e)%obj) + CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4)) + + END SELECT + + DEALLOCATE(p) + + CASE("2DCyl") + READ(10,*) n, elemType + BACKSPACE(10) + + SELECT CASE(elemType) + CASE (2) + !Triangular element + ALLOCATE(p(1:3)) + READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3) + ALLOCATE(meshVol2DCylTria:: self%vols(e)%obj) + CALL self%vols(e)%obj%init(n - self%numEdges, p(1:3)) + + CASE (3) + !Quadrilateral element + ALLOCATE(p(1:4)) + READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4) + ALLOCATE(meshVol2DCylQuad:: self%vols(e)%obj) + CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4)) + + END SELECT + + DEALLOCATE(p) + + CASE("2DCart") + READ(10,*) n, elemType + BACKSPACE(10) + + SELECT CASE(elemType) + CASE (2) + !Triangular element + ALLOCATE(p(1:3)) + READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3) + ALLOCATE(meshVol2DCartTria:: self%vols(e)%obj) + CALL self%vols(e)%obj%init(n - self%numEdges, p(1:3)) + + CASE (3) + !Quadrilateral element + ALLOCATE(p(1:4)) + READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4) + ALLOCATE(meshVol2DCartQuad:: self%vols(e)%obj) + CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4)) + + END SELECT + + DEALLOCATE(p) + + CASE("1DRad") + ALLOCATE(p(1:2)) + + READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2) + ALLOCATE(meshVol1DRadSegm:: self%vols(e)%obj) + CALL self%vols(e)%obj%init(n - self%numEdges, p(1:2)) + + DEALLOCATE(p) + + CASE("1DCart") + ALLOCATE(p(1:2)) + + READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2) + ALLOCATE(meshVol1DCartSegm:: self%vols(e)%obj) + CALL self%vols(e)%obj%init(n - self%numEdges, p(1:2)) + + DEALLOCATE(p) + + END SELECT + + END DO + + END SUBROUTINE readGmsh2 + +END MODULE moduleMeshInputGmsh2 diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 new file mode 100644 index 0000000..902206d --- /dev/null +++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 @@ -0,0 +1,197 @@ +MODULE moduleMeshOutputGmsh2 + + CONTAINS + !Prints the scattered properties of particles into the nodes + SUBROUTINE printOutputGmsh2(self, t) + USE moduleMesh + USE moduleRefParam + USE moduleSpecies + USE moduleOutput + IMPLICIT NONE + + CLASS(meshParticle), INTENT(in):: self + INTEGER, INTENT(in):: t + INTEGER:: n, i + TYPE(outputFormat):: output(1:self%numNodes) + REAL(8):: time + CHARACTER(:), ALLOCATABLE:: fileName + CHARACTER (LEN=iterationDigits):: tstring + + time = DBLE(t)*tauMin*ti_ref + + DO i = 1, nSpecies + WRITE(tstring, iterationFormat) t + fileName='OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh' + WRITE(*, "(6X,A15,A)") "Creating file: ", fileName + OPEN (60, file = path // folder // '/' // fileName) + WRITE(60, "(A)") '$MeshFormat' + WRITE(60, "(A)") '2.2 0 8' + WRITE(60, "(A)") '$EndMeshFormat' + WRITE(60, "(A)") '$NodeData' + WRITE(60, "(A)") '1' + WRITE(60, "(A)") '"Density ' // species(i)%obj%name // ' (m^-3)"' + WRITE(60, *) 1 + WRITE(60, *) time + WRITE(60, *) 3 + WRITE(60, *) t + WRITE(60, *) 1 + WRITE(60, *) self%numNodes + DO n=1, self%numNodes + CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj) + WRITE(60, "(I6,ES20.6E3)") n, output(n)%density + END DO + WRITE(60, "(A)") '$EndNodeData' + WRITE(60, "(A)") '$NodeData' + WRITE(60, "(A)") '1' + WRITE(60, "(A)") '"Velocity ' // species(i)%obj%name // ' (m/s)"' + WRITE(60, *) 1 + WRITE(60, *) time + WRITE(60, *) 3 + WRITE(60, *) t + WRITE(60, *) 3 + WRITE(60, *) self%numNodes + DO n=1, self%numNodes + WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity + END DO + WRITE(60, "(A)") '$EndNodeData' + WRITE(60, "(A)") '$NodeData' + WRITE(60, "(A)") '1' + WRITE(60, "(A)") '"Pressure ' // species(i)%obj%name // ' (Pa)"' + WRITE(60, *) 1 + WRITE(60, *) time + WRITE(60, *) 3 + WRITE(60, *) t + WRITE(60, *) 1 + WRITE(60, *) self%numNodes + DO n=1, self%numNodes + WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure + END DO + WRITE(60, "(A)") '$EndNodeData' + WRITE(60, "(A)") '$NodeData' + WRITE(60, "(A)") '1' + WRITE(60, "(A)") '"Temperature ' // species(i)%obj%name // ' (K)"' + WRITE(60, *) 1 + WRITE(60, *) time + WRITE(60, *) 3 + WRITE(60, *) t + WRITE(60, *) 1 + WRITE(60, *) self%numNodes + DO n=1, self%numNodes + WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature + END DO + WRITE(60, "(A)") '$EndNodeData' + CLOSE (60) + + END DO + + END SUBROUTINE printOutputGmsh2 + + !Prints the number of collisions into the volumes + SUBROUTINE printCollGmsh2(self, t) + USE moduleMesh + USE moduleRefParam + USE moduleCaseParam + USE moduleCollisions + USE moduleOutput + IMPLICIT NONE + + CLASS(meshParticle), INTENT(in):: self + INTEGER, INTENT(in):: t + INTEGER:: n + REAL(8):: time + CHARACTER(:), ALLOCATABLE:: fileName + CHARACTER (LEN=iterationDigits):: tstring + + + IF (collOutput) THEN + time = DBLE(t)*tauMin*ti_ref + WRITE(tstring, iterationFormat) t + + fileName='OUTPUT_' // tstring// '_Collisions.msh' + WRITE(*, "(6X,A15,A)") "Creating file: ", fileName + OPEN (60, file = path // folder // '/' // fileName) + WRITE(60, "(A)") '$MeshFormat' + WRITE(60, "(A)") '2.2 0 8' + WRITE(60, "(A)") '$EndMeshFormat' + WRITE(60, "(A)") '$ElementData' + WRITE(60, "(A)") '1' + WRITE(60, "(A)") '"Collisions"' + WRITE(60, *) 1 + WRITE(60, *) time + WRITE(60, *) 3 + WRITE(60, *) t + WRITE(60, *) 1 + WRITE(60, *) self%numVols + DO n=1, self%numVols + WRITE(60, "(I6,I10)") n + self%numEdges, self%vols(n)%obj%nColl + END DO + WRITE(60, "(A)") '$EndElementData' + + CLOSE(60) + + END IF + + END SUBROUTINE printCollGmsh2 + + !Prints the electrostatic EM properties into the nodes and volumes + SUBROUTINE printEMGmsh2(self, t) + USE moduleMesh + USE moduleRefParam + USE moduleCaseParam + USE moduleOutput + IMPLICIT NONE + + CLASS(meshParticle), INTENT(in):: self + INTEGER, INTENT(in):: t + INTEGER:: n, e + REAL(8):: time + CHARACTER(:), ALLOCATABLE:: fileName + CHARACTER (LEN=iterationDigits):: tstring + REAL(8):: xi(1:3) + + xi = (/ 0.D0, 0.D0, 0.D0 /) + + IF (emOutput) THEN + time = DBLE(t)*tauMin*ti_ref + WRITE(tstring, iterationFormat) t + + fileName='OUTPUT_' // tstring// '_EMField.msh' + WRITE(*, "(6X,A15,A)") "Creating file: ", fileName + OPEN (20, file = path // folder // '/' // fileName) + WRITE(20, "(A)") '$MeshFormat' + WRITE(20, "(A)") '2.2 0 8' + WRITE(20, "(A)") '$EndMeshFormat' + WRITE(20, "(A)") '$NodeData' + WRITE(20, "(A)") '1' + WRITE(20, "(A)") '"Potential (V)"' + WRITE(20, *) 1 + WRITE(20, *) time + WRITE(20, *) 3 + WRITE(20, *) t + WRITE(20, *) 1 + WRITE(20, *) self%numNodes + DO n=1, self%numNodes + WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref + END DO + WRITE(20, "(A)") '$EndNodeData' + + WRITE(20, "(A)") '$ElementData' + WRITE(20, "(A)") '1' + WRITE(20, "(A)") '"Electric Field (V/m)"' + WRITE(20, *) 1 + WRITE(20, *) time + WRITE(20, *) 3 + WRITE(20, *) t + WRITE(20, *) 3 + WRITE(20, *) self%numVols + DO e=1, self%numVols + WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF(xi)*EF_ref + END DO + WRITE(20, "(A)") '$EndElementData' + CLOSE(20) + + END IF + + END SUBROUTINE printEMGmsh2 + +END MODULE moduleMeshOutputGmsh2 diff --git a/src/modules/mesh/inout/makefile b/src/modules/mesh/inout/makefile new file mode 100644 index 0000000..7817608 --- /dev/null +++ b/src/modules/mesh/inout/makefile @@ -0,0 +1,4 @@ +all: gmsh2.o + +gmsh2.o: + $(MAKE) -C gmsh2 all diff --git a/src/modules/mesh/makefile b/src/modules/mesh/makefile index 9849ba9..0adc8f6 100644 --- a/src/modules/mesh/makefile +++ b/src/modules/mesh/makefile @@ -1,18 +1,18 @@ -all: moduleMesh.o moduleMeshBoundary.o 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o +all: moduleMesh.o moduleMeshBoundary.o inout.o 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o -3DCart.o: +3DCart.o: moduleMesh.o $(MAKE) -C 3DCart all -2DCyl.o: +2DCyl.o: moduleMesh.o $(MAKE) -C 2DCyl all -2DCart.o: +2DCart.o: moduleMesh.o $(MAKE) -C 2DCart all -1DCart.o: +1DCart.o: moduleMesh.o $(MAKE) -C 1DCart all -1DRad.o: +1DRad.o: moduleMesh.o $(MAKE) -C 1DRad all moduleMesh.o: moduleMesh.f90 @@ -20,3 +20,6 @@ moduleMesh.o: moduleMesh.f90 moduleMeshBoundary.o: moduleMesh.o moduleMeshBoundary.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ + +inout.o: 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o + $(MAKE) -C inout all diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 54ff883..e80b957 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -163,6 +163,7 @@ MODULE moduleMesh PROCEDURE(randPosVol_interface), DEFERRED, PASS:: randPos PROCEDURE(scatter_interface), DEFERRED, PASS:: scatter PROCEDURE(gatherEF_interface), DEFERRED, PASS:: gatherEF + PROCEDURE(elemK_interface), DEFERRED, PASS:: elemK PROCEDURE(elemF_interface), DEFERRED, PASS:: elemF PROCEDURE, PASS:: findCell PROCEDURE(phy2log_interface), DEFERRED, PASS:: phy2log @@ -206,6 +207,13 @@ MODULE moduleMesh END FUNCTION getNodesVol_interface + PURE FUNCTION elemK_interface(self) RESULT(localK) + IMPORT:: meshVol + CLASS(meshVol), INTENT(in):: self + REAL(8), ALLOCATABLE:: localK(:,:) + + END FUNCTION elemK_interface + PURE FUNCTION elemF_interface(self, source) RESULT(localF) IMPORT:: meshVol CLASS(meshVol), INTENT(in):: self @@ -252,8 +260,8 @@ MODULE moduleMesh END TYPE meshVolCont - !Abstract type of mesh - TYPE, PUBLIC, ABSTRACT:: meshGeneric + !Particle mesh + TYPE, PUBLIC:: meshParticle INTEGER:: numEdges, numNodes, numVols !Array of nodes TYPE(meshNodeCont), ALLOCATABLE:: nodes(:) @@ -261,69 +269,71 @@ MODULE moduleMesh TYPE(meshEdgeCont), ALLOCATABLE:: edges(:) !Array of volume elements TYPE(meshVolCont), ALLOCATABLE:: vols(:) + !Geometry of the mesh + CHARACTER(:), ALLOCATABLE:: geometry !Global stiffness matrix REAL(8), ALLOCATABLE, DIMENSION(:,:):: K !Permutation matrix for P L U factorization INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV - PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL() - PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL() - PROCEDURE(printEM_interface), POINTER, PASS:: printEM => NULL() + PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL() + PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL() + PROCEDURE(printEM_interface), POINTER, PASS:: printEM => NULL() + PROCEDURE(readMesh_interface), POINTER, PASS:: readMesh => NULL() + PROCEDURE(connectMesh_interface), POINTER, PASS:: connectMesh => NULL() CONTAINS - PROCEDURE(initMesh_interface), DEFERRED, PASS:: init - PROCEDURE(readMesh_interface), DEFERRED, PASS:: readMesh + PROCEDURE, PASS:: constructGlobalK - END TYPE meshGeneric + END TYPE meshParticle ABSTRACT INTERFACE - !Inits the mesh - SUBROUTINE initMesh_interface(self, meshFormat) - IMPORT meshGeneric - - CLASS(meshGeneric), INTENT(out):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat - - END SUBROUTINE initMesh_interface - - !Reads the mesh from a file - SUBROUTINE readMesh_interface(self, filename) - IMPORT meshGeneric - - CLASS(meshGeneric), INTENT(inout):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: filename - - END SUBROUTINE readMesh_interface - !Prints Species data SUBROUTINE printOutput_interface(self, t) - IMPORT meshGeneric + IMPORT meshParticle - CLASS(meshGeneric), INTENT(in):: self + CLASS(meshParticle), INTENT(in):: self INTEGER, INTENT(in):: t END SUBROUTINE printOutput_interface !Prints number of collisions SUBROUTINE printColl_interface(self, t) - IMPORT meshGeneric + IMPORT meshParticle - CLASS(meshGeneric), INTENT(in):: self + CLASS(meshParticle), INTENT(in):: self INTEGER, INTENT(in):: t END SUBROUTINE printColl_interface !Prints EM info SUBROUTINE printEM_interface(self, t) - IMPORT meshGeneric + IMPORT meshParticle - CLASS(meshGeneric), INTENT(in):: self + CLASS(meshParticle), INTENT(in):: self INTEGER, INTENT(in):: t END SUBROUTINE printEM_interface + !Reads the mesh from a file + SUBROUTINE readMesh_interface(self, filename) + IMPORT meshParticle + + CLASS(meshParticle), INTENT(inout):: self + CHARACTER(:), ALLOCATABLE, INTENT(in):: filename + + END SUBROUTINE readMesh_interface + + SUBROUTINE connectMesh_interface(self) + IMPORT meshParticle + + CLASS(meshParticle), INTENT(inout):: self + + END SUBROUTINE connectMesh_interface + + END INTERFACE - !Generic mesh - CLASS(meshGeneric), ALLOCATABLE, TARGET:: mesh + !Particle mesh + TYPE(meshParticle), TARGET:: mesh CONTAINS !Reset the output of node @@ -467,191 +477,31 @@ MODULE moduleMesh END SUBROUTINE collision - SUBROUTINE printOutputGmsh(self, t) - USE moduleRefParam - USE moduleSpecies - USE moduleOutput + !Constructs the global K matrix + SUBROUTINE constructGlobalK(self) IMPLICIT NONE - CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t - INTEGER:: n, i - TYPE(outputFormat):: output(1:self%numNodes) - REAL(8):: time - CHARACTER(:), ALLOCATABLE:: fileName - CHARACTER (LEN=iterationDigits):: tstring + CLASS(meshParticle), INTENT(inout):: self + INTEGER:: e + INTEGER, ALLOCATABLE:: n(:) + REAL(8), ALLOCATABLE:: localK(:,:) + INTEGER:: nNodes, i, j - time = DBLE(t)*tauMin*ti_ref + DO e = 1, self%numVols + n = self%vols(e)%obj%getNodes() + localK = self%vols(e)%obj%elemK() + nNodes = SIZE(n) + + DO i = 1, nNodes + DO j = 1, nNodes + self%K(n(i), n(j)) = self%K(n(i), n(j)) + localK(i, j) + + END DO - DO i = 1, nSpecies - WRITE(tstring, iterationFormat) t - fileName='OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh' - WRITE(*, "(6X,A15,A)") "Creating file: ", fileName - OPEN (60, file = path // folder // '/' // fileName) - WRITE(60, "(A)") '$MeshFormat' - WRITE(60, "(A)") '2.2 0 8' - WRITE(60, "(A)") '$EndMeshFormat' - WRITE(60, "(A)") '$NodeData' - WRITE(60, "(A)") '1' - WRITE(60, "(A)") '"Density ' // species(i)%obj%name // ' (m^-3)"' - WRITE(60, *) 1 - WRITE(60, *) time - WRITE(60, *) 3 - WRITE(60, *) t - WRITE(60, *) 1 - WRITE(60, *) self%numNodes - DO n=1, self%numNodes - CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj) - WRITE(60, "(I6,ES20.6E3)") n, output(n)%density END DO - WRITE(60, "(A)") '$EndNodeData' - WRITE(60, "(A)") '$NodeData' - WRITE(60, "(A)") '1' - WRITE(60, "(A)") '"Velocity ' // species(i)%obj%name // ' (m/s)"' - WRITE(60, *) 1 - WRITE(60, *) time - WRITE(60, *) 3 - WRITE(60, *) t - WRITE(60, *) 3 - WRITE(60, *) self%numNodes - DO n=1, self%numNodes - WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity - END DO - WRITE(60, "(A)") '$EndNodeData' - WRITE(60, "(A)") '$NodeData' - WRITE(60, "(A)") '1' - WRITE(60, "(A)") '"Pressure ' // species(i)%obj%name // ' (Pa)"' - WRITE(60, *) 1 - WRITE(60, *) time - WRITE(60, *) 3 - WRITE(60, *) t - WRITE(60, *) 1 - WRITE(60, *) self%numNodes - DO n=1, self%numNodes - WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure - END DO - WRITE(60, "(A)") '$EndNodeData' - WRITE(60, "(A)") '$NodeData' - WRITE(60, "(A)") '1' - WRITE(60, "(A)") '"Temperature ' // species(i)%obj%name // ' (K)"' - WRITE(60, *) 1 - WRITE(60, *) time - WRITE(60, *) 3 - WRITE(60, *) t - WRITE(60, *) 1 - WRITE(60, *) self%numNodes - DO n=1, self%numNodes - WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature - END DO - WRITE(60, "(A)") '$EndNodeData' - CLOSE (60) END DO - END SUBROUTINE printOutputGmsh + END SUBROUTINE constructGlobalK - SUBROUTINE printCollGmsh(self, t) - USE moduleRefParam - USE moduleCaseParam - USE moduleCollisions - USE moduleOutput - IMPLICIT NONE - - CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t - INTEGER:: n - REAL(8):: time - CHARACTER(:), ALLOCATABLE:: fileName - CHARACTER (LEN=iterationDigits):: tstring - - - IF (collOutput) THEN - time = DBLE(t)*tauMin*ti_ref - WRITE(tstring, iterationFormat) t - - fileName='OUTPUT_' // tstring// '_Collisions.msh' - WRITE(*, "(6X,A15,A)") "Creating file: ", fileName - OPEN (60, file = path // folder // '/' // fileName) - WRITE(60, "(A)") '$MeshFormat' - WRITE(60, "(A)") '2.2 0 8' - WRITE(60, "(A)") '$EndMeshFormat' - WRITE(60, "(A)") '$ElementData' - WRITE(60, "(A)") '1' - WRITE(60, "(A)") '"Collisions"' - WRITE(60, *) 1 - WRITE(60, *) time - WRITE(60, *) 3 - WRITE(60, *) t - WRITE(60, *) 1 - WRITE(60, *) self%numVols - DO n=1, self%numVols - WRITE(60, "(I6,I10)") n + self%numEdges, self%vols(n)%obj%nColl - END DO - WRITE(60, "(A)") '$EndElementData' - - CLOSE(60) - - END IF - - END SUBROUTINE printCollGmsh - - SUBROUTINE printEMGmsh(self, t) - USE moduleRefParam - USE moduleCaseParam - USE moduleOutput - IMPLICIT NONE - - CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t - INTEGER:: n, e - REAL(8):: time - CHARACTER(:), ALLOCATABLE:: fileName - CHARACTER (LEN=iterationDigits):: tstring - REAL(8):: xi(1:3) - - xi = (/ 0.D0, 0.D0, 0.D0 /) - - IF (emOutput) THEN - time = DBLE(t)*tauMin*ti_ref - WRITE(tstring, iterationFormat) t - - fileName='OUTPUT_' // tstring// '_EMField.msh' - WRITE(*, "(6X,A15,A)") "Creating file: ", fileName - OPEN (20, file = path // folder // '/' // fileName) - WRITE(20, "(A)") '$MeshFormat' - WRITE(20, "(A)") '2.2 0 8' - WRITE(20, "(A)") '$EndMeshFormat' - WRITE(20, "(A)") '$NodeData' - WRITE(20, "(A)") '1' - WRITE(20, "(A)") '"Potential (V)"' - WRITE(20, *) 1 - WRITE(20, *) time - WRITE(20, *) 3 - WRITE(20, *) t - WRITE(20, *) 1 - WRITE(20, *) self%numNodes - DO n=1, self%numNodes - WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref - END DO - WRITE(20, "(A)") '$EndNodeData' - - WRITE(20, "(A)") '$ElementData' - WRITE(20, "(A)") '1' - WRITE(20, "(A)") '"Electric Field (V/m)"' - WRITE(20, *) 1 - WRITE(20, *) time - WRITE(20, *) 3 - WRITE(20, *) t - WRITE(20, *) 3 - WRITE(20, *) self%numVols - DO e=1, self%numVols - WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF(xi)*EF_ref - END DO - WRITE(20, "(A)") '$EndElementData' - CLOSE(20) - - END IF - - END SUBROUTINE printEMGmsh - END MODULE moduleMesh diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index 98933de..395d031 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -689,11 +689,12 @@ MODULE moduleInput !Read the geometry (mesh) for the case SUBROUTINE readGeometry(config) USE moduleMesh - USE moduleMesh3DCartRead, ONLY: mesh3DCartGeneric - USE moduleMesh2DCylRead, ONLY: mesh2DCylGeneric - USE moduleMesh2DCartRead, ONLY: mesh2DCartGeneric - USE moduleMesh1DCartRead, ONLY: mesh1DCartGeneric - USE moduleMesh1DRadRead, ONLY: mesh1DRadGeneric + USE moduleMeshInputGmsh2, ONLY: initGmsh2 + USE moduleMesh3DCart, ONLY: connectMesh3DCart + USE moduleMesh2DCyl, ONLY: connectMesh2DCyl + USE moduleMesh2DCart, ONLY: connectMesh2DCart + USE moduleMesh1DRad, ONLY: connectMesh1DRad + USE moduleMesh1DCart, ONLY: connectMesh1DCart USE moduleErrors USE moduleOutput USE json_module @@ -701,45 +702,51 @@ MODULE moduleInput TYPE(json_file), INTENT(inout):: config LOGICAL:: found - CHARACTER(:), ALLOCATABLE:: geometryType, meshFormat, meshFile + CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile CHARACTER(:), ALLOCATABLE:: fullPath !Selects the type of geometry. - CALL config%get('geometry.type', geometryType, found) - SELECT CASE(geometryType) - CASE ("3DCart") - !Creates a 3D cylindrical mesh - ALLOCATE(mesh3DCartGeneric:: mesh) - - CASE ("2DCyl") - !Creates a 2D cylindrical mesh - ALLOCATE(mesh2DCylGeneric:: mesh) - - CASE ("2DCart") - !Creates a 2D cylindrical mesh - ALLOCATE(mesh2DCartGeneric:: mesh) - - CASE ("1DCart") - !Creates a 1D cartesian mesh - ALLOCATE(mesh1DCartGeneric:: mesh) - - CASE ("1DRad") - !Creates a 1D cartesian mesh - ALLOCATE(mesh1DRadGeneric:: mesh) - - CASE DEFAULT - CALL criticalError("Geometry type " // geometryType // " not supported.", "readGeometry") - - END SELECT + CALL config%get('geometry.type', mesh%geometry, found) !Gets the type of mesh CALL config%get('geometry.meshType', meshFormat, found) - CALL mesh%init(meshFormat) - !Reads the mesh + SELECT CASE(meshFormat) + CASE ("gmsh2") + CALL initGmsh2(mesh) + + CASE DEFAULT + CALL criticalError("Mesh format " // meshFormat // " not recogniced", "readGeometry") + + END SELECT + + !Reads the mesh file CALL config%get('geometry.meshFile', meshFile, found) fullpath = path // meshFile CALL mesh%readMesh(fullPath) + !Creates the connectivity between elements + SELECT CASE(mesh%geometry) + CASE("3DCart") + mesh%connectMesh => connectMesh3DCart + + CASE("2DCyl") + mesh%connectMesh => connectMesh2DCyl + + CASE("2DCart") + mesh%connectMesh => connectMesh2DCart + + CASE("1DRad") + mesh%connectMesh => connectMesh1DRad + + CASE("1DCart") + mesh%connectMesh => connectMesh1DCart + + END SELECT + CALL mesh%connectMesh + + !Builds the K matrix + CALL mesh%constructGlobalK() + END SUBROUTINE readGeometry SUBROUTINE readEMBoundary(config)