From 6b46a9c338aad8430c2e0d0d8f43ff56ec4675f5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Apr 2013 16:40:49 +0000 Subject: [PATCH] introduced output of initial geometry (vtk file) to mesh --- code/DAMASK_abaqus_exp.f | 1 + code/DAMASK_abaqus_std.f | 1 + code/DAMASK_marc.f90 | 1 + code/DAMASK_spectral_utilities.f90 | 13 +- code/libs.f90 | 6 +- code/mesh.f90 | 297 ++++++++++++++++------------- 6 files changed, 182 insertions(+), 137 deletions(-) diff --git a/code/DAMASK_abaqus_exp.f b/code/DAMASK_abaqus_exp.f index 77c01fa91..911c9cca7 100644 --- a/code/DAMASK_abaqus_exp.f +++ b/code/DAMASK_abaqus_exp.f @@ -92,6 +92,7 @@ end function getSolverJobName end module DAMASK_interface #include "IO.f90" +#include "libs.f90" #include "numerics.f90" #include "debug.f90" #include "math.f90" diff --git a/code/DAMASK_abaqus_std.f b/code/DAMASK_abaqus_std.f index e24f2df10..87a5fe3d1 100644 --- a/code/DAMASK_abaqus_std.f +++ b/code/DAMASK_abaqus_std.f @@ -92,6 +92,7 @@ end function getSolverJobName end module DAMASK_interface #include "IO.f90" +#include "libs.f90" #include "numerics.f90" #include "debug.f90" #include "math.f90" diff --git a/code/DAMASK_marc.f90 b/code/DAMASK_marc.f90 index e10daae61..75e8dca3c 100644 --- a/code/DAMASK_marc.f90 +++ b/code/DAMASK_marc.f90 @@ -124,6 +124,7 @@ end function getSolverJobName end module DAMASK_interface #include "IO.f90" +#include "libs.f90" #include "numerics.f90" #include "debug.f90" #include "math.f90" diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index 2f9da6ce5..3f980c51f 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -484,7 +484,7 @@ subroutine utilities_fourierConvolution(fieldAim) flush(6) !-------------------------------------------------------------------------------------------------- -! to the actual spectral method calculation (mechanical equilibrium) +! do the actual spectral method calculation (mechanical equilibrium) if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat do k = 1_pInt, res(3); do j = 1_pInt, res(2) ;do i = 1_pInt, res1_red if(any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 @@ -645,7 +645,7 @@ end function utilities_curlRMS !-------------------------------------------------------------------------------------------------- -!> @brief calculates mask compliance tensor +!> @brief calculates mask compliance tensor used to adjust F to fullfill stress BC !-------------------------------------------------------------------------------------------------- function utilities_maskedCompliance(rot_BC,mask_stress,C) use IO, only: & @@ -684,7 +684,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) if(debugGeneral) then write(6,'(/,a)') ' ... updating masked compliance ............................................' - write(6,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C rotated / GPa =',& + write(6,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',& transpose(temp99_Real)/1.e9_pReal flush(6) endif @@ -724,8 +724,9 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) if(debugGeneral .or. errmatinv) then ! report write(formatString, '(I16.16)') size_reduced formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' - write(6,trim(formatString),advance='no') ' C * S', transpose(matmul(c_reduced,s_reduced)) - write(6,trim(formatString),advance='no') ' S', transpose(s_reduced) + write(6,trim(formatString),advance='no') ' C * S (load) ', & + transpose(matmul(c_reduced,s_reduced)) + write(6,trim(formatString),advance='no') ' S (load) ', transpose(s_reduced) endif if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance') deallocate(c_reduced) @@ -735,7 +736,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) temp99_real = 0.0_pReal endif if(debugGeneral) & ! report - write(6,'(/,a,/,9(9(2x,f12.7,1x)/),/)',advance='no') ' Masked Compliance * GPa =', & + write(6,'(/,a,/,9(9(2x,f12.7,1x)/),/)',advance='no') ' Masked Compliance (load) * GPa =', & transpose(temp99_Real*1.e9_pReal) flush(6) utilities_maskedCompliance = math_Plain99to3333(temp99_Real) diff --git a/code/libs.f90 b/code/libs.f90 index 37a8bb94c..17c79d8d1 100644 --- a/code/libs.f90 +++ b/code/libs.f90 @@ -24,10 +24,10 @@ !-------------------------------------------------------------------------------------------------- #ifdef Spectral -#include "kdtree2.f90" -#include "IR_Precision.f90" -#include "Lib_VTK_IO.f90" +#include "../lib/kdtree2.f90" #endif +#include "../lib/IR_Precision.f90" +#include "../lib/Lib_VTK_IO.f90" module libs end module libs diff --git a/code/mesh.f90 b/code/mesh.f90 index 67cb7ffea..37de53f87 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -1,7 +1,7 @@ -! Copyright 2011-13 Max-Planck-Institut für Eisenforschung GmbH +! Copyright 2011-13 Max-Planck-Institut für Eisenforschung GmbH ! ! This file is part of DAMASK, -! the Düsseldorf Advanced Material Simulation Kit. +! the Düsseldorf Advanced Material Simulation Kit. ! ! DAMASK is free software: you can redistribute it and/or modify ! it under the terms of the GNU General Public License as published by @@ -19,11 +19,11 @@ !-------------------------------------------------------------------------------------------------- !* $Id$ !-------------------------------------------------------------------------------------------------- -!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH -!! Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH -!! Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH -!! Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!! Krishna Komerla, Max-Planck-Institut für Eisenforschung GmbH +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!> @author Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @author Krishna Komerla, Max-Planck-Institut für Eisenforschung GmbH !> @brief Sets up the mesh for the solvers MSC.Marc, Abaqus and the spectral solver !-------------------------------------------------------------------------------------------------- @@ -41,7 +41,7 @@ module mesh mesh_Nmaterials, & mesh_Nnodes, & !< total number of nodes in mesh mesh_Ncellnodes, & !< total number of cell nodes in mesh (including duplicates) - + mesh_Ncells, & !< total number of cells in mesh mesh_maxNnodes, & !< max number of nodes in any CP element mesh_maxNips, & !< max number of IPs in any CP element mesh_maxNipNeighbors, & !< max number of IP neighbors in any CP element @@ -81,10 +81,10 @@ module mesh initialcondTableStyle !< Table style (Marc only) #endif - type, private :: tCellnode !< set of parameters defining a cellnode - real(pReal), dimension(3) :: coords = 0.0_pReal - integer(pInt) :: elemParent = 0_pInt - integer(pInt) :: intraElemID = 0_pInt + type, private :: tCellnode !< set of parameters defining a cellnode + real(pReal), dimension(3) :: coords = 0.0_pReal !< coordinates of cell node + integer(pInt) :: elemParent = 0_pInt !< number of parent element + integer(pInt) :: intraElemID = 0_pInt !< internal number of cell node in element end type tCellnode integer(pInt), dimension(2), private :: & @@ -393,68 +393,72 @@ module mesh ],pInt) - public :: mesh_init, & - mesh_FEasCP, & - mesh_build_cells, & - mesh_build_ipVolumes, & - mesh_build_ipCoordinates, & - mesh_cellCenterCoordinates + public :: & + mesh_init, & + mesh_FEasCP, & + mesh_build_cells, & + mesh_build_ipVolumes, & + mesh_build_ipCoordinates, & + mesh_cellCenterCoordinates #ifdef Spectral - public :: mesh_regrid, & - mesh_nodesAroundCentres, & - mesh_deformedCoordsFFT, & - mesh_deformedCoordsLinear, & - mesh_volumeMismatch, & - mesh_shapeMismatch + public :: & + mesh_regrid, & + mesh_nodesAroundCentres, & + mesh_deformedCoordsFFT, & + mesh_deformedCoordsLinear, & + mesh_volumeMismatch, & + mesh_shapeMismatch #endif private :: & #ifdef Spectral - mesh_spectral_getGrid, & - mesh_spectral_getSize, & - mesh_spectral_getHomogenization, & - mesh_spectral_count_nodesAndElements, & - mesh_spectral_count_cpElements, & - mesh_spectral_map_elements, & - mesh_spectral_map_nodes, & - mesh_spectral_count_cpSizes, & - mesh_spectral_build_nodes, & - mesh_spectral_build_elements, & + mesh_spectral_getGrid, & + mesh_spectral_getSize, & + mesh_spectral_getHomogenization, & + mesh_spectral_count_nodesAndElements, & + mesh_spectral_count_cpElements, & + mesh_spectral_map_elements, & + mesh_spectral_map_nodes, & + mesh_spectral_count_cpSizes, & + mesh_spectral_build_nodes, & + mesh_spectral_build_elements, & #endif #ifdef Marc - mesh_marc_get_tableStyles, & - mesh_marc_count_nodesAndElements, & - mesh_marc_count_elementSets, & - mesh_marc_map_elementSets, & - mesh_marc_count_cpElements, & - mesh_marc_map_Elements, & - mesh_marc_map_nodes, & - mesh_marc_build_nodes, & - mesh_marc_count_cpSizes, & - mesh_marc_build_elements, & + mesh_marc_get_tableStyles, & + mesh_marc_count_nodesAndElements, & + mesh_marc_count_elementSets, & + mesh_marc_map_elementSets, & + mesh_marc_count_cpElements, & + mesh_marc_map_Elements, & + mesh_marc_map_nodes, & + mesh_marc_build_nodes, & + mesh_marc_count_cpSizes, & + mesh_marc_build_elements, & #endif #ifdef Abaqus - mesh_abaqus_count_nodesAndElements, & - mesh_abaqus_count_elementSets, & - mesh_abaqus_count_materials, & - mesh_abaqus_map_elementSets, & - mesh_abaqus_map_materials, & - mesh_abaqus_count_cpElements, & - mesh_abaqus_map_elements, & - mesh_abaqus_map_nodes, & - mesh_abaqus_build_nodes, & - mesh_abaqus_count_cpSizes, & - mesh_abaqus_build_elements, & + mesh_abaqus_count_nodesAndElements, & + mesh_abaqus_count_elementSets, & + mesh_abaqus_count_materials, & + mesh_abaqus_map_elementSets, & + mesh_abaqus_map_materials, & + mesh_abaqus_count_cpElements, & + mesh_abaqus_map_elements, & + mesh_abaqus_map_nodes, & + mesh_abaqus_build_nodes, & + mesh_abaqus_count_cpSizes, & + mesh_abaqus_build_elements, & #endif - mesh_get_damaskOptions, & - mesh_build_ipAreas, & - mesh_build_nodeTwins, & - mesh_build_sharedElems, & - mesh_build_ipNeighborhood, & - mesh_tell_statistics, & - FE_mapElemtype, & - mesh_faceMatch, & - mesh_build_FEdata + mesh_get_damaskOptions, & + mesh_build_ipAreas, & + mesh_build_nodeTwins, & + mesh_build_sharedElems, & + mesh_build_ipNeighborhood, & + mesh_tell_statistics, & + FE_mapElemtype, & + mesh_faceMatch, & + mesh_build_FEdata, & + mesh_writeGeom + contains @@ -591,7 +595,9 @@ subroutine mesh_init(ip,el) call mesh_build_ipNeighborhood call mesh_tell_statistics -if (usePingPong .and. (mesh_Nelems /= mesh_NcpElems)) call IO_error(600_pInt) ! ping-pong must be disabled when havin non-DAMASK-elements + call mesh_writeGeom + + if (usePingPong .and. (mesh_Nelems /= mesh_NcpElems)) call IO_error(600_pInt) ! ping-pong must be disabled when havin non-DAMASK-elements FEsolving_execElem = [ 1_pInt,mesh_NcpElems ] if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP) @@ -688,12 +694,14 @@ subroutine mesh_build_cells allocate(cellnodeParent(2_pInt,mesh_maxNcellnodes*mesh_NcpElems)) ; cellnodeParent = 0_pInt mesh_Ncellnodes = 0_pInt + mesh_Ncells = 0_pInt do e = 1_pInt,mesh_NcpElems ! loop over cpElems t = mesh_element(2_pInt,e) ! get element type g = FE_geomtype(t) ! get geometry type c = FE_celltype(g) ! get cell type localCellnode2globalCellnode = 0_pInt do i = 1_pInt,FE_Nips(g) ! loop over ips=cells in this element + mesh_Ncells = mesh_Ncells + FE_Nips(g) do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell localCellnodeID = FE_cell(n,i,g) if (localCellnodeID <= FE_NmatchingNodes(g)) then ! this cell node is a matching node @@ -1211,7 +1219,7 @@ subroutine mesh_spectral_build_nodes() / real(res(3),pReal) end forall - mesh_node = mesh_node0 ! why? + mesh_node = mesh_node0 end subroutine mesh_spectral_build_nodes @@ -1326,6 +1334,8 @@ function mesh_regrid(adaptive,resNewInput,minRes) use math, only: & math_periodicNearestNeighbor, & math_mul33x3 + + implicit none character(len=1024):: formatString, N_Digits logical, intent(in) :: adaptive ! if true, choose adaptive grid based on resNewInput, otherwise keep it constant integer(pInt), dimension(3), optional, intent(in) :: resNewInput ! f2py cannot handle optional arguments correctly (they are always present) @@ -1800,10 +1810,10 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes) do k = 0_pInt,iRes(3)+1_pInt do j = 0_pInt,iRes(2)+1_pInt do i = 0_pInt,iRes(1)+1_pInt - if (k==0_pInt .or. k==iRes(3)+1_pInt .or. & ! z skin - j==0_pInt .or. j==iRes(2)+1_pInt .or. & ! y skin - i==0_pInt .or. i==iRes(1)+1_pInt ) then ! x skin - me = [i,j,k] ! me on skin + if (k==0_pInt .or. k==iRes(3)+1_pInt .or. & ! z skin + j==0_pInt .or. j==iRes(2)+1_pInt .or. & ! y skin + i==0_pInt .or. i==iRes(1)+1_pInt ) then ! x skin + me = [i,j,k] ! me on skin shift = sign(abs(iRes+diag-2_pInt*me)/(iRes+diag),iRes+diag-2_pInt*me) lookup = me-diag+shift*iRes wrappedCentres(1:3,i+1_pInt, j+1_pInt, k+1_pInt) = & @@ -2883,11 +2893,10 @@ subroutine mesh_abaqus_count_nodesAndElements(myUnit) end subroutine mesh_abaqus_count_nodesAndElements -!******************************************************************** -! count overall number of element sets in mesh -! -! mesh_NelemSets, mesh_maxNelemInSet -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief count overall number of element sets in mesh and write 'mesh_NelemSets' and +!! 'mesh_maxNelemInSet' +!-------------------------------------------------------------------------------------------------- subroutine mesh_abaqus_count_elementSets(myUnit) use IO, only: IO_lc, & @@ -2904,7 +2913,7 @@ subroutine mesh_abaqus_count_elementSets(myUnit) logical :: inPart mesh_NelemSets = 0_pInt - mesh_maxNelemInSet = mesh_Nelems ! have to be conservative, since Abaqus allows for recursive definitons + mesh_maxNelemInSet = mesh_Nelems ! have to be conservative, since Abaqus allows for recursive definitons 610 FORMAT(A300) @@ -2927,11 +2936,13 @@ subroutine mesh_abaqus_count_elementSets(myUnit) end subroutine mesh_abaqus_count_elementSets +!-------------------------------------------------------------------------------------------------- !******************************************************************** ! count overall number of solid sections sets in mesh (Abaqus only) ! ! mesh_Nmaterials !******************************************************************** +!-------------------------------------------------------------------------------------------------- subroutine mesh_abaqus_count_materials(myUnit) use IO, only: IO_lc, & @@ -2971,11 +2982,13 @@ subroutine mesh_abaqus_count_materials(myUnit) end subroutine mesh_abaqus_count_materials +!-------------------------------------------------------------------------------------------------- !******************************************************************** ! Build element set mapping ! ! allocate globals: mesh_nameElemSet, mesh_mapElemSet !******************************************************************** +!-------------------------------------------------------------------------------------------------- subroutine mesh_abaqus_map_elementSets(myUnit) use IO, only: IO_lc, & @@ -3023,11 +3036,13 @@ subroutine mesh_abaqus_map_elementSets(myUnit) end subroutine mesh_abaqus_map_elementSets +!-------------------------------------------------------------------------------------------------- !******************************************************************** ! map solid section (Abaqus only) ! ! allocate globals: mesh_nameMaterial, mesh_mapMaterial !******************************************************************** +!-------------------------------------------------------------------------------------------------- subroutine mesh_abaqus_map_materials(myUnit) use IO, only: IO_lc, & @@ -3418,7 +3433,7 @@ subroutine mesh_abaqus_build_elements(myUnit) integer(pInt), parameter :: maxNchunks = 65_pInt integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - integer(pInt) :: i,j,k,c,e,t,homog,micro + integer(pInt) :: i,j,k,c,e,t,homog,micro, nNodesAlreadyRead logical inPart,materialFound character (len=64) :: materialName,elemSetName character(len=300) :: line @@ -3515,11 +3530,9 @@ subroutine mesh_abaqus_build_elements(myUnit) #endif -!******************************************************************** -! get any additional damask options from input file -! -! mesh_periodicSurface -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief get any additional damask options from input file, sets mesh_periodicSurface +!-------------------------------------------------------------------------------------------------- subroutine mesh_get_damaskOptions(myUnit) use IO, only: & @@ -3583,12 +3596,9 @@ use IO, only: & 620 end subroutine mesh_get_damaskOptions -!*********************************************************** -! calculation of IP interface areas -! -! allocate globals -! _ipArea, _ipAreaNormal -!*********************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal' +!-------------------------------------------------------------------------------------------------- subroutine mesh_build_ipAreas use math, only: & @@ -3663,12 +3673,9 @@ subroutine mesh_build_ipAreas end subroutine mesh_build_ipAreas -!*********************************************************** -! assignment of twin nodes for each cp node -! -! allocate globals -! _nodeTwins -!*********************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief assignment of twin nodes for each cp node, allocate globals '_nodeTwins' +!-------------------------------------------------------------------------------------------------- subroutine mesh_build_nodeTwins implicit none @@ -3736,14 +3743,10 @@ enddo end subroutine mesh_build_nodeTwins - -!******************************************************************** -! get maximum count of shared elements among cpElements and -! build list of elements shared by each node in mesh -! -! _maxNsharedElems -! _sharedElem -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief get maximum count of shared elements among cpElements and build list of elements shared +!! by each node in mesh. Allocate globals '_maxNsharedElems' and '_sharedElem' +!-------------------------------------------------------------------------------------------------- subroutine mesh_build_sharedElems implicit none @@ -3808,12 +3811,9 @@ deallocate(node_seen) end subroutine mesh_build_sharedElems -!*********************************************************** -! build up of IP neighborhood -! -! allocate globals -! _ipNeighborhood -!*********************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief build up of IP neighborhood, allocate globals '_ipNeighborhood' +!-------------------------------------------------------------------------------------------------- subroutine mesh_build_ipNeighborhood use math, only: math_mul3x3 @@ -3982,11 +3982,9 @@ enddo end subroutine mesh_build_ipNeighborhood -!*********************************************************** -! write statistics regarding input file parsing -! to the output file -! -!*********************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief write statistics regarding input file parsing to the output file +!-------------------------------------------------------------------------------------------------- subroutine mesh_tell_statistics use math, only: math_range @@ -4133,9 +4131,9 @@ enddo end subroutine mesh_tell_statistics -!*********************************************************** -! mapping of FE element types to internal representation -!*********************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief mapping of FE element types to internal representation +!-------------------------------------------------------------------------------------------------- integer(pInt) function FE_mapElemtype(what) use IO, only: IO_lc, IO_error @@ -4188,9 +4186,9 @@ integer(pInt) function FE_mapElemtype(what) end function FE_mapElemtype -!*********************************************************** -! find face-matching element of same type -!*********************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief find face-matching element of same type +!-------------------------------------------------------------------------------------------------- subroutine mesh_faceMatch(elem, face ,matchingElem, matchingFace) implicit none @@ -4292,11 +4290,11 @@ subroutine mesh_build_FEdata implicit none integer(pInt) :: me - allocate(FE_nodesAtIP(FE_maxmaxNnodesAtIP,FE_maxNips,FE_Ngeomtypes)) ; FE_nodesAtIP = 0_pInt - allocate(FE_ipNeighbor(FE_maxNipNeighbors,FE_maxNips,FE_Ngeomtypes)) ; FE_ipNeighbor = 0_pInt - allocate(FE_cell(FE_maxNcellnodesPerCell,FE_maxNips,FE_Ngeomtypes)) ; FE_cell = 0_pInt - allocate(FE_cellnodeParentnodeWeights(FE_maxNnodes,FE_maxNcellnodes,FE_Nelemtypes)) ; FE_cellnodeParentnodeWeights = 0.0_pReal - allocate(FE_cellface(FE_maxNcellnodesPerCellface,FE_maxNcellfaces,FE_Ncelltypes)) ; FE_cellface = 0.0_pReal + allocate(FE_nodesAtIP(FE_maxmaxNnodesAtIP,FE_maxNips,FE_Ngeomtypes)); FE_nodesAtIP = 0_pInt + allocate(FE_ipNeighbor(FE_maxNipNeighbors,FE_maxNips,FE_Ngeomtypes)); FE_ipNeighbor = 0_pInt + allocate(FE_cell(FE_maxNcellnodesPerCell,FE_maxNips,FE_Ngeomtypes)); FE_cell = 0_pInt + allocate(FE_cellnodeParentnodeWeights(FE_maxNnodes,FE_maxNcellnodes,FE_Nelemtypes)); FE_cellnodeParentnodeWeights = 0.0_pReal + allocate(FE_cellface(FE_maxNcellnodesPerCellface,FE_maxNcellfaces,FE_Ncelltypes)); FE_cellface = 0_pInt !*** fill FE_nodesAtIP with data *** @@ -4962,7 +4960,7 @@ subroutine mesh_build_FEdata me = 0_pInt me = me + 1_pInt - FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 2D 3node + FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 2D 3node, VTK_TRIANGLE (5) reshape(int([& 2,3, & 3,1, & @@ -4970,7 +4968,7 @@ subroutine mesh_build_FEdata ],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)]) me = me + 1_pInt - FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 2D 4node + FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 2D 4node, VTK_QUAD (9) reshape(int([& 2,3, & 4,1, & @@ -4979,7 +4977,7 @@ subroutine mesh_build_FEdata ],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)]) me = me + 1_pInt - FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 3D 4node + FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 3D 4node, VTK_TETRA (10) reshape(int([& 1,3,2, & 1,2,4, & @@ -4988,7 +4986,7 @@ subroutine mesh_build_FEdata ],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)]) me = me + 1_pInt - FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 3D 8node + FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 3D 8node, VTK_HEXAHEDRON (12) reshape(int([& 2,3,7,6, & 4,1,5,8, & @@ -5001,5 +4999,48 @@ subroutine mesh_build_FEdata end subroutine mesh_build_FEdata + +!-------------------------------------------------------------------------------------------------- +!> @brief writes out initial geometry +!-------------------------------------------------------------------------------------------------- +subroutine mesh_writeGeom + use DAMASK_interface, only: & + getSolverJobName + use Lib_VTK_IO, only: & + VTK_INI, & + VTK_GEO, & + VTK_CON, & + VTK_END + + implicit none + integer(pInt), dimension(1:mesh_Ncells) :: cell_type + integer(pInt), dimension(mesh_Ncells*(1_pInt+FE_maxNcellnodesPerCell)) :: connect + integer(pInt):: E_IO, g, c, e, CellID, i, NcellnodesSeen + integer(pInt), dimension(FE_Ncelltypes), parameter :: DAMASKTOVTK= [5,9,10,12] + + cellID = 0_pInt + NcellnodesSeen = 0_pInt + do e = 1_pInt, mesh_NcpElems ! loop over cpElems + g = FE_geomtype(mesh_element(2_pInt,e)) ! get element type) ! get geometry type + c = FE_celltype(g) ! get cell type + do i = 1_pInt,FE_Nips(g) ! loop over ips=cells in this element + cellID = cellID + 1_pInt + cell_type(cellID) = DAMASKTOVTK(c) + connect(NcellnodesSeen+1_pInt:NcellnodesSeen+FE_NcellnodesPerCell(c)+1_pInt) & + = [FE_NcellnodesPerCell(c),mesh_cell(1:FE_NcellnodesPerCell(c),i,e)-1_pInt] ! number of nodes per element, global element number (counting starts at 0) + NcellnodesSeen = NcellnodesSeen + FE_NcellnodesPerCell(c)+1_pInt + enddo + enddo + + E_IO = VTK_INI(output_format = 'ASCII', title=trim(getSolverJobName()), & + filename = trim(getSolverJobName())//'.vtk', mesh_topology = 'UNSTRUCTURED_GRID') + E_IO = VTK_GEO(NN = mesh_Ncellnodes, X = mesh_cellnode%coords(1), & + Y = mesh_cellnode%coords(2), Z = mesh_cellnode%coords(3)) + E_IO = VTK_CON(NC = mesh_Ncells, connect = connect(1:NcellnodesSeen), cell_type = cell_type) + E_IO = VTK_END() + +end subroutine mesh_writeGeom + + end module mesh