From aab19af13140a0bd3beec09a035ef890866a3feb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 15 Oct 2019 12:54:46 +0200 Subject: [PATCH 01/20] consistent use of subroutines, not functions. grouping file reading into meta function --- src/mesh_marc.f90 | 57 ++++++++++++++++++++++++++++++++++++----------- 1 file changed, 44 insertions(+), 13 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 482cde16c..695e8a797 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -102,10 +102,10 @@ subroutine mesh_init(ip,el) inputFile = IO_read_ASCII(trim(modelName)//trim(InputFileExtension)) ! parsing Marc input file - fileFormatVersion = inputRead_fileFormat(inputFile) + call inputRead_fileFormat(fileFormatVersion,inputFile) call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle,inputFile) if (fileFormatVersion > 12) & - marc_matNumber = inputRead_matNumber(hypoelasticTableStyle,inputFile) + call inputRead_matNumber(marc_matNumber,hypoelasticTableStyle,inputFile) call inputRead_NnodesAndElements(mesh_nNodes, mesh_nElems, inputFile) call IO_open_inputFile(FILEUNIT,modelName) @@ -222,12 +222,42 @@ subroutine writeGeometry(elemType, & end subroutine writeGeometry +subroutine inputRead() + + integer :: & + fileFormatVersion, & + hypoelasticTableStyle, & + initialcondTableStyle, & + nNodes, & + nElems + integer, parameter :: & + FILEUNIT = 222 + integer, dimension(:), allocatable :: & + matNumber !< material numbers for hypoelastic material + character(len=pStringLen), dimension(:), allocatable :: inputFile !< file content, separated per lines + + call inputRead_fileFormat(fileFormatVersion,inputFile) + call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle,inputFile) + if (fileFormatVersion > 12) & + call inputRead_matNumber(matNumber,hypoelasticTableStyle,inputFile) + call inputRead_NnodesAndElements(nNodes,nElems,inputFile) + + call IO_open_inputFile(FILEUNIT,modelName) + + !call inputRead_NelemSets(NelemSets,maxNelemInSet,FILEUNIT) + + close(FILEUNIT) + +end subroutine inputRead + + !-------------------------------------------------------------------------------------------------- !> @brief Figures out version of Marc input file format !-------------------------------------------------------------------------------------------------- -integer function inputRead_fileFormat(fileContent) +subroutine inputRead_fileFormat(fileFormat,fileContent) - character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines + integer, intent(out) :: fileFormat + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines integer, allocatable, dimension(:) :: chunkPos integer :: l @@ -235,12 +265,12 @@ integer function inputRead_fileFormat(fileContent) do l = 1, size(fileContent) chunkPos = IO_stringPos(fileContent(l)) if ( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'version') then - inputRead_fileFormat = IO_intValue(fileContent(l),chunkPos,2) + fileFormat = IO_intValue(fileContent(l),chunkPos,2) exit endif enddo -end function inputRead_fileFormat +end subroutine inputRead_fileFormat !-------------------------------------------------------------------------------------------------- @@ -272,12 +302,13 @@ end subroutine inputRead_tableStyles !-------------------------------------------------------------------------------------------------- !> @brief Figures out material number of hypoelastic material !-------------------------------------------------------------------------------------------------- -function inputRead_matNumber(tableStyle,fileContent) +subroutine inputRead_matNumber(matNumber, & + tableStyle,fileContent) - integer, intent(in) :: tableStyle - character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines + integer, allocatable, dimension(:), intent(out) :: matNumber + integer, intent(in) :: tableStyle + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines - integer, dimension(:), allocatable :: inputRead_matNumber integer, allocatable, dimension(:) :: chunkPos integer :: i, j, data_blocks, l @@ -291,17 +322,17 @@ function inputRead_matNumber(tableStyle,fileContent) else data_blocks = 1 endif - allocate(inputRead_matNumber(data_blocks), source = 0) + allocate(matNumber(data_blocks), source = 0) do i = 0, data_blocks - 1 j = i*(2+tableStyle) + 1 chunkPos = IO_stringPos(fileContent(l+1+j)) - inputRead_matNumber(i+1) = IO_intValue(fileContent(l+1+j),chunkPos,1) + matNumber(i+1) = IO_intValue(fileContent(l+1+j),chunkPos,1) enddo exit endif enddo -end function inputRead_matNumber +end subroutine inputRead_matNumber !-------------------------------------------------------------------------------------------------- From c4db2841ab66a888508a5dfcbafc2e3eb918961b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 15 Oct 2019 14:16:03 +0200 Subject: [PATCH 02/20] further cleaning --- src/mesh_marc.f90 | 109 +++++++++++++++++++++++++--------------------- 1 file changed, 59 insertions(+), 50 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 695e8a797..00f29a8c0 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -62,38 +62,37 @@ subroutine mesh_init(ip,el) integer, parameter :: FILEUNIT = 222 character(len=pStringLen), dimension(:), allocatable :: inputFile !< file content, separated per lines - integer :: & - mesh_NelemSets - real(pReal), dimension(:,:), allocatable :: & + + real(pReal), dimension(:,:), allocatable :: & node0_elem, & !< node x,y,z coordinates (initially!) node0_cell integer :: j, fileFormatVersion, elemType, & - mesh_maxNelemInSet, & - mesh_nElems, & + + nElems, & hypoelasticTableStyle, & initialcondTableStyle integer, dimension(:), allocatable :: & marc_matNumber !< array of material numbers for hypoelastic material (Marc only) - integer, dimension(:), allocatable :: & - microstructureAt, & - homogenizationAt - character(len=64), dimension(:), allocatable :: & - mesh_nameElemSet - integer, dimension(:,:), allocatable :: & - mesh_mapElemSet !< list of elements in elementSet - integer:: & - mesh_Nnodes !< total number of nodes in mesh + integer, dimension(:), allocatable :: & + microstructureAt, & + homogenizationAt + character(len=64), dimension(:), allocatable :: & + nameElemSet + integer, dimension(:,:), allocatable :: & + mapElemSet !< list of elements in elementSet + integer:: & + Nnodes !< total number of nodes in mesh - real(pReal), dimension(:,:), allocatable :: & - ip_reshaped + real(pReal), dimension(:,:), allocatable :: & + ip_reshaped integer,dimension(:,:,:), allocatable :: & - connectivity_cell !< cell connectivity for each element,ip/cell - integer, dimension(:,:), allocatable :: & - connectivity_elem + connectivity_cell !< cell connectivity for each element,ip/cell + integer, dimension(:,:), allocatable :: & + connectivity_elem - type(tMesh) :: theMesh + type(tMesh) :: theMesh write(6,'(/,a)') ' <<<+- mesh init -+>>>' @@ -106,41 +105,40 @@ subroutine mesh_init(ip,el) call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle,inputFile) if (fileFormatVersion > 12) & call inputRead_matNumber(marc_matNumber,hypoelasticTableStyle,inputFile) - call inputRead_NnodesAndElements(mesh_nNodes, mesh_nElems, inputFile) - + call inputRead_NnodesAndElements(nNodes, nElems, inputFile) + + allocate (mesh_mapFEtoCPelem(2,nElems), source = 0) + allocate (mesh_mapFEtoCPnode(2,Nnodes), source = 0) + call IO_open_inputFile(FILEUNIT,modelName) - call inputRead_NelemSets(mesh_NelemSets,mesh_maxNelemInSet,FILEUNIT) - allocate(mesh_nameElemSet(mesh_NelemSets)); mesh_nameElemSet = 'n/a' - allocate(mesh_mapElemSet(1+mesh_maxNelemInSet,mesh_NelemSets),source=0) - call inputRead_mapElemSets(mesh_nameElemSet,mesh_mapElemSet,FILEUNIT) - allocate (mesh_mapFEtoCPelem(2,mesh_nElems), source = 0) - call inputRead_mapElems(hypoelasticTableStyle,mesh_nameElemSet,mesh_mapElemSet,& - mesh_nElems,fileFormatVersion,marc_matNumber,FILEUNIT) - allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes),source=0) - call inputRead_mapNodes(mesh_Nnodes,inputFile) !ToDo: don't work on global variables + call inputRead_mapElemSets(nameElemSet,mapElemSet,FILEUNIT) + call inputRead_mapElems(hypoelasticTableStyle,nameElemSet,mapElemSet,& + nElems,fileFormatVersion,marc_matNumber,FILEUNIT) + + call inputRead_mapNodes(Nnodes,inputFile) - node0_elem = inputRead_elemNodes(mesh_Nnodes,inputFile) + node0_elem = inputRead_elemNodes(Nnodes,inputFile) - elemType = inputRead_elemType(mesh_nElems,FILEUNIT) + elemType = inputRead_elemType(nElems,FILEUNIT) call theMesh%init('mesh',elemType,node0_elem) - call theMesh%setNelems(mesh_nElems) + call theMesh%setNelems(nElems) - allocate(microstructureAt(theMesh%nElems), source=0) - allocate(homogenizationAt(theMesh%nElems), source=0) + allocate(microstructureAt(nElems), source=0) + allocate(homogenizationAt(nElems), source=0) - connectivity_elem = inputRead_connectivityElem(mesh_nElems,theMesh%elem%nNodes,FILEUNIT) + connectivity_elem = inputRead_connectivityElem(nElems,theMesh%elem%nNodes,FILEUNIT) call inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, & - mesh_nElems,theMesh%elem%nNodes,mesh_nameElemSet,mesh_mapElemSet,& + nElems,theMesh%elem%nNodes,nameElemSet,mapElemSet,& initialcondTableStyle,FILEUNIT) close (FILEUNIT) - allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal) + allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,nElems),source=0.0_pReal) allocate(cellNodeDefinition(theMesh%elem%nNodes-1)) - allocate(connectivity_cell(theMesh%elem%NcellNodesPerCell,theMesh%elem%nIPs,theMesh%nElems)) + allocate(connectivity_cell(theMesh%elem%NcellNodesPerCell,theMesh%elem%nIPs,nElems)) call buildCells(connectivity_cell,cellNodeDefinition,& theMesh%elem,connectivity_elem) allocate(node0_cell(3,maxval(connectivity_cell))) @@ -149,10 +147,7 @@ subroutine mesh_init(ip,el) allocate(ip_reshaped(3,theMesh%elem%nIPs*theMesh%nElems),source=0.0_pReal) call buildIPcoordinates(ip_reshaped,reshape(connectivity_cell,[theMesh%elem%NcellNodesPerCell,& theMesh%elem%nIPs*theMesh%nElems]),node0_cell) - - - if (usePingPong .and. (mesh_Nelems /= theMesh%nElems)) & - call IO_error(600) ! ping-pong must be disabled when having non-DAMASK elements + if (debug_e < 1 .or. debug_e > theMesh%nElems) & call IO_error(602,ext_msg='element') ! selected element does not exist if (debug_i < 1 .or. debug_i > theMesh%elem%nIPs) & @@ -235,6 +230,11 @@ subroutine inputRead() integer, dimension(:), allocatable :: & matNumber !< material numbers for hypoelastic material character(len=pStringLen), dimension(:), allocatable :: inputFile !< file content, separated per lines + + character(len=64), dimension(:), allocatable :: & + nameElemSet + integer, dimension(:,:), allocatable :: & + mapElemSet !< list of elements in elementSet call inputRead_fileFormat(fileFormatVersion,inputFile) call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle,inputFile) @@ -242,9 +242,14 @@ subroutine inputRead() call inputRead_matNumber(matNumber,hypoelasticTableStyle,inputFile) call inputRead_NnodesAndElements(nNodes,nElems,inputFile) + allocate (mesh_mapFEtoCPelem(2,nElems), source = 0) + allocate (mesh_mapFEtoCPnode(2,Nnodes), source = 0) + call IO_open_inputFile(FILEUNIT,modelName) - !call inputRead_NelemSets(NelemSets,maxNelemInSet,FILEUNIT) + call inputRead_mapElemSets(nameElemSet,mapElemSet,FILEUNIT) + call inputRead_mapElems(hypoelasticTableStyle,nameElemSet,mapElemSet,& + nElems,fileFormatVersion,matNumber,FILEUNIT) close(FILEUNIT) @@ -398,14 +403,18 @@ subroutine inputRead_NelemSets(nElemSets,maxNelemInSet,& !-------------------------------------------------------------------------------------------------- subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,fileUnit) - character(len=64), dimension(:), intent(out) :: nameElemSet - integer, dimension(:,:), intent(out) :: mapElemSet - integer, intent(in) :: fileUnit + character(len=64), dimension(:), allocatable, intent(out) :: nameElemSet + integer, dimension(:,:), allocatable, intent(out) :: mapElemSet + integer, intent(in) :: fileUnit integer, allocatable, dimension(:) :: chunkPos character(len=300) :: line - integer :: elemSet - + integer :: elemSet, NelemSets, maxNelemInSet + + + call inputRead_NelemSets(NelemSets,maxNelemInSet,fileUnit) + allocate(nameElemSet(NelemSets)); nameElemSet = 'n/a' + allocate(mapElemSet(1+maxNelemInSet,NelemSets),source=0) elemSet = 0 rewind(fileUnit) From 1d70d9b6aea81423fbe553ec4e08ace03f52ecf0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 15 Oct 2019 21:58:49 +0200 Subject: [PATCH 03/20] mesh type not very beneficial --- src/mesh_marc.f90 | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 00f29a8c0..9b8ff1022 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -44,7 +44,6 @@ module mesh mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) !-------------------------------------------------------------------------------------------------- - public :: & mesh_init, & mesh_FEasCP @@ -144,20 +143,20 @@ subroutine mesh_init(ip,el) allocate(node0_cell(3,maxval(connectivity_cell))) call buildCellNodes(node0_cell,& cellNodeDefinition,node0_elem) - allocate(ip_reshaped(3,theMesh%elem%nIPs*theMesh%nElems),source=0.0_pReal) + allocate(ip_reshaped(3,theMesh%elem%nIPs*nElems),source=0.0_pReal) call buildIPcoordinates(ip_reshaped,reshape(connectivity_cell,[theMesh%elem%NcellNodesPerCell,& - theMesh%elem%nIPs*theMesh%nElems]),node0_cell) + theMesh%elem%nIPs*nElems]),node0_cell) - if (debug_e < 1 .or. debug_e > theMesh%nElems) & + if (debug_e < 1 .or. debug_e > nElems) & call IO_error(602,ext_msg='element') ! selected element does not exist if (debug_i < 1 .or. debug_i > theMesh%elem%nIPs) & call IO_error(602,ext_msg='IP') ! selected element does not have requested IP - FEsolving_execElem = [ 1,theMesh%nElems ] ! parallel loop bounds set to comprise all DAMASK elements - allocate(FEsolving_execIP(2,theMesh%nElems), source=1) ! parallel loop bounds set to comprise from first IP... + FEsolving_execElem = [ 1,nElems ] ! parallel loop bounds set to comprise all DAMASK elements + allocate(FEsolving_execIP(2,nElems), source=1) ! parallel loop bounds set to comprise from first IP... FEsolving_execIP(2,:) = theMesh%elem%nIPs - allocate(calcMode(theMesh%elem%nIPs,theMesh%nElems)) + allocate(calcMode(theMesh%elem%nIPs,nElems)) calcMode = .false. ! pretend to have collected what first call is asking (F = I) calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc" From ffb4b2a455c11b6d007171c6330500fb059d48fd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 16 Oct 2019 03:57:51 +0200 Subject: [PATCH 04/20] theMesh is not useful maybe a shared mesh class becomes useful when Abaqus and PETSc mesh require similar functionality elem type can be reported directly --- src/mesh_marc.f90 | 59 ++++++++++++++++++++++------------------------- 1 file changed, 28 insertions(+), 31 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 9b8ff1022..b7d72b048 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -65,8 +65,9 @@ subroutine mesh_init(ip,el) real(pReal), dimension(:,:), allocatable :: & node0_elem, & !< node x,y,z coordinates (initially!) node0_cell + type(tElement) :: elem - integer :: j, fileFormatVersion, elemType, & + integer :: j, fileFormatVersion, & nElems, & hypoelasticTableStyle, & @@ -91,8 +92,6 @@ subroutine mesh_init(ip,el) connectivity_elem - type(tMesh) :: theMesh - write(6,'(/,a)') ' <<<+- mesh init -+>>>' @@ -119,44 +118,41 @@ subroutine mesh_init(ip,el) node0_elem = inputRead_elemNodes(Nnodes,inputFile) - elemType = inputRead_elemType(nElems,FILEUNIT) - - call theMesh%init('mesh',elemType,node0_elem) - call theMesh%setNelems(nElems) + call inputRead_elemType(elem,nElems,FILEUNIT) allocate(microstructureAt(nElems), source=0) allocate(homogenizationAt(nElems), source=0) - connectivity_elem = inputRead_connectivityElem(nElems,theMesh%elem%nNodes,FILEUNIT) + connectivity_elem = inputRead_connectivityElem(nElems,elem%nNodes,FILEUNIT) call inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, & - nElems,theMesh%elem%nNodes,nameElemSet,mapElemSet,& + nElems,elem%nNodes,nameElemSet,mapElemSet,& initialcondTableStyle,FILEUNIT) close (FILEUNIT) - allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,nElems),source=0.0_pReal) + allocate(mesh_ipCoordinates(3,elem%nIPs,nElems),source=0.0_pReal) - allocate(cellNodeDefinition(theMesh%elem%nNodes-1)) - allocate(connectivity_cell(theMesh%elem%NcellNodesPerCell,theMesh%elem%nIPs,nElems)) + allocate(cellNodeDefinition(elem%nNodes-1)) + allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems)) call buildCells(connectivity_cell,cellNodeDefinition,& - theMesh%elem,connectivity_elem) + elem,connectivity_elem) allocate(node0_cell(3,maxval(connectivity_cell))) call buildCellNodes(node0_cell,& cellNodeDefinition,node0_elem) - allocate(ip_reshaped(3,theMesh%elem%nIPs*nElems),source=0.0_pReal) - call buildIPcoordinates(ip_reshaped,reshape(connectivity_cell,[theMesh%elem%NcellNodesPerCell,& - theMesh%elem%nIPs*nElems]),node0_cell) + allocate(ip_reshaped(3,elem%nIPs*nElems),source=0.0_pReal) + call buildIPcoordinates(ip_reshaped,reshape(connectivity_cell,[elem%NcellNodesPerCell,& + elem%nIPs*nElems]),node0_cell) if (debug_e < 1 .or. debug_e > nElems) & call IO_error(602,ext_msg='element') ! selected element does not exist - if (debug_i < 1 .or. debug_i > theMesh%elem%nIPs) & + if (debug_i < 1 .or. debug_i > elem%nIPs) & call IO_error(602,ext_msg='IP') ! selected element does not have requested IP FEsolving_execElem = [ 1,nElems ] ! parallel loop bounds set to comprise all DAMASK elements allocate(FEsolving_execIP(2,nElems), source=1) ! parallel loop bounds set to comprise from first IP... - FEsolving_execIP(2,:) = theMesh%elem%nIPs + FEsolving_execIP(2,:) = elem%nIPs - allocate(calcMode(theMesh%elem%nIPs,nElems)) + allocate(calcMode(elem%nIPs,nElems)) calcMode = .false. ! pretend to have collected what first call is asking (F = I) calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc" @@ -165,7 +161,7 @@ subroutine mesh_init(ip,el) node0_elem) call writeGeometry(0,connectivity_elem,& - reshape(connectivity_cell,[theMesh%elem%NcellNodesPerCell,theMesh%elem%nIPs*theMesh%nElems]),& + reshape(connectivity_cell,[elem%NcellNodesPerCell,elem%nIPs*nElems]),& node0_cell,ip_reshaped) end subroutine mesh_init @@ -241,20 +237,22 @@ subroutine inputRead() call inputRead_matNumber(matNumber,hypoelasticTableStyle,inputFile) call inputRead_NnodesAndElements(nNodes,nElems,inputFile) - allocate (mesh_mapFEtoCPelem(2,nElems), source = 0) - allocate (mesh_mapFEtoCPnode(2,Nnodes), source = 0) - call IO_open_inputFile(FILEUNIT,modelName) call inputRead_mapElemSets(nameElemSet,mapElemSet,FILEUNIT) + allocate (mesh_mapFEtoCPelem(2,nElems), source = 0) call inputRead_mapElems(hypoelasticTableStyle,nameElemSet,mapElemSet,& nElems,fileFormatVersion,matNumber,FILEUNIT) + allocate (mesh_mapFEtoCPnode(2,Nnodes), source = 0) + call inputRead_mapNodes(Nnodes,inputFile) + close(FILEUNIT) end subroutine inputRead + !-------------------------------------------------------------------------------------------------- !> @brief Figures out version of Marc input file format !-------------------------------------------------------------------------------------------------- @@ -555,13 +553,14 @@ end function inputRead_elemNodes !-------------------------------------------------------------------------------------------------- !> @brief Gets element type (and checks if the whole mesh comprises of only one type) !-------------------------------------------------------------------------------------------------- -integer function inputRead_elemType(nElem,fileUnit) +subroutine inputRead_elemType(elem, & + nElem,fileUnit) - integer, intent(in) :: & + type(tElement), intent(out) :: elem + integer, intent(in) :: & nElem, & fileUnit - type(tElement) :: tempEl integer, allocatable, dimension(:) :: chunkPos character(len=300) :: line integer :: i,t @@ -579,12 +578,11 @@ integer function inputRead_elemType(nElem,fileUnit) chunkPos = IO_stringPos(line) if (t == -1) then t = mapElemtype(IO_stringValue(line,chunkPos,2)) - call tempEl%init(t) - inputRead_elemType = t + call elem%init(t) else if (t /= mapElemtype(IO_stringValue(line,chunkPos,2))) call IO_error(191,el=t,ip=i) endif - call IO_skipChunks(fileUnit,tempEl%nNodes-(chunkPos(1)-2)) + call IO_skipChunks(fileUnit,elem%nNodes-(chunkPos(1)-2)) enddo exit endif @@ -636,7 +634,7 @@ integer function inputRead_elemType(nElem,fileUnit) end function mapElemtype -620 end function inputRead_elemType +620 end subroutine inputRead_elemType !-------------------------------------------------------------------------------------------------- @@ -944,7 +942,6 @@ subroutine buildCellNodes(node_cell, & end subroutine buildCellNodes - !-------------------------------------------------------------------------------------------------- !> @brief Calculates IP coordinates as center of cell !-------------------------------------------------------------------------------------------------- From 0f77d2efdb389a5672fd49806b1f8cfe7a2d039f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 16 Oct 2019 04:22:21 +0200 Subject: [PATCH 05/20] consistent use of subroutines clearly reveal input(out) arguments --- src/mesh_marc.f90 | 61 ++++++++++++++++++++++++++--------------------- 1 file changed, 34 insertions(+), 27 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index b7d72b048..5d3572956 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -111,15 +111,16 @@ subroutine mesh_init(ip,el) call IO_open_inputFile(FILEUNIT,modelName) call inputRead_mapElemSets(nameElemSet,mapElemSet,FILEUNIT) call inputRead_mapElems(hypoelasticTableStyle,nameElemSet,mapElemSet,& - nElems,fileFormatVersion,marc_matNumber,FILEUNIT) + fileFormatVersion,marc_matNumber,FILEUNIT) + + call inputRead_mapNodes(inputFile) - call inputRead_mapNodes(Nnodes,inputFile) - - node0_elem = inputRead_elemNodes(Nnodes,inputFile) - - call inputRead_elemType(elem,nElems,FILEUNIT) + call inputRead_elemNodes(node0_elem,Nnodes,inputFile) + + + allocate(microstructureAt(nElems), source=0) allocate(homogenizationAt(nElems), source=0) @@ -231,21 +232,25 @@ subroutine inputRead() integer, dimension(:,:), allocatable :: & mapElemSet !< list of elements in elementSet - call inputRead_fileFormat(fileFormatVersion,inputFile) - call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle,inputFile) + call inputRead_fileFormat(fileFormatVersion, & + inputFile) + call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle, & + inputFile) if (fileFormatVersion > 12) & - call inputRead_matNumber(matNumber,hypoelasticTableStyle,inputFile) - call inputRead_NnodesAndElements(nNodes,nElems,inputFile) + call inputRead_matNumber(matNumber, & + hypoelasticTableStyle,inputFile) + call inputRead_NnodesAndElements(nNodes,nElems,& + inputFile) call IO_open_inputFile(FILEUNIT,modelName) - call inputRead_mapElemSets(nameElemSet,mapElemSet,FILEUNIT) + call inputRead_mapElemSets(nameElemSet,mapElemSet,& + FILEUNIT) allocate (mesh_mapFEtoCPelem(2,nElems), source = 0) - call inputRead_mapElems(hypoelasticTableStyle,nameElemSet,mapElemSet,& - nElems,fileFormatVersion,matNumber,FILEUNIT) + call inputRead_mapElems(hypoelasticTableStyle,nameElemSet,mapElemSet,fileFormatVersion,matNumber,FILEUNIT) allocate (mesh_mapFEtoCPnode(2,Nnodes), source = 0) - call inputRead_mapNodes(Nnodes,inputFile) + call inputRead_mapNodes(inputFile) close(FILEUNIT) @@ -433,9 +438,9 @@ subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,fileUnit) !-------------------------------------------------------------------------------------------------- !> @brief Maps elements from FE ID to internal (consecutive) representation. !-------------------------------------------------------------------------------------------------- -subroutine inputRead_mapElems(tableStyle,nameElemSet,mapElemSet,nElems,fileFormatVersion,matNumber,fileUnit) +subroutine inputRead_mapElems(tableStyle,nameElemSet,mapElemSet,fileFormatVersion,matNumber,fileUnit) - integer, intent(in) :: fileUnit,tableStyle,nElems,fileFormatVersion + integer, intent(in) :: fileUnit,tableStyle,fileFormatVersion integer, dimension(:), intent(in) :: matNumber character(len=64), intent(in), dimension(:) :: nameElemSet integer, dimension(:,:), intent(in) :: & @@ -445,8 +450,8 @@ subroutine inputRead_mapElems(tableStyle,nameElemSet,mapElemSet,nElems,fileForma character(len=300) :: line, & tmp - integer, dimension (1+nElems) :: contInts - integer :: i,cpElem + integer, dimension (1+size(mesh_mapFEtoCPelem,2)) :: contInts + integer :: i,cpElem cpElem = 0 contInts = 0 @@ -459,7 +464,7 @@ subroutine inputRead_mapElems(tableStyle,nameElemSet,mapElemSet,nElems,fileForma do i=1,3+TableStyle ! skip three (or four if new table style!) lines read (fileUnit,'(A300)') line enddo - contInts = IO_continuousIntValues(fileUnit,nElems,nameElemSet,& + contInts = IO_continuousIntValues(fileUnit,size(mesh_mapFEtoCPelem,2),nameElemSet,& mapElemSet,size(nameElemSet)) exit endif @@ -497,9 +502,8 @@ end subroutine inputRead_mapElems !-------------------------------------------------------------------------------------------------- !> @brief Maps node from FE ID to internal (consecutive) representation. !-------------------------------------------------------------------------------------------------- -subroutine inputRead_mapNodes(nNodes,fileContent) +subroutine inputRead_mapNodes(fileContent) - integer, intent(in) :: nNodes character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines integer, allocatable, dimension(:) :: chunkPos @@ -508,7 +512,7 @@ subroutine inputRead_mapNodes(nNodes,fileContent) do l = 1, size(fileContent) chunkPos = IO_stringPos(fileContent(l)) if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then - do i = 1,nNodes + do i = 1,size(mesh_mapFEtoCPnode,2) mesh_mapFEtoCPnode(1,i) = IO_fixedIntValue (fileContent(l+1+i),[0,10],1) mesh_mapFEtoCPnode(2,i) = i enddo @@ -524,16 +528,19 @@ end subroutine inputRead_mapNodes !-------------------------------------------------------------------------------------------------- !> @brief store x,y,z coordinates of all nodes in mesh. !-------------------------------------------------------------------------------------------------- -function inputRead_elemNodes(nNode,fileContent) result(nodes) +subroutine inputRead_elemNodes(nodes, & + nNode,fileContent) - integer, intent(in) :: nNode - character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines + real(pReal), allocatable, dimension(:,:), intent(out) :: nodes + integer, intent(in) :: nNode + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines - real(pReal), dimension(3,nNode) :: nodes integer, dimension(5), parameter :: node_ends = [0,10,30,50,70] integer, allocatable, dimension(:) :: chunkPos integer :: i,j,m,l + allocate(nodes(3,nNode)) + do l = 1, size(fileContent) chunkPos = IO_stringPos(fileContent(l)) if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then @@ -547,7 +554,7 @@ function inputRead_elemNodes(nNode,fileContent) result(nodes) endif enddo -end function inputRead_elemNodes +end subroutine inputRead_elemNodes !-------------------------------------------------------------------------------------------------- From 3abb549eab7a197aab0144cb5299da93fc23d2e4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 16 Oct 2019 04:33:26 +0200 Subject: [PATCH 06/20] better encapsulation inputRead fully parses the input file and gives back the required data --- src/mesh_marc.f90 | 103 ++++++++++++++++++---------------------------- 1 file changed, 40 insertions(+), 63 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 5d3572956..24284f889 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -59,28 +59,15 @@ subroutine mesh_init(ip,el) integer, intent(in) :: el, ip - integer, parameter :: FILEUNIT = 222 - character(len=pStringLen), dimension(:), allocatable :: inputFile !< file content, separated per lines - real(pReal), dimension(:,:), allocatable :: & node0_elem, & !< node x,y,z coordinates (initially!) node0_cell type(tElement) :: elem - integer :: j, fileFormatVersion, & - - nElems, & - hypoelasticTableStyle, & - initialcondTableStyle - integer, dimension(:), allocatable :: & - marc_matNumber !< array of material numbers for hypoelastic material (Marc only) + integer :: nElems integer, dimension(:), allocatable :: & microstructureAt, & homogenizationAt - character(len=64), dimension(:), allocatable :: & - nameElemSet - integer, dimension(:,:), allocatable :: & - mapElemSet !< list of elements in elementSet integer:: & Nnodes !< total number of nodes in mesh @@ -92,46 +79,13 @@ subroutine mesh_init(ip,el) connectivity_elem - write(6,'(/,a)') ' <<<+- mesh init -+>>>' - mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh - inputFile = IO_read_ASCII(trim(modelName)//trim(InputFileExtension)) - - ! parsing Marc input file - call inputRead_fileFormat(fileFormatVersion,inputFile) - call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle,inputFile) - if (fileFormatVersion > 12) & - call inputRead_matNumber(marc_matNumber,hypoelasticTableStyle,inputFile) - call inputRead_NnodesAndElements(nNodes, nElems, inputFile) - - allocate (mesh_mapFEtoCPelem(2,nElems), source = 0) - allocate (mesh_mapFEtoCPnode(2,Nnodes), source = 0) - - call IO_open_inputFile(FILEUNIT,modelName) - call inputRead_mapElemSets(nameElemSet,mapElemSet,FILEUNIT) - call inputRead_mapElems(hypoelasticTableStyle,nameElemSet,mapElemSet,& - fileFormatVersion,marc_matNumber,FILEUNIT) + mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh + call inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogenizationAt) + nElems = size(connectivity_elem,2) - call inputRead_mapNodes(inputFile) - - call inputRead_elemType(elem,nElems,FILEUNIT) - - call inputRead_elemNodes(node0_elem,Nnodes,inputFile) - - - - allocate(microstructureAt(nElems), source=0) - allocate(homogenizationAt(nElems), source=0) - - connectivity_elem = inputRead_connectivityElem(nElems,elem%nNodes,FILEUNIT) - call inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, & - nElems,elem%nNodes,nameElemSet,mapElemSet,& - initialcondTableStyle,FILEUNIT) - close (FILEUNIT) - - - allocate(mesh_ipCoordinates(3,elem%nIPs,nElems),source=0.0_pReal) + allocate(mesh_ipCoordinates(3,elem%nIPs,nElems),source=0.0_pReal) ! deprecated allocate(cellNodeDefinition(elem%nNodes-1)) allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems)) @@ -144,18 +98,15 @@ subroutine mesh_init(ip,el) call buildIPcoordinates(ip_reshaped,reshape(connectivity_cell,[elem%NcellNodesPerCell,& elem%nIPs*nElems]),node0_cell) - if (debug_e < 1 .or. debug_e > nElems) & - call IO_error(602,ext_msg='element') ! selected element does not exist - if (debug_i < 1 .or. debug_i > elem%nIPs) & - call IO_error(602,ext_msg='IP') ! selected element does not have requested IP + if (debug_e < 1 .or. debug_e > nElems) call IO_error(602,ext_msg='element') + if (debug_i < 1 .or. debug_i > elem%nIPs) call IO_error(602,ext_msg='IP') - FEsolving_execElem = [ 1,nElems ] ! parallel loop bounds set to comprise all DAMASK elements - allocate(FEsolving_execIP(2,nElems), source=1) ! parallel loop bounds set to comprise from first IP... + FEsolving_execElem = [ 1,nElems ] ! parallel loop bounds set to comprise all DAMASK elements + allocate(FEsolving_execIP(2,nElems), source=1) ! parallel loop bounds set to comprise from first IP... FEsolving_execIP(2,:) = elem%nIPs - allocate(calcMode(elem%nIPs,nElems)) - calcMode = .false. ! pretend to have collected what first call is asking (F = I) - calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc" + allocate(calcMode(elem%nIPs,nElems),source=.false.) ! pretend to have collected what first call is asking (F = I) + calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc" call discretization_init(microstructureAt,homogenizationAt,& ip_reshaped,& @@ -213,7 +164,18 @@ subroutine writeGeometry(elemType, & end subroutine writeGeometry -subroutine inputRead() +subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogenizationAt) + + type(tElement), intent(out) :: elem + real(pReal), dimension(:,:), allocatable, intent(out) :: & + node0_elem !< node x,y,z coordinates (initially!) + integer, dimension(:,:), allocatable, intent(out) :: & + connectivity_elem + + integer, dimension(:), allocatable, intent(out) :: & + microstructureAt, & + homogenizationAt + integer :: & fileFormatVersion, & @@ -242,16 +204,27 @@ subroutine inputRead() call inputRead_NnodesAndElements(nNodes,nElems,& inputFile) - call IO_open_inputFile(FILEUNIT,modelName) + call IO_open_inputFile(FILEUNIT,modelName) ! ToDo: It would be better to use fileContent call inputRead_mapElemSets(nameElemSet,mapElemSet,& FILEUNIT) + allocate (mesh_mapFEtoCPelem(2,nElems), source = 0) call inputRead_mapElems(hypoelasticTableStyle,nameElemSet,mapElemSet,fileFormatVersion,matNumber,FILEUNIT) allocate (mesh_mapFEtoCPnode(2,Nnodes), source = 0) call inputRead_mapNodes(inputFile) + call inputRead_elemType(elem, & + nElems,FILEUNIT) + call inputRead_elemNodes(node0_elem, & + Nnodes,inputFile) + + connectivity_elem = inputRead_connectivityElem(nElems,elem%nNodes,FILEUNIT) + + call inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, & + nElems,elem%nNodes,nameElemSet,mapElemSet,& + initialcondTableStyle,FILEUNIT) close(FILEUNIT) end subroutine inputRead @@ -704,7 +677,7 @@ function inputRead_connectivityElem(nElem,nNodes,fileUnit) subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, & nElem,nNodes,nameElemSet,mapElemSet,initialcondTableStyle,fileUnit) - integer, dimension(:), intent(out) :: & + integer, dimension(:), allocatable, intent(out) :: & microstructureAt, & homogenizationAt integer, intent(in) :: & @@ -723,6 +696,10 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza integer, dimension(1+nElem) :: contInts integer :: i,j,t,sv,myVal,e,nNodesAlreadyRead + + allocate(microstructureAt(nElem),source=0) + allocate(homogenizationAt(nElem),source=0) + rewind(fileUnit) read (fileUnit,'(A300)',END=630) line do From 3a1c4f95c38c494df5bc75678bcec3f9f46f3e1e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 16 Oct 2019 08:40:27 +0200 Subject: [PATCH 07/20] some polishing --- src/element.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/element.f90 b/src/element.f90 index d62b4fd93..3252d2813 100644 --- a/src/element.f90 +++ b/src/element.f90 @@ -146,11 +146,11 @@ module element 8 & ! 3D 8node ] !< number of cell nodes in a specific cell type - ! *** FE_ipNeighbor *** - ! is a list of the neighborhood of each IP. + ! *** IPneighbor *** + ! list of the neighborhood of each IP. ! It is sorted in (local) +x,-x, +y,-y, +z,-z direction. - ! Positive integers denote an intra-FE IP identifier. - ! Negative integers denote the interface behind which the neighboring (extra-FE) IP will be located. + ! Positive integers denote an intra-element IP identifier. + ! Negative integers denote the interface behind which the neighboring (extra-element) IP will be located. integer, dimension(nIPneighbor(cellType(1)),nIP(1)), parameter, private :: IPneighbor1 = & reshape([& From 4e213514bd40a8a1698d4693d0208877cd9708f2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 16 Oct 2019 13:04:20 +0200 Subject: [PATCH 08/20] cleaning --- src/mesh_FEM.f90 | 27 ++++----------------------- 1 file changed, 4 insertions(+), 23 deletions(-) diff --git a/src/mesh_FEM.f90 b/src/mesh_FEM.f90 index 8c0b69ca3..4b4ff18aa 100644 --- a/src/mesh_FEM.f90 +++ b/src/mesh_FEM.f90 @@ -28,8 +28,7 @@ module mesh integer, public, protected :: & mesh_Nboundaries, & mesh_NcpElems, & !< total number of CP elements in mesh - mesh_NcpElemsGlobal, & - mesh_Nnodes !< total number of nodes in mesh + mesh_NcpElemsGlobal !!!! BEGIN DEPRECATED !!!!! integer, public, protected :: & @@ -69,8 +68,7 @@ module mesh public :: & mesh_init, & mesh_FEM_build_ipVolumes, & - mesh_FEM_build_ipCoordinates, & - mesh_cellCenterCoordinates + mesh_FEM_build_ipCoordinates contains @@ -107,7 +105,8 @@ subroutine mesh_init integer, parameter :: FILEUNIT = 222 integer :: j integer, allocatable, dimension(:) :: chunkPos - integer :: dimPlex + integer :: dimPlex, & + mesh_Nnodes !< total number of nodes in mesh integer, parameter :: & mesh_ElemType=1 !< Element type of the mesh (only support homogeneous meshes) character(len=512) :: & @@ -228,26 +227,8 @@ subroutine mesh_init end subroutine mesh_init -!-------------------------------------------------------------------------------------------------- -!> @brief Calculates cell center coordinates. -!-------------------------------------------------------------------------------------------------- -pure function mesh_cellCenterCoordinates(ip,el) - - integer, intent(in) :: el, & !< element number - ip !< integration point number - real(pReal), dimension(3) :: mesh_cellCenterCoordinates !< x,y,z coordinates of the cell center of the requested IP cell - -end function mesh_cellCenterCoordinates - - !-------------------------------------------------------------------------------------------------- !> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume' -!> @details The IP volume is calculated differently depending on the cell type. -!> 2D cells assume an element depth of one in order to calculate the volume. -!> For the hexahedral cell we subdivide the cell into subvolumes of pyramidal -!> shape with a cell face as basis and the central ip at the tip. This subvolume is -!> calculated as an average of four tetrahedals with three corners on the cell face -!> and one corner at the central ip. !-------------------------------------------------------------------------------------------------- subroutine mesh_FEM_build_ipVolumes(dimPlex) From dbe15f88f2191c39ba7fca15edbc58de1bfc2390 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 16 Oct 2019 17:19:19 +0200 Subject: [PATCH 09/20] bugfix: forgot to read file first draft of nonlocal functions for area normal, area, and volume --- src/mesh_marc.f90 | 151 +++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 144 insertions(+), 7 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 24284f889..517865718 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -176,7 +176,6 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni microstructureAt, & homogenizationAt - integer :: & fileFormatVersion, & hypoelasticTableStyle, & @@ -194,6 +193,7 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni integer, dimension(:,:), allocatable :: & mapElemSet !< list of elements in elementSet + inputFile = IO_read_ASCII(trim(modelName)//trim(InputFileExtension)) call inputRead_fileFormat(fileFormatVersion, & inputFile) call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle, & @@ -209,17 +209,17 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni call inputRead_mapElemSets(nameElemSet,mapElemSet,& FILEUNIT) - allocate (mesh_mapFEtoCPelem(2,nElems), source = 0) + allocate (mesh_mapFEtoCPelem(2,nElems), source=0) call inputRead_mapElems(hypoelasticTableStyle,nameElemSet,mapElemSet,fileFormatVersion,matNumber,FILEUNIT) - allocate (mesh_mapFEtoCPnode(2,Nnodes), source = 0) + allocate (mesh_mapFEtoCPnode(2,Nnodes), source=0) call inputRead_mapNodes(inputFile) call inputRead_elemType(elem, & nElems,FILEUNIT) call inputRead_elemNodes(node0_elem, & Nnodes,inputFile) - + connectivity_elem = inputRead_connectivityElem(nElems,elem%nNodes,FILEUNIT) call inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, & @@ -423,22 +423,30 @@ subroutine inputRead_mapElems(tableStyle,nameElemSet,mapElemSet,fileFormatVersio character(len=300) :: line, & tmp - integer, dimension (1+size(mesh_mapFEtoCPelem,2)) :: contInts + integer, dimension(:), allocatable :: contInts integer :: i,cpElem + allocate(contInts(size(mesh_mapFEtoCPelem,2)+1)) + write(6,*) 'hallo';flush(6) + write(6,*) nameElemSet;flush(6) + write(6,*) mapElemSet;flush(6) + write(6,*) 'map', size(mesh_mapFEtoCPelem,2);flush(6) cpElem = 0 contInts = 0 rewind(fileUnit) do read (fileUnit,'(A300)',END=620) line + write(6,*) trim(line);flush(6) chunkPos = IO_stringPos(line) if (fileFormatVersion < 13) then ! Marc 2016 or earlier if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'hypoelastic' ) then do i=1,3+TableStyle ! skip three (or four if new table style!) lines read (fileUnit,'(A300)') line + write(6,*) i;flush(6) enddo contInts = IO_continuousIntValues(fileUnit,size(mesh_mapFEtoCPelem,2),nameElemSet,& - mapElemSet,size(nameElemSet)) + mapElemSet,size(nameElemSet)) + exit endif else ! Marc2017 and later @@ -739,7 +747,6 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza 630 end subroutine inputRead_microstructureAndHomogenization - !-------------------------------------------------------------------------------------------------- !> @brief Calculates cell node coordinates from element node coordinates !-------------------------------------------------------------------------------------------------- @@ -926,6 +933,7 @@ subroutine buildCellNodes(node_cell, & end subroutine buildCellNodes + !-------------------------------------------------------------------------------------------------- !> @brief Calculates IP coordinates as center of cell !-------------------------------------------------------------------------------------------------- @@ -950,6 +958,135 @@ subroutine buildIPcoordinates(IPcoordinates, & end subroutine buildIPcoordinates +!--------------------------------------------------------------------------------------------------- +!> @brief Calculates IP volume. +!> @details The IP volume is calculated differently depending on the cell type. +!> 2D cells assume an element depth of one in order to calculate the volume. +!--------------------------------------------------------------------------------------------------- +function IPvolume(elem,node,connectivity) + + type(tElement) :: elem + real(pReal), dimension(:,:), intent(in) :: node + integer, dimension(:,:,:), intent(in) :: connectivity + real(pReal), dimension(elem%nIPs,size(connectivity,3)) :: IPvolume + + integer :: e,i + + do e = 1,size(connectivity,3) + do i = 1,elem%nIPs + + select case (elem%cellType) + case (1) ! 2D 3node + IPvolume(i,e) = math_areaTriangle(node(1:3,connectivity(1,i,e)), & + node(1:3,connectivity(2,i,e)), & + node(1:3,connectivity(3,i,e))) + + case (2) ! 2D 4node + IPvolume(i,e) = math_areaTriangle(node(1:3,connectivity(1,i,e)), & ! assume planar shape, division in two triangles suffices + node(1:3,connectivity(2,i,e)), & + node(1:3,connectivity(3,i,e))) & + + math_areaTriangle(node(1:3,connectivity(3,i,e)), & + node(1:3,connectivity(4,i,e)), & + node(1:3,connectivity(1,i,e))) + case (3) ! 3D 4node + IPvolume(i,e) = math_volTetrahedron(node(1:3,connectivity(1,i,e)), & + node(1:3,connectivity(2,i,e)), & + node(1:3,connectivity(3,i,e)), & + node(1:3,connectivity(4,i,e))) + case (4) ! 3D 8node + IPvolume(i,e) = math_volTetrahedron(node(1:3,connectivity(1,i,e)), & + node(1:3,connectivity(2,i,e)), & + node(1:3,connectivity(3,i,e)), & + node(1:3,connectivity(6,i,e))) & + + math_volTetrahedron(node(1:3,connectivity(3,i,e)), & + node(1:3,connectivity(4,i,e)), & + node(1:3,connectivity(1,i,e)), & + node(1:3,connectivity(8,i,e))) & + + math_volTetrahedron(node(1:3,connectivity(6,i,e)), & + node(1:3,connectivity(7,i,e)), & + node(1:3,connectivity(8,i,e)), & + node(1:3,connectivity(3,i,e))) & + + math_volTetrahedron(node(1:3,connectivity(8,i,e)), & + node(1:3,connectivity(5,i,e)), & + node(1:3,connectivity(6,i,e)), & + node(1:3,connectivity(1,i,e))) ! first triangulation + IPvolume(i,e) = IPvolume(i,e) & + + math_volTetrahedron(node(1:3,connectivity(2,i,e)), & + node(1:3,connectivity(3,i,e)), & + node(1:3,connectivity(4,i,e)), & + node(1:3,connectivity(7,i,e))) & + + math_volTetrahedron(node(1:3,connectivity(4,i,e)), & + node(1:3,connectivity(1,i,e)), & + node(1:3,connectivity(2,i,e)), & + node(1:3,connectivity(5,i,e))) & + + math_volTetrahedron(node(1:3,connectivity(7,i,e)), & + node(1:3,connectivity(8,i,e)), & + node(1:3,connectivity(5,i,e)), & + node(1:3,connectivity(4,i,e))) & + + math_volTetrahedron(node(1:3,connectivity(5,i,e)), & + node(1:3,connectivity(6,i,e)), & + node(1:3,connectivity(7,i,e)), & + node(1:3,connectivity(2,i,e))) ! second triangulation + IPvolume(i,e) = IPvolume(i,e) * 0.5_pReal + end select + enddo + enddo + +end function IPvolume + +!-------------------------------------------------------------------------------------------------- +!> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal' +!-------------------------------------------------------------------------------------------------- +subroutine mesh_build_ipAreas(ipArea,ipAreaNormal,elem,nElem,connectivity,node) + + type(tElement) :: elem + integer :: nElem + integer, dimension(:,:,:), intent(in) :: connectivity + real(pReal), dimension(:,:), intent(in) :: node + real(pReal), dimension(elem%nIPneighbors,elem%nIPs,nElem), intent(out) :: ipArea + real(pReal), dimension(3,elem%nIPneighbors,elem%nIPs,nElem), intent(out) :: ipAreaNormal + + real(pReal), dimension (3,size(elem%cellFace,2)) :: nodePos, normals + real(pReal), dimension(3) :: normal + integer :: e,i,f,n,m + + m = size(elem%cellFace,1) + + do e = 1,nElem + do i = 1,elem%nIPs + do f = 1,size(elem%cellFace,2) !???? + + forall(n = 1: size(elem%cellface,1)) & + nodePos(1:3,n) = node(1:3,connectivity(elem%cellface(n,f),i,e)) + + select case (elem%cellType) + case (1,2) ! 2D 3 or 4 node + normal(1) = nodePos(2,2) - nodePos(2,1) ! x_normal = y_connectingVector + normal(2) = -(nodePos(1,2) - nodePos(1,1)) ! y_normal = -x_connectingVector + normal(3) = 0.0_pReal + case (3) ! 3D 4node + normal = math_cross(nodePos(1:3,2) - nodePos(1:3,1), & + nodePos(1:3,3) - nodePos(1:3,1)) + case (4) ! 3D 8node + ! for this cell type we get the normal of the quadrilateral face as an average of + ! four normals of triangular subfaces; since the face consists only of two triangles, + ! the sum has to be divided by two; this whole prcedure tries to compensate for + ! probable non-planar cell surfaces + normals(1:3,n) = 0.5_pReal & + * math_cross(nodePos(1:3,1+mod(n ,m)) - nodePos(1:3,n), & + nodePos(1:3,1+mod(n+1,m)) - nodePos(1:3,n)) + normal = 0.5_pReal * sum(normals,2) + end select + + ipArea(f,i,e) = norm2(normal) + ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal + + enddo + enddo + enddo + +end subroutine mesh_build_ipAreas + !-------------------------------------------------------------------------------------------------- !> @brief Gives the FE to CP ID mapping by binary search through lookup array !! valid questions (what) are 'elem', 'node' From 0ae0e2332583976c29ca858f63f2889e6bf3b2ab Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 16 Oct 2019 20:56:31 +0200 Subject: [PATCH 10/20] volume calculation was wrong could result in negative volumes --- src/math.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 26e63388c..453145700 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1307,10 +1307,10 @@ real(pReal) pure function math_volTetrahedron(v1,v2,v3,v4) real(pReal), dimension (3,3) :: m m(1:3,1) = v1-v2 - m(1:3,2) = v2-v3 - m(1:3,3) = v3-v4 + m(1:3,2) = v1-v3 + m(1:3,3) = v1-v4 - math_volTetrahedron = math_det33(m)/6.0_pReal + math_volTetrahedron = abs(math_det33(m))/6.0_pReal end function math_volTetrahedron From 95ecc05cb047432cb0a472e0ab223fd53d63d91e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 16 Oct 2019 21:13:26 +0200 Subject: [PATCH 11/20] better test ... --- src/math.f90 | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/math.f90 b/src/math.f90 index 453145700..5a3a5f467 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1404,6 +1404,7 @@ subroutine unitTest integer, dimension(5) :: range_out_ = [1,2,3,4,5] real(pReal) :: det + real(pReal), dimension(3) :: v3_1,v3_2,v3_3,v3_4 real(pReal), dimension(6) :: v6 real(pReal), dimension(9) :: v9 real(pReal), dimension(3,3) :: t33,t33_2 @@ -1452,6 +1453,15 @@ subroutine unitTest if(any(dNeq0(math_6toSym33(v6) - math_symmetric33(math_6toSym33(v6))))) & call IO_error(401,ext_msg='math_symmetric33') + call random_number(v3_1) + call random_number(v3_2) + call random_number(v3_3) + call random_number(v3_4) + + if(dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0, & + math_volTetrahedron(v3_1,v3_2,v3_3,v3_4),tol=1.0e-12_pReal)) & + call IO_error(401,ext_msg='math_volTetrahedron') + call random_number(t33) if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) & call IO_error(401,ext_msg='math_det33/math_detSym33') From 9b5545229f1a1deb77e432657d77d9e8086666d1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 16 Oct 2019 21:13:43 +0200 Subject: [PATCH 12/20] leftover comment ... --- src/element.f90 | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/element.f90 b/src/element.f90 index 3252d2813..8fb7d4c8b 100644 --- a/src/element.f90 +++ b/src/element.f90 @@ -256,10 +256,6 @@ module element -3,26,-4,24,-6,18 & ],[nIPneighbor(cellType(10)),nIP(10)]) -! MD: probably not needed END -! -------------------------------------------------------------------------------------------------- - - integer, dimension(nNode(1),NcellNode(geomType(1))), parameter :: cellNodeParentNodeWeights1 = & reshape([& From 3f481e1ceaa8a5fe8113b68054c8c3e1f7be6913 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 16 Oct 2019 22:00:25 +0200 Subject: [PATCH 13/20] corrected volume calculation and write to DADF5. follows https://www.osti.gov/servlets/purl/632793/ --- src/commercialFEM_fileList.f90 | 2 +- src/geometry_plastic_nonlocal.f90 | 19 ++++ src/mesh_marc.f90 | 140 +++++++++++++----------------- 3 files changed, 80 insertions(+), 81 deletions(-) diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index ac4aa009a..79a385818 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -14,11 +14,11 @@ #include "Lambert.f90" #include "rotations.f90" #include "FEsolving.f90" -#include "geometry_plastic_nonlocal.f90" #include "element.f90" #include "mesh_base.f90" #include "HDF5_utilities.f90" #include "results.f90" +#include "geometry_plastic_nonlocal.f90" #include "discretization.f90" #ifdef Abaqus #include "mesh_abaqus.f90" diff --git a/src/geometry_plastic_nonlocal.f90 b/src/geometry_plastic_nonlocal.f90 index 37909a4a5..07bb9a03b 100644 --- a/src/geometry_plastic_nonlocal.f90 +++ b/src/geometry_plastic_nonlocal.f90 @@ -7,6 +7,7 @@ !-------------------------------------------------------------------------------------------------- module geometry_plastic_nonlocal use prec + use results implicit none private @@ -32,6 +33,7 @@ module geometry_plastic_nonlocal geometry_plastic_nonlocal_setIPvolume, & geometry_plastic_nonlocal_setIParea, & geometry_plastic_nonlocal_setIPareaNormal, & + geometry_plastic_nonlocal_results, & geometry_plastic_nonlocal_disable contains @@ -112,4 +114,21 @@ subroutine geometry_plastic_nonlocal_disable end subroutine geometry_plastic_nonlocal_disable + +!--------------------------------------------------------------------------------------------------- +!> @brief Frees memory used by variables only needed by plastic_nonlocal +!--------------------------------------------------------------------------------------------------- +subroutine geometry_plastic_nonlocal_results + +#if defined(DAMASK_HDF5) + call results_openJobFile + call results_writeDataset('geometry',geometry_plastic_nonlocal_IPvolume0,'v_0',& + 'initial cell volume','m³') + call results_writeDataset('geometry',geometry_plastic_nonlocal_IParea0,'a_0',& + 'initial cell face area','m²') + call results_closeJobFile +#endif + +end subroutine geometry_plastic_nonlocal_results + end module geometry_plastic_nonlocal diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 517865718..f1e044d69 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -77,7 +77,8 @@ subroutine mesh_init(ip,el) connectivity_cell !< cell connectivity for each element,ip/cell integer, dimension(:,:), allocatable :: & connectivity_elem - + + real(pReal), dimension(:,:,:,:),allocatable :: x write(6,'(/,a)') ' <<<+- mesh init -+>>>' @@ -116,6 +117,11 @@ subroutine mesh_init(ip,el) reshape(connectivity_cell,[elem%NcellNodesPerCell,elem%nIPs*nElems]),& node0_cell,ip_reshaped) + call geometry_plastic_nonlocal_setIPvolume(IPvolume(elem,node0_cell,connectivity_cell)) + x = IPareaNormal(elem,nElems,connectivity_cell,node0_cell) + call geometry_plastic_nonlocal_setIParea(norm2(x,1)) + call geometry_plastic_nonlocal_results + end subroutine mesh_init !-------------------------------------------------------------------------------------------------- @@ -158,7 +164,7 @@ subroutine writeGeometry(elemType, & call results_writeDataset('geometry',coordinates_temp,'x_p', & 'coordinates of the material points','m') - call results_closeJobFile() + call results_closeJobFile #endif end subroutine writeGeometry @@ -427,22 +433,17 @@ subroutine inputRead_mapElems(tableStyle,nameElemSet,mapElemSet,fileFormatVersio integer :: i,cpElem allocate(contInts(size(mesh_mapFEtoCPelem,2)+1)) - write(6,*) 'hallo';flush(6) - write(6,*) nameElemSet;flush(6) - write(6,*) mapElemSet;flush(6) - write(6,*) 'map', size(mesh_mapFEtoCPelem,2);flush(6) + cpElem = 0 contInts = 0 rewind(fileUnit) do read (fileUnit,'(A300)',END=620) line - write(6,*) trim(line);flush(6) chunkPos = IO_stringPos(line) if (fileFormatVersion < 13) then ! Marc 2016 or earlier if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'hypoelastic' ) then do i=1,3+TableStyle ! skip three (or four if new table style!) lines read (fileUnit,'(A300)') line - write(6,*) i;flush(6) enddo contInts = IO_continuousIntValues(fileUnit,size(mesh_mapFEtoCPelem,2),nameElemSet,& mapElemSet,size(nameElemSet)) @@ -961,14 +962,16 @@ end subroutine buildIPcoordinates !--------------------------------------------------------------------------------------------------- !> @brief Calculates IP volume. !> @details The IP volume is calculated differently depending on the cell type. -!> 2D cells assume an element depth of one in order to calculate the volume. +!> 2D cells assume an element depth of 1.0 !--------------------------------------------------------------------------------------------------- function IPvolume(elem,node,connectivity) - type(tElement) :: elem - real(pReal), dimension(:,:), intent(in) :: node - integer, dimension(:,:,:), intent(in) :: connectivity + type(tElement), intent(in) :: elem + real(pReal), dimension(:,:), intent(in) :: node + integer, dimension(:,:,:), intent(in) :: connectivity + real(pReal), dimension(elem%nIPs,size(connectivity,3)) :: IPvolume + real(pReal), dimension(3) :: x0,x1,x2,x3,x4,x5,x6,x7 integer :: e,i @@ -994,98 +997,75 @@ function IPvolume(elem,node,connectivity) node(1:3,connectivity(3,i,e)), & node(1:3,connectivity(4,i,e))) case (4) ! 3D 8node - IPvolume(i,e) = math_volTetrahedron(node(1:3,connectivity(1,i,e)), & - node(1:3,connectivity(2,i,e)), & - node(1:3,connectivity(3,i,e)), & - node(1:3,connectivity(6,i,e))) & - + math_volTetrahedron(node(1:3,connectivity(3,i,e)), & - node(1:3,connectivity(4,i,e)), & - node(1:3,connectivity(1,i,e)), & - node(1:3,connectivity(8,i,e))) & - + math_volTetrahedron(node(1:3,connectivity(6,i,e)), & - node(1:3,connectivity(7,i,e)), & - node(1:3,connectivity(8,i,e)), & - node(1:3,connectivity(3,i,e))) & - + math_volTetrahedron(node(1:3,connectivity(8,i,e)), & - node(1:3,connectivity(5,i,e)), & - node(1:3,connectivity(6,i,e)), & - node(1:3,connectivity(1,i,e))) ! first triangulation - IPvolume(i,e) = IPvolume(i,e) & - + math_volTetrahedron(node(1:3,connectivity(2,i,e)), & - node(1:3,connectivity(3,i,e)), & - node(1:3,connectivity(4,i,e)), & - node(1:3,connectivity(7,i,e))) & - + math_volTetrahedron(node(1:3,connectivity(4,i,e)), & - node(1:3,connectivity(1,i,e)), & - node(1:3,connectivity(2,i,e)), & - node(1:3,connectivity(5,i,e))) & - + math_volTetrahedron(node(1:3,connectivity(7,i,e)), & - node(1:3,connectivity(8,i,e)), & - node(1:3,connectivity(5,i,e)), & - node(1:3,connectivity(4,i,e))) & - + math_volTetrahedron(node(1:3,connectivity(5,i,e)), & - node(1:3,connectivity(6,i,e)), & - node(1:3,connectivity(7,i,e)), & - node(1:3,connectivity(2,i,e))) ! second triangulation - IPvolume(i,e) = IPvolume(i,e) * 0.5_pReal + x0 = node(1:3,connectivity(1,i,e)) + x1 = node(1:3,connectivity(2,i,e)) + x2 = node(1:3,connectivity(4,i,e)) + x3 = node(1:3,connectivity(3,i,e)) + x4 = node(1:3,connectivity(5,i,e)) + x5 = node(1:3,connectivity(6,i,e)) + x6 = node(1:3,connectivity(8,i,e)) + x7 = node(1:3,connectivity(7,i,e)) + IPvolume(i,e) = dot_product((x7-x1)+(x6-x0),math_cross((x7-x2), (x3-x0))) & + + dot_product((x6-x0), math_cross((x7-x2)+(x5-x0),(x7-x4))) & + + dot_product((x7-x1), math_cross((x5-x0), (x7-x4)+(x3-x0))) + IPvolume(i,e) = IPvolume(i,e)/12.0_pReal end select enddo enddo end function IPvolume -!-------------------------------------------------------------------------------------------------- -!> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipAreas(ipArea,ipAreaNormal,elem,nElem,connectivity,node) - type(tElement) :: elem - integer :: nElem - integer, dimension(:,:,:), intent(in) :: connectivity - real(pReal), dimension(:,:), intent(in) :: node - real(pReal), dimension(elem%nIPneighbors,elem%nIPs,nElem), intent(out) :: ipArea - real(pReal), dimension(3,elem%nIPneighbors,elem%nIPs,nElem), intent(out) :: ipAreaNormal +!-------------------------------------------------------------------------------------------------- +!> @brief calculation of IP interface areas +!-------------------------------------------------------------------------------------------------- +function IPareaNormal(elem,nElem,connectivity,node) - real(pReal), dimension (3,size(elem%cellFace,2)) :: nodePos, normals - real(pReal), dimension(3) :: normal - integer :: e,i,f,n,m + type(tElement), intent(in) :: elem + integer, intent(in) :: nElem + integer, dimension(:,:,:), intent(in) :: connectivity + real(pReal), dimension(:,:), intent(in) :: node + + real(pReal), dimension(3,elem%nIPneighbors,elem%nIPs,nElem) :: ipAreaNormal + + real(pReal), dimension (3,size(elem%cellFace,2)) :: nodePos + integer :: e,i,f,n,m - m = size(elem%cellFace,1) + m = size(elem%cellFace,1) do e = 1,nElem do i = 1,elem%nIPs - do f = 1,size(elem%cellFace,2) !???? - - forall(n = 1: size(elem%cellface,1)) & + do f = 1,size(elem%cellFace,2) + do n = 1, m nodePos(1:3,n) = node(1:3,connectivity(elem%cellface(n,f),i,e)) + write(6,*) e,i,f,n,nodePos(1:3,n) + enddo select case (elem%cellType) - case (1,2) ! 2D 3 or 4 node - normal(1) = nodePos(2,2) - nodePos(2,1) ! x_normal = y_connectingVector - normal(2) = -(nodePos(1,2) - nodePos(1,1)) ! y_normal = -x_connectingVector - normal(3) = 0.0_pReal - case (3) ! 3D 4node - normal = math_cross(nodePos(1:3,2) - nodePos(1:3,1), & - nodePos(1:3,3) - nodePos(1:3,1)) - case (4) ! 3D 8node + case (1,2) ! 2D 3 or 4 node + IPareaNormal(1,n,i,e) = nodePos(2,2) - nodePos(2,1) ! x_normal = y_connectingVector + IPareaNormal(2,n,i,e) = -(nodePos(1,2) - nodePos(1,1)) ! y_normal = -x_connectingVector + IPareaNormal(3,n,i,e) = 0.0_pReal + case (3) ! 3D 4node + IPareaNormal(1:3,n,i,e) = math_cross(nodePos(1:3,2) - nodePos(1:3,1), & + nodePos(1:3,3) - nodePos(1:3,1)) + case (4) ! 3D 8node ! for this cell type we get the normal of the quadrilateral face as an average of ! four normals of triangular subfaces; since the face consists only of two triangles, ! the sum has to be divided by two; this whole prcedure tries to compensate for ! probable non-planar cell surfaces - normals(1:3,n) = 0.5_pReal & - * math_cross(nodePos(1:3,1+mod(n ,m)) - nodePos(1:3,n), & - nodePos(1:3,1+mod(n+1,m)) - nodePos(1:3,n)) - normal = 0.5_pReal * sum(normals,2) + do n = 1, m + IPareaNormal(1:3,n,i,e) = 0.5_pReal & + * math_cross(nodePos(1:3,1+mod(n ,m)) - nodePos(1:3,n), & + nodePos(1:3,1+mod(n+1,m)) - nodePos(1:3,n)) + enddo end select - - ipArea(f,i,e) = norm2(normal) - ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal - enddo enddo enddo -end subroutine mesh_build_ipAreas +end function IPareaNormal + !-------------------------------------------------------------------------------------------------- !> @brief Gives the FE to CP ID mapping by binary search through lookup array From f33a99d125408a423cf86c4844d6b454310d523e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 17 Oct 2019 00:21:48 +0200 Subject: [PATCH 14/20] polishing --- src/geometry_plastic_nonlocal.f90 | 23 +++++++-- src/mesh_marc.f90 | 85 ++++++++++++++++--------------- 2 files changed, 63 insertions(+), 45 deletions(-) diff --git a/src/geometry_plastic_nonlocal.f90 b/src/geometry_plastic_nonlocal.f90 index 07bb9a03b..b042943d0 100644 --- a/src/geometry_plastic_nonlocal.f90 +++ b/src/geometry_plastic_nonlocal.f90 @@ -119,13 +119,28 @@ end subroutine geometry_plastic_nonlocal_disable !> @brief Frees memory used by variables only needed by plastic_nonlocal !--------------------------------------------------------------------------------------------------- subroutine geometry_plastic_nonlocal_results + + integer, dimension(:), allocatable :: s #if defined(DAMASK_HDF5) call results_openJobFile - call results_writeDataset('geometry',geometry_plastic_nonlocal_IPvolume0,'v_0',& - 'initial cell volume','m³') - call results_writeDataset('geometry',geometry_plastic_nonlocal_IParea0,'a_0',& - 'initial cell face area','m²') + + writeVolume: block + real(pReal), dimension(:), allocatable :: temp + s = shape(geometry_plastic_nonlocal_IPvolume0) + temp = reshape(geometry_plastic_nonlocal_IPvolume0,[s(1)*s(2)]) + call results_writeDataset('geometry',temp,'v_0',& + 'initial cell volume','m³') + end block writeVolume + + writeArea: block + real(pReal), dimension(:,:), allocatable :: temp + s = shape(geometry_plastic_nonlocal_IParea0) + temp = reshape(geometry_plastic_nonlocal_IParea0,[s(1),s(2)*s(3)]) + call results_writeDataset('geometry',temp,'a_0',& + 'initial cell face area','m²') + end block writeArea + call results_closeJobFile #endif diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index f1e044d69..eae1e8ac7 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -83,9 +83,21 @@ subroutine mesh_init(ip,el) write(6,'(/,a)') ' <<<+- mesh init -+>>>' mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh + call inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogenizationAt) nElems = size(connectivity_elem,2) + if (debug_e < 1 .or. debug_e > nElems) call IO_error(602,ext_msg='element') + if (debug_i < 1 .or. debug_i > elem%nIPs) call IO_error(602,ext_msg='IP') + + FEsolving_execElem = [ 1,nElems ] ! parallel loop bounds set to comprise all DAMASK elements + allocate(FEsolving_execIP(2,nElems), source=1) ! parallel loop bounds set to comprise from first IP... + FEsolving_execIP(2,:) = elem%nIPs + + allocate(calcMode(elem%nIPs,nElems),source=.false.) ! pretend to have collected what first call is asking (F = I) + calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc" + + allocate(mesh_ipCoordinates(3,elem%nIPs,nElems),source=0.0_pReal) ! deprecated allocate(cellNodeDefinition(elem%nNodes-1)) @@ -97,25 +109,15 @@ subroutine mesh_init(ip,el) cellNodeDefinition,node0_elem) allocate(ip_reshaped(3,elem%nIPs*nElems),source=0.0_pReal) call buildIPcoordinates(ip_reshaped,reshape(connectivity_cell,[elem%NcellNodesPerCell,& - elem%nIPs*nElems]),node0_cell) + elem%nIPs*nElems]),node0_cell) - if (debug_e < 1 .or. debug_e > nElems) call IO_error(602,ext_msg='element') - if (debug_i < 1 .or. debug_i > elem%nIPs) call IO_error(602,ext_msg='IP') - - FEsolving_execElem = [ 1,nElems ] ! parallel loop bounds set to comprise all DAMASK elements - allocate(FEsolving_execIP(2,nElems), source=1) ! parallel loop bounds set to comprise from first IP... - FEsolving_execIP(2,:) = elem%nIPs - - allocate(calcMode(elem%nIPs,nElems),source=.false.) ! pretend to have collected what first call is asking (F = I) - calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc" - call discretization_init(microstructureAt,homogenizationAt,& ip_reshaped,& node0_elem) call writeGeometry(0,connectivity_elem,& - reshape(connectivity_cell,[elem%NcellNodesPerCell,elem%nIPs*nElems]),& - node0_cell,ip_reshaped) + reshape(connectivity_cell,[elem%NcellNodesPerCell,elem%nIPs*nElems]),& + node0_cell,ip_reshaped) call geometry_plastic_nonlocal_setIPvolume(IPvolume(elem,node0_cell,connectivity_cell)) x = IPareaNormal(elem,nElems,connectivity_cell,node0_cell) @@ -997,6 +999,9 @@ function IPvolume(elem,node,connectivity) node(1:3,connectivity(3,i,e)), & node(1:3,connectivity(4,i,e))) case (4) ! 3D 8node + ! J. Grandy, Efficient Calculation of Volume of Hexahedral Cells + ! Lawrence Livermore National Laboratory + ! https://www.osti.gov/servlets/purl/632793 x0 = node(1:3,connectivity(1,i,e)) x1 = node(1:3,connectivity(2,i,e)) x2 = node(1:3,connectivity(4,i,e)) @@ -1028,42 +1033,40 @@ function IPareaNormal(elem,nElem,connectivity,node) real(pReal), dimension(3,elem%nIPneighbors,elem%nIPs,nElem) :: ipAreaNormal - real(pReal), dimension (3,size(elem%cellFace,2)) :: nodePos + real(pReal), dimension (3,size(elem%cellFace,1)) :: nodePos integer :: e,i,f,n,m m = size(elem%cellFace,1) - do e = 1,nElem - do i = 1,elem%nIPs - do f = 1,size(elem%cellFace,2) - do n = 1, m - nodePos(1:3,n) = node(1:3,connectivity(elem%cellface(n,f),i,e)) - write(6,*) e,i,f,n,nodePos(1:3,n) - enddo + do e = 1,nElem + do i = 1,elem%nIPs + do f = 1,size(elem%cellFace,2) + nodePos = node(1:3,connectivity(elem%cellface(1:m,f),i,e)) - select case (elem%cellType) - case (1,2) ! 2D 3 or 4 node - IPareaNormal(1,n,i,e) = nodePos(2,2) - nodePos(2,1) ! x_normal = y_connectingVector - IPareaNormal(2,n,i,e) = -(nodePos(1,2) - nodePos(1,1)) ! y_normal = -x_connectingVector - IPareaNormal(3,n,i,e) = 0.0_pReal - case (3) ! 3D 4node - IPareaNormal(1:3,n,i,e) = math_cross(nodePos(1:3,2) - nodePos(1:3,1), & - nodePos(1:3,3) - nodePos(1:3,1)) - case (4) ! 3D 8node - ! for this cell type we get the normal of the quadrilateral face as an average of - ! four normals of triangular subfaces; since the face consists only of two triangles, - ! the sum has to be divided by two; this whole prcedure tries to compensate for - ! probable non-planar cell surfaces - do n = 1, m - IPareaNormal(1:3,n,i,e) = 0.5_pReal & - * math_cross(nodePos(1:3,1+mod(n ,m)) - nodePos(1:3,n), & - nodePos(1:3,1+mod(n+1,m)) - nodePos(1:3,n)) - enddo + select case (elem%cellType) + case (1,2) ! 2D 3 or 4 node + IPareaNormal(1,f,i,e) = nodePos(2,2) - nodePos(2,1) ! x_normal = y_connectingVector + IPareaNormal(2,f,i,e) = -(nodePos(1,2) - nodePos(1,1)) ! y_normal = -x_connectingVector + IPareaNormal(3,f,i,e) = 0.0_pReal + case (3) ! 3D 4node + IPareaNormal(1:3,f,i,e) = math_cross(nodePos(1:3,2) - nodePos(1:3,1), & + nodePos(1:3,3) - nodePos(1:3,1)) + case (4) ! 3D 8node + ! for this cell type we get the normal of the quadrilateral face as an average of + ! four normals of triangular subfaces; since the face consists only of two triangles, + ! the sum has to be divided by two; this whole prcedure tries to compensate for + ! probable non-planar cell surfaces + IPareaNormal(1:3,f,i,e) = 0.0_pReal + do n = 1, m + IPareaNormal(1:3,f,i,e) = IPareaNormal(1:3,f,i,e) & + + math_cross(nodePos(1:3,mod(n+0,m)+1) - nodePos(1:3,n), & + nodePos(1:3,mod(n+1,m)+1) - nodePos(1:3,n)) * 0.5_pReal + enddo end select enddo enddo - enddo - + enddo + end function IPareaNormal From 2b68c108f08f596f1b34ec86f4600c59430b98e2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 17 Oct 2019 05:54:08 +0200 Subject: [PATCH 15/20] read from memory, not from file --- src/IO.f90 | 46 ++++---------- src/mesh_marc.f90 | 155 +++++++++++++++++++++++----------------------- 2 files changed, 89 insertions(+), 112 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index 64ce6fc9a..5d96e67d2 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -44,7 +44,6 @@ module IO IO_extractValue, & IO_countDataLines #elif defined(Marc4DAMASK) - IO_skipChunks, & IO_fixedNoEFloatValue, & IO_fixedIntValue, & IO_countNumericalDataLines @@ -189,17 +188,17 @@ integer function IO_open_binary(fileName,mode) m = 'r' endif - if (m == 'w') then - open(newunit=IO_open_binary, file=trim(fileName),& - status='replace',access='stream',action='write',iostat=ierr) - if (ierr /= 0) call IO_error(100,ext_msg='could not open file (w): '//trim(fileName)) - elseif(m == 'r') then - open(newunit=IO_open_binary, file=trim(fileName),& - status='old', access='stream',action='read', iostat=ierr) - if (ierr /= 0) call IO_error(100,ext_msg='could not open file (r): '//trim(fileName)) - else - call IO_error(100,ext_msg='unknown access mode: '//m) - endif + if (m == 'w') then + open(newunit=IO_open_binary, file=trim(fileName),& + status='replace',access='stream',action='write',iostat=ierr) + if (ierr /= 0) call IO_error(100,ext_msg='could not open file (w): '//trim(fileName)) + elseif(m == 'r') then + open(newunit=IO_open_binary, file=trim(fileName),& + status='old', access='stream',action='read', iostat=ierr) + if (ierr /= 0) call IO_error(100,ext_msg='could not open file (r): '//trim(fileName)) + else + call IO_error(100,ext_msg='unknown access mode: '//m) + endif end function IO_open_binary @@ -403,7 +402,7 @@ pure function IO_stringPos(string) left = right + verify(string(right+1:),SEP) right = left + scan(string(left:),SEP) - 2 if ( string(left:left) == '#' ) exit - IO_stringPos = [IO_stringPos,left, right] + IO_stringPos = [IO_stringPos,left,right] IO_stringPos(1) = IO_stringPos(1)+1 endOfString: if (right < left) then IO_stringPos(IO_stringPos(1)*2+1) = len_trim(string) @@ -1023,27 +1022,6 @@ integer function IO_countNumericalDataLines(fileUnit) backspace(fileUnit) end function IO_countNumericalDataLines - - -!-------------------------------------------------------------------------------------------------- -!> @brief reads file to skip (at least) N chunks (may be over multiple lines) -!-------------------------------------------------------------------------------------------------- -subroutine IO_skipChunks(fileUnit,N) - - integer, intent(in) :: fileUnit, & !< file handle - N !< minimum number of chunks to skip - - integer :: remainingChunks - character(len=65536) :: line - - line = '' - remainingChunks = N - - do while (trim(line) /= IO_EOF .and. remainingChunks > 0) - line = IO_read(fileUnit) - remainingChunks = remainingChunks - (size(IO_stringPos(line))-1)/2 - enddo -end subroutine IO_skipChunks #endif diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index eae1e8ac7..82719f0d0 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -123,6 +123,7 @@ subroutine mesh_init(ip,el) x = IPareaNormal(elem,nElems,connectivity_cell,node0_cell) call geometry_plastic_nonlocal_setIParea(norm2(x,1)) call geometry_plastic_nonlocal_results + end subroutine mesh_init @@ -224,7 +225,7 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni call inputRead_mapNodes(inputFile) call inputRead_elemType(elem, & - nElems,FILEUNIT) + nElems,inputFile) call inputRead_elemNodes(node0_elem, & Nnodes,inputFile) @@ -355,7 +356,7 @@ end subroutine inputRead_NnodesAndElements !> @brief Count overall number of element sets in mesh. !-------------------------------------------------------------------------------------------------- subroutine inputRead_NelemSets(nElemSets,maxNelemInSet,& - fileUnit) + fileUnit) integer, intent(out) :: nElemSets, maxNelemInSet integer, intent(in) :: fileUnit @@ -421,62 +422,61 @@ subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,fileUnit) !-------------------------------------------------------------------------------------------------- subroutine inputRead_mapElems(tableStyle,nameElemSet,mapElemSet,fileFormatVersion,matNumber,fileUnit) - integer, intent(in) :: fileUnit,tableStyle,fileFormatVersion - integer, dimension(:), intent(in) :: matNumber - character(len=64), intent(in), dimension(:) :: nameElemSet - integer, dimension(:,:), intent(in) :: & - mapElemSet - - integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line, & - tmp - - integer, dimension(:), allocatable :: contInts - integer :: i,cpElem - - allocate(contInts(size(mesh_mapFEtoCPelem,2)+1)) - - cpElem = 0 - contInts = 0 - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - if (fileFormatVersion < 13) then ! Marc 2016 or earlier - if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'hypoelastic' ) then - do i=1,3+TableStyle ! skip three (or four if new table style!) lines - read (fileUnit,'(A300)') line - enddo - contInts = IO_continuousIntValues(fileUnit,size(mesh_mapFEtoCPelem,2),nameElemSet,& - mapElemSet,size(nameElemSet)) - - exit - endif - else ! Marc2017 and later - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'connectivity') then - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - if(any(matNumber==IO_intValue(line,chunkPos,6))) then - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - tmp = IO_lc(IO_stringValue(line,chunkPos,1)) - if (verify(trim(tmp),"0123456789")/=0) then ! found keyword - exit - else - contInts(1) = contInts(1) + 1 - read (tmp,*) contInts(contInts(1)+1) - endif - enddo - endif - endif - endif - enddo + integer, intent(in) :: fileUnit,tableStyle,fileFormatVersion + integer, dimension(:), intent(in) :: matNumber + character(len=64), intent(in), dimension(:) :: nameElemSet + integer, dimension(:,:), intent(in) :: & + mapElemSet + + integer, allocatable, dimension(:) :: chunkPos + character(len=300) :: line, & + tmp + + integer, dimension(:), allocatable :: contInts + integer :: i,cpElem + + allocate(contInts(size(mesh_mapFEtoCPelem,2)+1)) + + cpElem = 0 + contInts = 0 + rewind(fileUnit) + do + read (fileUnit,'(A300)',END=620) line + chunkPos = IO_stringPos(line) + Marc2016andEarlier: if (fileFormatVersion < 13) then + if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'hypoelastic' ) then + skipLines: do i=1,3+TableStyle + read (fileUnit,'(A300)') line + enddo skipLines + contInts = IO_continuousIntValues(fileUnit,size(mesh_mapFEtoCPelem,2),nameElemSet,& + mapElemSet,size(nameElemSet)) + exit + endif + else Marc2016andEarlier + if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'connectivity') then + read (fileUnit,'(A300)',END=620) line + chunkPos = IO_stringPos(line) + if(any(matNumber==IO_intValue(line,chunkPos,6))) then + do + read (fileUnit,'(A300)',END=620) line + chunkPos = IO_stringPos(line) + tmp = IO_lc(IO_stringValue(line,chunkPos,1)) + if (verify(trim(tmp),"0123456789")/=0) then ! found keyword + exit + else + contInts(1) = contInts(1) + 1 + read (tmp,*) contInts(contInts(1)+1) + endif + enddo + endif + endif + endif Marc2016andEarlier + enddo 620 do i = 1,contInts(1) - cpElem = cpElem+1 - mesh_mapFEtoCPelem(1,cpElem) = contInts(1+i) - mesh_mapFEtoCPelem(2,cpElem) = cpElem - enddo + cpElem = cpElem+1 + mesh_mapFEtoCPelem(1,cpElem) = contInts(1+i) + mesh_mapFEtoCPelem(2,cpElem) = cpElem + enddo call math_sort(mesh_mapFEtoCPelem) @@ -545,35 +545,34 @@ end subroutine inputRead_elemNodes !> @brief Gets element type (and checks if the whole mesh comprises of only one type) !-------------------------------------------------------------------------------------------------- subroutine inputRead_elemType(elem, & - nElem,fileUnit) + nElem,fileContent) - type(tElement), intent(out) :: elem - integer, intent(in) :: & - nElem, & - fileUnit + type(tElement), intent(out) :: elem + integer, intent(in) :: nElem + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line - integer :: i,t + integer :: i,j,t,l,remainingChunks t = -1 - - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'connectivity' ) then - read (fileUnit,'(A300)') line ! Garbage line + j = 0 + do l = 1, size(fileContent) + chunkPos = IO_stringPos(fileContent(l)) + if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'connectivity' ) then do i=1,nElem ! read all elements - read (fileUnit,'(A300)') line - chunkPos = IO_stringPos(line) + chunkPos = IO_stringPos(fileContent(l+1+i+j)) if (t == -1) then - t = mapElemtype(IO_stringValue(line,chunkPos,2)) + t = mapElemtype(IO_stringValue(fileContent(l+1+i+j),chunkPos,2)) call elem%init(t) else - if (t /= mapElemtype(IO_stringValue(line,chunkPos,2))) call IO_error(191,el=t,ip=i) + if (t /= mapElemtype(IO_stringValue(fileContent(l+1+i+j),chunkPos,2))) call IO_error(191,el=t,ip=i) endif - call IO_skipChunks(fileUnit,elem%nNodes-(chunkPos(1)-2)) + remainingChunks = elem%nNodes - (chunkPos(1) - 2) + do while(remainingChunks > 0) + j = j + 1 + chunkPos = IO_stringPos(fileContent(l+1+i+j)) + remainingChunks = remainingChunks - chunkPos(1) + enddo enddo exit endif @@ -622,10 +621,10 @@ subroutine inputRead_elemType(elem, & call IO_error(error_ID=190,ext_msg=IO_lc(what)) end select -end function mapElemtype + end function mapElemtype -620 end subroutine inputRead_elemType +end subroutine inputRead_elemType !-------------------------------------------------------------------------------------------------- From 008f717c08b53097cf8ddd5a9403ad147ef932cd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 17 Oct 2019 06:10:00 +0200 Subject: [PATCH 16/20] avoid reading from file --- src/mesh_marc.f90 | 47 +++++++++++++++++++++-------------------------- 1 file changed, 21 insertions(+), 26 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 82719f0d0..4154e75d2 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -229,7 +229,7 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni call inputRead_elemNodes(node0_elem, & Nnodes,inputFile) - connectivity_elem = inputRead_connectivityElem(nElems,elem%nNodes,FILEUNIT) + connectivity_elem = inputRead_connectivityElem(nElems,elem%nNodes,inputFile) call inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, & nElems,elem%nNodes,nameElemSet,mapElemSet,& @@ -555,10 +555,10 @@ subroutine inputRead_elemType(elem, & integer :: i,j,t,l,remainingChunks t = -1 - j = 0 do l = 1, size(fileContent) chunkPos = IO_stringPos(fileContent(l)) if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'connectivity' ) then + j = 0 do i=1,nElem ! read all elements chunkPos = IO_stringPos(fileContent(l+1+i+j)) if (t == -1) then @@ -630,45 +630,40 @@ end subroutine inputRead_elemType !-------------------------------------------------------------------------------------------------- !> @brief Stores node IDs !-------------------------------------------------------------------------------------------------- -function inputRead_connectivityElem(nElem,nNodes,fileUnit) +function inputRead_connectivityElem(nElem,nNodes,fileContent) integer, intent(in) :: & nElem, & - nNodes, & !< number of nodes per element - fileUnit + nNodes !< number of nodes per element + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines integer, dimension(nNodes,nElem) :: & inputRead_connectivityElem integer, allocatable, dimension(:) :: chunkPos - character(len=300) line integer, dimension(1+nElem) :: contInts - integer :: i,j,t,sv,myVal,e,nNodesAlreadyRead + integer :: i,k,j,t,e,l,nNodesAlreadyRead - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'connectivity' ) then - read (fileUnit,'(A300)',END=620) line ! garbage line + do l = 1, size(fileContent) + chunkPos = IO_stringPos(fileContent(l)) + if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'connectivity' ) then + j = 0 do i = 1,nElem - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - e = mesh_FEasCP('elem',IO_intValue(line,chunkPos,1)) + chunkPos = IO_stringPos(fileContent(l+1+i+j)) + e = mesh_FEasCP('elem',IO_intValue(fileContent(l+1+i+j),chunkPos,1)) if (e /= 0) then ! disregard non CP elems - nNodesAlreadyRead = 0 - do j = 1,chunkPos(1)-2 - inputRead_connectivityElem(j,e) = & - mesh_FEasCP('node',IO_IntValue(line,chunkPos,j+2)) + do k = 1,chunkPos(1)-2 + inputRead_connectivityElem(k,e) = & + mesh_FEasCP('node',IO_IntValue(fileContent(l+1+i+j),chunkPos,k+2)) enddo nNodesAlreadyRead = chunkPos(1) - 2 do while(nNodesAlreadyRead < nNodes) ! read on if not all nodes in one line - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - do j = 1,chunkPos(1) - inputRead_connectivityElem(nNodesAlreadyRead+j,e) = & - mesh_FEasCP('node',IO_IntValue(line,chunkPos,j)) + j = j + 1 + chunkPos = IO_stringPos(fileContent(l+1+i+j)) + do k = 1,chunkPos(1) + inputRead_connectivityElem(nNodesAlreadyRead+k,e) = & + mesh_FEasCP('node',IO_IntValue(fileContent(l+1+i+j),chunkPos,k)) enddo nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1) enddo @@ -678,7 +673,7 @@ function inputRead_connectivityElem(nElem,nNodes,fileUnit) endif enddo -620 end function inputRead_connectivityElem +end function inputRead_connectivityElem !-------------------------------------------------------------------------------------------------- From b386dc73b22814aae891809a4fa3232d5d6eea0d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 17 Oct 2019 07:46:20 +0200 Subject: [PATCH 17/20] enable to non-transposed tensor data usually, we store per data per cell, i.e. len(shape(x)) == 3 means x is a tensor. Due to the use of transposed tensors (due to column-major storage order in Fortran), we usually want to store the transpose of (3x3) tensors. Now the default can be changed --- src/element.f90 | 2 +- src/results.f90 | 24 +++++++++++++++++++----- 2 files changed, 20 insertions(+), 6 deletions(-) diff --git a/src/element.f90 b/src/element.f90 index 8fb7d4c8b..02f5fb762 100644 --- a/src/element.f90 +++ b/src/element.f90 @@ -753,7 +753,7 @@ subroutine tElement_init(self,elemType) self%cell = CELL10 end select - self%NcellNodesPerCell = NCELLNODEPERCELL(self%cellType) + self%NcellnodesPerCell = NCELLNODEPERCELL(self%cellType) select case(self%cellType) case(1) diff --git a/src/results.f90 b/src/results.f90 index 2498bf957..025f35602 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -296,21 +296,35 @@ end subroutine results_writeVectorDataset_real !-------------------------------------------------------------------------------------------------- !> @brief stores a tensor dataset in a group !-------------------------------------------------------------------------------------------------- -subroutine results_writeTensorDataset_real(group,dataset,label,description,SIunit) +subroutine results_writeTensorDataset_real(group,dataset,label,description,SIunit,transposed) character(len=*), intent(in) :: label,group,description character(len=*), intent(in), optional :: SIunit + logical, intent(in), optional :: transposed real(pReal), intent(in), dimension(:,:,:) :: dataset integer :: i + logical :: T integer(HID_T) :: groupHandle real(pReal), dimension(:,:,:), allocatable :: dataset_transposed + + if(present(transposed)) then + T = transposed + else + T = .true. + endif - allocate(dataset_transposed,mold=dataset) - do i=1,size(dataset,3) - dataset_transposed(1:3,1:3,i) = transpose(dataset(1:3,1:3,i)) - enddo + if(T) then + if(size(dataset_transposed,1) /= size(dataset_transposed,2)) & + call IO_error(0,ext_msg='transpose non-symmetric tensor') + allocate(dataset_transposed,mold=dataset) + do i=1,size(dataset_transposed,3) + dataset_transposed(:,:,i) = transpose(dataset(:,:,i)) + enddo + else + allocate(dataset_transposed,source=dataset) + endif groupHandle = results_openGroup(group) From 88df7f2957f793f5f384360d0d9e5f30d48a7d64 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 17 Oct 2019 07:48:57 +0200 Subject: [PATCH 18/20] store cell normal directions --- src/geometry_plastic_nonlocal.f90 | 23 ++++++++++++++++------- src/mesh_marc.f90 | 11 ++++++----- 2 files changed, 22 insertions(+), 12 deletions(-) diff --git a/src/geometry_plastic_nonlocal.f90 b/src/geometry_plastic_nonlocal.f90 index b042943d0..056c7cf57 100644 --- a/src/geometry_plastic_nonlocal.f90 +++ b/src/geometry_plastic_nonlocal.f90 @@ -120,26 +120,35 @@ end subroutine geometry_plastic_nonlocal_disable !--------------------------------------------------------------------------------------------------- subroutine geometry_plastic_nonlocal_results - integer, dimension(:), allocatable :: s + integer, dimension(:), allocatable :: shp #if defined(DAMASK_HDF5) call results_openJobFile writeVolume: block real(pReal), dimension(:), allocatable :: temp - s = shape(geometry_plastic_nonlocal_IPvolume0) - temp = reshape(geometry_plastic_nonlocal_IPvolume0,[s(1)*s(2)]) + shp = shape(geometry_plastic_nonlocal_IPvolume0) + temp = reshape(geometry_plastic_nonlocal_IPvolume0,[shp(1)*shp(2)]) call results_writeDataset('geometry',temp,'v_0',& 'initial cell volume','m³') end block writeVolume - writeArea: block + writeAreas: block real(pReal), dimension(:,:), allocatable :: temp - s = shape(geometry_plastic_nonlocal_IParea0) - temp = reshape(geometry_plastic_nonlocal_IParea0,[s(1),s(2)*s(3)]) + shp = shape(geometry_plastic_nonlocal_IParea0) + temp = reshape(geometry_plastic_nonlocal_IParea0,[shp(1),shp(2)*shp(3)]) call results_writeDataset('geometry',temp,'a_0',& 'initial cell face area','m²') - end block writeArea + end block writeAreas + + writeNormals: block + real(pReal), dimension(:,:,:), allocatable :: temp + shp = shape(geometry_plastic_nonlocal_IPareaNormal0) + temp = reshape(geometry_plastic_nonlocal_IPareaNormal0,[shp(1),shp(2),shp(3)*shp(4)]) + call results_writeDataset('geometry',temp,'n_0',& + 'initial cell face normals','-',transposed=.false.) + end block writeNormals + call results_closeJobFile #endif diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 4154e75d2..66906207a 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -77,8 +77,8 @@ subroutine mesh_init(ip,el) connectivity_cell !< cell connectivity for each element,ip/cell integer, dimension(:,:), allocatable :: & connectivity_elem - - real(pReal), dimension(:,:,:,:),allocatable :: x + real(pReal), dimension(:,:,:,:),allocatable :: & + unscaledNormals write(6,'(/,a)') ' <<<+- mesh init -+>>>' @@ -120,8 +120,9 @@ subroutine mesh_init(ip,el) node0_cell,ip_reshaped) call geometry_plastic_nonlocal_setIPvolume(IPvolume(elem,node0_cell,connectivity_cell)) - x = IPareaNormal(elem,nElems,connectivity_cell,node0_cell) - call geometry_plastic_nonlocal_setIParea(norm2(x,1)) + unscaledNormals = IPareaNormal(elem,nElems,connectivity_cell,node0_cell) + call geometry_plastic_nonlocal_setIParea(norm2(unscaledNormals,1)) + call geometry_plastic_nonlocal_setIPareaNormal(unscaledNormals/spread(norm2(unscaledNormals,1),1,3)) call geometry_plastic_nonlocal_results @@ -1034,7 +1035,7 @@ function IPareaNormal(elem,nElem,connectivity,node) do e = 1,nElem do i = 1,elem%nIPs - do f = 1,size(elem%cellFace,2) + do f = 1,elem%nIPneighbors nodePos = node(1:3,connectivity(elem%cellface(1:m,f),i,e)) select case (elem%cellType) From b9027e32578f5380d3982056d673efa538ffb5c7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 17 Oct 2019 08:40:05 +0200 Subject: [PATCH 19/20] checking size of unallocated array does not work --- src/results.f90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/results.f90 b/src/results.f90 index 025f35602..93351c5f2 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -316,8 +316,7 @@ subroutine results_writeTensorDataset_real(group,dataset,label,description,SIuni endif if(T) then - if(size(dataset_transposed,1) /= size(dataset_transposed,2)) & - call IO_error(0,ext_msg='transpose non-symmetric tensor') + if(size(dataset,1) /= size(dataset,2)) call IO_error(0,ext_msg='transpose non-symmetric tensor') allocate(dataset_transposed,mold=dataset) do i=1,size(dataset_transposed,3) dataset_transposed(:,:,i) = transpose(dataset(:,:,i)) From 7ddb724861fc2d0e4ac97d60448e813f9a581b18 Mon Sep 17 00:00:00 2001 From: Franz Roters Date: Fri, 18 Oct 2019 11:27:37 +0200 Subject: [PATCH 20/20] [skip ci] corrected wrongly copied comment --- src/geometry_plastic_nonlocal.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/geometry_plastic_nonlocal.f90 b/src/geometry_plastic_nonlocal.f90 index 056c7cf57..88634c245 100644 --- a/src/geometry_plastic_nonlocal.f90 +++ b/src/geometry_plastic_nonlocal.f90 @@ -116,7 +116,7 @@ end subroutine geometry_plastic_nonlocal_disable !--------------------------------------------------------------------------------------------------- -!> @brief Frees memory used by variables only needed by plastic_nonlocal +!> @brief Writes geometry data to results file !--------------------------------------------------------------------------------------------------- subroutine geometry_plastic_nonlocal_results