! ! * ! * Elmer, A Finite Element Software for Multiphysical Problems ! * ! * Copyright 1st April 1995 - , CSC - Scientific Computing 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. ! * ! *****************************************************************************/ ! * ! * $Id: ResultOutputSolve.src,v 1.21 2007/08/15 11:42:44 jpr Exp $ ! ***************************************************************************** SUBROUTINE ResultOutputSolver( Model,Solver,dt,TransientSimulation ) !DEC$ATTRIBUTES DLLEXPORT :: ResultOutputSolver !------------------------------------------------------------------------------ !****************************************************************************** ! ! Exports data to other suitable postprocessing software. ! Currently supported formats are GiD, VTK legacy, and Open DX. ! ! ARGUMENTS: ! ! TYPE(Model_t) :: Model, ! INPUT: All model information (mesh, materials, BCs, etc...) ! ! TYPE(Solver_t) :: Solver ! INPUT: Linear & nonlinear equation solver options ! ! REAL(KIND=dp) :: dt, ! INPUT: Timestep size for time dependent simulations ! ! LOGICAL :: TransientSimulation ! INPUT: Steady state or transient simulation ! !****************************************************************************** USE DefUtils IMPLICIT NONE TYPE(Solver_t) :: Solver TYPE(Model_t) :: Model REAL(KIND=dp) :: dt LOGICAL :: TransientSimulation, SaveGid, SaveVTK, SaveOpenDx, SaveGmsh CHARACTER(10) :: OutputFormat LOGICAL :: Found SaveGid = ListGetLogical(Solver % Values,'Gid Format',Found) SaveGmsh = ListGetLogical(Solver % Values,'Gmsh Format',Found) SaveVTK = ListGetLogical(Solver % Values,'VTK Format',Found) SaveOpenDx = ListGetLogical(Solver % Values,'Dx Format',Found) OutputFormat = GetString( Solver % Values, 'Output Format', Found ) IF(Found) THEN IF( OutputFormat == "gid" )THEN SaveGid = .TRUE. ELSE IF( OutputFormat == "vtk" )THEN SaveVTK = .TRUE. ELSE IF( OutputFormat == "dx" )THEN SaveOpenDx = .TRUE. ELSE IF( OutputFormat == "gmsh" )THEN SaveGmsh = .TRUE. ELSE CALL Warn( 'ResultOutputSolver', & 'Unknown output format "' // TRIM(OutputFormat) // '"' ) CALL Warn( 'ResultOutputSolver', & 'Available formats are "GiD", "VTK" and "DX"' ) END IF END IF IF(.NOT. (SaveGid .OR. SaveGmsh .OR. SaveVTK .OR. SaveOpenDx)) THEN CALL Warn('ResultOutputSolver','No output format was defined') RETURN END IF IF( SaveGid ) CALL GiDOutputSolver( Model,Solver,dt,TransientSimulation ) IF( SaveGmsh ) CALL GmshOutputSolver( Model,Solver,dt,TransientSimulation ) IF( SaveVTK ) CALL VtkOutputSolver( Model,Solver,dt,TransientSimulation ) IF( SaveOpenDx ) CALL DXOutputSolver( Model,Solver,dt,TransientSimulation ) END SUBROUTINE ResultOutputSolver !------------------------------------------------------------------------------ SUBROUTINE GiDOutputSolver( Model,Solver,dt,TransientSimulation ) !------------------------------------------------------------------------------ USE DefUtils IMPLICIT NONE !------------------------------------------------------------------------------ TYPE(Solver_t) :: Solver TYPE(Model_t) :: Model REAL(KIND=dp) :: dt LOGICAL :: TransientSimulation !------------------------------------------------------------------------------ ! Local variables !------------------------------------------------------------------------------ TYPE(Element_t),POINTER :: Element TYPE(Variable_t), POINTER :: Solution INTEGER, POINTER :: Perm(:) REAL(KIND=dp), POINTER :: Values(:) COMPLEX(KIND=dp), POINTER :: CValues(:) TYPE(Variable_t), POINTER :: TimeVariable TYPE(ValueList_t), POINTER :: SolverParams LOGICAL :: Found, CoordinatesWritten = .FALSE. LOGICAL :: FirstTimeStep = .TRUE., EigenAnalysis = .FALSE. INTEGER :: i,j,k,m,n,dim, Code, body_id, ElementCounter, Nloop, Loop INTEGER :: tensorComponents INTEGER, PARAMETER :: MaxElemCode = 1000 INTEGER :: ListElemTypes(MaxElemCode) CHARACTER(LEN=1024) :: OutputFile, ResFile, MshFile, Txt, Family, & ScalarFieldName, VectorFieldName, TensorFieldName, CompName CHARACTER(LEN=1024) :: Txt2, Txt3 INTEGER :: PyramidMap613(14,4), PyramidMap605(2,4) INTEGER :: WedgeMap706(3,4), WedgeMap715(21,4) SAVE FirstTimeStep !------------------------------------------------------------------------------ PyramidMap605(1,:) = (/ 3, 5, 4, 1 /) PyramidMap605(2,:) = (/ 3, 5, 2, 1 /) PyramidMap613(1,:) = (/ 7, 8, 3, 12 /) PyramidMap613(2,:) = (/ 10, 11, 12, 5 /) PyramidMap613(3,:) = (/ 10, 13, 12, 5 /) PyramidMap613(4,:) = (/ 9, 10, 13, 11 /) PyramidMap613(5,:) = (/ 9, 13, 11, 12 /) PyramidMap613(6,:) = (/ 9, 10, 11, 1 /) PyramidMap613(7,:) = (/ 9, 6, 11, 1 /) PyramidMap613(8,:) = (/ 9, 6, 11, 12 /) PyramidMap613(9,:) = (/ 9, 8, 12, 4 /) PyramidMap613(10,:) = (/ 9, 13, 12, 4 /) PyramidMap613(11,:) = (/ 7, 9, 8, 12 /) PyramidMap613(12,:) = (/ 7, 9, 6, 12 /) PyramidMap613(13,:) = (/ 7, 6, 11, 12 /) PyramidMap613(14,:) = (/ 7, 6, 11, 2 /) WedgeMap706(1,:) = (/ 5, 4, 3, 1 /) WedgeMap706(2,:) = (/ 5, 3, 2, 1 /) WedgeMap706(3,:) = (/ 5, 6, 4, 3 /) WedgeMap715(1,:) = (/ 10, 11, 5, 2 /) WedgeMap715(2,:) = (/ 12, 11, 6, 3 /) WedgeMap715(3,:) = (/ 12, 10, 4, 1 /) WedgeMap715(4,:) = (/ 7, 8, 11, 2 /) WedgeMap715(5,:) = (/ 7, 10, 11, 2 /) WedgeMap715(6,:) = (/ 13, 14, 11, 5 /) WedgeMap715(7,:) = (/ 13, 10, 11, 5 /) WedgeMap715(8,:) = (/ 9, 10, 8, 11 /) WedgeMap715(9,:) = (/ 9, 7, 10, 8 /) WedgeMap715(10,:) = (/ 9, 12, 10, 11 /) WedgeMap715(11,:) = (/ 9, 12, 10, 1 /) WedgeMap715(12,:) = (/ 9, 7, 10, 1 /) WedgeMap715(13,:) = (/ 9, 8, 11, 3 /) WedgeMap715(14,:) = (/ 9, 12, 11, 3 /) WedgeMap715(15,:) = (/ 15, 12, 10, 11 /) WedgeMap715(16,:) = (/ 15, 12, 10, 4 /) WedgeMap715(17,:) = (/ 15, 10, 14, 11 /) WedgeMap715(18,:) = (/ 15, 13, 10, 14 /) WedgeMap715(19,:) = (/ 15, 13, 10, 4 /) WedgeMap715(20,:) = (/ 15, 14, 11, 6 /) WedgeMap715(21,:) = (/ 15, 12, 11, 6 /) SolverParams => GetSolverParams() EigenAnalysis = GetLogical( SolverParams, 'Eigen Analysis', Found ) OutputFile = GetString( Solver % Values, 'Output File Name', Found ) IF( .NOT.Found ) THEN PRINT *,'Output File Name undefined.' OutputFile = 'Output' END IF PRINT *,'Output File Name = ',TRIM(OutputFile) WRITE(ResFile,'(A,A)') TRIM(OutputFile),'.flavia.res' WRITE(MshFile,'(A,A)') TRIM(OutputFile),'.flavia.msh' PRINT *,'res-file = ',TRIM(ResFile) PRINT *,'msh-file = ',TRIM(MshFile) ! Write the GiD msh-file: !------------------------ dim = CoordinateSystemDimension() IF( CoordinatesWritten ) GOTO 10 OPEN( UNIT=10, FILE=MshFile ) ! First check how many element types are involved in the analysis: !----------------------------------------------------------------- ListElemTypes = 0 ElementCounter = 0 ! DO i = 1, Solver % NumberOfActiveElements ! Element => GetActiveElement(i) DO i = 1, Model % NumberOfBulkElements !+ Model % NumberOfBoundaryElements Element => Model % Mesh % Elements(i) Code = Element % TYPE % ElementCode ListElemTypes(Code) = ListElemTypes(Code)+1 END DO PRINT *,'Total number of elements =',SUM(ListElemTypes) ! Write the different element types in different blocks: !------------------------------------------------------- DO i = 1,MaxElemCode IF(ListElemTypes(i) == 0) CYCLE PRINT *,ListElemTypes(i),'elements of type',i n = MOD(i,100) IF( INT(i/100) == 1 ) Family = 'Point' IF( INT(i/100) == 2 ) Family = 'Linear' IF( INT(i/100) == 3 ) Family = 'Triangle' IF( INT(i/100) == 4 ) Family = 'Quadrilateral' IF( INT(i/100) == 5 ) Family = 'Tetrahedra' IF( INT(i/100) == 6 ) THEN Family = 'Tetrahedra' ! PYRAMIDS WILL BE SPLITTED n = 4 ! INTO LINEAR TETRAHEDRA END IF IF( INT(i/100) == 7 ) THEN Family = 'Tetrahedra' ! WEDGES WILL BE SPLITTED n = 4 ! INTO LINEAR TETRAHEDRA END IF IF( INT(i/100) == 8 ) Family = 'Hexahedra' IF( n < 10 ) THEN WRITE(Txt,'(A,I1,A,A,A,I2)') 'MESH "Elmer Mesh" dimension ',& dim,' ElemType ', TRIM(Family),' Nnode', n ELSE WRITE(Txt,'(A,I1,A,A,A,I3)') 'MESH "Elmer Mesh" dimension ',& dim,' ElemType ', TRIM(Family),' Nnode', n END IF WRITE(10,'(A)') TRIM(Txt) ! Write all node coordinates in the first block: !----------------------------------------------- IF( .NOT.CoordinatesWritten ) THEN WRITE(10,'(A)') 'Coordinates' DO j = 1, Model % Mesh % NumberOfNodes WRITE(10,'(I6,3E16.6E3)') j, & Model % Mesh % Nodes % x(j), & Model % Mesh % Nodes % y(j), & Model % Mesh % Nodes % z(j) END DO WRITE(10,'(A)') 'end coordinates' WRITE(10,'(A)') ' ' CoordinatesWritten = .TRUE. END IF ! Write the element connectivity tables: !--------------------------------------- WRITE(10,'(A)') 'Elements' ! DO j = 1, Solver % NumberOfActiveElements ! Element => GetActiveElement( j ) DO j = 1, Model % NumberOfBulkElements !+ Model % NumberOfBoundaryElements Element => Model % Mesh % Elements(j) Code = Element % TYPE % ElementCode IF( Code /= i ) CYCLE body_id = Element % BodyId IF( Code == 613 ) THEN ! 13 noded pyramids will be splitted into 14 linear tetraheda !------------------------------------------------------------ DO m = 1,14 ElementCounter = ElementCounter + 1 WRITE(10,'(100I10)') ElementCounter, & Element % NodeIndexes(PyramidMap613(m,:)), body_id END DO ELSEIF( Code == 605 ) THEN ! 5 noded pyramids will be splitted into 2 linear tetraheda !---------------------------------------------------------- DO m = 1,2 ElementCounter = ElementCounter + 1 WRITE(10,'(100I10)') ElementCounter, & Element % NodeIndexes(PyramidMap605(m,:)), body_id END DO ELSEIF( Code == 706 ) THEN ! 6 noded wedges will be splitted into 3 linear tetraheda !--------------------------------------------------------- DO m = 1,3 ElementCounter = ElementCounter + 1 WRITE(10,'(100I10)') ElementCounter, & Element % NodeIndexes(WedgeMap706(m,:)), body_id END DO ELSEIF( Code == 715 ) THEN ! 15 noded wedges will be splitted into 21 linear tetraheda !---------------------------------------------------------- DO m = 1,21 ElementCounter = ElementCounter + 1 WRITE(10,'(100I10)') ElementCounter, & Element % NodeIndexes(WedgeMap715(m,:)), body_id END DO ELSE ! Standard elements are understood by GiD as such !------------------------------------------------ ElementCounter = ElementCounter + 1 WRITE(10,'(100I10)') ElementCounter, Element % NodeIndexes, body_id END IF END DO WRITE(10,'(A)') 'end elements' WRITE(10,'(A)') ' ' END DO 10 CONTINUE ! Write the GiD res-file: !------------------------ IF( TransientSimulation .AND. FirstTimeStep ) THEN OPEN(UNIT=10, FILE=ResFile ) WRITE(10,'(A)') 'GiD Post Result File 1.0' FirstTimeStep = .FALSE. ELSEIF( TransientSimulation .AND. (.NOT.FirstTimeStep ) ) THEN OPEN(UNIT=10, FILE=ResFile, POSITION='APPEND' ) ELSE OPEN(UNIT=10, FILE=ResFile ) WRITE(10,'(A)') 'GiD Post Result File 1.0' END IF Nloop = 1 IF( EigenAnalysis ) THEN Nloop = GetInteger( Solver % Values, 'Eigen System Values', Found ) IF( .NOT.Found ) Nloop = 1 END IF DO Loop = 1, Nloop PRINT *,'------------' ! First scalar fields: !---------------------- DO i = 1, 99 IF( i<10 ) THEN WRITE(Txt,'(A,I2)') 'Scalar Field',i ELSE WRITE(Txt,'(A,I3)') 'Scalar Field',i END IF ScalarFieldName = GetString( Solver % Values, TRIM(Txt), Found ) IF(.NOT. Found) EXIT Solution => VariableGet( Solver % Mesh % Variables, ScalarFieldName ) IF( .NOT.ASSOCIATED( Solution ) ) THEN PRINT *,'Scalar field "',TRIM(ScalarFieldName),'" not found' ELSE PRINT *,'Scarar field',i,'= "',TRIM(ScalarFieldName),'"' Perm => Solution % Perm IF( .NOT.EigenAnalysis ) THEN Values => Solution % Values ELSE Cvalues => Solution % EigenVectors(Loop,:) END IF IF( TransientSimulation ) THEN TimeVariable => VariableGet( Solver % Mesh % Variables, 'Time' ) PRINT *,'Current time=',TimeVariable % Values(1) WRITE(10,'(A,A,A,E16.6E3,A)') 'Result "',& TRIM(ScalarFieldName),'" "Transient analysis" ', & TimeVariable % Values(1) ,' Scalar OnNodes' ELSE IF( .NOT.EigenAnalysis ) THEN WRITE(10,'(A,A,A,I2,A)') 'Result "',& TRIM(ScalarFieldName),'" "Steady analysis"',Loop,' Scalar OnNodes' ELSE WRITE(10,'(A,A,A,I2,A)') 'Result "',& TRIM(ScalarFieldName),'" "Eigen analysis"',Loop,' Scalar OnNodes' END IF END IF WRITE(10,'(A,A,A)') 'ComponentNames "',TRIM(ScalarFieldName),'"' WRITE(10,'(A)') 'Values' DO j = 1, Model % Mesh % NumberOfNodes k = Perm(j) IF( .NOT.EigenAnalysis ) THEN WRITE(10,'(I6,E16.6E3)') j, Values(k) ELSE WRITE(10,'(I6,E16.6E3)') j, REAL(CValues(k)) END IF END DO WRITE(10,'(A)') 'end values' WRITE(10,'(A)') ' ' END IF END DO ! Then vector fields: !-------------------- DO i = 1, 99 IF( i<10 ) THEN WRITE(Txt,'(A,I2)') 'Vector Field',i ELSE WRITE(Txt,'(A,I3)') 'Vector Field',i END IF VectorFieldName = GetString( Solver % Values, TRIM(Txt), Found ) IF(.NOT. Found) EXIT PRINT *,'Vector field',i,'= "',TRIM(VectorFieldName),'"' IF( TransientSimulation ) THEN TimeVariable => VariableGet( Solver % Mesh % Variables, 'Time' ) PRINT *,'Current time=',TimeVariable % Values(1) WRITE(10,'(A,A,A,E16.6E3,A)') 'Result "',& TRIM(VectorFieldName),'" "Transient analysis" ', & TimeVariable % Values(1) ,' Vector OnNodes' ELSE IF( .NOT.EigenAnalysis ) THEN WRITE(10,'(A,A,A,I2,A)') 'Result "',& TRIM(VectorFieldName),'" "Steady analysis"',Loop,' Vector OnNodes' ELSE WRITE(10,'(A,A,A,I2,A)') 'Result "',& TRIM(VectorFieldName),'" "Eigen analysis"',Loop,' Vector OnNodes' END IF END IF WRITE(Txt,'(A)') 'ComponentNames ' DO j = 1, dim IF(j VariableGet( Solver % Mesh % Variables, TRIM(Txt) ) IF( .NOT.ASSOCIATED( Solution ) ) THEN PRINT *,'Vector field component',k,' not found' ELSE Perm => Solution % Perm IF( .NOT.EigenAnalysis ) THEN Values => Solution % Values WRITE(Txt2,'(A,E16.6E3)') TRIM(Txt2), Values( Perm(j) ) ELSE CValues => Solution % Eigenvectors(Loop,:) WRITE(Txt2,'(A,E16.6E3)') TRIM(Txt2), REAL(CValues( Perm(j) ) ) END IF END IF END DO WRITE(10,'(A)') TRIM(Txt2) END DO WRITE(10,'(A)') 'end values' END DO ! Finally tensor fields: !----------------------- DO i = 1, 99 IF( i<10 ) THEN WRITE(Txt,'(A,I2)') 'Tensor Field',i ELSE WRITE(Txt,'(A,I3)') 'Tensor Field',i END IF TensorFieldName = GetString( Solver % Values, TRIM(Txt), Found ) IF(.NOT. Found) EXIT PRINT *,'Tensor field',i,'= "',TRIM(TensorFieldName),'"' IF( TransientSimulation ) THEN TimeVariable => VariableGet( Solver % Mesh % Variables, 'Time' ) PRINT *,'Current time=',TimeVariable % Values(1) WRITE(10,'(A,A,A,E16.6E3,A)') 'Result "',& TRIM(TensorFieldName),'" "Transient analysis" ', & TimeVariable % Values(1) ,' Matrix OnNodes' ELSE IF( .NOT.EigenAnalysis ) THEN WRITE(10,'(A,A,A,I2,A)') 'Result "',& TRIM(TensorFieldName),'" "Steady analysis"',Loop,' Matrix OnNodes' ELSE WRITE(10,'(A,A,A,I2,A)') 'Result "',& TRIM(TensorFieldName),'" "Eigen analysis"',Loop,' Matrix OnNodes' END IF END IF WRITE(Txt,'(A)') 'ComponentNames ' IF( dim == 2 ) THEN TensorComponents = 3 ELSE TensorComponents = 6 END IF DO j = 1, TensorComponents IF(j VariableGet( Solver % Mesh % Variables, TRIM(Txt) ) IF( .NOT.ASSOCIATED( Solution ) ) THEN PRINT *,'Tensor field component',k,' not found' ELSE Perm => Solution % Perm !Values => Solution % Values !WRITE(Txt2,'(A,E16.6E3)') TRIM(Txt2), Values( Perm(j) ) IF( .NOT.EigenAnalysis ) THEN Values => Solution % Values WRITE(Txt2,'(A,E16.6E3)') TRIM(Txt2), Values( Perm(j) ) ELSE CValues => Solution % Eigenvectors(Loop,:) WRITE(Txt2,'(A,E16.6E3)') TRIM(Txt2), REAL(CValues( Perm(j) ) ) END IF END IF END DO WRITE(10,'(A)') TRIM(Txt2) END DO WRITE(10,'(A)') 'end values' END DO END DO ! Nloop CLOSE(10) PRINT *,'Output complete.' !------------------------------------------------------------------------------ END SUBROUTINE GiDOutputSolver !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ SUBROUTINE GmshOutputSolver( Model,Solver,dt,TransientSimulation ) !------------------------------------------------------------------------------ USE DefUtils IMPLICIT NONE !------------------------------------------------------------------------------ TYPE(Solver_t) :: Solver TYPE(Model_t) :: Model REAL(KIND=dp) :: dt LOGICAL :: TransientSimulation !------------------------------------------------------------------------------ ! Local variables !------------------------------------------------------------------------------ TYPE(Element_t),POINTER :: Element TYPE(Variable_t), POINTER :: Solution INTEGER, POINTER :: Perm(:) REAL(KIND=dp), POINTER :: Values(:) COMPLEX(KIND=dp), POINTER :: CValues(:) TYPE(Variable_t), POINTER :: TimeVariable TYPE(ValueList_t), POINTER :: SolverParams LOGICAL :: Found, GotField, FileAppend, AlterTopology LOGICAL :: EigenAnalysis = .FALSE., EigenActive INTEGER :: VisitedTimes = 0 INTEGER :: i,j,k,l,m,n,dim,dofs,Code, body_id, Nloop, Loop, Vari, Rank INTEGER :: ElemNo, NumberOfAllElements INTEGER, PARAMETER :: MaxElemCode = 827 INTEGER, TARGET :: PermIndexes(27) INTEGER :: ListElemTypes(MaxElemCode) INTEGER, ALLOCATABLE :: ElementOrder(:) INTEGER, POINTER :: NodeIndexes(:) INTEGER, PARAMETER :: LENGTH = 1024 CHARACTER(LEN=LENGTH) :: OutputFile, Txt, FieldName, CompName CHARACTER(LEN=6), PARAMETER :: BodyIdField = 'bodyid' INTEGER :: ElemTypeOrder(15) SAVE VisitedTimes, ListElemTypes, ElementOrder, ElemNo !------------------------------------------------------------------------------ ElemTypeOrder(:) = (/ 101, 202, 303, 404, 504, 808, 706, 605, & 203, 306, 408, 510, 820, 715, 613 /) CALL Info('GmshOutputSolver','Saving results in Gmsh format') SolverParams => GetSolverParams() EigenAnalysis = GetLogical( SolverParams, 'Eigen Analysis', Found ) FileAppend = GetLogical( SolverParams,'File Append',Found) IF(.NOT. Found) FileAppend = .TRUE. AlterTopology = GetLogical( SolverParams,'Alter Topology',Found) OutputFile = GetString( Solver % Values, 'Output File Name', Found ) IF( .NOT.Found ) THEN CALL Warn('GmshOutputSolver','Output File Name undefined') OutputFile = 'Output.gm' END IF dim = CoordinateSystemDimension() IF( VisitedTimes > 0 ) THEN IF ( AlterTopology ) THEN OPEN(UNIT=10, FILE=OutputFile, POSITION='APPEND' ) GOTO 5 END IF IF ( FileAppend ) THEN OPEN(UNIT=10, FILE=OutputFile, POSITION='APPEND' ) GOTO 10 END IF END IF ! Save the one-time information !------------------------------------------------- PRINT *,'Output File Name = ',TRIM(OutputFile) OPEN(UNIT=10, FILE=OutputFile ) WRITE(10,'(A)') '$PostFormat' WRITE(10,'(A)') '1.4 0 8' WRITE(10,'(A)') '$EndPostFormat' 5 CONTINUE IF( VisitedTimes > 0) DEALLOCATE(ElementOrder) NumberOfAllElements = Model % NumberOfBulkElements + Model % NumberOfBoundaryElements ALLOCATE(ElementOrder(NumberOfAllElements)) ElementOrder = 0 ListElemTypes = 0 DO i = 1, NumberOfAllElements Element => Model % Mesh % Elements(i) Code = Element % TYPE % ElementCode ListElemTypes(Code) = ListElemTypes(Code)+1 END DO ! PRINT *,'Total number of elements =',SUM(ListElemTypes) ! Order the elements so that they are saved in the order defined by the ! ElemTypeOrder vector. !---------------------------------------------------------------------- ElemNo = 0 DO j = 1,SIZE(ElemTypeOrder) Code = ElemTypeOrder(j) IF( ListElemTypes(code) == 0) CYCLE DO i = 1, NumberOfAllElements Element => Model % Mesh % Elements(i) IF(Code /= Element % TYPE % ElementCode) CYCLE ElemNo = ElemNo + 1 ElementOrder(ElemNo) = i END DO END DO ! Check that there is a field to work with !------------------------------------------------- GotField = .FALSE. DO Rank = 0,2 DO Vari = 1, 99 IF( Vari < 10 ) THEN IF(Rank==0) WRITE(Txt,'(A,I2)') 'Scalar Field',Vari IF(Rank==1) WRITE(Txt,'(A,I2)') 'Vector Field',Vari IF(Rank==2) WRITE(Txt,'(A,I2)') 'Tensor Field',Vari ELSE IF(Rank==0) WRITE(Txt,'(A,I3)') 'Scalar Field',Vari IF(Rank==1) WRITE(Txt,'(A,I3)') 'Vector Field',Vari IF(Rank==2) WRITE(Txt,'(A,I3)') 'Tensor Field',Vari END IF FieldName = GetString( Solver % Values, TRIM(Txt), Found ) IF(.NOT. Found) EXIT IF( Rank == 2) THEN CALL Warn('GmshOutputSolver','Not implemented yet for tensors!') CYCLE END IF Solution => VariableGet( Solver % Mesh % Variables, FieldName ) IF(ASSOCIATED(Solution)) THEN GotField = .TRUE. ELSE IF (FieldName==BodyIdField .AND. Rank==0) THEN GotField = .TRUE. END IF END DO END DO IF(.NOT. GotField) THEN CALL Warn('GmshOutputSolver','No variables defined for saving') GOTO 20 END IF 10 CONTINUE DO Rank = 0,2 DO Vari = 1, 99 IF( Vari < 10 ) THEN IF(Rank==0) WRITE(Txt,'(A,I2)') 'Scalar Field',Vari IF(Rank==1) WRITE(Txt,'(A,I2)') 'Vector Field',Vari IF(Rank==2) WRITE(Txt,'(A,I2)') 'Tensor Field',Vari ELSE IF(Rank==0) WRITE(Txt,'(A,I3)') 'Scalar Field',Vari IF(Rank==1) WRITE(Txt,'(A,I3)') 'Vector Field',Vari IF(Rank==2) WRITE(Txt,'(A,I3)') 'Tensor Field',Vari END IF FieldName = GetString( Solver % Values, TRIM(Txt), Found ) IF(.NOT. Found) EXIT IF(Rank == 2) THEN CALL Warn('GmshOutputSolver','Do the tensors') EXIT END IF Solution => VariableGet( Solver % Mesh % Variables, FieldName ) IF(.NOT. ASSOCIATED(Solution)) THEN IF(FieldName == BodyIdField .AND. Rank==0) THEN ! Gmsh Doesn't seem to be capable of extracting a physical ! volume. One workaround is to create a scalar field with ! nodal values equal to the body id of the element. CALL Warn('GmshOutputSolver', 'Implement Body Id output!') ELSE WRITE(Txt, '(A,A)') 'Nonexistent variable: ',TRIM(FieldName) CALL Warn('GmshOutputSolver', Txt) END IF CYCLE END IF ! Eliminate possible blancos that are problematic for gmsh !--------------------------------------------------------- CALL BlancoSubst(FieldName) Perm => Solution % Perm dofs = Solution % DOFs EigenActive = .FALSE. IF( EigenAnalysis ) EigenActive = ASSOCIATED(Solution % EigenVectors) IF(EigenActive) THEN PRINT *,'Do the eigen values' Cvalues => Solution % EigenVectors(Loop,:) Nloop = GetInteger( Solver % Values, 'Eigen System Values', Found ) IF( .NOT.Found ) Nloop = 1 ELSE Values => Solution % Values END IF ! Compute the number of active elements for this particular dof !-------------------------------------------------------------- ListElemTypes = 0 DO i = 1, ElemNo Element => Model % Mesh % Elements(ElementOrder(i)) Code = Element % TYPE % ElementCode NodeIndexes => Element % NodeIndexes IF( ANY(Perm(NodeIndexes) == 0)) CYCLE ListElemTypes(Code) = ListElemTypes(Code) + 1 END DO IF( SUM(ListElemTypes) == 0) CYCLE ! Now write the view related to the current dof !-------------------------------------------------------------- PRINT *,'Saving View for variable: ',TRIM(FieldName) WRITE(10,'(A)') '$View' WRITE(10,'(A,I8)') TRIM(FieldName),1 DO i = 1,SIZE(ElemTypeOrder) Code = ElemTypeOrder(i) n = ListElemTypes(code) IF(Rank == 0) THEN WRITE(10,'(I8,I8,I8)') n,0,0 ELSE IF(Rank == 1) THEN WRITE(10,'(I8,I8,I8)') 0,n,0 ELSE WRITE(10,'(I8,I8,I8)') 0,0,n END IF ! IF(n > 0) PRINT *,'there are',n,'elements of TYPE',Code END DO ! Saving of names data currently omitted WRITE(10,'(I8,I8,I8,I8)') 0,0,0,0 IF( TransientSimulation ) THEN TimeVariable => VariableGet( Solver % Mesh % Variables, 'Time' ) PRINT *,'Current time=',TimeVariable % Values(1) WRITE(10,'(E16.6E3)') TimeVariable % Values(1) ELSE IF(.NOT. EigenAnalysis) THEN WRITE(10,'(E16.6E3)') 1.0*VisitedTimes ELSE DO i=1,Nloop WRITE(10,'(E16.6E3)') 1.0*i END DO END IF DO i = 1, ElemNo Element => Model % Mesh % Elements(ElementOrder(i)) Code = Element % TYPE % ElementCode NodeIndexes => Element % NodeIndexes IF( ANY(Perm(NodeIndexes) == 0)) CYCLE n = element % TYPE % NumberOfNodes CALL ElmerToGmshIndex(Code,n,NodeIndexes) WRITE(10,'(100E18.8)') Model % Mesh % Nodes % x(NodeIndexes(1:n)) WRITE(10,'(100E18.8)') Model % Mesh % Nodes % y(NodeIndexes(1:n)) WRITE(10,'(100E18.8)') Model % Mesh % Nodes % z(NodeIndexes(1:n)) IF(Rank == 0) THEN DO k = 1, MOD(Code,100) ! NOTE: FMT currently allows writing of three digit exponents ! This limitiation is also present below for ! 2 and 3 dimensional fields. WRITE(10,'(E18.8E3)') Values(Perm(NodeIndexes(k))) END DO ELSE DO k = 1, MOD(Code,100) ! Always save three components for vectors! DO l = 1, 2 WRITE(10,'(E18.8E3)', ADVANCE='NO') Values(dofs*(Perm(NodeIndexes(k))-1)+l) END DO ! The third component usually only makes sense if it is related to the third coordinate IF(dofs >= 3 .AND. dim >= 3) THEN WRITE(10,'(E18.8E3)') Values(dofs*(Perm(NodeIndexes(k))-1)+3) ELSE WRITE(10,'(A)') ' 0.0' END IF END DO END IF END DO WRITE(10,'(A)') '$EndView' END DO END DO 20 CONTINUE CLOSE(10) VisitedTimes = VisitedTimes + 1 CALL Info('GmshOutputSolver','Gmsh output complete') CONTAINS SUBROUTINE BlancoSubst(variable_name) ! Subroutine by Kurt Baarman IMPLICIT NONE CHARACTER(LEN=LENGTH) :: variable_name CHARACTER(LEN=LENGTH), PARAMETER :: & NAMESET = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789_' CHARACTER, PARAMETER :: GOOD = '_', BAD = ' ' INTEGER :: i, last = LENGTH last = SCAN(variable_name, TRIM(NAMESET), .TRUE.) !WRITE (*,*) 'Original name: ',TRIM(variable_name) i = SCAN(variable_name(1:last), BAD) DO WHILE (i /= 0)! .AND. i < last) variable_name(i:i) = GOOD i = SCAN(variable_name(1:last), BAD) END DO !WRITE (*,*) 'Altered name: ',TRIM(variable_name) END SUBROUTINE BlancoSubst SUBROUTINE ElmerToGmshIndex(Code,n,NodeIndexes) INTEGER :: Code INTEGER, POINTER :: NodeIndexes(:) INTEGER i,n SELECT CASE( Code ) CASE(510) PermIndexes(1:n) = NodeIndexes(1:n) PermIndexes(10) = NodeIndexes(9) PermIndexes(9) = NodeIndexes(10) NodeIndexes => PermIndexes CASE(306) ! There seems to be conflicting data whether this is needed or not PermIndexes(1:n) = NodeIndexes(1:n) PermIndexes(5) = NodeIndexes(6) PermIndexes(6) = NodeIndexes(5) NodeIndexes => PermIndexes CASE DEFAULT END SELECT END SUBROUTINE ElmerToGmshIndex !------------------------------------------------------------------------------ END SUBROUTINE GmshOutputSolver !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ ! The Vtk output solver begins here !------------------------------------------------------------------------------ MODULE VtkLegacyFile USE MeshUtils USE ElementDescription IMPLICIT NONE ! PRIVATE SAVE PUBLIC :: WriteVtkLegacyFile TYPE :: VtkCell_t INTEGER :: nNodes ! FIXME: POINTER -> ALLOCATABLE when more compilers support it INTEGER, POINTER :: NodeIndex(:) END TYPE VtkCell_t INTEGER, PARAMETER :: VtkUnit = 58 CONTAINS SUBROUTINE WriteVtkLegacyFile( VtkFile, Model, SubtractDisp ) CHARACTER(LEN=*), INTENT(IN) :: VtkFile TYPE(Model_t) :: Model LOGICAL, INTENT(IN) :: SubtractDisp TYPE(Variable_t), POINTER :: Var,Var1 CHARACTER(LEN=512) :: str !LOGICAL :: FreeSurfaceFlag INTEGER :: i,j,k OPEN( UNIT=VtkUnit, FILE=VtkFile, STATUS='UNKNOWN' ) !FreeSurfaceFlag = FreeSurface( Model ) CALL WriteGrid( VtkUnit, Model, SubtractDisp ) WRITE( VtkUnit,'("POINT_DATA ",I0)' ) Model % NumberOfNodes Var => Model % Variables DO WHILE( ASSOCIATED( Var ) ) IF ( .NOT.Var % Output ) THEN Var => Var % Next CYCLE END IF SELECT CASE( Var % Name ) CASE( 'time','timestep','timestep size', 'timestep interval' ) CASE( 'mesh update' ) Var1 => Model % Variables DO WHILE( ASSOCIATED( Var1 ) ) IF ( TRIM( Var1 % Name ) == 'displacement' ) EXIT Var1 => Var1 % Next END DO IF ( .NOT.ASSOCIATED( Var1 ) ) THEN CALL WriteVector("Mesh.Update", Var, Model % NumberOfNodes,& 3, VtkUnit) END IF CASE( 'mesh update 1','mesh update 2', 'mesh update 3' ) CASE( 'displacement' ) WRITE( VtkUnit,'("VECTORS ",A," double")' ) "Displacement" DO i = 1, Model % NumberOfNodes k = i IF ( ASSOCIATED( Var % Perm ) ) k = Var % Perm(k) IF ( k > 0 ) THEN DO j=1,Var % DOFs WRITE( VtkUnit,'(E20.11E3)',ADVANCE='NO' ) & Var % Values(Var % DOFs*(k-1)+j) END DO IF ( Var % DOFs < 3 ) THEN WRITE( VtkUnit,'(" 0.0")',ADVANCE='NO' ) END IF WRITE( VtkUnit, * ) ELSE Var1 => Model % Variables DO WHILE( ASSOCIATED( Var1 ) ) IF ( TRIM( Var1 % Name ) == 'mesh update' ) EXIT Var1 => Var1 % Next END DO IF ( ASSOCIATED( Var1 ) ) THEN k = i IF ( ASSOCIATED(Var1 % Perm ) ) k = Var1 % Perm(k) IF ( k > 0 ) THEN DO j=1,Var1 % DOFs WRITE( VtkUnit,'(E20.11E3)',ADVANCE='NO' ) & Var1 % Values(Var1 % DOFs*(k-1)+j) END DO IF ( Var1 % DOFs < 3 ) THEN WRITE( VtkUnit,'(" 0.0")', ADVANCE='NO' ) END IF WRITE( VtkUnit, * ) ELSE WRITE( VtkUnit,'(" 0.0 0.0 0.0")' ) END IF ELSE WRITE( VtkUnit,'(" 0.0 0.0 0.0")' ) END IF END IF END DO CASE( 'displacement 1','displacement 2','displacement 3' ) CASE( 'flow solution' ) CALL WriteVector( "Velocity", Var, Model % NumberOfNodes, 4, & VtkUnit ) WRITE( VtkUnit,'("SCALARS ",A," double")' ) "Pressure" WRITE( VtkUnit,'("LOOKUP_TABLE default")' ) DO i = 1, Model % NumberOfNodes k = i IF ( ASSOCIATED( Var % Perm ) ) k = Var % Perm(k) WRITE( VtkUnit,'(E20.11E3)' ) & Var % Values(Var % DOFs*(k-1)+Var % DOFs) END DO CASE( 'velocity 1','velocity 2','velocity 3','pressure' ) CASE( 'magnetic field' ) CALL WriteVector( "MagField", Var, Model % NumberOfNodes, 3, & VtkUnit) CASE( 'magnetic field 1','magnetic field 2', 'magnetic field 3' ) CASE( 'electric current' ) CALL WriteVector( "Current", Var, Model % NumberOfNodes, 3, & VtkUnit) CASE('electric current 1','electric current 2','electric current 3') CASE( 'coordinate 1','coordinate 2','coordinate 3' ) CASE( 'magnetic flux density' ) CALL WriteVector("MagneticFlux", Var, Model % NumberOfNodes, 3,& VtkUnit) CASE( 'magnetic flux density 1','magnetic flux density 2', & 'magnetic flux density 3' ) CASE DEFAULT IF ( Var % DOFs == 1 ) THEN DO i=1,Var % NameLen str(i:i) = Var % Name(i:i) IF (str(i:i) == ' ') str(i:i) = '.' END DO str(1:1) = CHAR(ICHAR(str(1:1))-ICHAR('a')+ICHAR('A')) WRITE(VtkUnit,'("SCALARS ",A," double")') str(1:Var%NameLen) WRITE( VtkUnit,'("LOOKUP_TABLE default")' ) DO i = 1, Model % NumberOfNodes k = i IF ( ASSOCIATED( Var % Perm ) ) k = Var % Perm(k) IF ( k > 0 ) THEN WRITE( VtkUnit,'(E20.11E3)' ) Var % Values(k) ELSE WRITE( VtkUnit,'(" 0.0")' ) END IF END DO END IF END SELECT Var => Var % Next END DO ! IF ( FreeSurfaceFlag ) THEN ! WRITE( VtkUnit,'("VECTORS ",A," double")' ) "Coordinates" ! ! DO i = 1, Model % NumberOfNodes ! Var => VariableGet( Model % Variables,'Coordinate 1' ) ! WRITE( VtkUnit,'(E20.11E3)',ADVANCE='NO' ) Var % Values(i) ! ! Var => VariableGet( Model % Variables, 'Coordinate 2' ) ! WRITE( VtkUnit,'(E20.11E3)',ADVANCE='NO' ) Var % Values(i) ! ! Var => VariableGet( Model % Variables, 'Coordinate 3' ) ! WRITE( VtkUnit,'(E20.11E3)' ) Var % Values(i) ! END DO ! END IF CLOSE( VtkUnit ) END SUBROUTINE WriteVtkLegacyFile SUBROUTINE WriteVector( VarName, Var, nNodes, MaxDOF, IOUnit ) CHARACTER(*), INTENT(IN) :: VarName TYPE(Variable_t), INTENT(IN) :: Var INTEGER, INTENT(IN) :: nNodes, MaxDOF, IOUnit INTEGER :: i, j, k, n n = Var % DOFs - (MaxDOF - 3) WRITE( IOUnit, '("VECTORS ",A," double")' ) TRIM( VarName ) DO i = 1, nNodes k = i IF ( ASSOCIATED( Var % Perm ) ) k = Var % Perm(k) IF ( k > 0 ) THEN DO j=1, n WRITE( IOUnit,'(E20.11E3)',ADVANCE='NO' ) & Var % Values(Var % DOFs*(k-1)+j) END DO IF ( n < 3 ) THEN WRITE( IOUnit,'(" 0.0")',ADVANCE='NO' ) END IF WRITE( IOUnit, * ) ELSE WRITE( IOUnit,'(" 0.0 0.0 0.0")' ) END IF END DO END SUBROUTINE WriteVector LOGICAL FUNCTION FreeSurface( Model ) TYPE(Model_t), INTENT(IN) :: Model LOGICAL :: MoveBoundary, GotIt INTEGER :: i FreeSurface = .FALSE. MoveBoundary = .FALSE. DO i = 1, Model % NumberOfBCs FreeSurface = FreeSurface & .OR. ListGetLogical( Model % BCs(i) % Values, & 'Free Surface', GotIt ) IF ( FreeSurface ) THEN MoveBoundary = ListGetLogical( Model % BCs(i) % Values, & 'Internal Move Boundary', GotIt ) IF ( .NOT.GotIt ) MoveBoundary = .TRUE. FreeSurface = FreeSurface .AND. MoveBoundary END IF IF ( FreeSurface ) EXIT END DO END FUNCTION FreeSurface SUBROUTINE WriteGrid ( IOUnit, Model, SubtractDisp ) INTEGER, INTENT(IN) :: IOUnit TYPE(Model_t), INTENT(IN) :: Model LOGICAL, INTENT(IN) :: SubtractDisp ! Subtract Displacement from Coords. TYPE(Variable_t), POINTER :: Displacement, MeshUpdate, Var1, Var2 INTEGER :: i, k, l, nVtkCells, nVtkCellNum ! FIXME: POINTER -> ALLOCATABLE when more compilers support it TYPE(VtkCell_t), POINTER :: VtkCells(:) WRITE ( IOUnit, '("# vtk DataFile Version 3.0")' ) WRITE ( IOUnit, '("ElmerSolver output; started at ", A)' ) TRIM( FormatDate() ) WRITE ( IOUnit, '("ASCII")' ) WRITE ( IOUnit, '("DATASET UNSTRUCTURED_GRID")' ) ! ! Coordinates: ! WRITE ( IOUnit, '("POINTS ", I0, " double")' ) Model % NumberOfNodes ! First, look for displacements Displacement => NULL() MeshUpdate => NULL() IF ( SubtractDisp ) THEN Var1 => Model % Variables DO WHILE( ASSOCIATED( Var1 ) ) IF ( .NOT.Var1 % Output ) THEN Var1 => Var1 % Next CYCLE END IF SELECT CASE( Var1 % Name ) CASE( 'mesh update' ) Var2 => Model % Variables DO WHILE( ASSOCIATED( Var2 ) ) IF ( TRIM( Var2 % Name ) == 'displacement' ) EXIT Var2 => Var2 % Next END DO IF ( .NOT. ASSOCIATED( Var2 ) ) THEN Displacement => Var1 ELSE MeshUpdate => Var1 END IF CASE( 'displacement' ) Displacement => Var1 END SELECT Var1 => Var1 % Next END DO END IF DO i = 1, Model % NumberOfNodes IF ( .NOT.ASSOCIATED( Displacement ) ) THEN WRITE( IOUnit,'(3E20.11E3)' ) Model % Nodes % x(i), & Model % Nodes % y(i), & Model % Nodes % z(i) CYCLE END IF k = Displacement % Perm(i) l = 0 IF ( ASSOCIATED( MeshUpdate ) ) l = MeshUpdate % Perm(i) IF ( k > 0 ) THEN k = Displacement % DOFs * (k-1) SELECT CASE( Displacement % DOFs ) CASE(1) WRITE( IOUnit,'(3E20.11E3)' ) & Model % Nodes % x(i) - Displacement % Values(k+1), & Model % Nodes % y(i), & Model % Nodes % z(i) CASE( 2 ) WRITE( IOUnit,'(3E20.11E3)' ) & Model % Nodes % x(i) - Displacement % Values(k+1), & Model % Nodes % y(i) - Displacement % Values(k+2), & Model % Nodes % z(i) CASE(3) WRITE( IOUnit,'(3E20.11E3)' ) & Model % Nodes % x(i) - Displacement % Values(k+1), & Model % Nodes % y(i) - Displacement % Values(k+2), & Model % Nodes % z(i) - Displacement % Values(k+3) END SELECT ELSE IF ( l > 0 ) THEN k = MeshUpdate % DOFs * (l-1) SELECT CASE( MeshUpdate % DOFs ) CASE(1) WRITE( IOUnit,'(3E20.11E3)' ) & Model % Nodes % x(i) - MeshUpdate % Values(k+1), & Model % Nodes % y(i), & Model % Nodes % z(i) CASE(2) WRITE( IOUnit,'(3E20.11E3)' ) & Model % Nodes % x(i) - MeshUpdate % Values(k+1), & Model % Nodes % y(i) - MeshUpdate % Values(k+2), & Model % Nodes % z(i) CASE(3) WRITE( IOUnit,'(3E20.11E3)' ) & Model % Nodes % x(i) - MeshUpdate % Values(k+1), & Model % Nodes % y(i) - MeshUpdate % Values(k+2), & Model % Nodes % z(i) - MeshUpdate % Values(k+3) END SELECT ELSE WRITE( IOUnit,'(3E20.11E3)' ) Model % Nodes % x(i), & Model % Nodes % y(i), & Model % Nodes % z(i) ENDIF END DO WRITE ( IOUnit, * ) ! ! CELLS ! CALL GetNVtkCells ( Model, nVtkCells, nVtkCellNum ) WRITE( IOUnit,'("CELLS ",I0," ",I0)' ) nVtkCells, nVtkCellNum DO i = 1, Model%NumberOfBulkElements + Model%NumberOfBoundaryElements CALL Elements2Cells ( Model % Elements(i), VtkCells ) DO k = 1, SIZE( VtkCells ) WRITE( IOUnit, '(I0)', ADVANCE='NO' ) VtkCells(k) % nNodes WRITE( IOUnit, '(50(" ",I0))') VtkCells(k) % NodeIndex DEALLOCATE( VtkCells(k) % NodeIndex ) END DO ! FIXME: These explicit deallocations of VtkCells(j) % NodeIndex and ! VtkCells can be removed if/when they are changed to ALLOCATABLEs DEALLOCATE ( VtkCells ) END DO WRITE( IOUnit, * ) ! ! CELL_TYPES ! WRITE( IOUnit,'("CELL_TYPES ",I0)' ) nVtkCells DO i = 1, Model%NumberOfBulkElements + Model%NumberOfBoundaryElements CALL WriteCellType ( IOUnit, Model%Elements(i)%TYPE%ElementCode ) END DO WRITE( IOUnit, * ) CONTAINS ! Here's a bunch of routines to deal with ELMER Element -> VTK Cell ! translations. Most element types has a directly corresponding cell ! type, but some have to be expanded to a compound of simplier cell ! types. ! nVtkCells is the number of VtkCells, and VtkCellNum is the total ! number of integer numbers we have to write to the 'CELLS' section. SUBROUTINE GetNVtkCells( Model,nVtkCells,nVtkCellNum ) TYPE(Model_t), INTENT(IN) :: Model INTEGER, INTENT(OUT) :: nVtkCells, nVtkCellNum nVtkCells = 0 nVtkCellNum = 0 DO i = 1, Model%NumberOfBulkElements+Model%NumberOfBoundaryElements SELECT CASE ( Model % Elements(i) % TYPE % ElementCode ) CASE( 409 ) ! Translate to 4 VTK_QUAD cells nVtkCells = nVtkCells + 4 nVtkCellNum = nVtkCellNum + 16 + 4 CASE DEFAULT nVtkCells = nVtkCells + 1 nVtkCellNum = nVtkCellNum & + Model % Elements(i) % TYPE % NumberOfNodes + 1 END SELECT END DO END SUBROUTINE GetNVtkCells ! Translate an Element_t to an array of VtkCell_t. SUBROUTINE Elements2Cells( Elem, VtkCells ) TYPE(Element_t), INTENT(IN) :: Elem ! FIXME: POINTER -> ALLOCATABLE, INTENT(OUT) when more compilers ! support it TYPE(VtkCell_t), POINTER :: VtkCells(:) INTEGER :: k SELECT CASE ( Elem % TYPE % ElementCode ) CASE( 409 ) ! Translate to 4 VTK_QUAD cells ALLOCATE( VtkCells(4) ) VtkCells % nNodes = 4 DO k = 1, 4 ALLOCATE( VtkCells(k) % NodeIndex(4) ) END DO VtkCells(1) % NodeIndex = (/ Elem % NodeIndexes(1), & Elem % NodeIndexes(5), & Elem % NodeIndexes(9), & Elem % NodeIndexes(8) /) - 1 VtkCells(2) % NodeIndex = (/ Elem % NodeIndexes(5), & Elem % NodeIndexes(2), & Elem % NodeIndexes(6), & Elem % NodeIndexes(9) /) - 1 VtkCells(3) % NodeIndex = (/ Elem % NodeIndexes(9), & Elem % NodeIndexes(6), & Elem % NodeIndexes(3), & Elem % NodeIndexes(7) /) - 1 VtkCells(4) % NodeIndex = (/ Elem % NodeIndexes(8), & Elem % NodeIndexes(9), & Elem % NodeIndexes(7), & Elem % NodeIndexes(4) /) - 1 CASE DEFAULT ALLOCATE( VtkCells(1) ) VtkCells % nNodes = Elem % TYPE % NumberOfNodes ALLOCATE( VtkCells(1) % NodeIndex(VtkCells(1) % nNodes) ) VtkCells(1) % NodeIndex = Elem % NodeIndexes-1 END SELECT END SUBROUTINE Elements2Cells SUBROUTINE WriteCellType( IOUnit,Code ) INTEGER, INTENT(IN) :: IOUnit, Code SELECT CASE (Code) CASE( 101 ) WRITE( IOUnit,'(I0)' ) 1 CASE( 202 ) WRITE( IOUnit,'(I0)' ) 3 CASE( 203 ) WRITE( IOUnit,'(I0)' ) 21 CASE( 303 ) WRITE( IOUnit,'(I0)' ) 5 CASE( 306 ) WRITE( IOUnit,'(I0)' ) 22 CASE( 404 ) WRITE( IOUnit,'(I0)' ) 9 CASE( 408 ) WRITE( IOUnit,'(I0)' ) 23 CASE( 409 ) ! Translate to 4 VTK_QUAD cells WRITE( IOUnit,'(I0)' ) (/ 9, 9, 9, 9 /) CASE( 504 ) WRITE( IOUnit,'(I0)' ) 10 CASE( 510 ) WRITE( IOUnit,'(I0)' ) 24 CASE( 605 ) WRITE( IOUnit,'(I0)' ) 14 CASE( 706 ) WRITE( IOUnit,'(I0)' ) 13 CASE( 808 ) WRITE( IOUnit,'(I0)' ) 12 CASE( 820 ) WRITE( IOUnit,'(I0)' ) 25 CASE DEFAULT WRITE( IOUnit,'(I0)' ) -Code END SELECT END SUBROUTINE WriteCellType END SUBROUTINE WriteGrid END MODULE VtkLegacyFile !------------------------------------------------------------------------------ SUBROUTINE VtkOutputSolver( Model,Solver,dt,TransientSimulation ) !DEC$ATTRIBUTES DLLEXPORT :: VtkOutputSolver !------------------------------------------------------------------------------ USE DefUtils USE VtkLegacyFile IMPLICIT NONE TYPE(Solver_t) :: Solver TYPE(Model_t) :: Model REAL(dp) :: dt LOGICAL :: TransientSimulation INTEGER, SAVE :: nTime = 0 LOGICAL :: GotIt CHARACTER(MAX_NAME_LEN), SAVE :: FilePrefix ! Avoid compiler warings about unused variables IF ( TransientSimulation ) THEN; ENDIF IF ( dt > 0.0 ) THEN; ENDIF IF ( nTime == 0 ) THEN FilePrefix = GetString( Solver % Values,'Output File Name',GotIt ) IF ( .NOT.GotIt ) FilePrefix = "Output" END IF nTime = nTime + 1 CALL WriteData( TRIM(FilePrefix), Model, nTime ) CONTAINS SUBROUTINE WriteData( Prefix, Model, nTime ) CHARACTER(*), INTENT(IN) :: Prefix TYPE(Model_t) :: Model INTEGER, INTENT(IN) :: nTime CHARACTER(MAX_NAME_LEN) :: VtkFile TYPE(Mesh_t), POINTER :: Mesh TYPE(Variable_t), POINTER :: Var ! 27/04/09, SCAI: Parallel I/O for ResultOutputSolve INTEGER :: i, j, k, myrank, ierr LOGICAL :: EigAnal REAL(dp), POINTER :: OrigValues(:) INTEGER :: OrigDOFs CHARACTER(MAX_NAME_LEN) :: Dir ! 27/04/09, SCAI: Parallel I/O for ResultOutputSolve CALL mpi_comm_rank(mpi_comm_world,myrank,ierr) Mesh => Model % Meshes DO WHILE( ASSOCIATED(Mesh) ) IF ( .NOT.Mesh % OutputActive ) THEN Mesh => Mesh % next CYCLE END IF IF (LEN_TRIM(Mesh % Name) > 0 ) THEN Dir = TRIM(Mesh % Name) // "/" ELSE Dir = "./" END IF CALL SetCurrentMesh( Model, Mesh ) EigAnal = .FALSE. Solvers: DO i = 1, Model % NumberOfSolvers EigAnal = ListGetLogical( Model % Solvers(i) % Values, & "Eigen Analysis", GotIt ) Var => Model % Solvers(i) % Variable IF ( EigAnal .AND. ASSOCIATED(Var % EigenValues) ) THEN DO j = 1, Model % Solvers(i) % NOfEigenValues OrigValues => Var % Values OrigDOFs = Var % DOFs IF ( Model % Solvers(i) % Matrix % COMPLEX ) THEN Var % DOFs = Var % DOFs*2 ALLOCATE( Var % Values(2*SIZE(Var%EigenVectors,2)) ) FORALL ( k = 1:SIZE(Var % Values)/2 ) Var%Values(2*k-1) = REAL(Var%EigenVectors(j,k)) Var%Values(2*k) = AIMAG(Var%EigenVectors(j,k)) END FORALL ELSE ALLOCATE( Var % Values(SIZE(Var % EigenVectors,2)) ) Var % Values = Var % EigenVectors(j,:) END IF WRITE( VtkFile, '(A,A,I4.4,"_",I2.2,".vtk")' ) & TRIM(Dir), Prefix, nTime, j CALL WriteVtkLegacyFile( VtkFile, Model, .FALSE. ) DEALLOCATE( Var % Values ) Var % Values => OrigValues Var % DOFs = OrigDOFs END DO EXIT Solvers END IF END DO Solvers IF ( .NOT.EigAnal ) THEN ! 27/04/09, SCAI: Parallel I/O for ResultOutputSolve WRITE( VtkFile,'(A,A,".",I2.2,".",I4.4,".vtk")' ) TRIM(Dir),Prefix,myrank,nTime CALL WriteVtkLegacyFile( VtkFile, Model, .TRUE. ) END IF Mesh => Mesh % next END DO END SUBROUTINE WriteData END SUBROUTINE VtkOutputSolver !------------------------------------------------------------------------------ ! The DX output solver begins here !------------------------------------------------------------------------------ MODULE DXFile USE MeshUtils USE ElementDescription IMPLICIT NONE ! PRIVATE SAVE PUBLIC :: WriteDXFiles INTEGER, PARAMETER :: MAX_VERTIX = 4, MAX_PART_ELEM = 21 CONTAINS SUBROUTINE WriteDXFiles( Prefix, Model, SubtractDisp, nTime ) CHARACTER(LEN=*), INTENT(IN) :: Prefix TYPE(Model_t) :: Model LOGICAL, INTENT(IN) :: SubtractDisp INTEGER, INTENT(IN) :: nTime TYPE(Variable_t), POINTER :: Var,Var1 CHARACTER(LEN=512) :: str INTEGER :: i INTEGER, PARAMETER :: MasterUnit = 58 IF ( nTime == 1) THEN CALL WriteGrid( Prefix, Model, SubtractDisp ) OPEN( MasterUnit, FILE = Prefix // "Master.dx", STATUS="unknown" ) WRITE( MasterUnit, '(A)') 'object "group" class group' END IF Var => Model % Variables DO WHILE( ASSOCIATED( Var ) ) IF ( .NOT.Var % Output ) THEN Var => Var % Next CYCLE END IF SELECT CASE( Var % Name ) CASE( 'time','timestep','timestep size', 'timestep interval' ) CASE( 'mesh update' ) Var1 => Model % Variables DO WHILE( ASSOCIATED( Var1 ) ) IF ( TRIM( Var1 % Name ) == 'displacement' ) EXIT Var1 => Var1 % Next END DO IF ( .NOT.ASSOCIATED( Var1 ) ) THEN CALL WriteVariable( "MeshUpdate", Var, & Model % NumberOfNodes, Var % DOFs, 0, & nTime, MasterUnit, Prefix ) END IF CASE( 'mesh update 1','mesh update 2', 'mesh update 3' ) CASE( 'displacement' ) CALL WriteDisplacement( Var, Model, nTime, MasterUnit, Prefix ) CASE( 'displacement 1','displacement 2','displacement 3' ) CASE( 'flow solution' ) CALL WriteVariable( "Velocity", Var, Model % NumberOfNodes, & Var % DOFs-1, 0, nTime, MasterUnit, Prefix ) CALL WriteVariable( "Pressure", Var, Model % NumberOfNodes, 1, & Var % DOFs-1, nTime, MasterUnit, Prefix ) CASE( 'velocity 1','velocity 2','velocity 3','pressure' ) CASE( 'magnetic field' ) CALL WriteVariable( "MagField", Var, Model % NumberOfNodes, & Var % DOFs, 0, nTime, MasterUnit, Prefix ) CASE( 'magnetic field 1','magnetic field 2', 'magnetic field 3' ) CASE( 'electric current' ) CALL WriteVariable( "Current", Var, Model % NumberOfNodes, & Var % DOFs, 0, nTime, MasterUnit, Prefix ) CASE('electric current 1','electric current 2','electric current 3') CASE( 'coordinate 1','coordinate 2','coordinate 3' ) CASE( 'magnetic flux density' ) CALL WriteVariable( "MagneticFlux", Var, Model % NumberOfNodes,& Var % DOFs, 0, nTime, MasterUnit, Prefix ) CASE( 'magnetic flux density 1','magnetic flux density 2', & 'magnetic flux density 3' ) CASE DEFAULT DO i=1,Var % NameLen str(i:i) = Var % Name(i:i) IF( str(i:i) == ' ' ) str(i:i) = '_' END DO str(1:1) = CHAR(ICHAR(str(1:1))-ICHAR('a')+ICHAR('A')) CALL WriteVariable( TRIM(str), Var, Model % NumberOfNodes, & Var % DOFs, 0, nTime, MasterUnit, Prefix ) END SELECT Var => Var % Next END DO IF( nTime == 1) THEN CLOSE( MasterUnit ) END IF END SUBROUTINE WriteDXFiles SUBROUTINE WriteVariable( VarName, Var, nNodes, SelfDOF, Offset, nTime, & MasterUnit, Prefix ) CHARACTER(*), INTENT(IN) :: VarName TYPE(Variable_t), INTENT(IN) :: Var INTEGER, INTENT(IN) :: nNodes, SelfDOF, Offset, nTime, MasterUnit CHARACTER(*), INTENT(IN) :: Prefix INTEGER :: i, j, k INTEGER :: FUnit CHARACTER(MAX_NAME_LEN) :: FName, RelativeFName, MeshFile CHARACTER(7) :: VType FUnit = MasterUnit + 1 FName = Prefix // VarName // ".dx" i = INDEX(FName, '/', BACK=.TRUE.) RelativeFName = FName(i+1:) i = INDEX(Prefix, '/', BACK=.TRUE.) MeshFile = Prefix(i+1:) // "Mesh.dx" IF( nTime == 1 ) THEN IF ( SelfDOF == 1 ) THEN VType = "-scalar" ELSE VType = "-vector" END IF WRITE( MasterUnit, '(A)' ) 'member "' // VarName // VType & // '" value file "' // TRIM(RelativeFName) // '","' & // VarName // 'series' // '"' OPEN( FUnit, FILE=FName, STATUS="unknown" ) ELSE OPEN( FUnit, FILE=FName, STATUS="old", POSITION="append" ) DO i = 1, 2+(nTime-1)*7 BACKSPACE FUnit ! Yuck! END DO END IF IF( SelfDOF == 1 )THEN WRITE( FUnit,'("object ",I0," class array type double rank 0 & &items ",I0," data follows")') nTime, nNodes ELSE WRITE( FUnit,'("object ",I0," class array type double rank 1 shape & &",I0," items ",I0," data follows")') nTime, & SelfDOF, nNodes END IF DO i = 1, nNodes k = i IF( ASSOCIATED( Var % Perm ) ) k = Var % Perm(k) IF( k > 0 )THEN DO j=1, SelfDOF WRITE( FUnit,'(E20.11E3)',ADVANCE='NO' ) & Var % Values(Var % DOFs*(k-1)+j+Offset) END DO WRITE( FUnit, * ) ELSE WRITE( FUnit, '(9F4.1)' ) (/ (0.0, k=1,SelfDOF) /) END IF END DO WRITE( FUnit, '(A)' ) 'attribute "dep" string "positions"' WRITE( FUnit, * ) DO i = nTime + 1, 2*nTime WRITE( FUnit,'("object ",I0," class field")' ) i WRITE( FUnit,'(A,I0)' ) 'component "data" value ', i - nTime WRITE( FUnit,'(A,A,A)' ) 'component "positions" value file "',& TRIM( MeshFile ), '",1' WRITE( FUnit,'(A,A,A)' ) 'component "connections" value file "',& TRIM( MeshFile ), '",2' WRITE( FUnit, '(A,A,A)' ) 'attribute "name" string "', VarName,'"' WRITE( FUnit, * ) END DO WRITE( FUnit,'(A)' ) 'object "'//VarName//'series'//'" class series' DO i = 1, nTime WRITE( FUnit, '("member ",I0," value ",I0," position ",I0)' ) & i-1, i+nTime, i END DO WRITE( FUnit, '("end")' ) CLOSE( FUnit ) END SUBROUTINE WriteVariable ! WriteDisplacement is like WriteVariable, but specialized for ! Displacements; displacements need special treatment. SUBROUTINE WriteDisplacement( Var, Model, nTime, MasterUnit, Prefix ) TYPE(Variable_t), INTENT(IN) :: Var TYPE(Model_t), INTENT(IN) :: Model INTEGER, INTENT(IN) :: nTime, MasterUnit CHARACTER(*), INTENT(IN) :: Prefix INTEGER :: i, j, k INTEGER :: FUnit CHARACTER(MAX_NAME_LEN) :: FName, RelativeFName, MeshFile TYPE(Variable_t), POINTER :: Var1 FUnit = MasterUnit + 1 FName = Prefix // "Displacement.dx" i = INDEX(FName, '/', BACK=.TRUE.) RelativeFName = FName(i+1:) i = INDEX(Prefix, '/', BACK=.TRUE.) MeshFile = Prefix(i+1:) // "Mesh.dx" IF( nTime == 1 ) THEN WRITE( MasterUnit,'(A)' ) & 'member "Displacement-vector" value file "' & // TRIM(RelativeFName) // '", "Displacementseries"' OPEN( FUnit, FILE=FName, STATUS="unknown" ) ELSE OPEN( FUnit, FILE=FName, STATUS="old", POSITION="append" ) DO i = 1, 2+(nTime-1)*7 BACKSPACE FUnit ! Yuck! END DO END IF WRITE( FUnit,'("object ",I0," class array type double rank 1 shape ", & &I0," items ",I0," data follows")') nTime, Var % DOFs, & Model % NumberOfNodes DO i = 1, Model % NumberOfNodes k = i IF( ASSOCIATED( Var % Perm ) ) k = Var % Perm(k) IF( k > 0 )THEN DO j=1, Var % DOFs WRITE( FUnit,'(E20.11E3)',ADVANCE='NO' ) & Var % Values(Var % DOFs*(k-1)+j) END DO WRITE( FUnit, * ) ELSE Var1 => Model % Variables DO WHILE( ASSOCIATED( Var1 ) ) IF ( TRIM( Var1 % Name ) == 'mesh update' ) EXIT Var1 => Var1 % Next END DO IF( ASSOCIATED( Var1 ) )THEN k = i IF( ASSOCIATED(Var1 % Perm ) ) k = Var1 % Perm(k) IF( k > 0 )THEN DO j=1,Var1 % DOFs WRITE( FUnit,'(E20.11E3)',ADVANCE='NO' ) & Var1 % Values(Var1 % DOFs*(k-1)+j) END DO WRITE( FUnit, * ) ELSE WRITE( FUnit, '(9F4.1)' ) (/ (0.0, k=1,Var % DOFs) /) END IF ELSE WRITE( FUnit, '(9F4.1)' ) (/ (0.0, k=1,Var % DOFs) /) END IF END IF END DO WRITE( FUnit, '(A)' ) 'attribute "dep" string "positions"' WRITE( FUnit, * ) DO i = nTime + 1, 2*nTime WRITE( FUnit,'("object ",I0," class field")' ) i WRITE( FUnit,'(A,I0)' ) 'component "data" value ', i - nTime WRITE( FUnit,'(A,A,A)' ) 'component "positions" value file "',& TRIM( MeshFile ), '",1' WRITE( FUnit,'(A,A,A)' ) 'component "connections" value file "',& TRIM( MeshFile ), '",2' WRITE( FUnit,'(A,A,A)' ) 'attribute "name" string "Displacement"' WRITE( FUnit,* ) END DO WRITE( FUnit,'(A)' ) 'object "Displacementseries" class series' DO i = 1, nTime WRITE( FUnit,'("member ",I0," value ",I0," position ",I0)' ) & i-1, i+nTime, i END DO WRITE( FUnit, '("end")' ) CLOSE( FUnit ) END SUBROUTINE WriteDisplacement SUBROUTINE WriteGrid ( PRefix, Model, SubtractDisp ) CHARACTER(*), INTENT(IN) :: Prefix TYPE(Model_t), INTENT(IN) :: Model LOGICAL, INTENT(IN) :: SubtractDisp ! Subtract Displacement from Coords. TYPE(Variable_t), POINTER :: Displacement, MeshUpdate, Var1, Var2 INTEGER :: i, k, l, nElem, nVertix CHARACTER(MAX_NAME_LEN) :: FName INTEGER, PARAMETER :: FUnit = 58 INTEGER :: NodeIndex(MAX_PART_ELEM, MAX_VERTIX) FName = Prefix // "Mesh.dx" OPEN( UNIT=FUnit, FILE=FName, STATUS="unknown", ACTION="write" ) WRITE ( FUnit,'("# ElmerSolver output; started at ",A)' ) & TRIM( FormatDate() ) ! ! Coordinates: ! WRITE ( FUnit, '("# Node Coordinates")' ) WRITE ( FUnit, '("object 1 class array type float rank 1 shape 3 & &items ",I0," data follows")' ) Model % NumberOfNodes ! First, look for displacements Displacement => NULL() MeshUpdate => NULL() IF( SubtractDisp )THEN Var1 => Model % Variables DO WHILE( ASSOCIATED( Var1 ) ) IF( .NOT.Var1 % Output )THEN Var1 => Var1 % Next CYCLE END IF SELECT CASE( Var1 % Name ) CASE( 'mesh update' ) Var2 => Model % Variables DO WHILE( ASSOCIATED( Var2 ) ) IF ( TRIM( Var2 % Name ) == 'displacement' ) EXIT Var2 => Var2 % Next END DO IF( .NOT. ASSOCIATED( Var2 ) )THEN Displacement => Var1 ELSE MeshUpdate => Var1 END IF CASE( 'displacement' ) Displacement => Var1 END SELECT Var1 => Var1 % Next END DO END IF DO i = 1, Model % NumberOfNodes IF( .NOT.ASSOCIATED( Displacement ) )THEN WRITE( FUnit,'(3E20.11E3)' ) Model % Nodes % x(i), & Model % Nodes % y(i), & Model % Nodes % z(i) CYCLE END IF k = Displacement % Perm(i) l = 0 IF ( ASSOCIATED( MeshUpdate ) ) l = MeshUpdate % Perm(i) IF( k > 0 )THEN k = Displacement % DOFs * (k-1) SELECT CASE( Displacement % DOFs ) CASE(1) WRITE( FUnit,'(3E20.11E3)' ) & Model % Nodes % x(i) - Displacement % Values(k+1), & Model % Nodes % y(i), & Model % Nodes % z(i) CASE( 2 ) WRITE( FUnit,'(3E20.11E3)' ) & Model % Nodes % x(i) - Displacement % Values(k+1), & Model % Nodes % y(i) - Displacement % Values(k+2), & Model % Nodes % z(i) CASE(3) WRITE( FUnit,'(3E20.11E3)' ) & Model % Nodes % x(i) - Displacement % Values(k+1), & Model % Nodes % y(i) - Displacement % Values(k+2), & Model % Nodes % z(i) - Displacement % Values(k+3) END SELECT ELSE IF ( l > 0 ) THEN k = MeshUpdate % DOFs * (l-1) SELECT CASE( MeshUpdate % DOFs ) CASE(1) WRITE( FUnit,'(3E20.11E3)' ) & Model % Nodes % x(i) - MeshUpdate % Values(k+1), & Model % Nodes % y(i), & Model % Nodes % z(i) CASE(2) WRITE( FUnit,'(3E20.11E3)' ) & Model % Nodes % x(i) - MeshUpdate % Values(k+1), & Model % Nodes % y(i) - MeshUpdate % Values(k+2), & Model % Nodes % z(i) CASE(3) WRITE( FUnit,'(3E20.11E3)' ) & Model % Nodes % x(i) - MeshUpdate % Values(k+1), & Model % Nodes % y(i) - MeshUpdate % Values(k+2), & Model % Nodes % z(i) - MeshUpdate % Values(k+3) END SELECT ELSE WRITE( FUnit,'(3E20.11E3)' ) Model % Nodes % x(i), & Model % Nodes % y(i), & Model % Nodes % z(i) ENDIF END DO WRITE ( FUnit, * ) ! Elements. ! ! The only DX elements we use are triangles (for 2D) and tetrahedra ! (3D). All other element types are translated into compounds of either ! of these. ! ! Note that, at the moment, either all elements have to be 2D, or they ! all have to be 3D; therefore, boundary elements are not written. WRITE ( FUnit, '("# Element definitions")' ) CALL GetNElem( Model, nElem, nVertix ) WRITE( FUnit,'("object 2 class array type int rank 1 shape ", & &I0, " items ", I0, " data follows")' ) nVertix, nElem DO i = 1, Model % NumberOfBulkElements CALL TranslateElem( Model % Elements(i), NodeIndex, nElem ) DO k = 1, nElem WRITE( FUnit, '(4(" ",I0))') NodeIndex(k,1:nVertix) END DO END DO IF ( nVertix == 3 ) THEN WRITE( FUnit,'(A)') 'attribute "element type" string "triangles"' ELSE WRITE( FUnit,'(A)') 'attribute "element type" string "tetrahedra"' END IF WRITE( FUnit, '("end")' ) CONTAINS ! Here's a bunch of routines to deal with ELMER -> DX element ! translations. ! nElem is the number of DX elements, and nVertix is the numbe of ! vertices per element (either 3 for triangle or 4 for tetrahedra). SUBROUTINE GetNElem( Model, nElem, nVertix ) TYPE(Model_t), INTENT(IN) :: Model INTEGER, INTENT(OUT) :: nElem, nVertix IF (Model % Elements(1) % TYPE % ElementCode < 500) THEN nVertix = 3 ELSE nVertix = 4 END IF nElem = 0 DO i = 1, Model%NumberOfBulkElements SELECT CASE ( Model % Elements(i) % TYPE % ElementCode ) CASE( 303 ) nElem = nElem + 1 CASE( 404 ) nElem = nElem + 2 CASE( 409 ) nElem = nElem + 8 CASE( 504 ) nElem = nElem + 1 CASE( 605 ) nElem = nElem + 2 CASE( 613 ) nElem = nElem + 14 CASE( 706 ) nElem = nElem + 3 CASE( 715 ) nElem = nElem + 21 CASE( 808 ) nElem = nElem + 4 CASE DEFAULT WRITE (0, *) 'Sorry! Element type ', & Model % Elements(i) % TYPE % ElementCode, & ' not supported! ' STOP END SELECT END DO END SUBROUTINE GetNElem ! Translate an ELMER element to a (compound of) DX element(s). SUBROUTINE TranslateElem( Elem, NodeIndex, nElem ) TYPE(Element_t), INTENT(IN) :: Elem INTEGER, INTENT(OUT) :: NodeIndex(MAX_PART_ELEM, MAX_VERTIX) INTEGER, INTENT(OUT) :: nElem SELECT CASE ( Elem % TYPE % ElementCode ) CASE( 303 ) ! Triangle nElem = 1 NodeIndex(1,:3) = Elem % NodeIndexes(1:3) CASE( 404 ) ! Quad; translate into 2 triangles nElem = 2 NodeIndex(1,:3) = Elem % NodeIndexes(1:3) NodeIndex(2,:3) = Elem % NodeIndexes((/ 1,3,4 /)) CASE( 409 ) nElem = 8 NodeIndex(1,:3) = Elem % NodeIndexes((/ 1,5,9 /)) NodeIndex(2,:3) = Elem % NodeIndexes((/ 5,2,9 /)) NodeIndex(3,:3) = Elem % NodeIndexes((/ 2,6,9 /)) NodeIndex(4,:3) = Elem % NodeIndexes((/ 6,3,9 /)) NodeIndex(5,:3) = Elem % NodeIndexes((/ 3,7,9 /)) NodeIndex(6,:3) = Elem % NodeIndexes((/ 7,4,9 /)) NodeIndex(7,:3) = Elem % NodeIndexes((/ 4,8,9 /)) NodeIndex(8,:3) = Elem % NodeIndexes((/ 8,1,9 /)) CASE( 504 ) ! Tetrahedron nElem = 1 NodeIndex(1,:) = Elem % NodeIndexes(1:4) CASE( 605 ) ! Pyramid nElem = 2 NodeIndex(1,:) = Elem % NodeIndexes((/ 3,5,4,1 /)) NodeIndex(2,:) = Elem % NodeIndexes((/ 3,5,2,1 /)) CASE( 613 ) nElem = 14 NodeIndex(1,:) = Elem % NodeIndexes((/ 7,8,3,12 /)) NodeIndex(2,:) = Elem % NodeIndexes((/ 10,11,12,5 /)) NodeIndex(3,:) = Elem % NodeIndexes((/ 10,13,12,5 /)) NodeIndex(4,:) = Elem % NodeIndexes((/ 9,10,13,11 /)) NodeIndex(5,:) = Elem % NodeIndexes((/ 9,13,11,12 /)) NodeIndex(6,:) = Elem % NodeIndexes((/ 9,10,11,1 /)) NodeIndex(7,:) = Elem % NodeIndexes((/ 9,6,11,1 /)) NodeIndex(8,:) = Elem % NodeIndexes((/ 9,6,11,12 /)) NodeIndex(9,:) = Elem % NodeIndexes((/ 9,8,12,4 /)) NodeIndex(10,:) = Elem % NodeIndexes((/ 9,13,12,4 /)) NodeIndex(11,:) = Elem % NodeIndexes((/ 7,9,8,12 /)) NodeIndex(12,:) = Elem % NodeIndexes((/ 7,9,6,12 /)) NodeIndex(13,:) = Elem % NodeIndexes((/ 7,6,11,12 /)) NodeIndex(14,:) = Elem % NodeIndexes((/ 7,6,11,2 /)) CASE( 706 ) ! Wedge nElem = 3 NodeIndex(1,:) = Elem % NodeIndexes((/ 5,4,3,1 /)) NodeIndex(2,:) = Elem % NodeIndexes((/ 5,3,2,1 /)) NodeIndex(3,:) = Elem % NodeIndexes((/ 5,6,4,3 /)) CASE( 715 ) nElem = 21 NodeIndex(1,:) = Elem % NodeIndexes((/ 10,11,5,2 /)) NodeIndex(2,:) = Elem % NodeIndexes((/ 12,11,6,3 /)) NodeIndex(3,:) = Elem % NodeIndexes((/ 12,10,4,1 /)) NodeIndex(4,:) = Elem % NodeIndexes((/ 7,8,11,2 /)) NodeIndex(5,:) = Elem % NodeIndexes((/ 7,10,11,2 /)) NodeIndex(6,:) = Elem % NodeIndexes((/ 13,14,11,5 /)) NodeIndex(7,:) = Elem % NodeIndexes((/ 13,10,11,5 /)) NodeIndex(8,:) = Elem % NodeIndexes((/ 9,10,8,11 /)) NodeIndex(9,:) = Elem % NodeIndexes((/ 9,7,10,8 /)) NodeIndex(10,:) = Elem % NodeIndexes((/ 9,12,10,11 /)) NodeIndex(11,:) = Elem % NodeIndexes((/ 9,12,10,1 /)) NodeIndex(12,:) = Elem % NodeIndexes((/ 9,7,10,1 /)) NodeIndex(13,:) = Elem % NodeIndexes((/ 9,8,11,3 /)) NodeIndex(14,:) = Elem % NodeIndexes((/ 9,12,11,3 /)) NodeIndex(15,:) = Elem % NodeIndexes((/ 15,12,10,11 /)) NodeIndex(16,:) = Elem % NodeIndexes((/ 15,12,10,4 /)) NodeIndex(17,:) = Elem % NodeIndexes((/ 15,10,14,11 /)) NodeIndex(18,:) = Elem % NodeIndexes((/ 15,13,10,14 /)) NodeIndex(19,:) = Elem % NodeIndexes((/ 15,13,10,4 /)) NodeIndex(20,:) = Elem % NodeIndexes((/ 15,14,11,6 /)) NodeIndex(21,:) = Elem % NodeIndexes((/ 15,12,11,6 /)) CASE( 808 ) ! Hexahedron; translate into 4 tetrahedra nElem = 4 NodeIndex(1,:) = Elem % NodeIndexes((/ 1,2,4,5 /)) NodeIndex(2,:) = Elem % NodeIndexes((/ 2,3,4,7 /)) NodeIndex(3,:) = Elem % NodeIndexes((/ 2,7,5,6 /)) NodeIndex(4,:) = Elem % NodeIndexes((/ 4,5,7,8 /)) END SELECT NodeIndex = NodeIndex - 1 END SUBROUTINE TranslateElem END SUBROUTINE WriteGrid END MODULE DXFile !------------------------------------------------------------------------------ SUBROUTINE DXOutputSolver( Model,Solver,dt,TransientSimulation ) !DEC$ATTRIBUTES DLLEXPORT :: DXOutputSolver !------------------------------------------------------------------------------ USE DefUtils USE DXFile IMPLICIT NONE TYPE(Solver_t) :: Solver TYPE(Model_t) :: Model REAL(dp) :: dt LOGICAL :: TransientSimulation INTEGER, SAVE :: nTime = 0 LOGICAL :: GotIt CHARACTER(MAX_NAME_LEN), SAVE :: FilePrefix ! Avoid compiler warings about unused variables IF ( TransientSimulation ) THEN; ENDIF IF ( dt > 0.0 ) THEN; ENDIF IF ( nTime == 0 ) THEN FilePrefix = GetString( Solver % Values,'Output File Name',GotIt ) IF ( .NOT.GotIt ) FilePrefix = "Output" END IF nTime = nTime + 1 CALL WriteData( TRIM(FilePrefix), Model, nTime ) CONTAINS SUBROUTINE WriteData( Prefix, Model, nTime ) CHARACTER(*), INTENT(IN) :: Prefix TYPE(Model_t) :: Model INTEGER, INTENT(IN) :: nTime TYPE(Mesh_t), POINTER :: Mesh TYPE(Variable_t), POINTER :: Var INTEGER :: i, j, k LOGICAL :: EigAnal REAL(dp), POINTER :: OldValues(:) CHARACTER(MAX_NAME_LEN) :: Dir Mesh => Model % Meshes DO WHILE( ASSOCIATED(Mesh) ) IF ( .NOT.Mesh % OutputActive ) THEN Mesh => Mesh % next CYCLE END IF IF (LEN_TRIM(Mesh % Name) > 0 ) THEN Dir = TRIM(Mesh % Name) // "/" ELSE Dir = "./" END IF CALL SetCurrentMesh( Model, Mesh ) EigAnal = .FALSE. Solvers: DO i = 1, Model % NumberOfSolvers EigAnal = ListGetLogical( Model % Solvers(i) % Values, & "Eigen Analysis", GotIt ) Var => Model % Solvers(i) % Variable IF ( EigAnal .AND. ASSOCIATED(Var % EigenValues) ) THEN DO j = 1, Model % Solvers(i) % NOfEigenValues OldValues => Var % Values IF ( Model % Solvers(i) % Matrix % COMPLEX ) THEN ALLOCATE( Var % Values(2*SIZE(Var%EigenVectors,2)) ) FORALL ( k = 1:SIZE(Var % Values)/2 ) Var%Values(2*k-1) = REAL(Var%EigenVectors(j,k)) Var%Values(2*k) = AIMAG(Var%EigenVectors(j,k)) END FORALL ELSE ALLOCATE( Var % Values(SIZE(Var % EigenVectors,2)) ) Var % Values = Var % EigenVectors(j,:) END IF CALL WriteDXFiles( TRIM(Dir)//Prefix,Model,.FALSE.,j ) DEALLOCATE( Var % Values ) Var % Values => OldValues END DO EXIT Solvers END IF END DO Solvers IF ( .NOT.EigAnal ) THEN CALL WriteDXFiles( TRIM(Dir)//Prefix, Model, .TRUE., nTime ) END IF Mesh => Mesh % next END DO END SUBROUTINE WriteData END SUBROUTINE DXOutputSolver