diff --git a/fem/src/Adaptive.F90 b/fem/src/Adaptive.F90 index f853ee4a18..7f90f0ad9a 100644 --- a/fem/src/Adaptive.F90 +++ b/fem/src/Adaptive.F90 @@ -50,7 +50,8 @@ MODULE Adaptive USE LoadMod USE MeshUtils USE MeshRemeshing - + USE SaveUtils + IMPLICIT NONE @@ -1279,114 +1280,154 @@ FUNCTION External_ReMesh( RefMesh, ErrorLimit, HValue, NodalError, & TYPE(Mesh_t), POINTER :: NewMesh, RefMesh !------------------------------------------------------------------------------ TYPE(Mesh_t), POINTER :: Mesh - INTEGER :: i,j,k,n + INTEGER :: i,j,k,n,dim REAL(KIND=dp) :: Lambda + REAL(KIND=dp), POINTER :: HValueF(:) CHARACTER(:), ALLOCATABLE :: MeshCommand, Name, MeshInputFile + LOGICAL :: GmshFormat !------------------------------------------------------------------------------ - OPEN( 11, STATUS='UNKNOWN', FILE='bgmesh' ) - WRITE( 11,* ) COUNT( NodalError > 100*AEPS ) + dim = CoordinateSystemDimension() + ! Create a temporal field that includes the desired mesh density where Hvalue has been computed. + ! This is in terms of desired mesh density, currently -1 value is given if the nodal error is zero + ! implying that nothing is computed here. + ALLOCATE(HvalueF(SIZE(HValue))) + HValueF = -1.0_dp DO i=1,RefMesh % NumberOfNodes - IF ( NodalError(i) > 100*AEPS ) THEN - Lambda = ( ErrorLimit / NodalError(i) ) ** ( 1.0d0 / hConvergence(i) ) - - IF ( RefMesh % AdaptiveDepth < 1 ) THEN - Lambda = HValue(i) * MAX( MIN( Lambda, 1.33d0), 0.75d0) - ELSE - Lambda = HValue(i) * MAX(MIN(Lambda, MaxChange), 1.0d0/MaxChange) - END IF - - IF( .NOT.Coarsening ) Lambda = MIN( Lambda, Hvalue(i) ) - - IF ( maxH > 0 ) Lambda = MIN( Lambda, maxH ) - IF ( minH > 0 ) Lambda = MAX( Lambda, minH ) + IF ( NodalError(i) > 100*AEPS ) THEN + Lambda = ( ErrorLimit / NodalError(i) ) ** ( 1.0d0 / hConvergence(i) ) + IF ( RefMesh % AdaptiveDepth < 1 ) THEN + Lambda = HValue(i) * MAX( MIN( Lambda, 1.33d0), 0.75d0) + ELSE + Lambda = HValue(i) * MAX(MIN(Lambda, MaxChange), 1.0d0/MaxChange) + END IF + IF( .NOT.Coarsening ) Lambda = MIN( Lambda, Hvalue(i) ) - IF ( CoordinateSystemDimension() == 2 ) THEN - WRITE(11,'(3e23.15)') RefMesh % Nodes % x(i), & - RefMesh % Nodes % y(i), Lambda - ELSE - WRITE(11,'(4e23.15)') RefMesh % Nodes % x(i), & - RefMesh % Nodes % y(i), & - RefMesh % Nodes % z(i), Lambda - END IF - ELSE - IF ( CoordinateSystemDimension() == 2 ) THEN - WRITE(11,'(3e23.15)') RefMesh % Nodes % x(i), & - RefMesh % Nodes % y(i), HValue(i) - ELSE - WRITE(11,'(4e23.15)') RefMesh % Nodes % x(i), & - RefMesh % Nodes % y(i), & - RefMesh % Nodes % z(i), HValue(i) - END IF - END IF + IF ( maxH > 0 ) Lambda = MIN( Lambda, maxH ) + IF ( minH > 0 ) Lambda = MAX( Lambda, minH ) + HValueF(i) = Lambda + END IF END DO - - WRITE(11,*) 0 - CLOSE(11) + ! Save the current mesh in Elmer mesh format Path = ListGetString( Params, 'Adaptive Mesh Name', Found ) - IF ( .NOT. Found ) Path = 'RefinedMesh' - - i = RefMesh % AdaptiveDepth + 1 nLen = LEN_TRIM(Path) - Path = Path(1:nlen) // I2S(i) + + IF ( .NOT. Found ) THEN + i = RefMesh % AdaptiveDepth + 1 + Path = 'RefinedMesh'//I2S(i) + END IF nLen = LEN_TRIM(OutputPath) IF ( nlen > 0 ) THEN - Path = OutputPath(1:nlen) // '/' // TRIM(Path) + Path = OutputPath(1:nlen) // '/' // TRIM(Path) ELSE - Path = TRIM(Path) + Path = TRIM(Path) END IF + + GmshFormat = ListGetLogical( Params,'Adaptive Remesh Use Gmsh', Found ) + + IF( GmshFormat ) THEN + CALL Info( Caller,'Saving background mesh density in gmsh 2.0 (.msh) format' ) + + ! A cludge to change the pointer and save results in Gmsh format. + BLOCK + REAL(KIND=dp), POINTER :: PtoHvalue(:) + TYPE(Variable_t), POINTER :: HVar + HVar => VariableGet( RefMesh % Variables,'Hvalue') + pToHvalue => HVar % Values + HVar % Values => HvalueF + CALL ListAddString(Solver % Values,'Scalar Field 1','Hvalue') + CALL ListAddLogical(Solver % Values,'File Append',.FALSE.) + CALL ListAddLogical(Solver % Values,'Alter Topology',.TRUE.) + CALL ListAddNewString(Solver % Values, 'Output File Name', 'gmsh_bgmesh.msh') + CALL SaveGmshOutput( Model,Solver,0.0_dp,.FALSE.) + HVar % Values => PtoHvalue + END BLOCK + + ELSE + CALL Info( Caller,'Saving background mesh density in point cloud format' ) + + OPEN( 11, STATUS='UNKNOWN', FILE='bgmesh.nodes' ) + WRITE( 11,* ) COUNT( HValueF > 0.0_dp ) + DO i=1,RefMesh % NumberOfNodes + IF(.NOT. (HValueF(i) > 0.0_dp )) CYCLE + IF (dim == 2 ) THEN + WRITE(11,'(3e23.15)') RefMesh % Nodes % x(i), & + RefMesh % Nodes % y(i), HValueF(i) + ELSE + WRITE(11,'(4e23.15)') RefMesh % Nodes % x(i), & + RefMesh % Nodes % y(i), & + RefMesh % Nodes % z(i), HValueF(i) + END IF + END DO + WRITE(11,*) 0 + CLOSE(11) - CALL MakeDirectory( TRIM(Path) // CHAR(0) ) - CALL WriteMeshToDisk( RefMesh, Path ) + CALL MakeDirectory( TRIM(Path) // CHAR(0) ) + CALL WriteMeshToDisk( RefMesh, Path ) + END IF Mesh => RefMesh DO WHILE( ASSOCIATED( Mesh ) ) - IF ( Mesh % AdaptiveDepth == 0 ) EXIT - Mesh => Mesh % Parent + IF ( Mesh % AdaptiveDepth == 0 ) EXIT + Mesh => Mesh % Parent END DO - MeshInputFile = ListGetString( Params, 'Mesh Input File', Found ) - - IF ( .NOT. Found ) THEN - MeshInputFile = ListGetString( Model % Simulation, 'Mesh Input File' ) - END IF + MeshCommand = ListGetString( Solver % Values,'Mesh Command',Found) + IF(.NOT. Found ) THEN + IF( GmshFormat ) THEN + CALL Fatal('ReMesh','For now, provide "Mesh Command" for Gmsh meshing!') + END IF - MeshCommand = TRIM(OutputPath) // '/' // TRIM(Mesh % Name) // '/' // & - TRIM( MeshInputFile ) + MeshInputFile = ListGetString( Params, 'Mesh Input File', Found ) + IF ( .NOT. Found ) THEN + MeshInputFile = ListGetString( Model % Simulation, 'Mesh Input File' ) + END IF - SELECT CASE( CoordinateSystemDimension() ) - CASE(2) - MeshCommand = 'Mesh2D ' // TRIM(MeshCommand) // ' ' // & - TRIM(Path) // ' --bgmesh=bgmesh' + MeshCommand = TRIM(OutputPath) // '/' // TRIM(Mesh % Name) // '/' // & + TRIM( MeshInputFile ) - CASE(3) - MeshCommand = 'Mesh3D ' // TRIM(MeshCommand) // ' ' // & - TRIM(Path) // ' bgmesh' - END SELECT + SELECT CASE( dim ) + CASE(2) + ! Legacy mesh generator from the Elmer suite. + MeshCommand = 'Mesh2D '//TRIM(MeshCommand)//' '//TRIM(Path)// ' --bgmesh=bgmesh.nodes' + CASE(3) + MeshCommand = 'Mesh3D '//TRIM(MeshCommand)//' '//TRIM(Path)//' bgmesh.nodes' + END SELECT + END IF - CALL Info('ReMesh','System command: '//TRIM(MeshCommand),Level=10) + ! Remeshing command. + CALL Info('ReMesh','Meshing command: '//TRIM(MeshCommand),Level=10) CALL SystemCommand( MeshCommand ) + ! Check if also conversion command is given. + MeshCommand = ListGetString( Solver % Values,'Mesh Conversion Command',Found) + IF( Found ) THEN + CALL Info('ReMesh','Conversion command: '//TRIM(MeshCommand),Level=10) + CALL SystemCommand( MeshCommand ) + END IF + + ! Read the new mesh. NewMesh => LoadMesh2( Model, OutPutPath, Path, .FALSE., 1, 0 ) + ! Loading Gebhart factors is more or less obsolite. IF ( Solver % Variable % Name == 'temperature' ) THEN Name = ListGetString( Model % Simulation, 'Gebhart Factors', Found ) IF ( Found ) THEN MeshCommand = 'View ' // TRIM(OutputPath) // & '/' // TRIM(Mesh % Name) // ' ' // TRIM(Path) - CALL SystemCommand( MeshCommand ) Name = TRIM(OutputPath) // '/' // & - TRIM(Mesh % Name) // '/' // TRIM(Name) - + TRIM(Mesh % Name) // '/' // TRIM(Name) CALL LoadGebhartFactors( NewMesh, TRIM(Name) ) END IF END IF + DEALLOCATE(HvalueF) + !------------------------------------------------------------------------------ END FUNCTION External_ReMesh !------------------------------------------------------------------------------ diff --git a/fem/src/ParticleUtils.F90 b/fem/src/ParticleUtils.F90 index e7f7eae351..e186a4c1f3 100644 --- a/fem/src/ParticleUtils.F90 +++ b/fem/src/ParticleUtils.F90 @@ -49,6 +49,7 @@ MODULE ParticleUtils USE Lists USE MeshUtils USE GeneralUtils + USE SaveUtils IMPLICIT NONE diff --git a/fem/src/SaveUtils.F90 b/fem/src/SaveUtils.F90 index 4e788c3c22..52e787ebc4 100644 --- a/fem/src/SaveUtils.F90 +++ b/fem/src/SaveUtils.F90 @@ -774,6 +774,651 @@ SUBROUTINE GenerateSavePermutation(Mesh,DG,DN,LagN,SaveLinear,ActiveElem,NumberO END IF END SUBROUTINE GenerateSavePermutation - + + + !> Defines and potentially creates output directory. + !> The output directory may given in different ways, and even be part of the + !> filename, or be relative to home directory. We try to parse every possible + !> scenario here that user might have in mind. + !----------------------------------------------------------------------------- + SUBROUTINE SolverOutputDirectory( Solver, Filename, OutputDirectory, & + MakeDir, UseMeshDir ) + + USE ModelDescription + + TYPE(Solver_t) :: Solver + LOGICAL, OPTIONAL :: MakeDir, UseMeshDir + CHARACTER(*) :: Filename + CHARACTER(:), ALLOCATABLE :: OutputDirectory + + LOGICAL :: Found, AbsPathInName, DoDir, PartitioningSubDir + INTEGER :: nd, nf, n + CHARACTER(LEN=MAX_NAME_LEN) :: Str + + IF( PRESENT( MakeDir ) ) THEN + DoDir = MakeDir + ELSE + DoDir = ( Solver % TimesVisited == 0 ) .AND. ( ParEnv % MyPe == 0 ) + END IF + + ! Output directory is obtained in order + ! 1) solver section + ! 2) simulation section + ! 3) header section + OutputDirectory = ListGetString( Solver % Values,'Output Directory',Found) + IF(.NOT. Found) OutputDirectory = ListGetString( CurrentModel % Simulation,& + 'Output Directory',Found) + + IF(.NOT. Found) OutputDirectory = TRIM(OutputPath) + nd = LEN_TRIM(OutputDirectory) + + ! If the path is just working directory then that is not an excude + ! to not use the mesh name, or directory that comes with the filename + IF(.NOT. Found .AND. nd == 1 .AND. OutputDirectory(1:1)=='.') nd = 0 + + ! If requested by the optional parameter use the mesh directory when + ! no results directory given. This is an old convection used in some solvers. + IF( nd == 0 .AND. PRESENT( UseMeshDir ) ) THEN + IF( UseMeshDir ) THEN + OutputDirectory = TRIM(CurrentModel % Mesh % Name) + nd = LEN_TRIM(OutputDirectory) + END IF + END IF + + ! Use may have given part or all of the path in the filename. + ! This is not preferred, but we cannot trust the user. + nf = LEN_TRIM(Filename) + n = INDEX(Filename(1:nf),'/') + AbsPathInName = INDEX(FileName,':')>0 .OR. (Filename(1:1)=='/') & + .OR. (Filename(1:1)==Backslash) + + IF( nd > 0 .AND. .NOT. AbsPathInName ) THEN + ! Check that we have not given the path relative to home directory + ! because the code does not understand the meaning of tilde. + IF(nd>=2) THEN + IF( OutputDirectory(1:2) == '~/') THEN + CALL get_environment_variable('HOME',Str) + OutputDirectory = TRIM(Str)//'/'//OutputDirectory(3:nd) + nd = LEN_TRIM(OutputDirectory) + END IF + END IF + ! To be on the safe side create the directory. If it already exists no harm done. + ! Note that only one directory may be created. Hence if there is a path with many subdirectories + ! that will be a problem. Fortran does not have a standard ENQUIRE for directories hence + ! we just try to make it. + IF( DoDir ) THEN + CALL Info('SolverOutputDirectory','Creating directory: '//TRIM(OutputDirectory(1:nd)),Level=8) + CALL MakeDirectory( OutputDirectory(1:nd) // CHAR(0) ) + END IF + END IF + + ! In this case the filename includes also path and we remove it from there and + ! add it to the directory. + IF( n > 2 ) THEN + CALL Info('SolverOutputDirectory','Parcing path from filename: '//TRIM(Filename(1:n)),Level=10) + IF( AbsPathInName .OR. nd == 0) THEN + ! If the path is absolute then it overruns the given path! + OutputDirectory = Filename(1:n-1) + nd = n-1 + ELSE + ! If path is relative we add it to the OutputDirectory and take it away from Filename + OutputDirectory = OutputDirectory(1:nd)//'/'//Filename(1:n-1) + nd = nd + n + END IF + Filename = Filename(n+1:nf) + + IF( DoDir ) THEN + CALL Info('SolverOutputDirectory','Creating directory: '//TRIM(OutputDirectory(1:nd)),Level=8) + CALL MakeDirectory( OutputDirectory(1:nd) // CHAR(0) ) + END IF + END IF + + ! Finally, on request save each partitioning to different directory. + PartitioningSubDir = ListGetLogical( Solver % Values,'Output Partitioning Directory',Found) + IF(.NOT. Found ) THEN + PartitioningSubDir = ListGetLogical( CurrentModel % Simulation,'Output Partitioning Directory',Found) + END IF + IF( PartitioningSubDir ) THEN + OutputDirectory = TRIM(OutputDirectory)//'/np'//I2S(ParEnv % PEs) + nd = LEN_TRIM(OutputDirectory) + IF( DoDir ) THEN + CALL Info('SolverOutputDirectory','Creating directory: '//TRIM(OutputDirectory(1:nd)),Level=8) + CALL MakeDirectory( OutputDirectory(1:nd) // CHAR(0) ) + END IF + END IF + + END SUBROUTINE SolverOutputDirectory + !----------------------------------------------------------------------------- + + + !------------------------------------------------------------------------------ + !> Saves results in ascii format understood by the pre-/postprocessing software Gmsh. + !------------------------------------------------------------------------------ + SUBROUTINE SaveGmshOutput( Model,Solver,dt,Transient ) + !------------------------------------------------------------------------------ + USE Types + USE Lists + IMPLICIT NONE + !------------------------------------------------------------------------------ + TYPE(Solver_t) :: Solver + TYPE(Model_t) :: Model + REAL(KIND=dp) :: dt + LOGICAL :: Transient + !------------------------------------------------------------------------------ + ! Local variables + !------------------------------------------------------------------------------ + TYPE(Element_t),POINTER :: Element + INTEGER, POINTER :: Perm(:) + REAL(KIND=dp), POINTER :: Values(:),Values2(:),Values3(:) + REAL(KIND=dp) :: Vector(3), Time + COMPLEX(KIND=dp), POINTER :: CValues(:) + TYPE(Variable_t), POINTER :: Solution, TimeVariable + TYPE(ValueList_t), POINTER :: Params + TYPE(Mesh_t), POINTER :: Mesh + + LOGICAL :: Found, GotField, FileAppend, AlterTopology, MaskExists + LOGICAL :: EigenAnalysis = .FALSE., EigenActive, ComponentVector, Parallel + + INTEGER :: VisitedTimes = 0, ExtCount + INTEGER :: i,j,k,l,m,n,nsize,dim,dofs,ElmerCode, GmshCode,body_id, Vari, Rank, truedim + INTEGER :: Tag, NumberOfAllElements, BCOffSet + INTEGER, PARAMETER :: MaxElemCode = 827 + INTEGER :: ElmerToGmshType(MaxElemCode), GmshToElmerType(21), & + ElmerIndexes(27), GmshIndexes(27) + INTEGER, POINTER :: NodeIndexes(:) + + INTEGER, ALLOCATABLE :: NodePerm(:),DgPerm(:) + INTEGER, ALLOCATABLE, TARGET :: InvDgPerm(:), InvNodePerm(:) + LOGICAL, ALLOCATABLE :: ActiveElem(:) + LOGICAL :: NoPermutation, Numbering + INTEGER :: NumberOfGeomNodes, NumberOfDofNodes,NumberOfElements, ElemFirst, ElemLast,bc_id + INTEGER, POINTER :: InvFieldPerm(:), DGInvPerm(:) + + INTEGER, PARAMETER :: LENGTH = 1024 + CHARACTER(LEN=LENGTH) :: Txt, FieldName, CompName + CHARACTER(MAX_NAME_LEN) :: OutputFile + CHARACTER(:), ALLOCATABLE :: OutputDirectory + INTEGER :: GmshUnit + CHARACTER(*), PARAMETER :: Caller = 'SaveGmshOutput' + + SAVE VisitedTimes + + !------------------------------------------------------------------------------ + + CALL Info(Caller,'Saving results in Gmsh format') + + Mesh => Model % Mesh + Params => Solver % Values + Parallel = ( ParEnv % PEs > 1 ) + + ExtCount = ListGetInteger( Params,'Output Count',Found) + IF( Found ) THEN + VisitedTimes = ExtCount + ELSE + VisitedTimes = VisitedTimes + 1 + END IF + + Numbering = ListGetLogical( Params,'Filename Numbering',Found ) + IF(.NOT. Found) Numbering = .TRUE. + + GmshToElmerType = (/ 202, 303, 404, 504, 808, 706, 605, 203, 306, 409, & + 510, 827, 0, 0, 101, 408, 820, 715, 613, 0, 310 /) + ElmerToGmshType = 0 + + DO i=1,SIZE(GmshToElmerType) + j = GmshToElmerType(i) + IF( j > 0 ) ElmerToGmshType(j) = i + END DO + + EigenAnalysis = ListGetLogical( Params, 'Eigen Analysis', Found ) + FileAppend = ListGetLogical( Params,'File Append',Found) + IF(.NOT. Found) FileAppend = .TRUE. + AlterTopology = ListGetLogical( Params,'Alter Topology',Found) + + OutputFile = ListGetString( Params, 'Output File Name', Found ) + IF( Found ) THEN + IF(INDEX(OutputFile,'.') == 0) WRITE( OutputFile,'(A,A)') TRIM(OutputFile),".msh" + ELSE + OutputFile = 'Output.msh' + END IF + + CALL SolverOutputDirectory( Solver, OutputFile, OutputDirectory, UseMeshDir = .TRUE. ) + OutputFile = TRIM(OutputDirectory)// '/' //TRIM(OutputFile) + + !------------------------------------------------------------------------------ + ! Initialize stuff for masked saving + !------------------------------------------------------------------------------ + CALL GenerateSaveMask(Mesh,Params,Parallel,0,.FALSE.,& + NodePerm,ActiveElem,NumberOfGeomNodes,NumberOfElements,& + ElemFirst,ElemLast) + + IF( ParEnv % PEs > 1 ) THEN + IF( NumberOfElements == 0 ) THEN + CALL Info(Caller,'Nothing to save in partition: '//TRIM(I2S(ParEnv % MyPe)),Level=8) + RETURN + ELSE + OutputFile = TRIM(OutputFile)//'_'//I2S(ParEnv % PEs)//'np'//I2S(ParEnv % MyPe+1) + END IF + ELSE + IF( NumberOfElements == 0 ) THEN + CALL Warn(Caller,'Notging to save, this is suspicious') + RETURN + END IF + END IF + + CALL GenerateSavePermutation(Mesh,.FALSE.,.FALSE.,0,.FALSE.,& + ActiveElem,NumberOfGeomNodes,NoPermutation,NumberOfDofNodes,& + DgPerm,InvDgPerm,NodePerm,InvNodePerm) + + InvFieldPerm => InvNodePerm + + dim = CoordinateSystemDimension() + IF( VisitedTimes > 1 ) THEN + IF( AlterTopology ) THEN + IF( Numbering ) THEN + OutputFile = NextFreeFilename( OutputFile ) + END IF + CALL Info(Caller,'Writing mesh and data to a new file: '//TRIM(OutputFile)) + ELSE IF( FileAppend ) THEN + CALL Info(Caller,'Appending data to the same file: '//TRIM(OutputFile)) + OPEN(NEWUNIT=GmshUnit, FILE=OutputFile, POSITION='APPEND' ) + GOTO 10 + ELSE + IF( Numbering ) THEN + OutputFile = NextFreeFilename( OutputFile ) + END IF + CALL Info(Caller,'Writing data to a new file: '//TRIM(OutputFile)) + OPEN(NEWUNIT=GmshUnit, FILE=OutputFile ) + WRITE(GmshUnit,'(A)') '$MeshFormat' + WRITE(GmshUnit,'(A)') '2.0 0 8' + WRITE(GmshUnit,'(A)') '$EndMeshFormat' + GOTO 10 + END IF + END IF + + + ! Save the header + !------------------------------------------------- + CALL Info(Caller,'Saving results to file: '//TRIM(OutputFile)) + OPEN(NEWUNIT=GmshUnit, FILE=OutputFile ) + + WRITE(GmshUnit,'(A)') '$MeshFormat' + WRITE(GmshUnit,'(A)') '2.0 0 8' + WRITE(GmshUnit,'(A)') '$EndMeshFormat' + + + ! Save the mesh nodes + !------------------------------------------------- + CALL Info(Caller,'Writing the mesh nodes') + CALL WriteGmshNodes() + + ! Save the mesh elements + !------------------------------------------------- + CALL Info(Caller,'Writing the mesh elements') + CALL WriteGmshElements() + + ! With a mask the list of physical entities should be checked + !------------------------------------------------------------- + IF(.NOT. MaskExists ) THEN + ! CALL WritePhysicalNames() + END IF + +10 CONTINUE + + CALL Info(Caller,'Writing the nodal data') + CALL WriteGmshData() + + IF(.FALSE.) THEN + WRITE(GmshUnit,'(A)') '$ElementData' + WRITE(GmshUnit,'(A)') '$EndElementData' + END IF + + IF(.FALSE.) THEN + WRITE(GmshUnit,'(A)') '$ElementNodeData' + WRITE(GmshUnit,'(A)') '$EndElementNodeData' + END IF + + CLOSE(GmshUnit) + + IF(ALLOCATED(DgPerm)) DEALLOCATE(DgPerm) + IF(ALLOCATED(InvDgPerm)) DEALLOCATE(InvDgPerm) + IF(ALLOCATED(NodePerm)) DEALLOCATE(NodePerm) + IF(ALLOCATED(InvNodePerm)) DEALLOCATE(InvNodePerm) + IF(ALLOCATED(ActiveElem)) DEALLOCATE(ActiveElem) + + + CALL Info(Caller,'Gmsh output complete') + + CONTAINS + + SUBROUTINE WriteGmshNodes() + + nsize = NumberOfGeomNodes + + WRITE(GmshUnit,'(A)') '$Nodes' + WRITE(GmshUnit,'(I8)') nsize + DO i = 1, nsize + IF( NoPermutation ) THEN + j = i + ELSE + j = InvNodePerm(i) + END IF + + IF( dim == 3 ) THEN + WRITE(GmshUnit,'(I8,3ES16.7E3)') i,Mesh % Nodes % x(j),Mesh % Nodes % y(j), Mesh % Nodes % z(j) + ELSE + WRITE(GmshUnit,'(I8,2ES16.7E3,A)') i,Mesh % Nodes % x(j),Mesh % Nodes % y(j),' 0.0' + END IF + END DO + WRITE(GmshUnit,'(A)') '$EndNodes' + END SUBROUTINE WriteGmshNodes + + + SUBROUTINE WriteGmshElements() + + nsize = NumberOfElements + + BCOffSet = 100 + DO WHILE( BCOffset <= Model % NumberOfBodies ) + BCOffset = 10 * BCOffset + END DO + + WRITE(GmshUnit,'(A)') '$Elements' + WRITE(GmshUnit,'(I8)') nsize + + l = 0 + DO i = ElemFirst, ElemLast + IF(.NOT. ActiveElem(i) ) CYCLE + + l = l + 1 + Element => Mesh % Elements(i) + ElmerCode = Element % TYPE % ElementCode + + n = Element % Type % NumberOfNodes + IF( NoPermutation ) THEN + ElmerIndexes(1:n) = Element % NodeIndexes(1:n) + ELSE + ElmerIndexes(1:n) = NodePerm(Element % NodeIndexes(1:n)) + END IF + + GmshCode = ElmerToGmshType(ElmerCode) + IF( GmshCode == 0 ) THEN + CALL Warn(Caller,'Gmsh element index not found!') + CYCLE + END IF + + IF( i <= Model % NumberOfBulkElements ) THEN + Tag = Element % BodyId + ELSE + DO bc_id=1,Model % NumberOfBCs + IF ( Element % BoundaryInfo % Constraint == Model % BCs(bc_id) % Tag ) EXIT + END DO + Tag = bc_id + BCOffset + END IF + + WRITE(GmshUnit,'(I8,I3,I3,I5,I5)',ADVANCE='NO') l,GmshCode,2,Tag,Tag + k = MOD(ElmerCode,100) + + CALL ElmerToGmshIndex(ElmerCode,ElmerIndexes,GmshIndexes) + + DO j=1,k-1 + WRITE(GmshUnit,'(I8)',ADVANCE='NO') GmshIndexes(j) + END DO + WRITE(GmshUnit,'(I8)') GmshIndexes(k) + END DO + WRITE(GmshUnit,'(A)') '$EndElements' + END SUBROUTINE WriteGmshElements + + + SUBROUTINE WritePhysicalNames() + CALL Info(Caller,'Writing the physical entity names') + nsize = Model % NumberOfBodies + Model % NumberOfBCs + WRITE(GmshUnit,'(A)') '$PhysicalNames' + WRITE(GmshUnit,'(I8)') nsize + DO i=1,Model % NumberOfBodies + Txt = ListGetString( Model % Bodies(i) % Values,'Name',Found) + IF( Found ) THEN + WRITE(GmshUnit,'(I8,A)') i,'"'//TRIM(Txt)//'"' + ELSE + WRITE(GmshUnit,'(I8,A,I0,A)') i,'"Body ',i,'"' + END IF + END DO + DO i=1,Model % NumberOfBCs + Txt = ListGetString( Model % BCs(i) % Values,'Name',Found) + IF( Found ) THEN + WRITE(GmshUnit,'(I8,A)') i+BCOffset,'"'//TRIM(Txt)//'"' + ELSE + WRITE(GmshUnit,'(I8,A,I0,A)') i+BCOffset,'"Boundary Condition ',i,'"' + END IF + END DO + WRITE(GmshUnit,'(A)') '$EndPhysicalNames' + END SUBROUTINE WritePhysicalNames + + + ! In case of DG fields we average the fields on-the-fly to nodes. + !---------------------------------------------------------------- + SUBROUTINE CreateTemporalNodalField(Mesh,Solution,Revert) + TYPE(Mesh_t) :: Mesh + TYPE(Variable_t) :: Solution + LOGICAL, OPTIONAL :: revert + + REAL(KIND=dp), POINTER :: NodalVals(:), TmpVals(:) + INTEGER, POINTER :: NodalPerm(:), TmpPerm(:), NodalCnt(:) + INTEGER :: i,j,k,l,n,t,dofs,ElemFam + + SAVE NodalPerm, NodalVals, NodalCnt, TmpPerm, TmpVals + + IF( PRESENT( Revert ) ) THEN + IF( Revert ) THEN + DEALLOCATE( NodalVals, NodalPerm, NodalCnt ) + Solution % Perm => TmpPerm + Solution % Values => TmpVals + RETURN + END IF + END IF + + dofs = Solution % dofs + + n = Mesh % NumberOfNodes + ALLOCATE( NodalPerm(n), NodalCnt(n), NodalVals(n*dofs) ) + NodalPerm = 0 + NodalCnt = 0 + NodalVals = 0.0_dp + + DO t=1,Mesh % NumberOfBulkElements + Element => Mesh % Elements(t) + + ! This is just a quick hack to not consider those element in averaging that don't + ! even have one face on the active set of nodes. + IF( ALLOCATED(NodePerm) ) THEN + ElemFam = Element % TYPE % ElementCode / 100 + l = COUNT( NodePerm(Element % NodeIndexes ) > 0 ) + SELECT CASE(ElemFam) + CASE(3,4) + IF(l<2) CYCLE + CASE(5,6,7) + IF(l<3) CYCLE + CASE(8) + IF(l<4) CYCLE + END SELECT + END IF + + DO i=1,Element % TYPE % NumberOfNodes + j = Element % DGIndexes(i) + k = Element % NodeIndexes(i) + + NodalCnt(k) = NodalCnt(k) + 1 + NodalPerm(k) = k + + j = Solution % Perm(j) + IF(j==0) CYCLE + + DO l=1,dofs + NodalVals(dofs*(k-1)+l) = NodalVals(dofs*(k-1)+l) + Solution % Values(dofs*(j-1)+l) + END DO + END DO + END DO + + DO i=1,dofs + WHERE ( NodalCnt > 0 ) + NodalVals(i::dofs) = NodalVals(i::dofs) / NodalCnt + END WHERE + END DO + + TmpVals => Solution % Values + TmpPerm => Solution % Perm + + Solution % Perm => NodalPerm + Solution % Values => NodalVals + + + END SUBROUTINE CreateTemporalNodalField + + + + + SUBROUTINE WriteGmshData() + INTEGER :: ii + LOGICAL :: DgVar + + + ! Time is needed + !------------------------------------------------- + TimeVariable => VariableGet( Model % Variables, 'Time' ) + Time = TimeVariable % Values(1) + + ! Loop over different type of variables + !------------------------------------------------- + DO Rank = 0,2 + DO Vari = 1, 999 + IF(Rank==0) WRITE(Txt,'(A,I0)') 'Scalar Field ',Vari + IF(Rank==1) WRITE(Txt,'(A,I0)') 'Vector Field ',Vari + IF(Rank==2) WRITE(Txt,'(A,I0)') 'Tensor Field ',Vari + + FieldName = ListGetString( Params, TRIM(Txt), Found ) + IF(.NOT. Found) EXIT + IF( Rank == 2) THEN + CALL Warn(Caller,'Not implemented yet for tensors!') + CYCLE + END IF + + ComponentVector = .FALSE. + Solution => VariableGet( Mesh % Variables, FieldName ) + DGVar = .FALSE. + + IF(ASSOCIATED(Solution)) THEN + DGVar = ( Solution % TYPE == Variable_on_nodes_on_elements ) + IF(DgVar) CALL CreateTemporalNodalField(Mesh,Solution) + + Values => Solution % Values + Perm => Solution % Perm + dofs = Solution % DOFs + ELSE + IF( Rank == 1 ) THEN + Solution => VariableGet( Mesh % Variables, FieldName//' 1' ) + IF( ASSOCIATED( Solution ) ) THEN + ComponentVector = .TRUE. + Values => Solution % Values + Perm => Solution % Perm + dofs = 1 + Solution => VariableGet( Mesh % Variables, FieldName//' 2' ) + IF( ASSOCIATED(Solution)) THEN + Values2 => Solution % Values + dofs = 2 + END IF + Solution => VariableGet( Mesh % Variables, FieldName//' 3' ) + IF( ASSOCIATED(Solution)) THEN + Values3 => Solution % Values + dofs = 3 + END IF + END IF + END IF + IF( .NOT. ASSOCIATED(Solution)) THEN + CALL Warn(Caller,'Variable not present: '//TRIM(FieldName)) + CYCLE + END IF + END IF + + CALL Info(Caller,'Saving nodal variable: '//TRIM(FieldName),Level=12) + + IF( ASSOCIATED(Solution % EigenVectors) ) THEN + CALL Warn(Caller,'Eigenvectors related to field: '//TRIM(FieldName)) + CALL Warn(Caller,'Eigenvectors saving yet not supported') + END IF + + truedim = MIN(dofs, dim) + nsize = NumberOfGeomNodes + + WRITE(GmshUnit,'(A)') '$NodeData' + WRITE(GmshUnit,'(A)') '1' + WRITE(GmshUnit,'(A)') '"'//TRIM(FieldName)//'"' + WRITE(GmshUnit,'(A)') '1' + + ! Gmsh starts steady state indexes from zero, hence deductions by one + IF( Transient ) THEN + WRITE(GmshUnit,'(ES16.7E3)') Time + ELSE + WRITE(GmshUnit,'(ES16.7E3)') Time - 1.0_dp + END IF + WRITE(GmshUnit,'(A)') '3' + WRITE(GmshUnit,'(I8)') VisitedTimes-1 + IF(Rank == 0) THEN + WRITE(GmshUnit,'(A)') '1' + ELSE IF(Rank == 1) THEN + WRITE(GmshUnit,'(A)') '3' + ELSE + WRITE(GmshUnit,'(A)') '9' + END IF + WRITE(GmshUnit,'(I8)') nsize + + DO ii = 1, NumberOfGeomNodes + IF( NoPermutation ) THEN + i = ii + ELSE + i = InvFieldPerm(ii) + END IF + + IF( ASSOCIATED( Perm ) ) THEN + j = Perm(i) + ELSE + j = i + END IF + + IF( Rank == 0 ) THEN + WRITE(GmshUnit,'(I8,ES16.7E3)') ii,Values(j) + ELSE IF(Rank == 1) THEN + IF( j == 0 ) THEN + WRITE(GmshUnit,'(I8,A)') ii,' 0.0 0.0 0.0' + ELSE IF( ComponentVector ) THEN + IF( truedim == 2 ) THEN + WRITE(GmshUnit,'(I8,2ES16.7E3,A)') ii,& + Values(j),Values2(j),' 0.0' + ELSE + WRITE(GmshUnit,'(I8,3ES16.7E3)') ii,& + Values(j),Values2(j),Values3(j) + END IF + ELSE + IF( truedim == 2 ) THEN + WRITE(GmshUnit,'(I8,2ES16.7E3,A)') ii,& + Values(dofs*(j-1)+1),Values(dofs*(j-1)+2),' 0.0' + ELSE + WRITE(GmshUnit,'(I8,3ES16.7E3)') ii,& + Values(dofs*(j-1)+1),Values(dofs*(j-1)+2),Values(dofs*(j-1)+3) + END IF + END IF + END IF + END DO + WRITE(GmshUnit,'(A)') '$EndNodeData' + + IF(DgVar) CALL CreateTemporalNodalField(Mesh,Solution,Revert=.TRUE.) + + END DO + END DO + END SUBROUTINE WriteGmshData + +!------------------------------------------------------------------------------ + END SUBROUTINE SaveGmshOutput +!------------------------------------------------------------------------------ + END MODULE SaveUtils diff --git a/fem/src/SolverUtils.F90 b/fem/src/SolverUtils.F90 index f9f97770b7..868e7ee6d5 100644 --- a/fem/src/SolverUtils.F90 +++ b/fem/src/SolverUtils.F90 @@ -24828,117 +24828,6 @@ SUBROUTINE ReleaseProjectors(Model, Solver) END SUBROUTINE ReleaseProjectors - !> Defines and potentially creates output directory. - !> The output directory may given in different ways, and even be part of the - !> filename, or be relative to home directory. We try to parse every possible - !> scenario here that user might have in mind. - !----------------------------------------------------------------------------- - SUBROUTINE SolverOutputDirectory( Solver, Filename, OutputDirectory, & - MakeDir, UseMeshDir ) - - TYPE(Solver_t) :: Solver - LOGICAL, OPTIONAL :: MakeDir, UseMeshDir - CHARACTER(*) :: Filename - CHARACTER(:), ALLOCATABLE :: OutputDirectory - - LOGICAL :: Found, AbsPathInName, DoDir, PartitioningSubDir - INTEGER :: nd, nf, n - CHARACTER(LEN=MAX_NAME_LEN) :: Str - - IF( PRESENT( MakeDir ) ) THEN - DoDir = MakeDir - ELSE - DoDir = ( Solver % TimesVisited == 0 ) .AND. ( ParEnv % MyPe == 0 ) - END IF - - ! Output directory is obtained in order - ! 1) solver section - ! 2) simulation section - ! 3) header section - OutputDirectory = ListGetString( Solver % Values,'Output Directory',Found) - IF(.NOT. Found) OutputDirectory = ListGetString( CurrentModel % Simulation,& - 'Output Directory',Found) - - IF(.NOT. Found) OutputDirectory = TRIM(OutputPath) - nd = LEN_TRIM(OutputDirectory) - - ! If the path is just working directory then that is not an excude - ! to not use the mesh name, or directory that comes with the filename - IF(.NOT. Found .AND. nd == 1 .AND. OutputDirectory(1:1)=='.') nd = 0 - - ! If requested by the optional parameter use the mesh directory when - ! no results directory given. This is an old convection used in some solvers. - IF( nd == 0 .AND. PRESENT( UseMeshDir ) ) THEN - IF( UseMeshDir ) THEN - OutputDirectory = TRIM(CurrentModel % Mesh % Name) - nd = LEN_TRIM(OutputDirectory) - END IF - END IF - - ! Use may have given part or all of the path in the filename. - ! This is not preferred, but we cannot trust the user. - nf = LEN_TRIM(Filename) - n = INDEX(Filename(1:nf),'/') - AbsPathInName = INDEX(FileName,':')>0 .OR. (Filename(1:1)=='/') & - .OR. (Filename(1:1)==Backslash) - - IF( nd > 0 .AND. .NOT. AbsPathInName ) THEN - ! Check that we have not given the path relative to home directory - ! because the code does not understand the meaning of tilde. - IF(nd>=2) THEN - IF( OutputDirectory(1:2) == '~/') THEN - CALL get_environment_variable('HOME',Str) - OutputDirectory = TRIM(Str)//'/'//OutputDirectory(3:nd) - nd = LEN_TRIM(OutputDirectory) - END IF - END IF - ! To be on the safe side create the directory. If it already exists no harm done. - ! Note that only one directory may be created. Hence if there is a path with many subdirectories - ! that will be a problem. Fortran does not have a standard ENQUIRE for directories hence - ! we just try to make it. - IF( DoDir ) THEN - CALL Info('SolverOutputDirectory','Creating directory: '//TRIM(OutputDirectory(1:nd)),Level=8) - CALL MakeDirectory( OutputDirectory(1:nd) // CHAR(0) ) - END IF - END IF - - ! In this case the filename includes also path and we remove it from there and - ! add it to the directory. - IF( n > 2 ) THEN - CALL Info('SolverOutputDirectory','Parcing path from filename: '//TRIM(Filename(1:n)),Level=10) - IF( AbsPathInName .OR. nd == 0) THEN - ! If the path is absolute then it overruns the given path! - OutputDirectory = Filename(1:n-1) - nd = n-1 - ELSE - ! If path is relative we add it to the OutputDirectory and take it away from Filename - OutputDirectory = OutputDirectory(1:nd)//'/'//Filename(1:n-1) - nd = nd + n - END IF - Filename = Filename(n+1:nf) - - IF( DoDir ) THEN - CALL Info('SolverOutputDirectory','Creating directory: '//TRIM(OutputDirectory(1:nd)),Level=8) - CALL MakeDirectory( OutputDirectory(1:nd) // CHAR(0) ) - END IF - END IF - - ! Finally, on request save each partitioning to different directory. - PartitioningSubDir = ListGetLogical( Solver % Values,'Output Partitioning Directory',Found) - IF(.NOT. Found ) THEN - PartitioningSubDir = ListGetLogical( CurrentModel % Simulation,'Output Partitioning Directory',Found) - END IF - IF( PartitioningSubDir ) THEN - OutputDirectory = TRIM(OutputDirectory)//'/np'//I2S(ParEnv % PEs) - nd = LEN_TRIM(OutputDirectory) - IF( DoDir ) THEN - CALL Info('SolverOutputDirectory','Creating directory: '//TRIM(OutputDirectory(1:nd)),Level=8) - CALL MakeDirectory( OutputDirectory(1:nd) // CHAR(0) ) - END IF - END IF - - END SUBROUTINE SolverOutputDirectory - !----------------------------------------------------------------------------- ! This routine changes the IP field to DG field just while the results are being written. !--------------------------------------------------------------------------------------- diff --git a/fem/src/modules/GmshOutputReader.F90 b/fem/src/modules/GmshOutputReader.F90 index 12ab5041f1..54a0879952 100644 --- a/fem/src/modules/GmshOutputReader.F90 +++ b/fem/src/modules/GmshOutputReader.F90 @@ -26,8 +26,9 @@ !> \ingroup Solvers !------------------------------------------------------------------------------ SUBROUTINE GmshOutputReader( Model,Solver,dt,TransientSimulation ) -!------------------------------------------------------------------------------ +!------------------------------------------------------------------------------ USE DefUtils + USE SaveUtils IMPLICIT NONE !------------------------------------------------------------------------------ TYPE(Solver_t) :: Solver diff --git a/fem/src/modules/ResultOutputSolve/GidOutputSolver.F90 b/fem/src/modules/ResultOutputSolve/GidOutputSolver.F90 index 370061f148..262b065834 100644 --- a/fem/src/modules/ResultOutputSolve/GidOutputSolver.F90 +++ b/fem/src/modules/ResultOutputSolve/GidOutputSolver.F90 @@ -28,6 +28,7 @@ SUBROUTINE GiDOutputSolver( Model,Solver,dt,TransientSimulation ) !------------------------------------------------------------------------------ USE DefUtils + USE SaveUtils, ONLY : SolverOutputDirectory IMPLICIT NONE !------------------------------------------------------------------------------ TYPE(Solver_t) :: Solver diff --git a/fem/src/modules/ResultOutputSolve/GmshOutputSolver.F90 b/fem/src/modules/ResultOutputSolve/GmshOutputSolver.F90 deleted file mode 100644 index 23d244640f..0000000000 --- a/fem/src/modules/ResultOutputSolve/GmshOutputSolver.F90 +++ /dev/null @@ -1,545 +0,0 @@ -!/*****************************************************************************/ -! * -! * Elmer, A Finite Element Software for Multiphysical Problems -! * -! * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland -! * -! * This program is free software; you can redistribute it and/or -! * modify it under the terms of the GNU General Public License -! * as published by the Free Software Foundation; either version 2 -! * of the License, or (at your option) any later version. -! * -! * This program is distributed in the hope that it will be useful, -! * but WITHOUT ANY WARRANTY; without even the implied warranty of -! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -! * GNU General Public License for more details. -! * -! * You should have received a copy of the GNU General Public License -! * along with this program (in file fem/GPL-2); if not, write to the -! * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, -! * Boston, MA 02110-1301, USA. -! * -! *****************************************************************************/ - -!------------------------------------------------------------------------------ -!> Saves results in ascii format understood by the pre-/postprocessing software Gmsh. -!> \ingroup Solvers -!------------------------------------------------------------------------------ -SUBROUTINE GmshOutputSolver( Model,Solver,dt,TransientSimulation ) -!------------------------------------------------------------------------------ - USE DefUtils - USE SaveUtils - IMPLICIT NONE -!------------------------------------------------------------------------------ - TYPE(Solver_t) :: Solver - TYPE(Model_t) :: Model - REAL(KIND=dp) :: dt - LOGICAL :: TransientSimulation -!------------------------------------------------------------------------------ -! Local variables -!------------------------------------------------------------------------------ - TYPE(Element_t),POINTER :: Element - INTEGER, POINTER :: Perm(:) - REAL(KIND=dp), POINTER :: Values(:),Values2(:),Values3(:) - REAL(KIND=dp) :: Vector(3), Time - COMPLEX(KIND=dp), POINTER :: CValues(:) - TYPE(Variable_t), POINTER :: Solution, TimeVariable - TYPE(ValueList_t), POINTER :: Params - TYPE(Mesh_t), POINTER :: Mesh - - LOGICAL :: Found, GotField, FileAppend, AlterTopology, MaskExists - LOGICAL :: EigenAnalysis = .FALSE., EigenActive, ComponentVector, Parallel - - INTEGER :: VisitedTimes = 0, ExtCount - INTEGER :: i,j,k,l,m,n,nsize,dim,dofs,ElmerCode, GmshCode,body_id, Vari, Rank, truedim - INTEGER :: Tag, NumberOfAllElements, BCOffSet - INTEGER, PARAMETER :: MaxElemCode = 827 - INTEGER :: ElmerToGmshType(MaxElemCode), GmshToElmerType(21), & - ElmerIndexes(27), GmshIndexes(27) - INTEGER, POINTER :: NodeIndexes(:) - - INTEGER, ALLOCATABLE :: NodePerm(:),DgPerm(:) - INTEGER, ALLOCATABLE, TARGET :: InvDgPerm(:), InvNodePerm(:) - LOGICAL, ALLOCATABLE :: ActiveElem(:) - LOGICAL :: NoPermutation - INTEGER :: NumberOfGeomNodes, NumberOfDofNodes,NumberOfElements, ElemFirst, ElemLast - INTEGER, POINTER :: InvFieldPerm(:), DGInvPerm(:) - - INTEGER, PARAMETER :: LENGTH = 1024 - CHARACTER(LEN=LENGTH) :: Txt, FieldName, CompName - CHARACTER(MAX_NAME_LEN) :: OutputFile - CHARACTER(:), ALLOCATABLE :: OutputDirectory - INTEGER :: GmshUnit - CHARACTER(*), PARAMETER :: Caller = 'GmshOutputSolver' - - SAVE VisitedTimes - -!------------------------------------------------------------------------------ - - CALL Info(Caller,'Saving results in Gmsh format') - - Mesh => Model % Mesh - Params => Solver % Values - Parallel = ( ParEnv % PEs > 1 ) - - ExtCount = ListGetInteger( Params,'Output Count',Found) - IF( Found ) THEN - VisitedTimes = ExtCount - ELSE - VisitedTimes = VisitedTimes + 1 - END IF - - GmshToElmerType = (/ 202, 303, 404, 504, 808, 706, 605, 203, 306, 409, & - 510, 827, 0, 0, 101, 408, 820, 715, 613, 0, 310 /) - ElmerToGmshType = 0 - - DO i=1,SIZE(GmshToElmerType) - j = GmshToElmerType(i) - IF( j > 0 ) ElmerToGmshType(j) = i - END DO - - EigenAnalysis = GetLogical( Params, 'Eigen Analysis', Found ) - FileAppend = GetLogical( Params,'File Append',Found) - IF(.NOT. Found) FileAppend = .TRUE. - AlterTopology = GetLogical( Params,'Alter Topology',Found) - - OutputFile = GetString( Solver % Values, 'Output File Name', Found ) - IF( Found ) THEN - IF(INDEX(OutputFile,'.') == 0) WRITE( OutputFile,'(A,A)') TRIM(OutputFile),".msh" - ELSE - OutputFile = 'Output.msh' - END IF - - CALL SolverOutputDirectory( Solver, OutputFile, OutputDirectory, UseMeshDir = .TRUE. ) - OutputFile = TRIM(OutputDirectory)// '/' //TRIM(OutputFile) - - !------------------------------------------------------------------------------ - ! Initialize stuff for masked saving - !------------------------------------------------------------------------------ - CALL GenerateSaveMask(Mesh,Params,Parallel,0,.FALSE.,& - NodePerm,ActiveElem,NumberOfGeomNodes,NumberOfElements,& - ElemFirst,ElemLast) - - IF( ParEnv % PEs > 1 ) THEN - IF( NumberOfElements == 0 ) THEN - CALL Info(Caller,'Nothing to save in partition: '//TRIM(I2S(ParEnv % MyPe)),Level=8) - RETURN - ELSE - OutputFile = TRIM(OutputFile)//'_'//I2S(ParEnv % PEs)//'np'//I2S(ParEnv % MyPe+1) - END IF - ELSE - IF( NumberOfElements == 0 ) THEN - CALL Warn(Caller,'Notging to save, this is suspicious') - RETURN - END IF - END IF - - CALL GenerateSavePermutation(Mesh,.FALSE.,.FALSE.,0,.FALSE.,& - ActiveElem,NumberOfGeomNodes,NoPermutation,NumberOfDofNodes,& - DgPerm,InvDgPerm,NodePerm,InvNodePerm) - - InvFieldPerm => InvNodePerm - - dim = CoordinateSystemDimension() - IF( VisitedTimes > 1 ) THEN - IF( AlterTopology ) THEN - OutputFile = NextFreeFilename( OutputFile ) - CALL Info(Caller,'Writing mesh and data to a new file: '//TRIM(OutputFile)) - ELSE IF( FileAppend ) THEN - CALL Info(Caller,'Appending data to the same file: '//TRIM(OutputFile)) - OPEN(NEWUNIT=GmshUnit, FILE=OutputFile, POSITION='APPEND' ) - GOTO 10 - ELSE - OutputFile = NextFreeFilename( OutputFile ) - CALL Info(Caller,'Writing data to a new file: '//TRIM(OutputFile)) - OPEN(NEWUNIT=GmshUnit, FILE=OutputFile ) - WRITE(GmshUnit,'(A)') '$MeshFormat' - WRITE(GmshUnit,'(A)') '2.0 0 8' - WRITE(GmshUnit,'(A)') '$EndMeshFormat' - GOTO 10 - END IF - END IF - - - ! Save the header - !------------------------------------------------- - CALL Info('GmshOutputSolver','Saving results to file: '//TRIM(OutputFile)) - OPEN(NEWUNIT=GmshUnit, FILE=OutputFile ) - - WRITE(GmshUnit,'(A)') '$MeshFormat' - WRITE(GmshUnit,'(A)') '2.0 0 8' - WRITE(GmshUnit,'(A)') '$EndMeshFormat' - - - ! Save the mesh nodes - !------------------------------------------------- - CALL Info(Caller,'Writing the mesh nodes') - CALL WriteGmshNodes() - - ! Save the mesh elements - !------------------------------------------------- - CALL Info(Caller,'Writing the mesh elements') - CALL WriteGmshElements() - - ! With a mask the list of physical entities should be checked - !------------------------------------------------------------- - IF(.NOT. MaskExists ) THEN -! CALL WritePhysicalNames() - END IF - -10 CONTINUE - - CALL Info(Caller,'Writing the nodal data') - CALL WriteGmshData() - - IF(.FALSE.) THEN - WRITE(GmshUnit,'(A)') '$ElementData' - WRITE(GmshUnit,'(A)') '$EndElementData' - END IF - - IF(.FALSE.) THEN - WRITE(GmshUnit,'(A)') '$ElementNodeData' - WRITE(GmshUnit,'(A)') '$EndElementNodeData' - END IF - - CLOSE(GmshUnit) - - - - IF(ALLOCATED(DgPerm)) DEALLOCATE(DgPerm) - IF(ALLOCATED(InvDgPerm)) DEALLOCATE(InvDgPerm) - IF(ALLOCATED(NodePerm)) DEALLOCATE(NodePerm) - IF(ALLOCATED(InvNodePerm)) DEALLOCATE(InvNodePerm) - IF(ALLOCATED(ActiveElem)) DEALLOCATE(ActiveElem) - - - CALL Info(Caller,'Gmsh output complete') - -CONTAINS - - SUBROUTINE WriteGmshNodes() - - nsize = NumberOfGeomNodes - - WRITE(GmshUnit,'(A)') '$Nodes' - WRITE(GmshUnit,'(I8)') nsize - DO i = 1, nsize - IF( NoPermutation ) THEN - j = i - ELSE - j = InvNodePerm(i) - END IF - - IF( dim == 3 ) THEN - WRITE(GmshUnit,'(I8,3ES16.7E3)') i,Mesh % Nodes % x(j),Mesh % Nodes % y(j), Mesh % Nodes % z(j) - ELSE - WRITE(GmshUnit,'(I8,2ES16.7E3,A)') i,Mesh % Nodes % x(j),Mesh % Nodes % y(j),' 0.0' - END IF - END DO - WRITE(GmshUnit,'(A)') '$EndNodes' - END SUBROUTINE WriteGmshNodes - - - SUBROUTINE WriteGmshElements() - - nsize = NumberOfElements - - BCOffSet = 100 - DO WHILE( BCOffset <= Model % NumberOfBodies ) - BCOffset = 10 * BCOffset - END DO - - WRITE(GmshUnit,'(A)') '$Elements' - WRITE(GmshUnit,'(I8)') nsize - - l = 0 - DO i = ElemFirst, ElemLast - IF(.NOT. ActiveElem(i) ) CYCLE - - l = l + 1 - Element => Mesh % Elements(i) - ElmerCode = Element % TYPE % ElementCode - - n = Element % Type % NumberOfNodes - IF( NoPermutation ) THEN - ElmerIndexes(1:n) = Element % NodeIndexes(1:n) - ELSE - ElmerIndexes(1:n) = NodePerm(Element % NodeIndexes(1:n)) - END IF - - GmshCode = ElmerToGmshType(ElmerCode) - IF( GmshCode == 0 ) THEN - CALL Warn(Caller,'Gmsh element index not found!') - CYCLE - END IF - - IF( i <= Model % NumberOfBulkElements ) THEN - Tag = Element % BodyId - ELSE - Tag = GetBCId( Element ) + BCOffset - END IF - - WRITE(GmshUnit,'(I8,I3,I3,I5,I5)',ADVANCE='NO') l,GmshCode,2,Tag,Tag - k = MOD(ElmerCode,100) - - CALL ElmerToGmshIndex(ElmerCode,ElmerIndexes,GmshIndexes) - - DO j=1,k-1 - WRITE(GmshUnit,'(I8)',ADVANCE='NO') GmshIndexes(j) - END DO - WRITE(GmshUnit,'(I8)') GmshIndexes(k) - END DO - WRITE(GmshUnit,'(A)') '$EndElements' - END SUBROUTINE WriteGmshElements - - - SUBROUTINE WritePhysicalNames() - CALL Info(Caller,'Writing the physical entity names') - nsize = Model % NumberOfBodies + Model % NumberOfBCs - WRITE(GmshUnit,'(A)') '$PhysicalNames' - WRITE(GmshUnit,'(I8)') nsize - DO i=1,Model % NumberOfBodies - Txt = ListGetString( Model % Bodies(i) % Values,'Name',Found) - IF( Found ) THEN - WRITE(GmshUnit,'(I8,A)') i,'"'//TRIM(Txt)//'"' - ELSE - WRITE(GmshUnit,'(I8,A,I0,A)') i,'"Body ',i,'"' - END IF - END DO - DO i=1,Model % NumberOfBCs - Txt = ListGetString( Model % BCs(i) % Values,'Name',Found) - IF( Found ) THEN - WRITE(GmshUnit,'(I8,A)') i+BCOffset,'"'//TRIM(Txt)//'"' - ELSE - WRITE(GmshUnit,'(I8,A,I0,A)') i+BCOffset,'"Boundary Condition ',i,'"' - END IF - END DO - WRITE(GmshUnit,'(A)') '$EndPhysicalNames' - END SUBROUTINE WritePhysicalNames - - - ! In case of DG fields we average the fields on-the-fly to nodes. - !---------------------------------------------------------------- - SUBROUTINE CreateTemporalNodalField(Mesh,Solution,Revert) - TYPE(Mesh_t) :: Mesh - TYPE(Variable_t) :: Solution - LOGICAL, OPTIONAL :: revert - - REAL(KIND=dp), POINTER :: NodalVals(:), TmpVals(:) - INTEGER, POINTER :: NodalPerm(:), TmpPerm(:), NodalCnt(:) - INTEGER :: i,j,k,l,n,t,dofs,ElemFam - - SAVE NodalPerm, NodalVals, NodalCnt, TmpPerm, TmpVals - - IF( PRESENT( Revert ) ) THEN - IF( Revert ) THEN - DEALLOCATE( NodalVals, NodalPerm, NodalCnt ) - Solution % Perm => TmpPerm - Solution % Values => TmpVals - RETURN - END IF - END IF - - dofs = Solution % dofs - - n = Mesh % NumberOfNodes - ALLOCATE( NodalPerm(n), NodalCnt(n), NodalVals(n*dofs) ) - NodalPerm = 0 - NodalCnt = 0 - NodalVals = 0.0_dp - - DO t=1,Mesh % NumberOfBulkElements - Element => Mesh % Elements(t) - - ! This is just a quick hack to not consider those element in averaging that don't - ! even have one face on the active set of nodes. - IF( ALLOCATED(NodePerm) ) THEN - ElemFam = Element % TYPE % ElementCode / 100 - l = COUNT( NodePerm(Element % NodeIndexes ) > 0 ) - SELECT CASE(ElemFam) - CASE(3,4) - IF(l<2) CYCLE - CASE(5,6,7) - IF(l<3) CYCLE - CASE(8) - IF(l<4) CYCLE - END SELECT - END IF - - DO i=1,Element % TYPE % NumberOfNodes - j = Element % DGIndexes(i) - k = Element % NodeIndexes(i) - - NodalCnt(k) = NodalCnt(k) + 1 - NodalPerm(k) = k - - j = Solution % Perm(j) - IF(j==0) CYCLE - - DO l=1,dofs - NodalVals(dofs*(k-1)+l) = NodalVals(dofs*(k-1)+l) + Solution % Values(dofs*(j-1)+l) - END DO - END DO - END DO - - DO i=1,dofs - WHERE ( NodalCnt > 0 ) - NodalVals(i::dofs) = NodalVals(i::dofs) / NodalCnt - END WHERE - END DO - - TmpVals => Solution % Values - TmpPerm => Solution % Perm - - Solution % Perm => NodalPerm - Solution % Values => NodalVals - - - END SUBROUTINE CreateTemporalNodalField - - - - - SUBROUTINE WriteGmshData() - INTEGER :: ii - LOGICAL :: DgVar - - - ! Time is needed - !------------------------------------------------- - TimeVariable => VariableGet( Model % Variables, 'Time' ) - Time = TimeVariable % Values(1) - - ! Loop over different type of variables - !------------------------------------------------- - CALL Info(Caller,'Writing the nodal data') - DO Rank = 0,2 - DO Vari = 1, 999 - IF(Rank==0) WRITE(Txt,'(A,I0)') 'Scalar Field ',Vari - IF(Rank==1) WRITE(Txt,'(A,I0)') 'Vector Field ',Vari - IF(Rank==2) WRITE(Txt,'(A,I0)') 'Tensor Field ',Vari - - FieldName = GetString( Solver % Values, TRIM(Txt), Found ) - IF(.NOT. Found) EXIT - IF( Rank == 2) THEN - CALL Warn(Caller,'Not implemented yet for tensors!') - CYCLE - END IF - - ComponentVector = .FALSE. - Solution => VariableGet( Mesh % Variables, FieldName ) - DGVar = .FALSE. - - IF(ASSOCIATED(Solution)) THEN - DGVar = ( Solution % TYPE == Variable_on_nodes_on_elements ) - IF(DgVar) CALL CreateTemporalNodalField(Mesh,Solution) - - Values => Solution % Values - Perm => Solution % Perm - dofs = Solution % DOFs - ELSE - IF( Rank == 1 ) THEN - Solution => VariableGet( Mesh % Variables, FieldName//' 1' ) - IF( ASSOCIATED( Solution ) ) THEN - ComponentVector = .TRUE. - Values => Solution % Values - Perm => Solution % Perm - dofs = 1 - Solution => VariableGet( Mesh % Variables, FieldName//' 2' ) - IF( ASSOCIATED(Solution)) THEN - Values2 => Solution % Values - dofs = 2 - END IF - Solution => VariableGet( Mesh % Variables, FieldName//' 3' ) - IF( ASSOCIATED(Solution)) THEN - Values3 => Solution % Values - dofs = 3 - END IF - END IF - END IF - IF( .NOT. ASSOCIATED(Solution)) THEN - CALL Warn('GmshOutputSolver','Variable not present: '//TRIM(FieldName)) - CYCLE - END IF - END IF - IF( ASSOCIATED(Solution % EigenVectors) ) THEN - CALL Warn(Caller,'Eigenvectors related to field: '//TRIM(FieldName)) - CALL Warn(Caller,'Eigenvectors saving yet not supported') - END IF - - truedim = MIN(dofs, dim) - nsize = NumberOfGeomNodes - - WRITE(GmshUnit,'(A)') '$NodeData' - WRITE(GmshUnit,'(A)') '1' - WRITE(GmshUnit,'(A)') '"'//TRIM(FieldName)//'"' - WRITE(GmshUnit,'(A)') '1' - - ! Gmsh starts steady state indexes from zero, hence deductions by one - IF( TransientSimulation ) THEN - WRITE(GmshUnit,'(ES16.7E3)') Time - ELSE - WRITE(GmshUnit,'(ES16.7E3)') Time - 1.0_dp - END IF - WRITE(GmshUnit,'(A)') '3' - WRITE(GmshUnit,'(I8)') VisitedTimes-1 - IF(Rank == 0) THEN - WRITE(GmshUnit,'(A)') '1' - ELSE IF(Rank == 1) THEN - WRITE(GmshUnit,'(A)') '3' - ELSE - WRITE(GmshUnit,'(A)') '9' - END IF - WRITE(GmshUnit,'(I8)') nsize - - DO ii = 1, NumberOfGeomNodes - IF( NoPermutation ) THEN - i = ii - ELSE - i = InvFieldPerm(ii) - END IF - - IF( ASSOCIATED( Perm ) ) THEN - j = Perm(i) - ELSE - j = i - END IF - - IF( Rank == 0 ) THEN - WRITE(GmshUnit,'(I8,ES16.7E3)') ii,Values(j) - ELSE IF(Rank == 1) THEN - IF( j == 0 ) THEN - WRITE(GmshUnit,'(I8,A)') ii,' 0.0 0.0 0.0' - ELSE IF( ComponentVector ) THEN - IF( truedim == 2 ) THEN - WRITE(GmshUnit,'(I8,2ES16.7E3,A)') ii,& - Values(j),Values2(j),' 0.0' - ELSE - WRITE(GmshUnit,'(I8,3ES16.7E3)') ii,& - Values(j),Values2(j),Values3(j) - END IF - ELSE - IF( truedim == 2 ) THEN - WRITE(GmshUnit,'(I8,2ES16.7E3,A)') ii,& - Values(dofs*(j-1)+1),Values(dofs*(j-1)+2),' 0.0' - ELSE - WRITE(GmshUnit,'(I8,3ES16.7E3)') ii,& - Values(dofs*(j-1)+1),Values(dofs*(j-1)+2),Values(dofs*(j-1)+3) - END IF - END IF - END IF - END DO - WRITE(GmshUnit,'(A)') '$EndNodeData' - - IF(DgVar) CALL CreateTemporalNodalField(Mesh,Solution,Revert=.TRUE.) - - END DO - END DO - END SUBROUTINE WriteGmshData - - - -!------------------------------------------------------------------------------ -END SUBROUTINE GmshOutputSolver -!------------------------------------------------------------------------------ - diff --git a/fem/src/modules/ResultOutputSolve/ResultOutputSolver.F90 b/fem/src/modules/ResultOutputSolve/ResultOutputSolver.F90 index 30b7a5f891..41ca2412bb 100644 --- a/fem/src/modules/ResultOutputSolve/ResultOutputSolver.F90 +++ b/fem/src/modules/ResultOutputSolve/ResultOutputSolver.F90 @@ -271,7 +271,8 @@ END SUBROUTINE ElmerPostOutputSolver IF( SaveGmsh ) THEN CALL Info( Caller,'Saving in gmsh 2.0 (.msh) format' ) CALL ListAddInteger( Params,'Output Count',OutputCount(2)) - CALL GmshOutputSolver( Model,Solver,dt,TransientSimulation ) + ! For other call uses this recides in SaveUtils. + CALL SaveGmshOutput( Model,Solver,dt,TransientSimulation ) END IF IF( SaveVTK ) THEN CALL Info( Caller,'Saving in legacy VTK (.vtk) format' ) diff --git a/fem/src/modules/ResultOutputSolve/VtkOutputSolver.F90 b/fem/src/modules/ResultOutputSolve/VtkOutputSolver.F90 index d69885bbf3..f1183b9d0d 100644 --- a/fem/src/modules/ResultOutputSolve/VtkOutputSolver.F90 +++ b/fem/src/modules/ResultOutputSolve/VtkOutputSolver.F90 @@ -28,6 +28,7 @@ MODULE VtkLegacyFile USE MeshUtils USE ElementDescription + USE SaveUtils, ONLY : SolverOutputDirectory IMPLICIT NONE ! PRIVATE