From 019f0556e68078023c2981925b194f08746b9912 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 5 Oct 2019 17:18:21 +0200 Subject: [PATCH 01/48] better have explicit arguments --- src/mesh_marc.f90 | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 4a5f17674..8f286c0bb 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -207,7 +207,7 @@ subroutine mesh_init(ip,el) mesh_element(1,:) = -1 ! DEPRECATED mesh_element(2,:) = elemType ! DEPRECATED - call mesh_marc_buildElements(mesh_nElems,initialcondTableStyle,FILEUNIT) + call mesh_marc_buildElements(mesh_nElems,theMesh%elem%nNodes,initialcondTableStyle,FILEUNIT) if (myDebug) write(6,'(a)') ' Built elements'; flush(6) close (FILEUNIT) @@ -607,7 +607,7 @@ integer function mesh_marc_getElemType(nElem,fileUnit) endif enddo -contains + contains !-------------------------------------------------------------------------------------------------- !> @brief mapping of Marc element types to internal representation @@ -659,10 +659,11 @@ end function mapElemtype !-------------------------------------------------------------------------------------------------- !> @brief Stores node IDs and homogenization and microstructure ID !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_buildElements(nElem,initialcondTableStyle,fileUnit) +subroutine mesh_marc_buildElements(nElem,nNodes,initialcondTableStyle,fileUnit) integer, intent(in) :: & nElem, & + nNodes, & !< number of nodes per element initialcondTableStyle, & fileUnit @@ -688,7 +689,7 @@ subroutine mesh_marc_buildElements(nElem,initialcondTableStyle,fileUnit) mesh_element(4+j,e) = mesh_FEasCP('node',IO_IntValue(line,chunkPos,j+2)) ! CP ids of nodes enddo nNodesAlreadyRead = chunkPos(1) - 2 - do while(nNodesAlreadyRead < theMesh%elem%nNodes) ! read on if not all nodes in one line + 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) From 076aa3f72bcb888bddb2cf19e9e01488e62627e5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 5 Oct 2019 17:37:14 +0200 Subject: [PATCH 02/48] not needed --- src/mesh_marc.f90 | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 8f286c0bb..ae1eba117 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -1191,18 +1191,17 @@ end subroutine mesh_build_ipCoordinates !-------------------------------------------------------------------------------------------------- subroutine mesh_build_ipAreas - integer :: e,t,g,c,i,f,n,m + integer :: e,c,i,f,n,m real(pReal), dimension (3,FE_maxNcellnodesPerCellface) :: nodePos, normals real(pReal), dimension(3) :: normal allocate(mesh_ipArea(theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) allocate(mesh_ipAreaNormal(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) + + c = theMesh%elem%cellType - !$OMP PARALLEL DO PRIVATE(t,g,c,nodePos,normal,normals) + !$OMP PARALLEL DO PRIVATE(nodePos,normal,normals) do e = 1,theMesh%nElems ! loop over cpElems - t = mesh_element(2,e) ! get element type - g = theMesh%elem%geomType - c = theMesh%elem%cellType select case (c) case (1,2) ! 2D 3 or 4 node From 1d35699884128ffb8eeb9f8a684d6b333ce73674 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 5 Oct 2019 19:46:08 +0200 Subject: [PATCH 03/48] ip volume is only needed by plastic nonlocal --- src/mesh_marc.f90 | 93 ++++++++++++++++++++++------------------------- 1 file changed, 44 insertions(+), 49 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index ae1eba117..7a0a0bd27 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -45,7 +45,6 @@ module mesh real(pReal), dimension(:,:), allocatable :: & mesh_node, & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!! - mesh_ipVolume, & !< volume associated with IP (initially!) mesh_node0 !< node x,y,z coordinates (initially!) real(pReal), dimension(:,:,:), allocatable:: & @@ -219,8 +218,7 @@ subroutine mesh_init(ip,el) allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal) call mesh_build_ipCoordinates if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6) - call mesh_build_ipVolumes - if (myDebug) write(6,'(a)') ' Built IP volumes'; flush(6) + call mesh_build_ipAreas if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) @@ -245,7 +243,8 @@ subroutine mesh_init(ip,el) call discretization_init(mesh_element(3,:),mesh_element(4,:),& reshape(mesh_ipCoordinates,[3,theMesh%elem%nIPs*theMesh%nElems]),& mesh_node0) - call geometry_plastic_nonlocal_setIPvolume(mesh_ipVolume) + + call geometry_plastic_nonlocal_setIPvolume(IPvolume()) call geometry_plastic_nonlocal_setIPneighborhood(mesh_ipNeighborhood2) call geometry_plastic_nonlocal_setIParea(mesh_IParea) call geometry_plastic_nonlocal_setIPareaNormal(mesh_IPareaNormal) @@ -1020,7 +1019,7 @@ end function mesh_build_cellnodes !-------------------------------------------------------------------------------------------------- -!> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume' +!> @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. !> For the hexahedral cell we subdivide the cell into subvolumes of pyramidal @@ -1028,58 +1027,57 @@ end function mesh_build_cellnodes !> calculated as an average of four tetrahedals with three corners on the cell face !> and one corner at the central ip. !-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipVolumes +function IPvolume() + real(pReal), dimension(theMesh%elem%nIPs,theMesh%nElems) :: IPvolume + integer :: e,i,c,m,f,n real(pReal), dimension(size(theMesh%elem%cellFace,1),size(theMesh%elem%cellFace,2)) :: subvolume - allocate(mesh_ipVolume(theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal) c = theMesh%elem%cellType m = FE_NcellnodesPerCellface(c) - !$OMP PARALLEL DO PRIVATE(f,n,subvolume) - do e = 1,theMesh%nElems - select case (c) + do e = 1,theMesh%nElems + select case (c) - case (1) ! 2D 3node - forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element - mesh_ipVolume(i,e) = math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & - theMesh%node_0(1:3,mesh_cell2(2,i,e)), & - theMesh%node_0(1:3,mesh_cell2(3,i,e))) + case (1) ! 2D 3node + forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element + IPvolume(i,e) = math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & + theMesh%node_0(1:3,mesh_cell2(2,i,e)), & + theMesh%node_0(1:3,mesh_cell2(3,i,e))) - case (2) ! 2D 4node - forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element - mesh_ipVolume(i,e) = math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & ! here we assume a planar shape, so division in two triangles suffices - theMesh%node_0(1:3,mesh_cell2(2,i,e)), & - theMesh%node_0(1:3,mesh_cell2(3,i,e))) & - + math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(3,i,e)), & - theMesh%node_0(1:3,mesh_cell2(4,i,e)), & - theMesh%node_0(1:3,mesh_cell2(1,i,e))) + case (2) ! 2D 4node + forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element + IPvolume(i,e) = math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & ! here we assume a planar shape, so division in two triangles suffices + theMesh%node_0(1:3,mesh_cell2(2,i,e)), & + theMesh%node_0(1:3,mesh_cell2(3,i,e))) & + + math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(3,i,e)), & + theMesh%node_0(1:3,mesh_cell2(4,i,e)), & + theMesh%node_0(1:3,mesh_cell2(1,i,e))) - case (3) ! 3D 4node - forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element - mesh_ipVolume(i,e) = math_volTetrahedron(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & - theMesh%node_0(1:3,mesh_cell2(2,i,e)), & - theMesh%node_0(1:3,mesh_cell2(3,i,e)), & - theMesh%node_0(1:3,mesh_cell2(4,i,e))) + case (3) ! 3D 4node + forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element + IPvolume(i,e) = math_volTetrahedron(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & + theMesh%node_0(1:3,mesh_cell2(2,i,e)), & + theMesh%node_0(1:3,mesh_cell2(3,i,e)), & + theMesh%node_0(1:3,mesh_cell2(4,i,e))) - case (4) ! 3D 8node - do i = 1,theMesh%elem%nIPs ! loop over ips=cells in this element - subvolume = 0.0_pReal - forall(f = 1:FE_NipNeighbors(c), n = 1:m) & - subvolume(n,f) = math_volTetrahedron(& - mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface( n ,f),i,e)), & - mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(1+mod(n ,m),f),i,e)), & - mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(1+mod(n+1,m),f),i,e)), & - mesh_ipCoordinates(1:3,i,e)) - mesh_ipVolume(i,e) = 0.5_pReal * sum(subvolume) ! each subvolume is based on four tetrahedrons, altough the face consists of only two triangles -> averaging factor two - enddo + case (4) ! 3D 8node + do i = 1,theMesh%elem%nIPs ! loop over ips=cells in this element + subvolume = 0.0_pReal + forall(f = 1:FE_NipNeighbors(c), n = 1:m) & + subvolume(n,f) = math_volTetrahedron(& + mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface( n ,f),i,e)), & + mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(1+mod(n ,m),f),i,e)), & + mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(1+mod(n+1,m),f),i,e)), & + mesh_ipCoordinates(1:3,i,e)) + IPvolume(i,e) = 0.5_pReal * sum(subvolume) ! each subvolume is based on four tetrahedrons, altough the face consists of only two triangles -> averaging factor two + enddo - end select - enddo - !$OMP END PARALLEL DO + end select + enddo -end subroutine mesh_build_ipVolumes +end function IPvolume subroutine IP_neighborhood2 @@ -1155,15 +1153,12 @@ subroutine IP_neighborhood2 end subroutine IP_neighborhood2 !-------------------------------------------------------------------------------------------------- -!> @brief Calculates IP Coordinates. Allocates global array 'mesh_ipCoordinates' -! Called by all solvers in mesh_init in order to initialize the ip coordinates. -! Later on the current ip coordinates are directly prvided by the spectral solver and by Abaqus, -! so no need to use this subroutine anymore; Marc however only provides nodal displacements, +!> @brief Calculates IP Coordinates. +! Marc however only provides nodal displacements, ! so in this case the ip coordinates are always calculated on the basis of this subroutine. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! FOR THE MOMENT THIS SUBROUTINE ACTUALLY CALCULATES THE CELL CENTER AND NOT THE IP COORDINATES, ! AS THE IP IS NOT (ALWAYS) LOCATED IN THE CENTER OF THE IP VOLUME. -! HAS TO BE CHANGED IN A LATER VERSION. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !-------------------------------------------------------------------------------------------------- subroutine mesh_build_ipCoordinates From 18b8e71f69f19091a55835a773095d39ac618b46 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 5 Oct 2019 20:09:01 +0200 Subject: [PATCH 04/48] polishing --- src/mesh_marc.f90 | 66 +++++++++++++++++++++++++++-------------------- 1 file changed, 38 insertions(+), 28 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 7a0a0bd27..134c908d4 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -219,19 +219,35 @@ subroutine mesh_init(ip,el) call mesh_build_ipCoordinates if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6) - call mesh_build_ipAreas - if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) call IP_neighborhood2 if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6) - 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) & - call IO_error(602,ext_msg='IP') ! selected element does not have requested IP + call discretization_init(mesh_element(3,:),mesh_element(4,:),& + reshape(mesh_ipCoordinates,[3,theMesh%elem%nIPs*theMesh%nElems]),& + mesh_node0) + + ! calculate and store information needed for nonlocal plasticity model + call geometry_plastic_nonlocal_setIPvolume(IPvolume()) + call geometry_plastic_nonlocal_setIPneighborhood(mesh_ipNeighborhood2) + + call mesh_build_ipAreas + if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) + call geometry_plastic_nonlocal_setIParea(mesh_IParea) + call geometry_plastic_nonlocal_setIPareaNormal(mesh_IPareaNormal) + deallocate(mesh_IParea) + deallocate(mesh_IPareaNormal) + + ! sanity checks + if (usePingPong .and. (mesh_Nelems /= theMesh%nElems)) & + call IO_error(600) + if (debug_e < 1 .or. debug_e > theMesh%nElems) & + call IO_error(602,ext_msg='element') + if (debug_i < 1 .or. debug_i > theMesh%elem%nIPs) & + call IO_error(602,ext_msg='IP') + + ! solving related 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_execIP(2,:) = theMesh%elem%nIPs @@ -239,15 +255,6 @@ subroutine mesh_init(ip,el) allocate(calcMode(theMesh%elem%nIPs,theMesh%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" - - call discretization_init(mesh_element(3,:),mesh_element(4,:),& - reshape(mesh_ipCoordinates,[3,theMesh%elem%nIPs*theMesh%nElems]),& - mesh_node0) - - call geometry_plastic_nonlocal_setIPvolume(IPvolume()) - call geometry_plastic_nonlocal_setIPneighborhood(mesh_ipNeighborhood2) - call geometry_plastic_nonlocal_setIParea(mesh_IParea) - call geometry_plastic_nonlocal_setIPareaNormal(mesh_IPareaNormal) end subroutine mesh_init @@ -1040,30 +1047,30 @@ function IPvolume() do e = 1,theMesh%nElems select case (c) - case (1) ! 2D 3node - forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element + case (1) ! 2D 3node + forall (i = 1:theMesh%elem%nIPs) & IPvolume(i,e) = math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & theMesh%node_0(1:3,mesh_cell2(2,i,e)), & theMesh%node_0(1:3,mesh_cell2(3,i,e))) - case (2) ! 2D 4node - forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element - IPvolume(i,e) = math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & ! here we assume a planar shape, so division in two triangles suffices + case (2) ! 2D 4node + forall (i = 1:theMesh%elem%nIPs) & + IPvolume(i,e) = math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & ! here we assume a planar shape, so division in two triangles suffices theMesh%node_0(1:3,mesh_cell2(2,i,e)), & theMesh%node_0(1:3,mesh_cell2(3,i,e))) & + math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(3,i,e)), & theMesh%node_0(1:3,mesh_cell2(4,i,e)), & theMesh%node_0(1:3,mesh_cell2(1,i,e))) - case (3) ! 3D 4node - forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element + case (3) ! 3D 4node + forall (i = 1:theMesh%elem%nIPs) & IPvolume(i,e) = math_volTetrahedron(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & theMesh%node_0(1:3,mesh_cell2(2,i,e)), & theMesh%node_0(1:3,mesh_cell2(3,i,e)), & theMesh%node_0(1:3,mesh_cell2(4,i,e))) - case (4) ! 3D 8node - do i = 1,theMesh%elem%nIPs ! loop over ips=cells in this element + case (4) ! 3D 8node + do i = 1,theMesh%elem%nIPs subvolume = 0.0_pReal forall(f = 1:FE_NipNeighbors(c), n = 1:m) & subvolume(n,f) = math_volTetrahedron(& @@ -1071,7 +1078,7 @@ function IPvolume() mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(1+mod(n ,m),f),i,e)), & mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(1+mod(n+1,m),f),i,e)), & mesh_ipCoordinates(1:3,i,e)) - IPvolume(i,e) = 0.5_pReal * sum(subvolume) ! each subvolume is based on four tetrahedrons, altough the face consists of only two triangles -> averaging factor two + IPvolume(i,e) = 0.5_pReal * sum(subvolume) ! each subvolume is based on four tetrahedrons, but the face consists of two triangles -> average by 2 enddo end select @@ -1080,6 +1087,9 @@ function IPvolume() end function IPvolume +!-------------------------------------------------------------------------------------------------- +!> @brief Calculates IP neighborhood +!-------------------------------------------------------------------------------------------------- subroutine IP_neighborhood2 integer, dimension(:,:), allocatable :: faces @@ -1190,7 +1200,7 @@ subroutine mesh_build_ipAreas real(pReal), dimension (3,FE_maxNcellnodesPerCellface) :: nodePos, normals real(pReal), dimension(3) :: normal - allocate(mesh_ipArea(theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) + allocate(mesh_ipArea (theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) allocate(mesh_ipAreaNormal(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) c = theMesh%elem%cellType From 7f403ad50e24999ec3e04d94a94cf1e0a97f80a8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 5 Oct 2019 20:23:33 +0200 Subject: [PATCH 05/48] avoid global variables --- src/mesh_marc.f90 | 37 ++++++++++++++++++------------------- 1 file changed, 18 insertions(+), 19 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 134c908d4..20f1b62ff 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -47,10 +47,6 @@ module mesh mesh_node, & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!! mesh_node0 !< node x,y,z coordinates (initially!) - real(pReal), dimension(:,:,:), allocatable:: & - mesh_ipArea !< area of interface to neighboring IP (initially!) - real(pReal),dimension(:,:,:,:), allocatable :: & - mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!) ! -------------------------------------------------------------------------------------------------- @@ -154,6 +150,10 @@ subroutine mesh_init(ip,el) hypoelasticTableStyle, & initialcondTableStyle logical :: myDebug + real(pReal), dimension(:,:,:), allocatable:: & + mesh_ipArea !< area of interface to neighboring IP (initially!) + real(pReal),dimension(:,:,:,:), allocatable :: & + mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!) write(6,'(/,a)') ' <<<+- mesh init -+>>>' @@ -232,12 +232,10 @@ subroutine mesh_init(ip,el) call geometry_plastic_nonlocal_setIPvolume(IPvolume()) call geometry_plastic_nonlocal_setIPneighborhood(mesh_ipNeighborhood2) - call mesh_build_ipAreas + call mesh_build_ipAreas(mesh_IParea,mesh_IPareaNormal) if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) call geometry_plastic_nonlocal_setIParea(mesh_IParea) call geometry_plastic_nonlocal_setIPareaNormal(mesh_IPareaNormal) - deallocate(mesh_IParea) - deallocate(mesh_IPareaNormal) ! sanity checks if (usePingPong .and. (mesh_Nelems /= theMesh%nElems)) & @@ -1194,19 +1192,21 @@ end subroutine mesh_build_ipCoordinates !-------------------------------------------------------------------------------------------------- !> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal' !-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipAreas +subroutine mesh_build_ipAreas(area,areaNormal) + + real(pReal), dimension(:,:,:), allocatable, intent(out) :: area + real(pReal), dimension(:,:,:,:), allocatable, intent(out) :: areaNormal integer :: e,c,i,f,n,m real(pReal), dimension (3,FE_maxNcellnodesPerCellface) :: nodePos, normals real(pReal), dimension(3) :: normal - allocate(mesh_ipArea (theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) - allocate(mesh_ipAreaNormal(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) + allocate(area (theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) + allocate(areaNormal(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) c = theMesh%elem%cellType - !$OMP PARALLEL DO PRIVATE(nodePos,normal,normals) - do e = 1,theMesh%nElems ! loop over cpElems + do e = 1,theMesh%nElems select case (c) case (1,2) ! 2D 3 or 4 node @@ -1217,8 +1217,8 @@ subroutine mesh_build_ipAreas 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 - mesh_ipArea(f,i,e) = norm2(normal) - mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal + area(f,i,e) = norm2(normal) + areaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal enddo enddo @@ -1229,8 +1229,8 @@ subroutine mesh_build_ipAreas nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(n,f),i,e)) normal = math_cross(nodePos(1:3,2) - nodePos(1:3,1), & nodePos(1:3,3) - nodePos(1:3,1)) - mesh_ipArea(f,i,e) = norm2(normal) - mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal + area(f,i,e) = norm2(normal) + areaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal enddo enddo @@ -1249,14 +1249,13 @@ subroutine mesh_build_ipAreas * 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) - mesh_ipArea(f,i,e) = norm2(normal) - mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) + area(f,i,e) = norm2(normal) + areaNormal(1:3,f,i,e) = normal / norm2(normal) enddo enddo end select enddo - !$OMP END PARALLEL DO end subroutine mesh_build_ipAreas From 04272f88d5f522eb590389094778604f71b12ded Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 8 Oct 2019 09:30:32 +0200 Subject: [PATCH 06/48] untangling - avoid public variables - openMP in initialization hardly useful - structure of init should reflect tasks: 1) reading 2) discretization 3) nonlocal stuff --- src/DAMASK_marc.f90 | 2 +- src/mesh_marc.f90 | 138 +++++++++++++++++++------------------------- 2 files changed, 61 insertions(+), 79 deletions(-) diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index 3f527daf5..0cc815e8e 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -265,7 +265,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & call debug_reset() ! resets debugging outdatedFFN1 = .false. cycleCounter = cycleCounter + 1 - mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates + !mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates call mesh_build_ipCoordinates() ! update ip coordinates endif if (outdatedByNewInc) then diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 20f1b62ff..0582163d9 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -47,6 +47,10 @@ module mesh mesh_node, & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!! mesh_node0 !< node x,y,z coordinates (initially!) + real(pReal), dimension(:,:,:), allocatable:: & + mesh_ipArea !< area of interface to neighboring IP (initially!) + real(pReal),dimension(:,:,:,:), allocatable :: & + mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!) ! -------------------------------------------------------------------------------------------------- @@ -108,7 +112,10 @@ integer, dimension(:,:), allocatable :: & 6 & ! (3D 8node) ],pInt) - + integer, dimension(:), allocatable :: & + microstructureAt, & + homogenizationAt + integer :: & mesh_NelemSets character(len=64), dimension(:), allocatable :: & @@ -150,10 +157,6 @@ subroutine mesh_init(ip,el) hypoelasticTableStyle, & initialcondTableStyle logical :: myDebug - real(pReal), dimension(:,:,:), allocatable:: & - mesh_ipArea !< area of interface to neighboring IP (initially!) - real(pReal),dimension(:,:,:,:), allocatable :: & - mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!) write(6,'(/,a)') ' <<<+- mesh init -+>>>' @@ -161,40 +164,24 @@ subroutine mesh_init(ip,el) myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0) + ! parsing Marc input file call IO_open_inputFile(FILEUNIT,modelName) fileFormatVersion = mesh_marc_get_fileFormat(FILEUNIT) - if (myDebug) write(6,'(a)') ' Got input file format'; flush(6) - call mesh_marc_get_tableStyles(initialcondTableStyle,hypoelasticTableStyle,FILEUNIT) - if (myDebug) write(6,'(a)') ' Got table styles'; flush(6) - - if (fileFormatVersion > 12) then + if (fileFormatVersion > 12) & Marc_matNumber = mesh_marc_get_matNumber(FILEUNIT,hypoelasticTableStyle) - if (myDebug) write(6,'(a)') ' Got hypoleastic material number'; flush(6) - endif - call mesh_marc_count_nodesAndElements(mesh_nNodes, mesh_nElems, FILEUNIT) - if (myDebug) write(6,'(a)') ' Counted nodes/elements'; flush(6) - - call mesh_marc_count_elementSets(mesh_NelemSets,mesh_maxNelemInSet,FILEUNIT) - if (myDebug) write(6,'(a)') ' Counted element sets'; flush(6) - + call mesh_marc_count_elementSets(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 mesh_marc_map_elementSets(mesh_nameElemSet,mesh_mapElemSet,FILEUNIT) - if (myDebug) write(6,'(a)') ' Mapped element sets'; flush(6) - allocate (mesh_mapFEtoCPelem(2,mesh_nElems), source = 0) - call mesh_marc_map_elements(hypoelasticTableStyle,mesh_nameElemSet,mesh_mapElemSet,mesh_nElems,fileFormatVersion,FILEUNIT) - if (myDebug) write(6,'(a)') ' Mapped elements'; flush(6) - + call mesh_marc_map_elements(hypoelasticTableStyle,mesh_nameElemSet,mesh_mapElemSet,mesh_nElems,fileFormatVersion,FILEUNIT) allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes),source=0) call mesh_marc_map_nodes(mesh_Nnodes,FILEUNIT) !ToDo: don't work on global variables - if (myDebug) write(6,'(a)') ' Mapped nodes'; flush(6) mesh_node0 = mesh_marc_build_nodes(mesh_Nnodes,FILEUNIT) mesh_node = mesh_node0 - if (myDebug) write(6,'(a)') ' Built nodes'; flush(6) elemType = mesh_marc_getElemType(mesh_nElems,FILEUNIT) if (myDebug) write(6,'(a)') ' Counted CP sizes'; flush(6) @@ -202,11 +189,14 @@ subroutine mesh_init(ip,el) call theMesh%init('mesh',elemType,mesh_node0) call theMesh%setNelems(mesh_nElems) + allocate(microstructureAt(theMesh%nElems), source=0) + allocate(homogenizationAt(theMesh%nElems), source=0) allocate(mesh_element(4+theMesh%elem%nNodes,theMesh%nElems), source=0) mesh_element(1,:) = -1 ! DEPRECATED mesh_element(2,:) = elemType ! DEPRECATED - call mesh_marc_buildElements(mesh_nElems,theMesh%elem%nNodes,initialcondTableStyle,FILEUNIT) + call mesh_marc_buildElements(microstructureAt,homogenizationAt, & + mesh_nElems,theMesh%elem%nNodes,initialcondTableStyle,FILEUNIT) if (myDebug) write(6,'(a)') ' Built elements'; flush(6) close (FILEUNIT) @@ -219,33 +209,19 @@ subroutine mesh_init(ip,el) call mesh_build_ipCoordinates if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6) + call mesh_build_ipAreas + if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) call IP_neighborhood2 if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6) - - call discretization_init(mesh_element(3,:),mesh_element(4,:),& - reshape(mesh_ipCoordinates,[3,theMesh%elem%nIPs*theMesh%nElems]),& - mesh_node0) - - ! calculate and store information needed for nonlocal plasticity model - call geometry_plastic_nonlocal_setIPvolume(IPvolume()) - call geometry_plastic_nonlocal_setIPneighborhood(mesh_ipNeighborhood2) - - call mesh_build_ipAreas(mesh_IParea,mesh_IPareaNormal) - if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) - call geometry_plastic_nonlocal_setIParea(mesh_IParea) - call geometry_plastic_nonlocal_setIPareaNormal(mesh_IPareaNormal) - - ! sanity checks if (usePingPong .and. (mesh_Nelems /= theMesh%nElems)) & - call IO_error(600) + 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') + 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') + call IO_error(602,ext_msg='IP') ! selected element does not have requested IP - ! solving related 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_execIP(2,:) = theMesh%elem%nIPs @@ -253,6 +229,15 @@ subroutine mesh_init(ip,el) allocate(calcMode(theMesh%elem%nIPs,theMesh%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" + + call discretization_init(microstructureAt,homogenizationAt,& + reshape(mesh_ipCoordinates,[3,theMesh%elem%nIPs*theMesh%nElems]),& + mesh_node0) + + call geometry_plastic_nonlocal_setIPvolume(IPvolume()) + call geometry_plastic_nonlocal_setIPneighborhood(mesh_ipNeighborhood2) + call geometry_plastic_nonlocal_setIParea(mesh_IParea) + call geometry_plastic_nonlocal_setIPareaNormal(mesh_IPareaNormal) end subroutine mesh_init @@ -663,8 +648,12 @@ end function mapElemtype !-------------------------------------------------------------------------------------------------- !> @brief Stores node IDs and homogenization and microstructure ID !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_buildElements(nElem,nNodes,initialcondTableStyle,fileUnit) +subroutine mesh_marc_buildElements(microstructureAt,homogenizationAt, & + nElem,nNodes,initialcondTableStyle,fileUnit) + integer, dimension(:), intent(out) :: & + microstructureAt, & + homogenizationAt integer, intent(in) :: & nElem, & nNodes, & !< number of nodes per element @@ -740,7 +729,8 @@ subroutine mesh_marc_buildElements(nElem,nNodes,initialcondTableStyle,fileUnit) (fileUnit,theMesh%nElems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) do i = 1,contInts(1) e = mesh_FEasCP('elem',contInts(1+i)) - mesh_element(1+sv,e) = myVal + if (sv == 2) microstructureAt(e) = myVal + if (sv == 3) homogenizationAt(e) = myVal enddo if (initialcondTableStyle == 0) read (fileUnit,'(A300)',END=630) line ! ignore IP range for old table style read (fileUnit,'(A300)',END=630) line @@ -1007,7 +997,6 @@ function mesh_build_cellnodes() myCoords mesh_build_cellnodes = 0.0_pReal -!$OMP PARALLEL DO PRIVATE(e,localCellnodeID,myCoords) do n = 1,mesh_Ncellnodes ! loop over cell nodes e = mesh_cellnodeParent(1,n) localCellnodeID = mesh_cellnodeParent(2,n) @@ -1018,7 +1007,6 @@ function mesh_build_cellnodes() enddo mesh_build_cellnodes(1:3,n) = myCoords / sum(theMesh%elem%cellNodeParentNodeWeights(:,localCellnodeID)) enddo -!$OMP END PARALLEL DO end function mesh_build_cellnodes @@ -1045,30 +1033,30 @@ function IPvolume() do e = 1,theMesh%nElems select case (c) - case (1) ! 2D 3node - forall (i = 1:theMesh%elem%nIPs) & + case (1) ! 2D 3node + forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element IPvolume(i,e) = math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & theMesh%node_0(1:3,mesh_cell2(2,i,e)), & theMesh%node_0(1:3,mesh_cell2(3,i,e))) - case (2) ! 2D 4node - forall (i = 1:theMesh%elem%nIPs) & - IPvolume(i,e) = math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & ! here we assume a planar shape, so division in two triangles suffices + case (2) ! 2D 4node + forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element + IPvolume(i,e) = math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & ! here we assume a planar shape, so division in two triangles suffices theMesh%node_0(1:3,mesh_cell2(2,i,e)), & theMesh%node_0(1:3,mesh_cell2(3,i,e))) & + math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(3,i,e)), & theMesh%node_0(1:3,mesh_cell2(4,i,e)), & theMesh%node_0(1:3,mesh_cell2(1,i,e))) - case (3) ! 3D 4node - forall (i = 1:theMesh%elem%nIPs) & + case (3) ! 3D 4node + forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element IPvolume(i,e) = math_volTetrahedron(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & theMesh%node_0(1:3,mesh_cell2(2,i,e)), & theMesh%node_0(1:3,mesh_cell2(3,i,e)), & theMesh%node_0(1:3,mesh_cell2(4,i,e))) - case (4) ! 3D 8node - do i = 1,theMesh%elem%nIPs + case (4) ! 3D 8node + do i = 1,theMesh%elem%nIPs ! loop over ips=cells in this element subvolume = 0.0_pReal forall(f = 1:FE_NipNeighbors(c), n = 1:m) & subvolume(n,f) = math_volTetrahedron(& @@ -1076,7 +1064,7 @@ function IPvolume() mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(1+mod(n ,m),f),i,e)), & mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(1+mod(n+1,m),f),i,e)), & mesh_ipCoordinates(1:3,i,e)) - IPvolume(i,e) = 0.5_pReal * sum(subvolume) ! each subvolume is based on four tetrahedrons, but the face consists of two triangles -> average by 2 + IPvolume(i,e) = 0.5_pReal * sum(subvolume) ! each subvolume is based on four tetrahedrons, altough the face consists of only two triangles -> averaging factor two enddo end select @@ -1085,9 +1073,6 @@ function IPvolume() end function IPvolume -!-------------------------------------------------------------------------------------------------- -!> @brief Calculates IP neighborhood -!-------------------------------------------------------------------------------------------------- subroutine IP_neighborhood2 integer, dimension(:,:), allocatable :: faces @@ -1174,7 +1159,7 @@ subroutine mesh_build_ipCoordinates integer :: e,i,n real(pReal), dimension(3) :: myCoords - !$OMP PARALLEL DO PRIVATE(myCoords) + do e = 1,theMesh%nElems do i = 1,theMesh%elem%nIPs myCoords = 0.0_pReal @@ -1184,7 +1169,6 @@ subroutine mesh_build_ipCoordinates mesh_ipCoordinates(1:3,i,e) = myCoords / real(theMesh%elem%nCellnodesPerCell,pReal) enddo enddo - !$OMP END PARALLEL DO end subroutine mesh_build_ipCoordinates @@ -1192,21 +1176,19 @@ end subroutine mesh_build_ipCoordinates !-------------------------------------------------------------------------------------------------- !> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal' !-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipAreas(area,areaNormal) - - real(pReal), dimension(:,:,:), allocatable, intent(out) :: area - real(pReal), dimension(:,:,:,:), allocatable, intent(out) :: areaNormal +subroutine mesh_build_ipAreas integer :: e,c,i,f,n,m real(pReal), dimension (3,FE_maxNcellnodesPerCellface) :: nodePos, normals real(pReal), dimension(3) :: normal - allocate(area (theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) - allocate(areaNormal(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) + allocate(mesh_ipArea(theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) + allocate(mesh_ipAreaNormal(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) c = theMesh%elem%cellType - do e = 1,theMesh%nElems + + do e = 1,theMesh%nElems ! loop over cpElems select case (c) case (1,2) ! 2D 3 or 4 node @@ -1217,8 +1199,8 @@ subroutine mesh_build_ipAreas(area,areaNormal) 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 - area(f,i,e) = norm2(normal) - areaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal + mesh_ipArea(f,i,e) = norm2(normal) + mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal enddo enddo @@ -1229,8 +1211,8 @@ subroutine mesh_build_ipAreas(area,areaNormal) nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(n,f),i,e)) normal = math_cross(nodePos(1:3,2) - nodePos(1:3,1), & nodePos(1:3,3) - nodePos(1:3,1)) - area(f,i,e) = norm2(normal) - areaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal + mesh_ipArea(f,i,e) = norm2(normal) + mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal enddo enddo @@ -1249,14 +1231,14 @@ subroutine mesh_build_ipAreas(area,areaNormal) * 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) - area(f,i,e) = norm2(normal) - areaNormal(1:3,f,i,e) = normal / norm2(normal) + mesh_ipArea(f,i,e) = norm2(normal) + mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) enddo enddo end select enddo - + end subroutine mesh_build_ipAreas From 5e79c360e96ec0fcb254673c02ef4ebda7d2b351 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 8 Oct 2019 10:03:03 +0200 Subject: [PATCH 07/48] no need for double definition --- src/mesh_marc.f90 | 20 ++++++-------------- 1 file changed, 6 insertions(+), 14 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 0582163d9..75ec282e6 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -96,14 +96,6 @@ integer, dimension(:,:), allocatable :: & ],pInt) - integer, dimension(FE_Ncelltypes), parameter :: FE_NcellnodesPerCellface = & !< number of cell nodes per cell face in a specific cell type - int([& - 2, & ! (2D 3node) - 2, & ! (2D 4node) - 3, & ! (3D 4node) - 4 & ! (3D 8node) - ],pInt) - integer, dimension(FE_Ncelltypes), parameter :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type int([& 3, & ! (2D 3node) @@ -1028,7 +1020,7 @@ function IPvolume() real(pReal), dimension(size(theMesh%elem%cellFace,1),size(theMesh%elem%cellFace,2)) :: subvolume c = theMesh%elem%cellType - m = FE_NcellnodesPerCellface(c) + m = size(theMesh%elem%cellFace,1) do e = 1,theMesh%nElems select case (c) @@ -1194,7 +1186,7 @@ subroutine mesh_build_ipAreas case (1,2) ! 2D 3 or 4 node do i = 1,theMesh%elem%nIPs do f = 1,FE_NipNeighbors(c) ! loop over cell faces - forall(n = 1:FE_NcellnodesPerCellface(c)) & + forall(n = 1: size(theMesh%elem%cellface,1)) & nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(n,f),i,e)) normal(1) = nodePos(2,2) - nodePos(2,1) ! x_normal = y_connectingVector normal(2) = -(nodePos(1,2) - nodePos(1,1)) ! y_normal = -x_connectingVector @@ -1207,7 +1199,7 @@ subroutine mesh_build_ipAreas case (3) ! 3D 4node do i = 1,theMesh%elem%nIPs do f = 1,FE_NipNeighbors(c) ! loop over cell faces - forall(n = 1:FE_NcellnodesPerCellface(c)) & + forall(n = 1: size(theMesh%elem%cellface,1)) & nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(n,f),i,e)) normal = math_cross(nodePos(1:3,2) - nodePos(1:3,1), & nodePos(1:3,3) - nodePos(1:3,1)) @@ -1221,12 +1213,12 @@ subroutine mesh_build_ipAreas ! 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 - m = FE_NcellnodesPerCellface(c) + m = size(theMesh%elem%cellFace,1) do i = 1,theMesh%elem%nIPs do f = 1,FE_NipNeighbors(c) ! loop over cell faces - forall(n = 1:FE_NcellnodesPerCellface(c)) & + forall(n = 1: size(theMesh%elem%cellface,1)) & nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(n,f),i,e)) - forall(n = 1:FE_NcellnodesPerCellface(c)) & + forall(n = 1: size(theMesh%elem%cellface,1)) & 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)) From 098f2903ea0703a4c460812072d51a3214e0475b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 8 Oct 2019 17:04:27 +0200 Subject: [PATCH 08/48] not used anymore --- src/mesh_marc.f90 | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 75ec282e6..a60d754d6 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -39,9 +39,6 @@ module mesh integer, dimension(:,:), allocatable :: & mesh_element - - integer, dimension(:,:,:,:), allocatable :: & - mesh_ipNeighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] real(pReal), dimension(:,:), allocatable :: & mesh_node, & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!! From 16fd608da61953c99650d7720ac9a450a12a3a58 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 8 Oct 2019 17:12:53 +0200 Subject: [PATCH 09/48] more reasonable name --- src/mesh_marc.f90 | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index a60d754d6..60808a80e 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -38,7 +38,7 @@ module mesh !-------------------------------------------------------------------------------------------------- integer, dimension(:,:), allocatable :: & - mesh_element + mesh_FEnodes real(pReal), dimension(:,:), allocatable :: & mesh_node, & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!! @@ -180,9 +180,7 @@ subroutine mesh_init(ip,el) allocate(microstructureAt(theMesh%nElems), source=0) allocate(homogenizationAt(theMesh%nElems), source=0) - allocate(mesh_element(4+theMesh%elem%nNodes,theMesh%nElems), source=0) - mesh_element(1,:) = -1 ! DEPRECATED - mesh_element(2,:) = elemType ! DEPRECATED + allocate(mesh_FEnodes(theMesh%elem%nNodes,theMesh%nElems), source=0) call mesh_marc_buildElements(microstructureAt,homogenizationAt, & mesh_nElems,theMesh%elem%nNodes,initialcondTableStyle,FILEUNIT) @@ -668,14 +666,14 @@ subroutine mesh_marc_buildElements(microstructureAt,homogenizationAt, & if (e /= 0) then ! disregard non CP elems nNodesAlreadyRead = 0 do j = 1,chunkPos(1)-2 - mesh_element(4+j,e) = mesh_FEasCP('node',IO_IntValue(line,chunkPos,j+2)) ! CP ids of nodes + mesh_FEnodes(j,e) = mesh_FEasCP('node',IO_IntValue(line,chunkPos,j+2)) ! CP ids of nodes 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) - mesh_element(4+nNodesAlreadyRead+j,e) = mesh_FEasCP('node',IO_IntValue(line,chunkPos,j)) ! CP ids of nodes + mesh_FEnodes(nNodesAlreadyRead+j,e) = mesh_FEasCP('node',IO_IntValue(line,chunkPos,j)) ! CP ids of nodes enddo nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1) enddo @@ -689,12 +687,12 @@ subroutine mesh_marc_buildElements(microstructureAt,homogenizationAt, & #if defined(DAMASK_HDF5) call results_openJobFile call HDF5_closeGroup(results_addGroup('geometry')) - call results_writeDataset('geometry',mesh_element(5:,:),'C',& + call results_writeDataset('geometry',mesh_FEnodes,'C',& 'connectivity of the elements','-') call results_closeJobFile #endif - call buildCells(theMesh,theMesh%elem,mesh_element(5:,:)) + call buildCells(theMesh,theMesh%elem,mesh_FEnodes) read (fileUnit,'(A300)',END=630) line do @@ -937,7 +935,7 @@ subroutine mesh_build_cellconnectivity do n = 1,theMesh%elem%NcellnodesPerCell localCellnodeID = theMesh%elem%cell(n,i) if (localCellnodeID <= FE_NmatchingNodes(theMesh%elem%geomType)) then ! this cell node is a matching node - matchingNodeID = mesh_element(4+localCellnodeID,e) + matchingNodeID = mesh_FEnodes(localCellnodeID,e) if (matchingNode2cellnode(matchingNodeID) == 0) then ! if this matching node does not yet exist in the glbal cell node list ... mesh_Ncellnodes = mesh_Ncellnodes + 1 ! ... count it as cell node ... matchingNode2cellnode(matchingNodeID) = mesh_Ncellnodes ! ... and remember its global ID @@ -991,7 +989,7 @@ function mesh_build_cellnodes() localCellnodeID = mesh_cellnodeParent(2,n) myCoords = 0.0_pReal do m = 1,theMesh%elem%nNodes - myCoords = myCoords + mesh_node(1:3,mesh_element(4+m,e)) & + myCoords = myCoords + mesh_node(1:3,mesh_FEnodes(m,e)) & * theMesh%elem%cellNodeParentNodeWeights(m,localCellnodeID) enddo mesh_build_cellnodes(1:3,n) = myCoords / sum(theMesh%elem%cellNodeParentNodeWeights(:,localCellnodeID)) From 040cd3e35d0982980fab841eda752a18137a14a6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 8 Oct 2019 17:26:02 +0200 Subject: [PATCH 10/48] no public variables --- src/mesh_marc.f90 | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 60808a80e..9e4c515ba 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -119,8 +119,7 @@ integer, dimension(:,:), allocatable :: & integer, dimension(:,:,:,:), allocatable :: & mesh_ipNeighborhood2 !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] - integer, dimension(:), allocatable :: & - Marc_matNumber !< array of material numbers for hypoelastic material (Marc only) + public :: & mesh_init, & @@ -145,6 +144,8 @@ subroutine mesh_init(ip,el) mesh_nElems, & hypoelasticTableStyle, & initialcondTableStyle + integer, dimension(:), allocatable :: & + marc_matNumber !< array of material numbers for hypoelastic material (Marc only) logical :: myDebug write(6,'(/,a)') ' <<<+- mesh init -+>>>' @@ -158,14 +159,15 @@ subroutine mesh_init(ip,el) fileFormatVersion = mesh_marc_get_fileFormat(FILEUNIT) call mesh_marc_get_tableStyles(initialcondTableStyle,hypoelasticTableStyle,FILEUNIT) if (fileFormatVersion > 12) & - Marc_matNumber = mesh_marc_get_matNumber(FILEUNIT,hypoelasticTableStyle) + marc_matNumber = mesh_marc_get_matNumber(FILEUNIT,hypoelasticTableStyle) call mesh_marc_count_nodesAndElements(mesh_nNodes, mesh_nElems, FILEUNIT) call mesh_marc_count_elementSets(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 mesh_marc_map_elementSets(mesh_nameElemSet,mesh_mapElemSet,FILEUNIT) allocate (mesh_mapFEtoCPelem(2,mesh_nElems), source = 0) - call mesh_marc_map_elements(hypoelasticTableStyle,mesh_nameElemSet,mesh_mapElemSet,mesh_nElems,fileFormatVersion,FILEUNIT) + call mesh_marc_map_elements(hypoelasticTableStyle,mesh_nameElemSet,mesh_mapElemSet,& + mesh_nElems,fileFormatVersion,marc_matNumber,FILEUNIT) allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes),source=0) call mesh_marc_map_nodes(mesh_Nnodes,FILEUNIT) !ToDo: don't work on global variables @@ -418,9 +420,10 @@ subroutine mesh_marc_map_elementSets(nameElemSet,mapElemSet,fileUnit) !-------------------------------------------------------------------------------------------------- !> @brief Maps elements from FE ID to internal (consecutive) representation. !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_map_elements(tableStyle,nameElemSet,mapElemSet,nElems,fileFormatVersion,fileUnit) +subroutine mesh_marc_map_elements(tableStyle,nameElemSet,mapElemSet,nElems,fileFormatVersion,matNumber,fileUnit) integer, intent(in) :: fileUnit,tableStyle,nElems,fileFormatVersion + integer, dimension(:), intent(in) :: matNumber character(len=64), intent(in), dimension(:) :: nameElemSet integer, dimension(:,:), intent(in) :: & mapElemSet @@ -451,7 +454,7 @@ subroutine mesh_marc_map_elements(tableStyle,nameElemSet,mapElemSet,nElems,fileF if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'connectivity') then read (fileUnit,'(A300)',END=620) line chunkPos = IO_stringPos(line) - if(any(Marc_matNumber==IO_intValue(line,chunkPos,6))) then + if(any(matNumber==IO_intValue(line,chunkPos,6))) then do read (fileUnit,'(A300)',END=620) line chunkPos = IO_stringPos(line) From 7d438d38686a1a8b24c8fc69e3bfb8e800a12fc7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 8 Oct 2019 17:37:30 +0200 Subject: [PATCH 11/48] not needed --- src/mesh_base.f90 | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/src/mesh_base.f90 b/src/mesh_base.f90 index 102f53e9e..dab7059ee 100644 --- a/src/mesh_base.f90 +++ b/src/mesh_base.f90 @@ -1,4 +1,3 @@ - !-------------------------------------------------------------------------------------------------- !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH @@ -8,14 +7,13 @@ !-------------------------------------------------------------------------------------------------- module mesh_base - use, intrinsic :: iso_c_binding use prec use element implicit none !--------------------------------------------------------------------------------------------------- -!> Properties of a the whole mesh (consisting of one type of elements) +!> Properties of a whole mesh (consisting of one type of elements) !--------------------------------------------------------------------------------------------------- type, public :: tMesh type(tElement) :: & @@ -33,11 +31,7 @@ module mesh_base elemType, & Ncells, & nIPneighbors, & - NcellNodes, & - maxElemsPerNode - integer(pInt), dimension(:), allocatable, public :: & - homogenizationAt, & - microstructureAt + NcellNodes integer(pInt), dimension(:,:), allocatable, public :: & connectivity contains @@ -47,6 +41,7 @@ module mesh_base end type tMesh contains + subroutine tMesh_base_init(self,meshType,elemType,nodes) class(tMesh) :: self From b647245e394cf897bc7fcc2d88be21eef2b9e271 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 8 Oct 2019 18:52:34 +0200 Subject: [PATCH 12/48] general polishing --- src/DAMASK_abaqus.f | 9 ++++----- src/DAMASK_marc.f90 | 3 --- src/mesh_marc.f90 | 49 +++++++++++++++++++++++---------------------- 3 files changed, 29 insertions(+), 32 deletions(-) diff --git a/src/DAMASK_abaqus.f b/src/DAMASK_abaqus.f index 514307fe8..7a43f688e 100644 --- a/src/DAMASK_abaqus.f +++ b/src/DAMASK_abaqus.f @@ -313,11 +313,10 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& call CPFEM_general(computationMode,usePingPong,dfgrd0,dfgrd1,temperature,dtime,noel,npt,stress_h,ddsdde_h) -! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 -! straight: 11, 22, 33, 12, 23, 13 -! ABAQUS explicit: 11, 22, 33, 12, 23, 13 -! ABAQUS implicit: 11, 22, 33, 12, 13, 23 -! ABAQUS implicit: 11, 22, 33, 12 +! DAMASK: 11, 22, 33, 12, 23, 13 +! ABAQUS explicit: 11, 22, 33, 12, 23, 13 +! ABAQUS implicit: 11, 22, 33, 12, 13, 23 +! ABAQUS implicit: 11, 22, 33, 12 ddsdde = ddsdde_h(1:ntens,1:ntens) stress = stress_h(1:ntens) diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index 0cc815e8e..bea952fca 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -315,9 +315,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & lastLovl = lovl ! record lovl call CPFEM_general(computationMode,usePingPong,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde) -! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 -! Marc: 11, 22, 33, 12, 23, 13 -! Marc: 11, 22, 33, 12 d = ddsdde(1:ngens,1:ngens) s = stress(1:ndi+nshear) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 9e4c515ba..73f684d85 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -44,11 +44,6 @@ module mesh mesh_node, & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!! mesh_node0 !< node x,y,z coordinates (initially!) - real(pReal), dimension(:,:,:), allocatable:: & - mesh_ipArea !< area of interface to neighboring IP (initially!) - real(pReal),dimension(:,:,:,:), allocatable :: & - mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!) - ! -------------------------------------------------------------------------------------------------- type(tMesh) :: theMesh @@ -107,9 +102,9 @@ integer, dimension(:,:), allocatable :: & integer :: & mesh_NelemSets + character(len=64), dimension(:), allocatable :: & mesh_nameElemSet - integer, dimension(:,:), allocatable :: & mesh_mapElemSet !< list of elements in elementSet integer, dimension(:,:), allocatable, target :: & @@ -147,6 +142,10 @@ subroutine mesh_init(ip,el) integer, dimension(:), allocatable :: & marc_matNumber !< array of material numbers for hypoelastic material (Marc only) logical :: myDebug + real(pReal), dimension(:,:,:), allocatable:: & + mesh_ipArea !< area of interface to neighboring IP (initially!) + real(pReal),dimension(:,:,:,:), allocatable :: & + mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!) write(6,'(/,a)') ' <<<+- mesh init -+>>>' @@ -198,7 +197,9 @@ subroutine mesh_init(ip,el) call mesh_build_ipCoordinates if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6) - call mesh_build_ipAreas + allocate(mesh_ipArea(theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems)) + allocate(mesh_ipAreaNormal(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems)) + call mesh_build_ipAreas(mesh_ipArea,mesh_ipAreaNormal) if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) call IP_neighborhood2 @@ -225,8 +226,8 @@ subroutine mesh_init(ip,el) call geometry_plastic_nonlocal_setIPvolume(IPvolume()) call geometry_plastic_nonlocal_setIPneighborhood(mesh_ipNeighborhood2) - call geometry_plastic_nonlocal_setIParea(mesh_IParea) - call geometry_plastic_nonlocal_setIPareaNormal(mesh_IPareaNormal) + call geometry_plastic_nonlocal_setIParea(mesh_ipArea) + call geometry_plastic_nonlocal_setIPareaNormal(mesh_ipAreaNormal) end subroutine mesh_init @@ -737,13 +738,14 @@ subroutine mesh_marc_buildElements(microstructureAt,homogenizationAt, & subroutine buildCells(thisMesh,elem,connectivity_elem) - class(tMesh) :: thisMesh - type(tElement) :: elem + class(tMesh) :: thisMesh + type(tElement), intent(in) :: elem integer,dimension(:,:), intent(in) :: connectivity_elem - integer,dimension(:,:), allocatable :: parentsAndWeights,candidates_global - integer,dimension(:), allocatable :: candidates_local + + integer,dimension(:), allocatable :: candidates_local + integer,dimension(:,:), allocatable :: parentsAndWeights,candidates_global, connectivity_cell_reshape integer,dimension(:,:,:), allocatable :: connectivity_cell - integer,dimension(:,:), allocatable :: connectivity_cell_reshape + real(pReal), dimension(:,:), allocatable :: nodes_new,nodes integer :: e, n, c, p, s,i,m,j,nParentNodes,nCellNode,Nelem,candidateID @@ -1166,14 +1168,13 @@ end subroutine mesh_build_ipCoordinates !-------------------------------------------------------------------------------------------------- !> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal' !-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipAreas +subroutine mesh_build_ipAreas(ipArea,ipAreaNormal) integer :: e,c,i,f,n,m + real(pReal), dimension(theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), intent(out) :: ipArea + real(pReal), dimension(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), intent(out) :: ipAreaNormal real(pReal), dimension (3,FE_maxNcellnodesPerCellface) :: nodePos, normals real(pReal), dimension(3) :: normal - - allocate(mesh_ipArea(theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) - allocate(mesh_ipAreaNormal(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) c = theMesh%elem%cellType @@ -1189,8 +1190,8 @@ subroutine mesh_build_ipAreas 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 - mesh_ipArea(f,i,e) = norm2(normal) - mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal + ipArea(f,i,e) = norm2(normal) + ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal enddo enddo @@ -1201,8 +1202,8 @@ subroutine mesh_build_ipAreas nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(n,f),i,e)) normal = math_cross(nodePos(1:3,2) - nodePos(1:3,1), & nodePos(1:3,3) - nodePos(1:3,1)) - mesh_ipArea(f,i,e) = norm2(normal) - mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal + ipArea(f,i,e) = norm2(normal) + ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal enddo enddo @@ -1221,8 +1222,8 @@ subroutine mesh_build_ipAreas * 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) - mesh_ipArea(f,i,e) = norm2(normal) - mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) + ipArea(f,i,e) = norm2(normal) + ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) enddo enddo From 90c03d94d15672ab18d8be6642a94f694b65a0a7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 8 Oct 2019 19:08:29 +0200 Subject: [PATCH 13/48] further separation for clearer structure --- src/mesh_marc.f90 | 31 +++++++++++++++---------------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 73f684d85..7105f54ac 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -64,8 +64,6 @@ integer, dimension(:,:), allocatable :: & mesh_cell2, & !< cell connectivity for each element,ip/cell mesh_cell !< cell connectivity for each element,ip/cell -! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS) -! Hence, I suggest to prefix with "FE_" integer, parameter :: & FE_Ngeomtypes = 10, & @@ -73,20 +71,6 @@ integer, dimension(:,:), allocatable :: & FE_maxNcellnodesPerCell = 8, & FE_maxNcellnodesPerCellface = 4 - integer, dimension(FE_Ngeomtypes), parameter :: FE_NmatchingNodes = & !< number of nodes that are needed for face matching in a specific type of element geometry - int([ & - 3, & ! element 6 (2D 3node 1ip) - 3, & ! element 125 (2D 6node 3ip) - 4, & ! element 11 (2D 4node 4ip) - 4, & ! element 27 (2D 8node 9ip) - 4, & ! element 134 (3D 4node 1ip) - 4, & ! element 127 (3D 10node 4ip) - 6, & ! element 136 (3D 6node 6ip) - 8, & ! element 117 (3D 8node 1ip) - 8, & ! element 7 (3D 8node 8ip) - 8 & ! element 21 (3D 20node 27ip) - ],pInt) - integer, dimension(FE_Ncelltypes), parameter :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type int([& @@ -916,6 +900,21 @@ subroutine mesh_build_cellconnectivity integer, dimension(:), allocatable :: & matchingNode2cellnode + + integer, dimension(FE_Ngeomtypes), parameter :: FE_NmatchingNodes = & !< number of nodes that are needed for face matching in a specific type of element geometry + int([ & + 3, & ! element 6 (2D 3node 1ip) + 3, & ! element 125 (2D 6node 3ip) + 4, & ! element 11 (2D 4node 4ip) + 4, & ! element 27 (2D 8node 9ip) + 4, & ! element 134 (3D 4node 1ip) + 4, & ! element 127 (3D 10node 4ip) + 6, & ! element 136 (3D 6node 6ip) + 8, & ! element 117 (3D 8node 1ip) + 8, & ! element 7 (3D 8node 8ip) + 8 & ! element 21 (3D 20node 27ip) + ],pInt) + integer, dimension(:,:), allocatable :: & cellnodeParent integer, dimension(theMesh%elem%Ncellnodes) :: & From df3bb7c5599ebd34294478dcd0674a1a7d22626d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 8 Oct 2019 23:53:16 +0200 Subject: [PATCH 14/48] cleaning --- src/mesh_marc.f90 | 15 ++++----------- 1 file changed, 4 insertions(+), 11 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 7105f54ac..6bd3453bc 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -65,14 +65,7 @@ integer, dimension(:,:), allocatable :: & mesh_cell !< cell connectivity for each element,ip/cell - integer, parameter :: & - FE_Ngeomtypes = 10, & - FE_Ncelltypes = 4, & - FE_maxNcellnodesPerCell = 8, & - FE_maxNcellnodesPerCellface = 4 - - - integer, dimension(FE_Ncelltypes), parameter :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type + integer, dimension(4), parameter :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type int([& 3, & ! (2D 3node) 4, & ! (2D 4node) @@ -901,7 +894,7 @@ subroutine mesh_build_cellconnectivity integer, dimension(:), allocatable :: & matchingNode2cellnode - integer, dimension(FE_Ngeomtypes), parameter :: FE_NmatchingNodes = & !< number of nodes that are needed for face matching in a specific type of element geometry + integer, dimension(10), parameter :: FE_NmatchingNodes = & !< number of nodes that are needed for face matching in a specific type of element geometry int([ & 3, & ! element 6 (2D 3node 1ip) 3, & ! element 125 (2D 6node 3ip) @@ -924,7 +917,7 @@ subroutine mesh_build_cellconnectivity matchingNodeID, & localCellnodeID - allocate(mesh_cell(FE_maxNcellnodesPerCell,theMesh%elem%nIPs,theMesh%nElems), source=0) + allocate(mesh_cell(theMesh%elem%NcellNodesPerCell,theMesh%elem%nIPs,theMesh%nElems), source=0) allocate(matchingNode2cellnode(theMesh%nNodes), source=0) allocate(cellnodeParent(2,theMesh%elem%Ncellnodes*theMesh%nElems), source=0) @@ -1172,7 +1165,7 @@ subroutine mesh_build_ipAreas(ipArea,ipAreaNormal) integer :: e,c,i,f,n,m real(pReal), dimension(theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), intent(out) :: ipArea real(pReal), dimension(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), intent(out) :: ipAreaNormal - real(pReal), dimension (3,FE_maxNcellnodesPerCellface) :: nodePos, normals + real(pReal), dimension (3,size(theMesh%elem%cellFace,2)) :: nodePos, normals real(pReal), dimension(3) :: normal c = theMesh%elem%cellType From a2a05158f2964df2348a96bbd2b610e451543e5b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 9 Oct 2019 01:32:19 +0200 Subject: [PATCH 15/48] not existing anymore --- src/mesh_FEM.f90 | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/mesh_FEM.f90 b/src/mesh_FEM.f90 index 2fd63c0b3..8c0b69ca3 100644 --- a/src/mesh_FEM.f90 +++ b/src/mesh_FEM.f90 @@ -221,9 +221,6 @@ subroutine mesh_init call theMesh%init(dimplex,integrationOrder,mesh_node0) call theMesh%setNelems(mesh_NcpElems) - theMesh%homogenizationAt = mesh_element(3,:) - theMesh%microstructureAt = mesh_element(4,:) - call discretization_init(mesh_element(3,:),mesh_element(4,:),& reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), & mesh_node0) From 513f1c6726cf93e4424e6adaccd5dc325eb4b436 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 12 Oct 2019 12:34:37 +0200 Subject: [PATCH 16/48] write out point data might be possible to use a rectilinear grid for this also in the case of grid solvers --- processing/post/DADF5_vtk_cells.py | 2 +- processing/post/DADF5_vtk_points.py | 126 ++++++++++++++++++++++++++++ 2 files changed, 127 insertions(+), 1 deletion(-) create mode 100755 processing/post/DADF5_vtk_points.py diff --git a/processing/post/DADF5_vtk_cells.py b/processing/post/DADF5_vtk_cells.py index 6d1e8e340..473e887f8 100755 --- a/processing/post/DADF5_vtk_cells.py +++ b/processing/post/DADF5_vtk_cells.py @@ -40,7 +40,7 @@ if options.con is None: options.con=[] for filename in options.filenames: results = damask.DADF5(filename) - if results.structured: # for grid solvers use rectilinear grid + if results.structured: # for grid solvers use rectilinear grid rGrid = vtk.vtkRectilinearGrid() coordArray = [vtk.vtkDoubleArray(), vtk.vtkDoubleArray(), diff --git a/processing/post/DADF5_vtk_points.py b/processing/post/DADF5_vtk_points.py new file mode 100755 index 000000000..ee71f787d --- /dev/null +++ b/processing/post/DADF5_vtk_points.py @@ -0,0 +1,126 @@ +#!/usr/bin/env python3 + +import os +import argparse + +import numpy as np +import vtk +from vtk.util import numpy_support + +import damask + +scriptName = os.path.splitext(os.path.basename(__file__))[0] +scriptID = ' '.join([scriptName,damask.version]) + +# -------------------------------------------------------------------- +# MAIN +# -------------------------------------------------------------------- +parser = argparse.ArgumentParser() + +#ToDo: We need to decide on a way of handling arguments of variable lentght +#https://stackoverflow.com/questions/15459997/passing-integer-lists-to-python + +#parser.add_argument('--version', action='version', version='%(prog)s {}'.format(scriptID)) +parser.add_argument('filenames', nargs='+', + help='DADF5 files') +parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string', + help='name of subdirectory relative to the location of the DADF5 file to hold output') +parser.add_argument('--mat', nargs='+', + help='labels for materialpoint',dest='mat') +parser.add_argument('--con', nargs='+', + help='labels for constituent',dest='con') + +options = parser.parse_args() + +if options.mat is None: options.mat=[] +if options.con is None: options.con=[] + +# --- loop over input files ------------------------------------------------------------------------ + +for filename in options.filenames: + results = damask.DADF5(filename) + + Polydata = vtk.vtkPolyData() + Points = vtk.vtkPoints() + Vertices = vtk.vtkCellArray() + + if results.structured: # for grid solvers calculate points + delta = results.size/results.grid*0.5 + for z in np.linspace(delta[2],results.size[2]-delta[2],results.grid[2]): + for y in np.linspace(delta[1],results.size[1]-delta[1],results.grid[1]): + for x in np.linspace(delta[0],results.size[0]-delta[0],results.grid[0]): + pointID = Points.InsertNextPoint([x,y,z]) + Vertices.InsertNextCell(1) + Vertices.InsertCellPoint(pointID) + + Polydata.SetPoints(Points) + Polydata.SetVerts(Vertices) + Polydata.Modified() + + for i,inc in enumerate(results.iter_visible('increments')): + print('Output step {}/{}'.format(i+1,len(results.increments))) + vtk_data = [] + + results.set_visible('materialpoints',False) + results.set_visible('constituents', True) + for label in options.con: + + for p in results.iter_visible('con_physics'): + if p != 'generic': + for c in results.iter_visible('constituents'): + x = results.get_dataset_location(label) + if len(x) == 0: + continue + array = results.read_dataset(x,0) + shape = [array.shape[0],np.product(array.shape[1:])] + vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) + vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) + Polydata.GetCellData().AddArray(vtk_data[-1]) + else: + x = results.get_dataset_location(label) + if len(x) == 0: + continue + array = results.read_dataset(x,0) + shape = [array.shape[0],np.product(array.shape[1:])] + vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) + vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) + Polydata.GetCellData().AddArray(vtk_data[-1]) + + results.set_visible('constituents', False) + results.set_visible('materialpoints',True) + for label in options.mat: + for p in results.iter_visible('mat_physics'): + if p != 'generic': + for m in results.iter_visible('materialpoints'): + x = results.get_dataset_location(label) + if len(x) == 0: + continue + array = results.read_dataset(x,0) + shape = [array.shape[0],np.product(array.shape[1:])] + vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) + vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) + Polydata.GetCellData().AddArray(vtk_data[-1]) + else: + x = results.get_dataset_location(label) + if len(x) == 0: + continue + array = results.read_dataset(x,0) + shape = [array.shape[0],np.product(array.shape[1:])] + vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) + vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) + Polydata.GetCellData().AddArray(vtk_data[-1]) + + writer = vtk.vtkXMLPolyDataWriter() + + + dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir)) + if not os.path.isdir(dirname): + os.mkdir(dirname,0o755) + file_out = '{}_{}.{}'.format(os.path.splitext(os.path.split(filename)[-1])[0],inc,writer.GetDefaultFileExtension()) + + writer.SetCompressorTypeToZLib() + writer.SetDataModeToBinary() + writer.SetFileName(os.path.join(dirname,file_out)) + writer.SetInputData(Polydata) + + writer.Write() From 734e6ef15f0cff7d796258a8566bfebde67bd2c7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 12 Oct 2019 15:03:26 +0200 Subject: [PATCH 17/48] writing initial coordinates to DADF5 --- src/mesh_marc.f90 | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 6bd3453bc..1b6cb5586 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -123,6 +123,9 @@ subroutine mesh_init(ip,el) mesh_ipArea !< area of interface to neighboring IP (initially!) real(pReal),dimension(:,:,:,:), allocatable :: & mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!) + + real(pReal), dimension(:,:), allocatable :: & + ip_reshaped write(6,'(/,a)') ' <<<+- mesh init -+>>>' @@ -197,10 +200,18 @@ subroutine mesh_init(ip,el) 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" + ip_reshaped = reshape(mesh_ipCoordinates,[3,theMesh%elem%nIPs*theMesh%nElems]) call discretization_init(microstructureAt,homogenizationAt,& - reshape(mesh_ipCoordinates,[3,theMesh%elem%nIPs*theMesh%nElems]),& + ip_reshaped,& mesh_node0) + call results_openJobFile + call HDF5_closeGroup(results_addGroup('geometry')) + call results_writeDataset('geometry',ip_reshaped,'x_c', & + 'cell center coordinates','m') + call results_writeDataset('geometry',mesh_node0,'x_n', & + 'nodal coordinates','m') + call results_closeJobFile() call geometry_plastic_nonlocal_setIPvolume(IPvolume()) call geometry_plastic_nonlocal_setIPneighborhood(mesh_ipNeighborhood2) call geometry_plastic_nonlocal_setIParea(mesh_ipArea) @@ -667,7 +678,6 @@ subroutine mesh_marc_buildElements(microstructureAt,homogenizationAt, & #if defined(DAMASK_HDF5) call results_openJobFile - call HDF5_closeGroup(results_addGroup('geometry')) call results_writeDataset('geometry',mesh_FEnodes,'C',& 'connectivity of the elements','-') call results_closeJobFile From 8cebf8a10d4bae84ca1e5f72ffcbd3365a9f1f2d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 12 Oct 2019 15:50:10 +0200 Subject: [PATCH 18/48] needed --- src/CPFEM.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 6649de542..32098b109 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -381,6 +381,7 @@ subroutine CPFEM_results(inc,time) call results_addIncrement(inc,time) call constitutive_results call crystallite_results + call homogenization_results call results_removeLink('current') ! ToDo: put this into closeJobFile call results_closeJobFile #endif From 708bbd3cb92b09810bbc995cf23c4f5891608ff6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 12 Oct 2019 16:15:04 +0200 Subject: [PATCH 19/48] mesh/grid type transparent handling of coordinates --- processing/post/DADF5_vtk_points.py | 17 ++++++----------- python/damask/dadf5.py | 14 ++++++++++++++ 2 files changed, 20 insertions(+), 11 deletions(-) diff --git a/processing/post/DADF5_vtk_points.py b/processing/post/DADF5_vtk_points.py index ee71f787d..87c1ad93e 100755 --- a/processing/post/DADF5_vtk_points.py +++ b/processing/post/DADF5_vtk_points.py @@ -40,19 +40,14 @@ if options.con is None: options.con=[] for filename in options.filenames: results = damask.DADF5(filename) - Polydata = vtk.vtkPolyData() Points = vtk.vtkPoints() - Vertices = vtk.vtkCellArray() - - if results.structured: # for grid solvers calculate points - delta = results.size/results.grid*0.5 - for z in np.linspace(delta[2],results.size[2]-delta[2],results.grid[2]): - for y in np.linspace(delta[1],results.size[1]-delta[1],results.grid[1]): - for x in np.linspace(delta[0],results.size[0]-delta[0],results.grid[0]): - pointID = Points.InsertNextPoint([x,y,z]) - Vertices.InsertNextCell(1) - Vertices.InsertCellPoint(pointID) + Vertices = vtk.vtkCellArray() + for c in results.cell_coordinates(): + pointID = Points.InsertNextPoint(c) + Vertices.InsertNextCell(1) + Vertices.InsertCellPoint(pointID) + Polydata = vtk.vtkPolyData() Polydata.SetPoints(Points) Polydata.SetVerts(Vertices) Polydata.Modified() diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 8444fa259..955765fcc 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -335,6 +335,20 @@ class DADF5(): return dataset + def cell_coordinates(self): + """Initial coordinates of the cell centers.""" + if self.structured: + delta = self.size/self.grid*0.5 + z, y, x = np.meshgrid(np.linspace(delta[2],self.size[2]-delta[2],self.grid[2]), + np.linspace(delta[1],self.size[1]-delta[1],self.grid[1]), + np.linspace(delta[0],self.size[0]-delta[0],self.grid[0]), + ) + return np.concatenate((x[:,:,:,None],y[:,:,:,None],y[:,:,:,None]),axis = 3).reshape([np.product(self.grid),3]) + else: + with h5py.File(self.filename,'r') as f: + return f['geometry/x_c'][()] + + def add_Cauchy(self,P='P',F='F'): """ Adds Cauchy stress calculated from 1st Piola-Kirchhoff stress and deformation gradient. From fabab089366bfa44327dc8e3f1c0305ad8d41d63 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 12 Oct 2019 16:25:59 +0200 Subject: [PATCH 20/48] did not run without HDF5 support --- src/mesh_marc.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 1b6cb5586..1f4738bec 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -204,7 +204,7 @@ subroutine mesh_init(ip,el) call discretization_init(microstructureAt,homogenizationAt,& ip_reshaped,& mesh_node0) - +#if defined(DAMASK_HDF5) call results_openJobFile call HDF5_closeGroup(results_addGroup('geometry')) call results_writeDataset('geometry',ip_reshaped,'x_c', & @@ -212,6 +212,7 @@ subroutine mesh_init(ip,el) call results_writeDataset('geometry',mesh_node0,'x_n', & 'nodal coordinates','m') call results_closeJobFile() +#endif call geometry_plastic_nonlocal_setIPvolume(IPvolume()) call geometry_plastic_nonlocal_setIPneighborhood(mesh_ipNeighborhood2) call geometry_plastic_nonlocal_setIParea(mesh_ipArea) From d8d99f3694ba25244875e82297173608897a8f4a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 12 Oct 2019 16:39:44 +0200 Subject: [PATCH 21/48] geometry folder needs to be created earlier --- src/mesh_marc.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 1f4738bec..efa4e43d8 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -206,7 +206,6 @@ subroutine mesh_init(ip,el) mesh_node0) #if defined(DAMASK_HDF5) call results_openJobFile - call HDF5_closeGroup(results_addGroup('geometry')) call results_writeDataset('geometry',ip_reshaped,'x_c', & 'cell center coordinates','m') call results_writeDataset('geometry',mesh_node0,'x_n', & @@ -679,6 +678,7 @@ subroutine mesh_marc_buildElements(microstructureAt,homogenizationAt, & #if defined(DAMASK_HDF5) call results_openJobFile + call HDF5_closeGroup(results_addGroup('geometry')) call results_writeDataset('geometry',mesh_FEnodes,'C',& 'connectivity of the elements','-') call results_closeJobFile From e6d25bfdab9e271bbda6db717418a1fc65c1dbd7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 12 Oct 2019 19:24:03 +0200 Subject: [PATCH 22/48] almost no shared functionality --- src/mesh_marc.f90 | 56 +++++++++++++++++++++++++++++++++-------------- 1 file changed, 40 insertions(+), 16 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index efa4e43d8..7b58d79d6 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -89,8 +89,7 @@ integer, dimension(:,:), allocatable :: & mesh_mapFEtoCPnode !< [sorted FEid, corresponding CPid] integer, dimension(:,:,:,:), allocatable :: & - mesh_ipNeighborhood2 !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] - + mesh_ipNeighborhood2 !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] public :: & @@ -163,7 +162,8 @@ subroutine mesh_init(ip,el) allocate(homogenizationAt(theMesh%nElems), source=0) allocate(mesh_FEnodes(theMesh%elem%nNodes,theMesh%nElems), source=0) - call mesh_marc_buildElements(microstructureAt,homogenizationAt, & + call mesh_marc_buildElements(mesh_nElems,theMesh%elem%nNodes,FILEUNIT) + call mesh_marc_buildElements2(microstructureAt,homogenizationAt, & mesh_nElems,theMesh%elem%nNodes,initialcondTableStyle,FILEUNIT) if (myDebug) write(6,'(a)') ' Built elements'; flush(6) close (FILEUNIT) @@ -625,18 +625,13 @@ end function mapElemtype !-------------------------------------------------------------------------------------------------- -!> @brief Stores node IDs and homogenization and microstructure ID +!> @brief Stores node IDs !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_buildElements(microstructureAt,homogenizationAt, & - nElem,nNodes,initialcondTableStyle,fileUnit) +subroutine mesh_marc_buildElements(nElem,nNodes,fileUnit) - integer, dimension(:), intent(out) :: & - microstructureAt, & - homogenizationAt integer, intent(in) :: & nElem, & nNodes, & !< number of nodes per element - initialcondTableStyle, & fileUnit integer, allocatable, dimension(:) :: chunkPos @@ -658,7 +653,7 @@ subroutine mesh_marc_buildElements(microstructureAt,homogenizationAt, & if (e /= 0) then ! disregard non CP elems nNodesAlreadyRead = 0 do j = 1,chunkPos(1)-2 - mesh_FEnodes(j,e) = mesh_FEasCP('node',IO_IntValue(line,chunkPos,j+2)) ! CP ids of nodes + mesh_FEnodes(j,e) = mesh_FEasCP('node',IO_IntValue(line,chunkPos,j+2)) ! CP ids of nodes enddo nNodesAlreadyRead = chunkPos(1) - 2 do while(nNodesAlreadyRead < nNodes) ! read on if not all nodes in one line @@ -674,7 +669,7 @@ subroutine mesh_marc_buildElements(microstructureAt,homogenizationAt, & exit endif enddo -620 rewind(fileUnit) ! just in case "initial state" appears before "connectivity" +620 rewind(fileUnit) #if defined(DAMASK_HDF5) call results_openJobFile @@ -686,6 +681,30 @@ subroutine mesh_marc_buildElements(microstructureAt,homogenizationAt, & call buildCells(theMesh,theMesh%elem,mesh_FEnodes) +end subroutine mesh_marc_buildElements + + +!-------------------------------------------------------------------------------------------------- +!> @brief Stores homogenization and microstructure ID +!-------------------------------------------------------------------------------------------------- +subroutine mesh_marc_buildElements2(microstructureAt,homogenizationAt, & + nElem,nNodes,initialcondTableStyle,fileUnit) + + integer, dimension(:), intent(out) :: & + microstructureAt, & + homogenizationAt + integer, intent(in) :: & + nElem, & + nNodes, & !< number of nodes per element + initialcondTableStyle, & + fileUnit + + integer, allocatable, dimension(:) :: chunkPos + character(len=300) line + + integer, dimension(1+nElem) :: contInts + integer :: i,j,t,sv,myVal,e,nNodesAlreadyRead + read (fileUnit,'(A300)',END=630) line do chunkPos = IO_stringPos(line) @@ -721,7 +740,7 @@ subroutine mesh_marc_buildElements(microstructureAt,homogenizationAt, & endif enddo -630 end subroutine mesh_marc_buildElements +630 end subroutine mesh_marc_buildElements2 subroutine buildCells(thisMesh,elem,connectivity_elem) @@ -1006,7 +1025,7 @@ function mesh_build_cellnodes() end function mesh_build_cellnodes -!-------------------------------------------------------------------------------------------------- +!--------------------------------------------------------------------------------------------------- !> @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. @@ -1014,7 +1033,7 @@ end function mesh_build_cellnodes !> 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. -!-------------------------------------------------------------------------------------------------- +!--------------------------------------------------------------------------------------------------- function IPvolume() real(pReal), dimension(theMesh%elem%nIPs,theMesh%nElems) :: IPvolume @@ -1068,12 +1087,16 @@ function IPvolume() end function IPvolume +!--------------------------------------------------------------------------------------------------- +!> @brief cell neighborhood +!--------------------------------------------------------------------------------------------------- subroutine IP_neighborhood2 integer, dimension(:,:), allocatable :: faces - integer, dimension(:), allocatable :: face + integer, dimension(:), allocatable :: face integer :: e,i,f,c,m,n,j,k,l,p, current, next,i2,e2,n2,k2 logical :: match + allocate(faces(size(theMesh%elem%cellface,1)+3,size(theMesh%elem%cellface,2)*theMesh%elem%nIPs*theMesh%Nelems)) ! store cell face definitions @@ -1140,6 +1163,7 @@ subroutine IP_neighborhood2 end subroutine IP_neighborhood2 + !-------------------------------------------------------------------------------------------------- !> @brief Calculates IP Coordinates. ! Marc however only provides nodal displacements, From 4c9bc326f642b3648d241a79d8b9519c6c03225a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 12 Oct 2019 21:58:26 +0200 Subject: [PATCH 23/48] better read file only once --- src/mesh_marc.f90 | 55 +++++++++++++++++++++-------------------------- 1 file changed, 25 insertions(+), 30 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 7b58d79d6..bb9526fd5 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -110,6 +110,8 @@ 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 + integer :: j, fileFormatVersion, elemType, & mesh_maxNelemInSet, & mesh_nElems, & @@ -129,13 +131,13 @@ subroutine mesh_init(ip,el) 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)) myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0) ! parsing Marc input file call IO_open_inputFile(FILEUNIT,modelName) - fileFormatVersion = mesh_marc_get_fileFormat(FILEUNIT) - call mesh_marc_get_tableStyles(initialcondTableStyle,hypoelasticTableStyle,FILEUNIT) + fileFormatVersion = mesh_marc_get_fileFormat(inputFile) + call mesh_marc_get_tableStyles(initialcondTableStyle,hypoelasticTableStyle,inputFile) if (fileFormatVersion > 12) & marc_matNumber = mesh_marc_get_matNumber(FILEUNIT,hypoelasticTableStyle) call mesh_marc_count_nodesAndElements(mesh_nNodes, mesh_nElems, FILEUNIT) @@ -223,55 +225,48 @@ end subroutine mesh_init !-------------------------------------------------------------------------------------------------- !> @brief Figures out version of Marc input file format !-------------------------------------------------------------------------------------------------- -integer function mesh_marc_get_fileFormat(fileUnit) +integer function mesh_marc_get_fileFormat(fileContent) - integer, intent(in) :: fileUnit - + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines + integer, allocatable, dimension(:) :: chunkPos - character(len=300) line - - - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'version') then - mesh_marc_get_fileFormat = IO_intValue(line,chunkPos,2) + integer :: l + + do l = 1, size(fileContent) + chunkPos = IO_stringPos(fileContent(l)) + if ( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'version') then + mesh_marc_get_fileFormat = IO_intValue(fileContent(l),chunkPos,2) exit endif enddo -620 end function mesh_marc_get_fileFormat +end function mesh_marc_get_fileFormat !-------------------------------------------------------------------------------------------------- !> @brief Figures out table styles for initial cond and hypoelastic !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_get_tableStyles(initialcond,hypoelastic,fileUnit) +subroutine mesh_marc_get_tableStyles(initialcond,hypoelastic,fileContent) integer, intent(out) :: initialcond, hypoelastic - integer, intent(in) :: fileUnit + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines integer, allocatable, dimension(:) :: chunkPos - character(len=300) line - + integer :: l + initialcond = 0 hypoelastic = 0 - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'table' .and. chunkPos(1) > 5) then - initialcond = IO_intValue(line,chunkPos,4) - hypoelastic = IO_intValue(line,chunkPos,5) + do l = 1, size(fileContent) + chunkPos = IO_stringPos(fileContent(l)) + if ( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'table' .and. chunkPos(1) > 5) then + initialcond = IO_intValue(fileContent(l),chunkPos,4) + hypoelastic = IO_intValue(fileContent(l),chunkPos,5) exit endif enddo -620 end subroutine mesh_marc_get_tableStyles +end subroutine mesh_marc_get_tableStyles !-------------------------------------------------------------------------------------------------- From 7611513bb87689e07bf8391d0e5d1593c466c7f4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 12 Oct 2019 22:18:55 +0200 Subject: [PATCH 24/48] strange indentation --- src/IO.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index 066ebbe8c..64ce6fc9a 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -1072,8 +1072,8 @@ integer function IO_countContinuousIntValues(fileUnit) if (chunkPos(1) < 1) then ! empty line exit elseif (IO_lc(IO_stringValue(line,chunkPos,2)) == 'to' ) then ! found range indicator - IO_countContinuousIntValues = 1 + abs( IO_intValue(line,chunkPos,3) & - - IO_intValue(line,chunkPos,1)) + IO_countContinuousIntValues = 1 + abs( IO_intValue(line,chunkPos,3) & + -IO_intValue(line,chunkPos,1)) exit ! only one single range indicator allowed else IO_countContinuousIntValues = IO_countContinuousIntValues+chunkPos(1)-1 ! add line's count when assuming 'c' From d6b6096007cc14d54ba797e9566d8adf3e0bb557 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 12 Oct 2019 22:40:00 +0200 Subject: [PATCH 25/48] avoid global variables --- src/mesh_marc.f90 | 35 +++++++++++++++++++---------------- 1 file changed, 19 insertions(+), 16 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index bb9526fd5..0aab1b106 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -162,13 +162,23 @@ subroutine mesh_init(ip,el) allocate(microstructureAt(theMesh%nElems), source=0) allocate(homogenizationAt(theMesh%nElems), source=0) - allocate(mesh_FEnodes(theMesh%elem%nNodes,theMesh%nElems), source=0) - call mesh_marc_buildElements(mesh_nElems,theMesh%elem%nNodes,FILEUNIT) + mesh_FEnodes = mesh_marc_buildElements(theMesh%nElems,theMesh%elem%nNodes,FILEUNIT) call mesh_marc_buildElements2(microstructureAt,homogenizationAt, & mesh_nElems,theMesh%elem%nNodes,initialcondTableStyle,FILEUNIT) if (myDebug) write(6,'(a)') ' Built elements'; flush(6) close (FILEUNIT) + + +#if defined(DAMASK_HDF5) + call results_openJobFile + call HDF5_closeGroup(results_addGroup('geometry')) + call results_writeDataset('geometry',mesh_FEnodes,'C',& + 'connectivity of the elements','-') + call results_closeJobFile +#endif + + call buildCells(theMesh,theMesh%elem,mesh_FEnodes) call mesh_build_cellconnectivity if (myDebug) write(6,'(a)') ' Built cell connectivity'; flush(6) @@ -622,12 +632,15 @@ end function mapElemtype !-------------------------------------------------------------------------------------------------- !> @brief Stores node IDs !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_buildElements(nElem,nNodes,fileUnit) +function mesh_marc_buildElements(nElem,nNodes,fileUnit) integer, intent(in) :: & nElem, & nNodes, & !< number of nodes per element fileUnit + + integer, dimension(nElem,nNodes) :: & + mesh_marc_buildElements integer, allocatable, dimension(:) :: chunkPos character(len=300) line @@ -648,14 +661,14 @@ subroutine mesh_marc_buildElements(nElem,nNodes,fileUnit) if (e /= 0) then ! disregard non CP elems nNodesAlreadyRead = 0 do j = 1,chunkPos(1)-2 - mesh_FEnodes(j,e) = mesh_FEasCP('node',IO_IntValue(line,chunkPos,j+2)) ! CP ids of nodes + mesh_marc_buildElements(j,e) = mesh_FEasCP('node',IO_IntValue(line,chunkPos,j+2)) ! CP ids of nodes 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) - mesh_FEnodes(nNodesAlreadyRead+j,e) = mesh_FEasCP('node',IO_IntValue(line,chunkPos,j)) ! CP ids of nodes + mesh_marc_buildElements(nNodesAlreadyRead+j,e) = mesh_FEasCP('node',IO_IntValue(line,chunkPos,j)) ! CP ids of nodes enddo nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1) enddo @@ -665,18 +678,8 @@ subroutine mesh_marc_buildElements(nElem,nNodes,fileUnit) endif enddo 620 rewind(fileUnit) - -#if defined(DAMASK_HDF5) - call results_openJobFile - call HDF5_closeGroup(results_addGroup('geometry')) - call results_writeDataset('geometry',mesh_FEnodes,'C',& - 'connectivity of the elements','-') - call results_closeJobFile -#endif - call buildCells(theMesh,theMesh%elem,mesh_FEnodes) - -end subroutine mesh_marc_buildElements +end function mesh_marc_buildElements !-------------------------------------------------------------------------------------------------- From 680ed535d79cc3b3449151e44406b8ac40a3a32c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Oct 2019 10:42:34 +0200 Subject: [PATCH 26/48] avoid file operations and line labels --- src/mesh_marc.f90 | 55 +++++++++++++++++++---------------------------- 1 file changed, 22 insertions(+), 33 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 0aab1b106..950ab671a 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -149,9 +149,9 @@ subroutine mesh_init(ip,el) call mesh_marc_map_elements(hypoelasticTableStyle,mesh_nameElemSet,mesh_mapElemSet,& mesh_nElems,fileFormatVersion,marc_matNumber,FILEUNIT) allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes),source=0) - call mesh_marc_map_nodes(mesh_Nnodes,FILEUNIT) !ToDo: don't work on global variables + call mesh_marc_map_nodes(mesh_Nnodes,inputFile) !ToDo: don't work on global variables - mesh_node0 = mesh_marc_build_nodes(mesh_Nnodes,FILEUNIT) + mesh_node0 = mesh_marc_build_nodes(mesh_Nnodes,inputFile) mesh_node = mesh_node0 elemType = mesh_marc_getElemType(mesh_nElems,FILEUNIT) @@ -478,34 +478,26 @@ end subroutine mesh_marc_map_elements !-------------------------------------------------------------------------------------------------- !> @brief Maps node from FE ID to internal (consecutive) representation. !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_map_nodes(nNodes,fileUnit) +subroutine mesh_marc_map_nodes(nNodes,fileContent) - integer, intent(in) :: fileUnit, nNodes + integer, intent(in) :: nNodes + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines integer, allocatable, dimension(:) :: chunkPos - character(len=300) line + integer :: i, l - integer, dimension (nNodes) :: node_count - integer :: i - - node_count = 0 - - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'coordinates' ) then - read (fileUnit,'(A300)') line ! skip crap line + 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 - read (fileUnit,'(A300)') line - mesh_mapFEtoCPnode(1,i) = IO_fixedIntValue (line,[0,10],1) + mesh_mapFEtoCPnode(1,i) = IO_fixedIntValue (fileContent(l+1+i),[0,10],1) mesh_mapFEtoCPnode(2,i) = i enddo exit endif enddo -620 call math_sort(mesh_mapFEtoCPnode) + call math_sort(mesh_mapFEtoCPnode) end subroutine mesh_marc_map_nodes @@ -513,33 +505,30 @@ end subroutine mesh_marc_map_nodes !-------------------------------------------------------------------------------------------------- !> @brief store x,y,z coordinates of all nodes in mesh. !-------------------------------------------------------------------------------------------------- -function mesh_marc_build_nodes(nNode,fileUnit) result(nodes) +function mesh_marc_build_nodes(nNode,fileContent) result(nodes) - integer, intent(in) :: nNode,fileUnit + 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 - character(len=300) :: line - integer :: i,j,m + integer :: i,j,m,l - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'coordinates' ) then - read (fileUnit,'(A300)') line ! skip crap line + do l = 1, size(fileContent) + chunkPos = IO_stringPos(fileContent(l)) + if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then do i=1,nNode - read (fileUnit,'(A300)') line - m = mesh_FEasCP('node',IO_fixedIntValue(line,node_ends,1)) + m = mesh_FEasCP('node',IO_fixedIntValue(fileContent(l+1+i),node_ends,1)) do j = 1,3 - nodes(j,m) = mesh_unitlength * IO_fixedNoEFloatValue(line,node_ends,j+1) + nodes(j,m) = mesh_unitlength * IO_fixedNoEFloatValue(fileContent(l+1+i),node_ends,j+1) enddo enddo exit endif enddo -620 end function mesh_marc_build_nodes +end function mesh_marc_build_nodes !-------------------------------------------------------------------------------------------------- From 9ea91b84e819dd576f2ba0276f726d0cddbf1a9c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Oct 2019 11:36:21 +0200 Subject: [PATCH 27/48] easier to read and without file access --- src/mesh_marc.f90 | 93 +++++++++++++++++++++-------------------------- 1 file changed, 42 insertions(+), 51 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 950ab671a..4c854a2f8 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -135,12 +135,13 @@ subroutine mesh_init(ip,el) myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0) ! parsing Marc input file - call IO_open_inputFile(FILEUNIT,modelName) fileFormatVersion = mesh_marc_get_fileFormat(inputFile) call mesh_marc_get_tableStyles(initialcondTableStyle,hypoelasticTableStyle,inputFile) if (fileFormatVersion > 12) & - marc_matNumber = mesh_marc_get_matNumber(FILEUNIT,hypoelasticTableStyle) - call mesh_marc_count_nodesAndElements(mesh_nNodes, mesh_nElems, FILEUNIT) + marc_matNumber = mesh_marc_get_matNumber(hypoelasticTableStyle,inputFile) + call mesh_marc_count_nodesAndElements(mesh_nNodes, mesh_nElems, inputFile) + + call IO_open_inputFile(FILEUNIT,modelName) call mesh_marc_count_elementSets(mesh_NelemSets,mesh_maxNelemInSet,FILEUNIT) allocate(mesh_nameElemSet(mesh_NelemSets)); mesh_nameElemSet = 'n/a' allocate(mesh_mapElemSet(1+mesh_maxNelemInSet,mesh_NelemSets),source=0) @@ -258,8 +259,8 @@ end function mesh_marc_get_fileFormat !-------------------------------------------------------------------------------------------------- subroutine mesh_marc_get_tableStyles(initialcond,hypoelastic,fileContent) - integer, intent(out) :: initialcond, hypoelastic - character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines + integer, intent(out) :: initialcond, hypoelastic + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines integer, allocatable, dimension(:) :: chunkPos integer :: l @@ -282,73 +283,63 @@ end subroutine mesh_marc_get_tableStyles !-------------------------------------------------------------------------------------------------- !> @brief Figures out material number of hypoelastic material !-------------------------------------------------------------------------------------------------- -function mesh_marc_get_matNumber(fileUnit,tableStyle) +function mesh_marc_get_matNumber(tableStyle,fileContent) - integer, intent(in) :: fileUnit, tableStyle + integer, intent(in) :: tableStyle + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines + integer, dimension(:), allocatable :: mesh_marc_get_matNumber integer, allocatable, dimension(:) :: chunkPos - integer :: i, j, data_blocks - character(len=300) :: line + integer :: i, j, data_blocks, l - data_blocks = 1 - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'hypoelastic') then - read (fileUnit,'(A300)',END=620) line - if (len(trim(line))/=0) then - chunkPos = IO_stringPos(line) - data_blocks = IO_intValue(line,chunkPos,1) - endif - allocate(mesh_marc_get_matNumber(data_blocks), source = 0) - do i=1,data_blocks ! read all data blocks - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - mesh_marc_get_matNumber(i) = IO_intValue(line,chunkPos,1) - do j=1,2 + tableStyle ! read 2 or 3 remaining lines of data block - read (fileUnit,'(A300)') line - enddo - enddo - exit - endif + do l = 1, size(fileContent) + chunkPos = IO_stringPos(fileContent(l)) + if ( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'hypoelastic') then + if (len(trim(fileContent(l+1)))/=0) then + chunkPos = IO_stringPos(fileContent(l+1)) + data_blocks = IO_intValue(fileContent(l+1),chunkPos,1) + else + data_blocks = 1 + endif + allocate(mesh_marc_get_matNumber(data_blocks), source = 0) + do i = 0, data_blocks - 1 + j = i*(2+tableStyle) + 1 + chunkPos = IO_stringPos(fileContent(l+1+j)) + mesh_marc_get_matNumber(i+1) = IO_intValue(fileContent(l+1+j),chunkPos,1) + enddo + exit + endif enddo -620 end function mesh_marc_get_matNumber +end function mesh_marc_get_matNumber !-------------------------------------------------------------------------------------------------- !> @brief Count overall number of nodes and elements !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_count_nodesAndElements(nNodes, nElems, fileUnit) - - integer, intent(in) :: fileUnit - integer, intent(out) :: nNodes, nElems +subroutine mesh_marc_count_nodesAndElements(nNodes, nElems, fileContent) + integer, intent(out) :: nNodes, nElems + character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines + integer, allocatable, dimension(:) :: chunkPos - character(len=300) :: line + integer :: l nNodes = 0 nElems = 0 - rewind(fileUnit) - do - read (fileUnit,'(A300)',END=620) line - chunkPos = IO_stringPos(line) - - if ( IO_lc(IO_StringValue(line,chunkPos,1)) == 'sizing') & - nElems = IO_IntValue (line,chunkPos,3) - if ( IO_lc(IO_StringValue(line,chunkPos,1)) == 'coordinates') then - read (fileUnit,'(A300)') line - chunkPos = IO_stringPos(line) - nNodes = IO_IntValue (line,chunkPos,2) - exit ! assumes that "coordinates" comes later in file + do l = 1, size(fileContent) + chunkPos = IO_stringPos(fileContent(l)) + if (IO_lc(IO_StringValue(fileContent(l),chunkPos,1)) == 'sizing') then + nElems = IO_IntValue (fileContent(l),chunkPos,3) + elseif (IO_lc(IO_StringValue(fileContent(l),chunkPos,1)) == 'coordinates') then + chunkPos = IO_stringPos(fileContent(l+1)) + nNodes = IO_IntValue (fileContent(l+1),chunkPos,2) endif enddo -620 end subroutine mesh_marc_count_nodesAndElements +end subroutine mesh_marc_count_nodesAndElements !-------------------------------------------------------------------------------------------------- From ed1d06d6f14a4a3cf85854556878c1c4b33b81db Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Oct 2019 11:53:35 +0200 Subject: [PATCH 28/48] make functions independent of file state --- src/mesh_marc.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 4c854a2f8..402731345 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -657,9 +657,8 @@ function mesh_marc_buildElements(nElem,nNodes,fileUnit) exit endif enddo -620 rewind(fileUnit) - -end function mesh_marc_buildElements + +620 end function mesh_marc_buildElements !-------------------------------------------------------------------------------------------------- @@ -683,6 +682,7 @@ subroutine mesh_marc_buildElements2(microstructureAt,homogenizationAt, & integer, dimension(1+nElem) :: contInts integer :: i,j,t,sv,myVal,e,nNodesAlreadyRead + rewind(fileUnit) read (fileUnit,'(A300)',END=630) line do chunkPos = IO_stringPos(line) From 004abc2d4e41ffb10f8c2519dc3fce765af1ffaa Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Oct 2019 12:19:37 +0200 Subject: [PATCH 29/48] cell node definition needs to be stored --- src/mesh_marc.f90 | 73 ++++++++++++++++++++++++++++------------------- 1 file changed, 44 insertions(+), 29 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 402731345..598af3921 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -23,7 +23,12 @@ module mesh implicit none private - + + type tCellNodeDefinition + integer, dimension(:,:), allocatable :: parents + integer, dimension(:,:), allocatable :: weights + end type tCellNodeDefinition + real(pReal), public, protected :: & mesh_unitlength !< physical length of one unit in mesh @@ -38,7 +43,7 @@ module mesh !-------------------------------------------------------------------------------------------------- integer, dimension(:,:), allocatable :: & - mesh_FEnodes + connectivity_elem real(pReal), dimension(:,:), allocatable :: & mesh_node, & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!! @@ -164,7 +169,7 @@ subroutine mesh_init(ip,el) allocate(microstructureAt(theMesh%nElems), source=0) allocate(homogenizationAt(theMesh%nElems), source=0) - mesh_FEnodes = mesh_marc_buildElements(theMesh%nElems,theMesh%elem%nNodes,FILEUNIT) + connectivity_elem = mesh_marc_buildElements(theMesh%nElems,theMesh%elem%nNodes,FILEUNIT) call mesh_marc_buildElements2(microstructureAt,homogenizationAt, & mesh_nElems,theMesh%elem%nNodes,initialcondTableStyle,FILEUNIT) if (myDebug) write(6,'(a)') ' Built elements'; flush(6) @@ -174,12 +179,12 @@ subroutine mesh_init(ip,el) #if defined(DAMASK_HDF5) call results_openJobFile call HDF5_closeGroup(results_addGroup('geometry')) - call results_writeDataset('geometry',mesh_FEnodes,'C',& + call results_writeDataset('geometry',connectivity_elem,'C',& 'connectivity of the elements','-') call results_closeJobFile #endif - call buildCells(theMesh,theMesh%elem,mesh_FEnodes) + call buildCells(theMesh,theMesh%elem,connectivity_elem) call mesh_build_cellconnectivity if (myDebug) write(6,'(a)') ' Built cell connectivity'; flush(6) @@ -641,14 +646,16 @@ function mesh_marc_buildElements(nElem,nNodes,fileUnit) if (e /= 0) then ! disregard non CP elems nNodesAlreadyRead = 0 do j = 1,chunkPos(1)-2 - mesh_marc_buildElements(j,e) = mesh_FEasCP('node',IO_IntValue(line,chunkPos,j+2)) ! CP ids of nodes + mesh_marc_buildElements(j,e) = & + mesh_FEasCP('node',IO_IntValue(line,chunkPos,j+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) - mesh_marc_buildElements(nNodesAlreadyRead+j,e) = mesh_FEasCP('node',IO_IntValue(line,chunkPos,j)) ! CP ids of nodes + mesh_marc_buildElements(nNodesAlreadyRead+j,e) = & + mesh_FEasCP('node',IO_IntValue(line,chunkPos,j)) enddo nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1) enddo @@ -731,6 +738,8 @@ subroutine buildCells(thisMesh,elem,connectivity_elem) integer,dimension(:,:), allocatable :: parentsAndWeights,candidates_global, connectivity_cell_reshape integer,dimension(:,:,:), allocatable :: connectivity_cell + type(tCellNodeDefinition), dimension(:), allocatable :: cellNodeDefinition + real(pReal), dimension(:,:), allocatable :: nodes_new,nodes integer :: e, n, c, p, s,i,m,j,nParentNodes,nCellNode,Nelem,candidateID @@ -755,6 +764,8 @@ subroutine buildCells(thisMesh,elem,connectivity_elem) enddo nCellNode = thisMesh%nNodes + + allocate(cellNodeDefinition(elem%nNodes-1)) !--------------------------------------------------------------------------------------------------- ! set connectivity of cell nodes that are defined by 2,...,nNodes real nodes @@ -815,7 +826,9 @@ subroutine buildCells(thisMesh,elem,connectivity_elem) enddo i = uniqueRows(candidates_global(1:2*nParentNodes,:)) - + + allocate(cellNodeDefinition(nParentNodes-1)%parents(i,nParentNodes)) + allocate(cellNodeDefinition(nParentNodes-1)%weights(i,nParentNodes)) ! calculate coordinates of cell nodes and insert their ID into the cell conectivity nodes_new = reshape([(0.0_pReal,j = 1, 3*i)], [3,i]) @@ -823,27 +836,29 @@ subroutine buildCells(thisMesh,elem,connectivity_elem) i = 1 n = 1 do while(n <= size(candidates_local)*Nelem) - j=0 - parentsAndWeights(:,1) = candidates_global(1:nParentNodes,n+j) - parentsAndWeights(:,2) = candidates_global(nParentNodes+1:nParentNodes*2,n+j) - e = candidates_global(nParentNodes*2+1,n+j) - c = candidates_global(nParentNodes*2+2,n+j) - do m = 1, nParentNodes - nodes_new(:,i) = nodes_new(:,i) & - + thisMesh%node_0(:,parentsAndWeights(m,1)) * real(parentsAndWeights(m,2),pReal) - enddo - nodes_new(:,i) = nodes_new(:,i)/real(sum(parentsAndWeights(:,2)),pReal) + j=0 + parentsAndWeights(:,1) = candidates_global(1:nParentNodes,n+j) + parentsAndWeights(:,2) = candidates_global(nParentNodes+1:nParentNodes*2,n+j) + e = candidates_global(nParentNodes*2+1,n+j) + c = candidates_global(nParentNodes*2+2,n+j) + do m = 1, nParentNodes + nodes_new(:,i) = nodes_new(:,i) & + + thisMesh%node_0(:,parentsAndWeights(m,1)) * real(parentsAndWeights(m,2),pReal) + enddo + nodes_new(:,i) = nodes_new(:,i)/real(sum(parentsAndWeights(:,2)),pReal) - do while (n+j<= size(candidates_local)*Nelem) - if (any(candidates_global(1:2*nParentNodes,n+j)/=candidates_global(1:2*nParentNodes,n))) exit - where (connectivity_cell(:,:,candidates_global(nParentNodes*2+1,n+j)) == -candidates_global(nParentNodes*2+2,n+j)) ! still locally defined - connectivity_cell(:,:,candidates_global(nParentNodes*2+1,n+j)) = nCellNode + i ! gets current new cell node id - end where + do while (n+j<= size(candidates_local)*Nelem) + if (any(candidates_global(1:2*nParentNodes,n+j)/=candidates_global(1:2*nParentNodes,n))) exit + where (connectivity_cell(:,:,candidates_global(nParentNodes*2+1,n+j)) == -candidates_global(nParentNodes*2+2,n+j)) ! still locally defined + connectivity_cell(:,:,candidates_global(nParentNodes*2+1,n+j)) = nCellNode + i ! gets current new cell node id + end where - j = j + 1 - enddo - i=i+1 - n = n+j + j = j+1 + enddo + cellNodeDefinition(nParentNodes-1)%parents(i,:) = parentsAndWeights(:,1) + cellNodeDefinition(nParentNodes-1)%weights(i,:) = parentsAndWeights(:,2) + i = i+1 + n = n+j enddo nCellNode = nCellNode + i @@ -940,7 +955,7 @@ subroutine mesh_build_cellconnectivity do n = 1,theMesh%elem%NcellnodesPerCell localCellnodeID = theMesh%elem%cell(n,i) if (localCellnodeID <= FE_NmatchingNodes(theMesh%elem%geomType)) then ! this cell node is a matching node - matchingNodeID = mesh_FEnodes(localCellnodeID,e) + matchingNodeID = connectivity_elem(localCellnodeID,e) if (matchingNode2cellnode(matchingNodeID) == 0) then ! if this matching node does not yet exist in the glbal cell node list ... mesh_Ncellnodes = mesh_Ncellnodes + 1 ! ... count it as cell node ... matchingNode2cellnode(matchingNodeID) = mesh_Ncellnodes ! ... and remember its global ID @@ -994,7 +1009,7 @@ function mesh_build_cellnodes() localCellnodeID = mesh_cellnodeParent(2,n) myCoords = 0.0_pReal do m = 1,theMesh%elem%nNodes - myCoords = myCoords + mesh_node(1:3,mesh_FEnodes(m,e)) & + myCoords = myCoords + mesh_node(1:3,connectivity_elem(m,e)) & * theMesh%elem%cellNodeParentNodeWeights(m,localCellnodeID) enddo mesh_build_cellnodes(1:3,n) = myCoords / sum(theMesh%elem%cellNodeParentNodeWeights(:,localCellnodeID)) From bbc2ed90a2e052fcaaed7eac8a2029bd8fd37bb6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Oct 2019 17:48:38 +0200 Subject: [PATCH 30/48] use own HDF5 libraries requires compilation with openMP (and Intel compiler) --- .../mods_MarcMentat/2019/Marc_tools/include_linux64 | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/installation/mods_MarcMentat/2019/Marc_tools/include_linux64 b/installation/mods_MarcMentat/2019/Marc_tools/include_linux64 index 0ef7c5733..290055ed3 100644 --- a/installation/mods_MarcMentat/2019/Marc_tools/include_linux64 +++ b/installation/mods_MarcMentat/2019/Marc_tools/include_linux64 @@ -63,6 +63,7 @@ else INTEGER_PATH=/$MARC_INTEGER_SIZE fi +FCOMP=ifort INTELPATH="/opt/intel/compilers_and_libraries_2017/linux" # find the root directory of the compiler installation: @@ -103,9 +104,6 @@ if test "$DAMASK_HDF5" = "ON";then H5FC="$(h5fc -shlib -show)" HDF5_LIB=${H5FC//ifort/} FCOMP="$H5FC -DDAMASK_HDF5" - echo $FCOMP -else - FCOMP=ifort fi # AEM @@ -531,7 +529,7 @@ else FORT_OPT=" $FORT_OPT -save -zero" fi if test "$MARCHDF" = "HDF"; then - FORT_OPT="$FORT_OPT -DMARCHDF=$MARCHDF $HDF_INCLUDE" + FORT_OPT="$FORT_OPT -DMARCHDF=$MARCHDF" fi FORTLOW="$FCOMP $FORT_OPT $PROFILE -O0 $I8FFLAGS -I$MARC_SOURCE/common \ @@ -757,7 +755,7 @@ SECLIBS="-L$MARC_LIB -llapi" SOLVERLIBS="${BCSSOLVERLIBS} ${VKISOLVERLIBS} ${CASISOLVERLIBS} ${MF2SOLVERLIBS} \ $MKLLIB -L$MARC_MKL -liomp5 \ - $MARC_LIB/blas_src.a ${ACSI_LIB}/ACSI_MarcLib.a $KDTREE2_LIB/kdtree2.a $HDF_LIBS $HDF_LIB" + $MARC_LIB/blas_src.a ${ACSI_LIB}/ACSI_MarcLib.a $KDTREE2_LIB/kdtree2.a $HDF5_LIB" SOLVERLIBS_DLL=${SOLVERLIBS} if test "$AEM_DLL" -eq 1 From e3b16639bf10e9f3ad2b67f2ea4808e2973664d9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Oct 2019 18:20:54 +0200 Subject: [PATCH 31/48] native integer needs to match otherwise, results are wrong. Therefore, we need to use our own HDF5 library since MSC provides one for 4 byte integers --- .gitlab-ci.yml | 2 +- src/HDF5_utilities.f90 | 4 ---- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index ad22c7998..df5a10f60 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -75,7 +75,7 @@ variables: MSC: "$MSC2019" IntelMarc: "$IntelCompiler17_8" IntelAbaqus: "$IntelCompiler16_4" - HDF5Marc: "HDF5/1.10.4/Intel-17.8" + HDF5Marc: "HDF5/1.10.5/Intel-17.8" # ++++++++++++ Documentation ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Doxygen1_8_15: "Documentation/Doxygen/1.8.15" # ------------ Defaults ---------------------------------------------- diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90 index 68281a6a0..e4819431e 100644 --- a/src/HDF5_utilities.f90 +++ b/src/HDF5_utilities.f90 @@ -112,14 +112,10 @@ subroutine HDF5_utilities_init call h5open_f(hdferr) if (hdferr < 0) call IO_error(1,ext_msg='HDF5_Utilities_init: h5open_f') -#ifndef Marc4DAMASK -! This test should ensure that integer size matches. For some reasons, the HDF5 libraries -! that come with MSC.Marc>=2019 seem to be of 4byte even though it is a 8byte Marc version call h5tget_size_f(H5T_NATIVE_INTEGER,typeSize, hdferr) if (hdferr < 0) call IO_error(1,ext_msg='HDF5_Utilities_init: h5tget_size_f (int)') if (int(bit_size(0),SIZE_T)/=typeSize*8) & call IO_error(0,ext_msg='Default integer size does not match H5T_NATIVE_INTEGER') -#endif call h5tget_size_f(H5T_NATIVE_DOUBLE,typeSize, hdferr) if (hdferr < 0) call IO_error(1,ext_msg='HDF5_Utilities_init: h5tget_size_f (double)') From e0cb1a87cd459e105f6ae0fd64aafe0f59e81b88 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Oct 2019 19:22:57 +0200 Subject: [PATCH 32/48] cleaning --- src/mesh_marc.f90 | 140 ---------------------------------------------- 1 file changed, 140 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 598af3921..567bc3f81 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -125,10 +125,6 @@ subroutine mesh_init(ip,el) integer, dimension(:), allocatable :: & marc_matNumber !< array of material numbers for hypoelastic material (Marc only) logical :: myDebug - real(pReal), dimension(:,:,:), allocatable:: & - mesh_ipArea !< area of interface to neighboring IP (initially!) - real(pReal),dimension(:,:,:,:), allocatable :: & - mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!) real(pReal), dimension(:,:), allocatable :: & ip_reshaped @@ -195,9 +191,6 @@ subroutine mesh_init(ip,el) call mesh_build_ipCoordinates if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6) - allocate(mesh_ipArea(theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems)) - allocate(mesh_ipAreaNormal(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems)) - call mesh_build_ipAreas(mesh_ipArea,mesh_ipAreaNormal) if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) call IP_neighborhood2 @@ -230,10 +223,7 @@ subroutine mesh_init(ip,el) 'nodal coordinates','m') call results_closeJobFile() #endif - call geometry_plastic_nonlocal_setIPvolume(IPvolume()) call geometry_plastic_nonlocal_setIPneighborhood(mesh_ipNeighborhood2) - call geometry_plastic_nonlocal_setIParea(mesh_ipArea) - call geometry_plastic_nonlocal_setIPareaNormal(mesh_ipAreaNormal) end subroutine mesh_init @@ -1018,68 +1008,6 @@ function mesh_build_cellnodes() end function mesh_build_cellnodes -!--------------------------------------------------------------------------------------------------- -!> @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. -!> 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. -!--------------------------------------------------------------------------------------------------- -function IPvolume() - - real(pReal), dimension(theMesh%elem%nIPs,theMesh%nElems) :: IPvolume - - integer :: e,i,c,m,f,n - real(pReal), dimension(size(theMesh%elem%cellFace,1),size(theMesh%elem%cellFace,2)) :: subvolume - - c = theMesh%elem%cellType - m = size(theMesh%elem%cellFace,1) - - do e = 1,theMesh%nElems - select case (c) - - case (1) ! 2D 3node - forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element - IPvolume(i,e) = math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & - theMesh%node_0(1:3,mesh_cell2(2,i,e)), & - theMesh%node_0(1:3,mesh_cell2(3,i,e))) - - case (2) ! 2D 4node - forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element - IPvolume(i,e) = math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & ! here we assume a planar shape, so division in two triangles suffices - theMesh%node_0(1:3,mesh_cell2(2,i,e)), & - theMesh%node_0(1:3,mesh_cell2(3,i,e))) & - + math_areaTriangle(theMesh%node_0(1:3,mesh_cell2(3,i,e)), & - theMesh%node_0(1:3,mesh_cell2(4,i,e)), & - theMesh%node_0(1:3,mesh_cell2(1,i,e))) - - case (3) ! 3D 4node - forall (i = 1:theMesh%elem%nIPs) & ! loop over ips=cells in this element - IPvolume(i,e) = math_volTetrahedron(theMesh%node_0(1:3,mesh_cell2(1,i,e)), & - theMesh%node_0(1:3,mesh_cell2(2,i,e)), & - theMesh%node_0(1:3,mesh_cell2(3,i,e)), & - theMesh%node_0(1:3,mesh_cell2(4,i,e))) - - case (4) ! 3D 8node - do i = 1,theMesh%elem%nIPs ! loop over ips=cells in this element - subvolume = 0.0_pReal - forall(f = 1:FE_NipNeighbors(c), n = 1:m) & - subvolume(n,f) = math_volTetrahedron(& - mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface( n ,f),i,e)), & - mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(1+mod(n ,m),f),i,e)), & - mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(1+mod(n+1,m),f),i,e)), & - mesh_ipCoordinates(1:3,i,e)) - IPvolume(i,e) = 0.5_pReal * sum(subvolume) ! each subvolume is based on four tetrahedrons, altough the face consists of only two triangles -> averaging factor two - enddo - - end select - enddo - -end function IPvolume - - !--------------------------------------------------------------------------------------------------- !> @brief cell neighborhood !--------------------------------------------------------------------------------------------------- @@ -1185,74 +1113,6 @@ subroutine mesh_build_ipCoordinates end subroutine mesh_build_ipCoordinates -!-------------------------------------------------------------------------------------------------- -!> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipAreas(ipArea,ipAreaNormal) - - integer :: e,c,i,f,n,m - real(pReal), dimension(theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), intent(out) :: ipArea - real(pReal), dimension(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), intent(out) :: ipAreaNormal - real(pReal), dimension (3,size(theMesh%elem%cellFace,2)) :: nodePos, normals - real(pReal), dimension(3) :: normal - - c = theMesh%elem%cellType - - - do e = 1,theMesh%nElems ! loop over cpElems - select case (c) - - case (1,2) ! 2D 3 or 4 node - do i = 1,theMesh%elem%nIPs - do f = 1,FE_NipNeighbors(c) ! loop over cell faces - forall(n = 1: size(theMesh%elem%cellface,1)) & - nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(n,f),i,e)) - 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 - ipArea(f,i,e) = norm2(normal) - ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal - enddo - enddo - - case (3) ! 3D 4node - do i = 1,theMesh%elem%nIPs - do f = 1,FE_NipNeighbors(c) ! loop over cell faces - forall(n = 1: size(theMesh%elem%cellface,1)) & - nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(n,f),i,e)) - normal = math_cross(nodePos(1:3,2) - nodePos(1:3,1), & - nodePos(1:3,3) - nodePos(1:3,1)) - ipArea(f,i,e) = norm2(normal) - ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal - enddo - enddo - - 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 - m = size(theMesh%elem%cellFace,1) - do i = 1,theMesh%elem%nIPs - do f = 1,FE_NipNeighbors(c) ! loop over cell faces - forall(n = 1: size(theMesh%elem%cellface,1)) & - nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(theMesh%elem%cellface(n,f),i,e)) - forall(n = 1: size(theMesh%elem%cellface,1)) & - 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) - ipArea(f,i,e) = norm2(normal) - ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) - enddo - enddo - - end select - 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 4999aa4e14cb640a4473e0b58dfd30acc91bf0b5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Oct 2019 19:36:30 +0200 Subject: [PATCH 33/48] more cleaning --- src/DAMASK_marc.f90 | 2 +- src/mesh_marc.f90 | 30 ------------------------------ 2 files changed, 1 insertion(+), 31 deletions(-) diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index bea952fca..09a1c83c8 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -266,7 +266,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & outdatedFFN1 = .false. cycleCounter = cycleCounter + 1 !mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates - call mesh_build_ipCoordinates() ! update ip coordinates + !call mesh_build_ipCoordinates() ! update ip coordinates endif if (outdatedByNewInc) then computationMode = ior(computationMode,CPFEM_AGERESULTS) ! calc and age results diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 567bc3f81..6a7fcf71a 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -100,7 +100,6 @@ integer, dimension(:,:), allocatable :: & public :: & mesh_init, & mesh_build_cellnodes, & - mesh_build_ipCoordinates, & mesh_FEasCP @@ -188,7 +187,6 @@ subroutine mesh_init(ip,el) if (myDebug) write(6,'(a)') ' Built cell nodes'; flush(6) allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal) - call mesh_build_ipCoordinates if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6) if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) @@ -1085,34 +1083,6 @@ subroutine IP_neighborhood2 end subroutine IP_neighborhood2 -!-------------------------------------------------------------------------------------------------- -!> @brief Calculates IP Coordinates. -! Marc however only provides nodal displacements, -! so in this case the ip coordinates are always calculated on the basis of this subroutine. -! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! FOR THE MOMENT THIS SUBROUTINE ACTUALLY CALCULATES THE CELL CENTER AND NOT THE IP COORDINATES, -! AS THE IP IS NOT (ALWAYS) LOCATED IN THE CENTER OF THE IP VOLUME. -! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipCoordinates - - integer :: e,i,n - real(pReal), dimension(3) :: myCoords - - - do e = 1,theMesh%nElems - do i = 1,theMesh%elem%nIPs - myCoords = 0.0_pReal - do n = 1,theMesh%elem%nCellnodesPerCell - myCoords = myCoords + mesh_cellnode(1:3,mesh_cell2(n,i,e)) - enddo - mesh_ipCoordinates(1:3,i,e) = myCoords / real(theMesh%elem%nCellnodesPerCell,pReal) - enddo - enddo - -end subroutine mesh_build_ipCoordinates - - !-------------------------------------------------------------------------------------------------- !> @brief Gives the FE to CP ID mapping by binary search through lookup array !! valid questions (what) are 'elem', 'node' From c2a87019400e7d8b612f662252f40eda4dcc7486 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Oct 2019 19:55:25 +0200 Subject: [PATCH 34/48] further removal of deprecated stuff --- src/mesh_marc.f90 | 147 +--------------------------------------------- 1 file changed, 1 insertion(+), 146 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 6a7fcf71a..ebc14bc13 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -37,9 +37,6 @@ module mesh real(pReal), dimension(:,:,:), allocatable, public :: & mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) - - real(pReal), dimension(:,:), allocatable, public :: & - mesh_cellnode !< cell node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) !-------------------------------------------------------------------------------------------------- integer, dimension(:,:), allocatable :: & @@ -55,29 +52,13 @@ module mesh integer:: & - mesh_Ncellnodes, & !< total number of cell nodes in mesh (including duplicates) - mesh_elemType, & !< Element type of the mesh (only support homogeneous meshes) - mesh_Nnodes, & !< total number of nodes in mesh - mesh_Ncells, & !< total number of cells in mesh - mesh_maxNsharedElems !< max number of CP elements sharing a node + mesh_Nnodes !< total number of nodes in mesh -integer, dimension(:,:), allocatable :: & - mesh_cellnodeParent !< cellnode's parent element ID, cellnode's intra-element ID - integer,dimension(:,:,:), allocatable :: & mesh_cell2, & !< cell connectivity for each element,ip/cell mesh_cell !< cell connectivity for each element,ip/cell - - integer, dimension(4), parameter :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type - int([& - 3, & ! (2D 3node) - 4, & ! (2D 4node) - 4, & ! (3D 4node) - 6 & ! (3D 8node) - ],pInt) - integer, dimension(:), allocatable :: & microstructureAt, & homogenizationAt @@ -99,7 +80,6 @@ integer, dimension(:,:), allocatable :: & public :: & mesh_init, & - mesh_build_cellnodes, & mesh_FEasCP @@ -123,7 +103,6 @@ subroutine mesh_init(ip,el) initialcondTableStyle integer, dimension(:), allocatable :: & marc_matNumber !< array of material numbers for hypoelastic material (Marc only) - logical :: myDebug real(pReal), dimension(:,:), allocatable :: & ip_reshaped @@ -132,7 +111,6 @@ subroutine mesh_init(ip,el) mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh inputFile = IO_read_ASCII(trim(modelName)//trim(InputFileExtension)) - myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0) ! parsing Marc input file fileFormatVersion = mesh_marc_get_fileFormat(inputFile) @@ -156,7 +134,6 @@ subroutine mesh_init(ip,el) mesh_node = mesh_node0 elemType = mesh_marc_getElemType(mesh_nElems,FILEUNIT) - if (myDebug) write(6,'(a)') ' Counted CP sizes'; flush(6) call theMesh%init('mesh',elemType,mesh_node0) call theMesh%setNelems(mesh_nElems) @@ -167,7 +144,6 @@ subroutine mesh_init(ip,el) connectivity_elem = mesh_marc_buildElements(theMesh%nElems,theMesh%elem%nNodes,FILEUNIT) call mesh_marc_buildElements2(microstructureAt,homogenizationAt, & mesh_nElems,theMesh%elem%nNodes,initialcondTableStyle,FILEUNIT) - if (myDebug) write(6,'(a)') ' Built elements'; flush(6) close (FILEUNIT) @@ -181,18 +157,9 @@ subroutine mesh_init(ip,el) call buildCells(theMesh,theMesh%elem,connectivity_elem) - call mesh_build_cellconnectivity - if (myDebug) write(6,'(a)') ' Built cell connectivity'; flush(6) - mesh_cellnode = mesh_build_cellnodes() - if (myDebug) write(6,'(a)') ' Built cell nodes'; flush(6) - allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal) - if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6) - - if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) call IP_neighborhood2 - if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6) if (usePingPong .and. (mesh_Nelems /= theMesh%nElems)) & call IO_error(600) ! ping-pong must be disabled when having non-DAMASK elements @@ -894,118 +861,6 @@ subroutine buildCells(thisMesh,elem,connectivity_elem) end subroutine buildCells -!-------------------------------------------------------------------------------------------------- -!> @brief Split CP elements into cells. -!> @details Build a mapping between cells and the corresponding cell nodes ('mesh_cell'). -!> Cell nodes that are also matching nodes are unique in the list of cell nodes, -!> all others (currently) might be stored more than once. -!-------------------------------------------------------------------------------------------------- -subroutine mesh_build_cellconnectivity - - integer, dimension(:), allocatable :: & - matchingNode2cellnode - - integer, dimension(10), parameter :: FE_NmatchingNodes = & !< number of nodes that are needed for face matching in a specific type of element geometry - int([ & - 3, & ! element 6 (2D 3node 1ip) - 3, & ! element 125 (2D 6node 3ip) - 4, & ! element 11 (2D 4node 4ip) - 4, & ! element 27 (2D 8node 9ip) - 4, & ! element 134 (3D 4node 1ip) - 4, & ! element 127 (3D 10node 4ip) - 6, & ! element 136 (3D 6node 6ip) - 8, & ! element 117 (3D 8node 1ip) - 8, & ! element 7 (3D 8node 8ip) - 8 & ! element 21 (3D 20node 27ip) - ],pInt) - - integer, dimension(:,:), allocatable :: & - cellnodeParent - integer, dimension(theMesh%elem%Ncellnodes) :: & - localCellnode2globalCellnode - integer :: & - e,n,i, & - matchingNodeID, & - localCellnodeID - - allocate(mesh_cell(theMesh%elem%NcellNodesPerCell,theMesh%elem%nIPs,theMesh%nElems), source=0) - allocate(matchingNode2cellnode(theMesh%nNodes), source=0) - allocate(cellnodeParent(2,theMesh%elem%Ncellnodes*theMesh%nElems), source=0) - - mesh_Ncells = theMesh%nElems*theMesh%elem%nIPs -!-------------------------------------------------------------------------------------------------- -! Count cell nodes (including duplicates) and generate cell connectivity list - mesh_Ncellnodes = 0 - - do e = 1,theMesh%nElems - localCellnode2globalCellnode = 0 - do i = 1,theMesh%elem%nIPs - do n = 1,theMesh%elem%NcellnodesPerCell - localCellnodeID = theMesh%elem%cell(n,i) - if (localCellnodeID <= FE_NmatchingNodes(theMesh%elem%geomType)) then ! this cell node is a matching node - matchingNodeID = connectivity_elem(localCellnodeID,e) - if (matchingNode2cellnode(matchingNodeID) == 0) then ! if this matching node does not yet exist in the glbal cell node list ... - mesh_Ncellnodes = mesh_Ncellnodes + 1 ! ... count it as cell node ... - matchingNode2cellnode(matchingNodeID) = mesh_Ncellnodes ! ... and remember its global ID - cellnodeParent(1,mesh_Ncellnodes) = e ! ... and where it belongs to - cellnodeParent(2,mesh_Ncellnodes) = localCellnodeID - endif - mesh_cell(n,i,e) = matchingNode2cellnode(matchingNodeID) - else ! this cell node is no matching node - if (localCellnode2globalCellnode(localCellnodeID) == 0) then ! if this local cell node does not yet exist in the global cell node list ... - mesh_Ncellnodes = mesh_Ncellnodes + 1 ! ... count it as cell node ... - localCellnode2globalCellnode(localCellnodeID) = mesh_Ncellnodes ! ... and remember its global ID ... - cellnodeParent(1,mesh_Ncellnodes) = e ! ... and it belongs to - cellnodeParent(2,mesh_Ncellnodes) = localCellnodeID - endif - mesh_cell(n,i,e) = localCellnode2globalCellnode(localCellnodeID) - endif - enddo - enddo - enddo - - allocate(mesh_cellnodeParent(2,mesh_Ncellnodes)) - allocate(mesh_cellnode(3,mesh_Ncellnodes)) - - forall(n = 1:mesh_Ncellnodes) - mesh_cellnodeParent(1,n) = cellnodeParent(1,n) - mesh_cellnodeParent(2,n) = cellnodeParent(2,n) - endforall - -end subroutine mesh_build_cellconnectivity - - -!-------------------------------------------------------------------------------------------------- -!> @brief Calculate position of cellnodes from the given position of nodes -!> Build list of cellnodes' coordinates. -!> Cellnode coordinates are calculated from a weighted sum of node coordinates. -!-------------------------------------------------------------------------------------------------- -function mesh_build_cellnodes() - - - real(pReal), dimension(3,mesh_Ncellnodes) :: mesh_build_cellnodes - - integer :: & - e,n,m, & - localCellnodeID - real(pReal), dimension(3) :: & - myCoords - - mesh_build_cellnodes = 0.0_pReal - do n = 1,mesh_Ncellnodes ! loop over cell nodes - e = mesh_cellnodeParent(1,n) - localCellnodeID = mesh_cellnodeParent(2,n) - myCoords = 0.0_pReal - do m = 1,theMesh%elem%nNodes - myCoords = myCoords + mesh_node(1:3,connectivity_elem(m,e)) & - * theMesh%elem%cellNodeParentNodeWeights(m,localCellnodeID) - enddo - mesh_build_cellnodes(1:3,n) = myCoords / sum(theMesh%elem%cellNodeParentNodeWeights(:,localCellnodeID)) - enddo - -end function mesh_build_cellnodes - - !--------------------------------------------------------------------------------------------------- !> @brief cell neighborhood !--------------------------------------------------------------------------------------------------- From 6456f9891b385815659092b7f40239018737e2e9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Oct 2019 20:34:20 +0200 Subject: [PATCH 35/48] avoid public variables --- src/mesh_marc.f90 | 62 +++++++++++++++++++++++------------------------ 1 file changed, 30 insertions(+), 32 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index ebc14bc13..a85924e2a 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -30,15 +30,20 @@ module mesh end type tCellNodeDefinition real(pReal), public, protected :: & - mesh_unitlength !< physical length of one unit in mesh + mesh_unitlength !< physical length of one unit in mesh + + integer, dimension(:,:), allocatable, target :: & + mesh_mapFEtoCPelem, & !< [sorted FEid, corresponding CPid] + mesh_mapFEtoCPnode !< [sorted FEid, corresponding CPid] !-------------------------------------------------------------------------------------------------- -! public variables (DEPRECATED) - - real(pReal), dimension(:,:,:), allocatable, public :: & - mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) +! DEPRECATED + real(pReal), dimension(:,:,:), allocatable, public :: & + mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) !-------------------------------------------------------------------------------------------------- + + integer, dimension(:,:), allocatable :: & connectivity_elem @@ -46,33 +51,10 @@ module mesh mesh_node, & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!! mesh_node0 !< node x,y,z coordinates (initially!) -! -------------------------------------------------------------------------------------------------- - type(tMesh) :: theMesh - - integer:: & - mesh_Nnodes !< total number of nodes in mesh - - integer,dimension(:,:,:), allocatable :: & - mesh_cell2, & !< cell connectivity for each element,ip/cell - mesh_cell !< cell connectivity for each element,ip/cell - - integer, dimension(:), allocatable :: & - microstructureAt, & - homogenizationAt - - integer :: & - mesh_NelemSets - - character(len=64), dimension(:), allocatable :: & - mesh_nameElemSet - integer, dimension(:,:), allocatable :: & - mesh_mapElemSet !< list of elements in elementSet - integer, dimension(:,:), allocatable, target :: & - mesh_mapFEtoCPelem, & !< [sorted FEid, corresponding CPid] - mesh_mapFEtoCPnode !< [sorted FEid, corresponding CPid] + mesh_cell2 !< cell connectivity for each element,ip/cell integer, dimension(:,:,:,:), allocatable :: & mesh_ipNeighborhood2 !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] @@ -95,6 +77,8 @@ subroutine mesh_init(ip,el) integer, parameter :: FILEUNIT = 222 character(len=pStringLen), dimension(:), allocatable :: inputFile !< file content, separated per lines + integer :: & + mesh_NelemSets integer :: j, fileFormatVersion, elemType, & mesh_maxNelemInSet, & @@ -103,6 +87,15 @@ subroutine mesh_init(ip,el) 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 real(pReal), dimension(:,:), allocatable :: & ip_reshaped @@ -143,7 +136,8 @@ subroutine mesh_init(ip,el) connectivity_elem = mesh_marc_buildElements(theMesh%nElems,theMesh%elem%nNodes,FILEUNIT) call mesh_marc_buildElements2(microstructureAt,homogenizationAt, & - mesh_nElems,theMesh%elem%nNodes,initialcondTableStyle,FILEUNIT) + mesh_nElems,theMesh%elem%nNodes,mesh_nameElemSet,mesh_mapElemSet,& + initialcondTableStyle,FILEUNIT) close (FILEUNIT) @@ -627,7 +621,7 @@ function mesh_marc_buildElements(nElem,nNodes,fileUnit) !> @brief Stores homogenization and microstructure ID !-------------------------------------------------------------------------------------------------- subroutine mesh_marc_buildElements2(microstructureAt,homogenizationAt, & - nElem,nNodes,initialcondTableStyle,fileUnit) + nElem,nNodes,nameElemSet,mapElemSet,initialcondTableStyle,fileUnit) integer, dimension(:), intent(out) :: & microstructureAt, & @@ -637,6 +631,10 @@ subroutine mesh_marc_buildElements2(microstructureAt,homogenizationAt, & nNodes, & !< number of nodes per element initialcondTableStyle, & fileUnit + character(len=64), dimension(:), intent(in) :: & + nameElemSet + integer, dimension(:,:), intent(in) :: & + mapElemSet !< list of elements in elementSet integer, allocatable, dimension(:) :: chunkPos character(len=300) line @@ -664,7 +662,7 @@ subroutine mesh_marc_buildElements2(microstructureAt,homogenizationAt, & read (fileUnit,'(A300)',END=630) line ! read extra line endif contInts = IO_continuousIntValues& ! get affected elements - (fileUnit,theMesh%nElems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) + (fileUnit,nElem,nameElemSet,mapElemSet,size(nameElemSet)) do i = 1,contInts(1) e = mesh_FEasCP('elem',contInts(1+i)) if (sv == 2) microstructureAt(e) = myVal From 369cae5332f71b1c7a3407575354d781375255b8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Oct 2019 20:52:49 +0200 Subject: [PATCH 36/48] no public variables --- src/mesh_marc.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index a85924e2a..d564674e8 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -47,9 +47,6 @@ module mesh integer, dimension(:,:), allocatable :: & connectivity_elem - real(pReal), dimension(:,:), allocatable :: & - mesh_node, & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!! - mesh_node0 !< node x,y,z coordinates (initially!) type(tMesh) :: theMesh @@ -79,6 +76,9 @@ subroutine mesh_init(ip,el) character(len=pStringLen), dimension(:), allocatable :: inputFile !< file content, separated per lines integer :: & mesh_NelemSets + real(pReal), dimension(:,:), allocatable :: & + mesh_node, & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!! + mesh_node0 !< node x,y,z coordinates (initially!) integer :: j, fileFormatVersion, elemType, & mesh_maxNelemInSet, & From d03efade06218826b611c2df51925efd84913088 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Oct 2019 21:09:49 +0200 Subject: [PATCH 37/48] better calculate nodes later --- src/mesh_marc.f90 | 24 +++--------------------- 1 file changed, 3 insertions(+), 21 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index d564674e8..7a5aaa698 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -43,7 +43,6 @@ module mesh !-------------------------------------------------------------------------------------------------- - integer, dimension(:,:), allocatable :: & connectivity_elem @@ -693,7 +692,6 @@ subroutine buildCells(thisMesh,elem,connectivity_elem) type(tCellNodeDefinition), dimension(:), allocatable :: cellNodeDefinition - real(pReal), dimension(:,:), allocatable :: nodes_new,nodes integer :: e, n, c, p, s,i,m,j,nParentNodes,nCellNode,Nelem,candidateID Nelem = thisMesh%Nelems @@ -783,8 +781,6 @@ subroutine buildCells(thisMesh,elem,connectivity_elem) allocate(cellNodeDefinition(nParentNodes-1)%parents(i,nParentNodes)) allocate(cellNodeDefinition(nParentNodes-1)%weights(i,nParentNodes)) - ! calculate coordinates of cell nodes and insert their ID into the cell conectivity - nodes_new = reshape([(0.0_pReal,j = 1, 3*i)], [3,i]) i = 1 n = 1 @@ -794,11 +790,6 @@ subroutine buildCells(thisMesh,elem,connectivity_elem) parentsAndWeights(:,2) = candidates_global(nParentNodes+1:nParentNodes*2,n+j) e = candidates_global(nParentNodes*2+1,n+j) c = candidates_global(nParentNodes*2+2,n+j) - do m = 1, nParentNodes - nodes_new(:,i) = nodes_new(:,i) & - + thisMesh%node_0(:,parentsAndWeights(m,1)) * real(parentsAndWeights(m,2),pReal) - enddo - nodes_new(:,i) = nodes_new(:,i)/real(sum(parentsAndWeights(:,2)),pReal) do while (n+j<= size(candidates_local)*Nelem) if (any(candidates_global(1:2*nParentNodes,n+j)/=candidates_global(1:2*nParentNodes,n))) exit @@ -815,21 +806,12 @@ subroutine buildCells(thisMesh,elem,connectivity_elem) enddo nCellNode = nCellNode + i - if (i/=0) nodes = reshape([nodes,nodes_new],[3,nCellNode]) + + mesh_cell2 = connectivity_cell + enddo - thisMesh%node_0 = nodes - mesh_cell2 = connectivity_cell - -#if defined(DAMASK_HDF5) - connectivity_cell_reshape = reshape(connectivity_cell,[elem%NcellNodesPerCell,elem%nIPs*Nelem]) - call results_openJobFile - call results_writeDataset('geometry',connectivity_cell_reshape,'c',& - 'connectivity of the cells','-') - call results_closeJobFile -#endif contains - !------------------------------------------------------------------------------------------------ !> @brief count unique rows (same rows need to be stored consecutively) !------------------------------------------------------------------------------------------------ From 33e639426af897d328c90020c32905a635aee200 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Oct 2019 21:45:08 +0200 Subject: [PATCH 38/48] polishing --- src/mesh_marc.f90 | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 7a5aaa698..85fa10a53 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -133,7 +133,7 @@ subroutine mesh_init(ip,el) allocate(microstructureAt(theMesh%nElems), source=0) allocate(homogenizationAt(theMesh%nElems), source=0) - connectivity_elem = mesh_marc_buildElements(theMesh%nElems,theMesh%elem%nNodes,FILEUNIT) + connectivity_elem = mesh_marc_buildElements(mesh_nElems,theMesh%elem%nNodes,FILEUNIT) call mesh_marc_buildElements2(microstructureAt,homogenizationAt, & mesh_nElems,theMesh%elem%nNodes,mesh_nameElemSet,mesh_mapElemSet,& initialcondTableStyle,FILEUNIT) @@ -148,7 +148,7 @@ subroutine mesh_init(ip,el) call results_closeJobFile #endif - call buildCells(theMesh,theMesh%elem,connectivity_elem) + call buildCells(mesh_Nnodes,theMesh%elem,connectivity_elem) allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal) @@ -680,9 +680,9 @@ subroutine mesh_marc_buildElements2(microstructureAt,homogenizationAt, & 630 end subroutine mesh_marc_buildElements2 -subroutine buildCells(thisMesh,elem,connectivity_elem) +subroutine buildCells(nNodes,elem,connectivity_elem) - class(tMesh) :: thisMesh + integer, intent(in) :: nNodes type(tElement), intent(in) :: elem integer,dimension(:,:), intent(in) :: connectivity_elem @@ -694,7 +694,7 @@ subroutine buildCells(thisMesh,elem,connectivity_elem) integer :: e, n, c, p, s,i,m,j,nParentNodes,nCellNode,Nelem,candidateID - Nelem = thisMesh%Nelems + Nelem = size(connectivity_elem,1) !--------------------------------------------------------------------------------------------------- ! initialize global connectivity to negative local connectivity @@ -713,7 +713,7 @@ subroutine buildCells(thisMesh,elem,connectivity_elem) endif realNode enddo enddo - nCellNode = thisMesh%nNodes + nCellNode = nNodes allocate(cellNodeDefinition(elem%nNodes-1)) @@ -777,11 +777,9 @@ subroutine buildCells(thisMesh,elem,connectivity_elem) enddo i = uniqueRows(candidates_global(1:2*nParentNodes,:)) - allocate(cellNodeDefinition(nParentNodes-1)%parents(i,nParentNodes)) allocate(cellNodeDefinition(nParentNodes-1)%weights(i,nParentNodes)) - i = 1 n = 1 do while(n <= size(candidates_local)*Nelem) @@ -807,9 +805,10 @@ subroutine buildCells(thisMesh,elem,connectivity_elem) enddo nCellNode = nCellNode + i - mesh_cell2 = connectivity_cell enddo + + mesh_cell2 = connectivity_cell contains !------------------------------------------------------------------------------------------------ From 2b65c888c45c507614f39e0f3a7fa43e6a2ae4e4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Oct 2019 22:16:42 +0200 Subject: [PATCH 39/48] avoid public variables --- src/mesh_marc.f90 | 136 +++++++++------------------------------------- 1 file changed, 27 insertions(+), 109 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 85fa10a53..a4d5c169e 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -29,6 +29,8 @@ module mesh integer, dimension(:,:), allocatable :: weights end type tCellNodeDefinition + type(tCellNodeDefinition), dimension(:), allocatable :: cellNodeDefinition + real(pReal), public, protected :: & mesh_unitlength !< physical length of one unit in mesh @@ -43,19 +45,6 @@ module mesh !-------------------------------------------------------------------------------------------------- - integer, dimension(:,:), allocatable :: & - connectivity_elem - - - type(tMesh) :: theMesh - - integer,dimension(:,:,:), allocatable :: & - mesh_cell2 !< cell connectivity for each element,ip/cell - - integer, dimension(:,:,:,:), allocatable :: & - mesh_ipNeighborhood2 !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] - - public :: & mesh_init, & mesh_FEasCP @@ -98,6 +87,14 @@ subroutine mesh_init(ip,el) real(pReal), dimension(:,:), allocatable :: & ip_reshaped + integer,dimension(:,:,:), allocatable :: & + connectivity_cell !< cell connectivity for each element,ip/cell + integer, dimension(:,:), allocatable :: & + connectivity_elem + + + type(tMesh) :: theMesh + write(6,'(/,a)') ' <<<+- mesh init -+>>>' @@ -129,6 +126,9 @@ subroutine mesh_init(ip,el) call theMesh%init('mesh',elemType,mesh_node0) call theMesh%setNelems(mesh_nElems) + allocate(cellNodeDefinition(theMesh%elem%nNodes-1)) + + allocate(microstructureAt(theMesh%nElems), source=0) allocate(homogenizationAt(theMesh%nElems), source=0) @@ -148,11 +148,12 @@ subroutine mesh_init(ip,el) call results_closeJobFile #endif - call buildCells(mesh_Nnodes,theMesh%elem,connectivity_elem) + allocate(connectivity_cell(theMesh%elem%NcellNodesPerCell,theMesh%elem%nIPs,theMesh%nElems)) + call buildCells(connectivity_cell,cellNodeDefinition,& + mesh_Nnodes,theMesh%elem,connectivity_elem) allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal) - call IP_neighborhood2 if (usePingPong .and. (mesh_Nelems /= theMesh%nElems)) & call IO_error(600) ! ping-pong must be disabled when having non-DAMASK elements @@ -181,7 +182,6 @@ subroutine mesh_init(ip,el) 'nodal coordinates','m') call results_closeJobFile() #endif - call geometry_plastic_nonlocal_setIPneighborhood(mesh_ipNeighborhood2) end subroutine mesh_init @@ -680,17 +680,20 @@ subroutine mesh_marc_buildElements2(microstructureAt,homogenizationAt, & 630 end subroutine mesh_marc_buildElements2 -subroutine buildCells(nNodes,elem,connectivity_elem) - integer, intent(in) :: nNodes - type(tElement), intent(in) :: elem - integer,dimension(:,:), intent(in) :: connectivity_elem + +subroutine buildCells(connectivity_cell,cellNodeDefinition, & + nNodes,elem,connectivity_elem) + + type(tCellNodeDefinition), dimension(:), intent(out) :: cellNodeDefinition + integer,dimension(:,:,:),intent(out):: connectivity_cell + + integer, intent(in) :: nNodes + type(tElement), intent(in) :: elem + integer,dimension(:,:), intent(in) :: connectivity_elem integer,dimension(:), allocatable :: candidates_local - integer,dimension(:,:), allocatable :: parentsAndWeights,candidates_global, connectivity_cell_reshape - integer,dimension(:,:,:), allocatable :: connectivity_cell - - type(tCellNodeDefinition), dimension(:), allocatable :: cellNodeDefinition + integer,dimension(:,:), allocatable :: parentsAndWeights,candidates_global integer :: e, n, c, p, s,i,m,j,nParentNodes,nCellNode,Nelem,candidateID @@ -698,7 +701,6 @@ subroutine buildCells(nNodes,elem,connectivity_elem) !--------------------------------------------------------------------------------------------------- ! initialize global connectivity to negative local connectivity - allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,Nelem)) connectivity_cell = -spread(elem%cell,3,Nelem) ! local cell node ID !--------------------------------------------------------------------------------------------------- @@ -714,9 +716,6 @@ subroutine buildCells(nNodes,elem,connectivity_elem) enddo enddo nCellNode = nNodes - - - allocate(cellNodeDefinition(elem%nNodes-1)) !--------------------------------------------------------------------------------------------------- ! set connectivity of cell nodes that are defined by 2,...,nNodes real nodes @@ -804,11 +803,7 @@ subroutine buildCells(nNodes,elem,connectivity_elem) enddo nCellNode = nCellNode + i - - enddo - - mesh_cell2 = connectivity_cell contains !------------------------------------------------------------------------------------------------ @@ -840,83 +835,6 @@ subroutine buildCells(nNodes,elem,connectivity_elem) end subroutine buildCells -!--------------------------------------------------------------------------------------------------- -!> @brief cell neighborhood -!--------------------------------------------------------------------------------------------------- -subroutine IP_neighborhood2 - - integer, dimension(:,:), allocatable :: faces - integer, dimension(:), allocatable :: face - integer :: e,i,f,c,m,n,j,k,l,p, current, next,i2,e2,n2,k2 - logical :: match - - allocate(faces(size(theMesh%elem%cellface,1)+3,size(theMesh%elem%cellface,2)*theMesh%elem%nIPs*theMesh%Nelems)) - - ! store cell face definitions - f = 0 - do e = 1,theMesh%nElems - do i = 1,theMesh%elem%nIPs - do n = 1, theMesh%elem%nIPneighbors - f = f + 1 - face = mesh_cell2(theMesh%elem%cellFace(:,n),i,e) - storeSorted: do j = 1, size(face) - faces(j,f) = maxval(face) - face(maxloc(face)) = -huge(1) - enddo storeSorted - faces(j:j+2,f) = [e,i,n] - enddo - enddo - enddo - - ! sort .. - call math_sort(faces,sortDim=1) - do p = 2, size(faces,1)-2 - n = 1 - do while(n <= size(faces,2)) - j=0 - do while (n+j<= size(faces,2)) - if (faces(p-1,n+j)/=faces(p-1,n)) exit - j = j + 1 - enddo - e = n+j-1 - if (any(faces(p,n:e)/=faces(p,n))) call math_sort(faces(:,n:e),sortDim=p) - n = e+1 - enddo - enddo - - allocate(mesh_ipNeighborhood2(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems),source=0) - - ! find IP neighbors - f = 1 - do while(f <= size(faces,2)) - e = faces(size(theMesh%elem%cellface,1)+1,f) - i = faces(size(theMesh%elem%cellface,1)+2,f) - n = faces(size(theMesh%elem%cellface,1)+3,f) - - if (f < size(faces,2)) then - match = all(faces(1:size(theMesh%elem%cellface,1),f) == faces(1:size(theMesh%elem%cellface,1),f+1)) - e2 = faces(size(theMesh%elem%cellface,1)+1,f+1) - i2 = faces(size(theMesh%elem%cellface,1)+2,f+1) - n2 = faces(size(theMesh%elem%cellface,1)+3,f+1) - else - match = .false. - endif - - if (match) then - if (e == e2) then ! same element. MD: I don't think that we need this (not even for other elements) - k = theMesh%elem%IPneighbor(n, i) - k2 = theMesh%elem%IPneighbor(n2,i2) - endif - mesh_ipNeighborhood2(1:3,n, i, e) = [e2,i2,n2] - mesh_ipNeighborhood2(1:3,n2,i2,e2) = [e, i, n] - f = f +1 - endif - f = f +1 - enddo - -end subroutine IP_neighborhood2 - - !-------------------------------------------------------------------------------------------------- !> @brief Gives the FE to CP ID mapping by binary search through lookup array !! valid questions (what) are 'elem', 'node' From c41d8eb25777c8cbbc72dc80f14b8b51253d6ae4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Oct 2019 23:34:03 +0200 Subject: [PATCH 40/48] use new data structure to calculate cell node coordinates --- src/mesh_marc.f90 | 51 +++++++++++++++++++++++++++++++++++++---------- 1 file changed, 41 insertions(+), 10 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index a4d5c169e..6baae8638 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -65,8 +65,8 @@ subroutine mesh_init(ip,el) integer :: & mesh_NelemSets real(pReal), dimension(:,:), allocatable :: & - mesh_node, & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!! - mesh_node0 !< node x,y,z coordinates (initially!) + node0_elem, & !< node x,y,z coordinates (initially!) + node0_cell integer :: j, fileFormatVersion, elemType, & mesh_maxNelemInSet, & @@ -88,7 +88,7 @@ subroutine mesh_init(ip,el) real(pReal), dimension(:,:), allocatable :: & ip_reshaped integer,dimension(:,:,:), allocatable :: & - connectivity_cell !< cell connectivity for each element,ip/cell + connectivity_cell !< cell connectivity for each element,ip/cell integer, dimension(:,:), allocatable :: & connectivity_elem @@ -119,14 +119,13 @@ subroutine mesh_init(ip,el) allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes),source=0) call mesh_marc_map_nodes(mesh_Nnodes,inputFile) !ToDo: don't work on global variables - mesh_node0 = mesh_marc_build_nodes(mesh_Nnodes,inputFile) - mesh_node = mesh_node0 + node0_elem = mesh_marc_build_nodes(mesh_Nnodes,inputFile) elemType = mesh_marc_getElemType(mesh_nElems,FILEUNIT) - call theMesh%init('mesh',elemType,mesh_node0) + call theMesh%init('mesh',elemType,node0_elem) call theMesh%setNelems(mesh_nElems) - allocate(cellNodeDefinition(theMesh%elem%nNodes-1)) + @@ -148,9 +147,13 @@ subroutine mesh_init(ip,el) call results_closeJobFile #endif - allocate(connectivity_cell(theMesh%elem%NcellNodesPerCell,theMesh%elem%nIPs,theMesh%nElems)) + allocate(cellNodeDefinition(theMesh%elem%nNodes-1)) + allocate(connectivity_cell(theMesh%elem%NcellNodesPerCell,theMesh%elem%nIPs,theMesh%nElems)) call buildCells(connectivity_cell,cellNodeDefinition,& mesh_Nnodes,theMesh%elem,connectivity_elem) + allocate(node0_cell(3,maxval(connectivity_cell))) + call buildCellNodes(node0_cell,& + cellNodeDefinition,node0_elem) allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal) @@ -173,12 +176,12 @@ subroutine mesh_init(ip,el) ip_reshaped = reshape(mesh_ipCoordinates,[3,theMesh%elem%nIPs*theMesh%nElems]) call discretization_init(microstructureAt,homogenizationAt,& ip_reshaped,& - mesh_node0) + node0_elem) #if defined(DAMASK_HDF5) call results_openJobFile call results_writeDataset('geometry',ip_reshaped,'x_c', & 'cell center coordinates','m') - call results_writeDataset('geometry',mesh_node0,'x_n', & + call results_writeDataset('geometry',node0_elem,'x_n', & 'nodal coordinates','m') call results_closeJobFile() #endif @@ -835,6 +838,34 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, & end subroutine buildCells +subroutine buildCellNodes(node_cell,& + definition,node_elem) + + real(pReal), dimension(:,:), intent(out) :: node_cell + type(tCellNodeDefinition), dimension(:), intent(in) :: definition + real(pReal), dimension(:,:), intent(in) :: node_elem + + integer :: i, j, k, n + + n = size(node_elem,2) + node_cell(:,1:n) = node_elem + + do i = 1, size(cellNodeDefinition,1) + do j = 1, size(cellNodeDefinition(i)%parents,1) + n = n+1 + node_cell(:,n) = 0.0_pReal + do k = 1, size(cellNodeDefinition(i)%parents,2) + node_cell(:,n) = node_cell(:,n) & + + node_cell(:,definition(i)%parents(j,k)) * real(definition(i)%weights(j,k),pReal) + enddo + node_cell(:,n) = node_cell(:,n)/real(sum(definition(i)%weights(j,:)),pReal) + enddo + enddo + +end subroutine buildCellNodes + + + !-------------------------------------------------------------------------------------------------- !> @brief Gives the FE to CP ID mapping by binary search through lookup array !! valid questions (what) are 'elem', 'node' From 513b1906f6d173e0fefae30ecc10ceea67316e66 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 14 Oct 2019 08:23:21 +0200 Subject: [PATCH 41/48] bugfixes: cell definition was not stored correctly due to wrong indexing --- src/mesh_marc.f90 | 47 ++++++++++++++++++++++++++++++++++------------- 1 file changed, 34 insertions(+), 13 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 6baae8638..33e2265ed 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -121,13 +121,11 @@ subroutine mesh_init(ip,el) node0_elem = mesh_marc_build_nodes(mesh_Nnodes,inputFile) + elemType = mesh_marc_getElemType(mesh_nElems,FILEUNIT) call theMesh%init('mesh',elemType,node0_elem) call theMesh%setNelems(mesh_nElems) - - - allocate(microstructureAt(theMesh%nElems), source=0) allocate(homogenizationAt(theMesh%nElems), source=0) @@ -146,7 +144,8 @@ subroutine mesh_init(ip,el) 'connectivity of the elements','-') call results_closeJobFile #endif - + allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal) + allocate(cellNodeDefinition(theMesh%elem%nNodes-1)) allocate(connectivity_cell(theMesh%elem%NcellNodesPerCell,theMesh%elem%nIPs,theMesh%nElems)) call buildCells(connectivity_cell,cellNodeDefinition,& @@ -154,8 +153,9 @@ subroutine mesh_init(ip,el) allocate(node0_cell(3,maxval(connectivity_cell))) call buildCellNodes(node0_cell,& cellNodeDefinition,node0_elem) - - allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal) + 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)) & @@ -173,7 +173,6 @@ subroutine mesh_init(ip,el) 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" - ip_reshaped = reshape(mesh_ipCoordinates,[3,theMesh%elem%nIPs*theMesh%nElems]) call discretization_init(microstructureAt,homogenizationAt,& ip_reshaped,& node0_elem) @@ -718,6 +717,7 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, & endif realNode enddo enddo + nCellNode = nNodes !--------------------------------------------------------------------------------------------------- @@ -787,25 +787,26 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, & do while(n <= size(candidates_local)*Nelem) j=0 parentsAndWeights(:,1) = candidates_global(1:nParentNodes,n+j) - parentsAndWeights(:,2) = candidates_global(nParentNodes+1:nParentNodes*2,n+j) + parentsAndWeights(:,2) = candidates_global(nParentNodes+1:nParentNodes*2,n+j) + e = candidates_global(nParentNodes*2+1,n+j) c = candidates_global(nParentNodes*2+2,n+j) do while (n+j<= size(candidates_local)*Nelem) if (any(candidates_global(1:2*nParentNodes,n+j)/=candidates_global(1:2*nParentNodes,n))) exit where (connectivity_cell(:,:,candidates_global(nParentNodes*2+1,n+j)) == -candidates_global(nParentNodes*2+2,n+j)) ! still locally defined - connectivity_cell(:,:,candidates_global(nParentNodes*2+1,n+j)) = nCellNode + i ! gets current new cell node id + connectivity_cell(:,:,candidates_global(nParentNodes*2+1,n+j)) = nCellNode + 1 ! gets current new cell node id end where j = j+1 enddo + nCellNode = nCellNode + 1 cellNodeDefinition(nParentNodes-1)%parents(i,:) = parentsAndWeights(:,1) cellNodeDefinition(nParentNodes-1)%weights(i,:) = parentsAndWeights(:,2) - i = i+1 + i = i + 1 n = n+j - enddo - nCellNode = nCellNode + i + enddo contains @@ -849,7 +850,7 @@ subroutine buildCellNodes(node_cell,& n = size(node_elem,2) node_cell(:,1:n) = node_elem - + do i = 1, size(cellNodeDefinition,1) do j = 1, size(cellNodeDefinition(i)%parents,1) n = n+1 @@ -865,6 +866,26 @@ subroutine buildCellNodes(node_cell,& end subroutine buildCellNodes +subroutine buildIPcoordinates(IPcoordinates,& + connectivity_cell,node_cell) + + real(pReal), dimension(:,:), intent(out) :: IPcoordinates + integer,dimension(:,:), intent(in) :: connectivity_cell !< cell connectivity for each element,ip/cell + real(pReal), dimension(:,:), intent(in) :: node_cell + + integer :: i, j + + do i = 1, size(connectivity_cell,2) + IPcoordinates(:,i) = 0.0_pReal + do j = 1, size(connectivity_cell,1) + IPcoordinates(:,i) = IPcoordinates(:,i) & + + node_cell(:,connectivity_cell(j,i)) + enddo + IPcoordinates(:,i) = IPcoordinates(:,i)/real(size(connectivity_cell,1),pReal) + enddo + +end subroutine buildIPcoordinates + !-------------------------------------------------------------------------------------------------- !> @brief Gives the FE to CP ID mapping by binary search through lookup array From fae4546cfd2c2504478ad70af1c9659648656866 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 14 Oct 2019 09:37:31 +0200 Subject: [PATCH 42/48] polishing --- src/mesh_marc.f90 | 49 +++++++++++++++++++++++++++-------------------- 1 file changed, 28 insertions(+), 21 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 33e2265ed..653d81fea 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -149,7 +149,7 @@ subroutine mesh_init(ip,el) allocate(cellNodeDefinition(theMesh%elem%nNodes-1)) allocate(connectivity_cell(theMesh%elem%NcellNodesPerCell,theMesh%elem%nIPs,theMesh%nElems)) call buildCells(connectivity_cell,cellNodeDefinition,& - mesh_Nnodes,theMesh%elem,connectivity_elem) + theMesh%elem,connectivity_elem) allocate(node0_cell(3,maxval(connectivity_cell))) call buildCellNodes(node0_cell,& cellNodeDefinition,node0_elem) @@ -683,16 +683,17 @@ subroutine mesh_marc_buildElements2(microstructureAt,homogenizationAt, & - +!-------------------------------------------------------------------------------------------------- +!> @brief Calculates cell node coordinates from element node coordinates +!-------------------------------------------------------------------------------------------------- subroutine buildCells(connectivity_cell,cellNodeDefinition, & - nNodes,elem,connectivity_elem) + elem,connectivity_elem) - type(tCellNodeDefinition), dimension(:), intent(out) :: cellNodeDefinition - integer,dimension(:,:,:),intent(out):: connectivity_cell + type(tCellNodeDefinition), dimension(:), intent(out) :: cellNodeDefinition ! definition of cell nodes for increasing number of parents + integer, dimension(:,:,:),intent(out) :: connectivity_cell - integer, intent(in) :: nNodes - type(tElement), intent(in) :: elem - integer,dimension(:,:), intent(in) :: connectivity_elem + type(tElement), intent(in) :: elem ! element definition + integer, dimension(:,:), intent(in) :: connectivity_elem ! connectivity of the elements integer,dimension(:), allocatable :: candidates_local integer,dimension(:,:), allocatable :: parentsAndWeights,candidates_global @@ -718,7 +719,7 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, & enddo enddo - nCellNode = nNodes + nCellNode = maxval(connectivity_elem) !--------------------------------------------------------------------------------------------------- ! set connectivity of cell nodes that are defined by 2,...,nNodes real nodes @@ -839,17 +840,20 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, & end subroutine buildCells -subroutine buildCellNodes(node_cell,& +!-------------------------------------------------------------------------------------------------- +!> @brief Calculates cell node coordinates from element node coordinates +!-------------------------------------------------------------------------------------------------- +subroutine buildCellNodes(node_cell, & definition,node_elem) - real(pReal), dimension(:,:), intent(out) :: node_cell - type(tCellNodeDefinition), dimension(:), intent(in) :: definition - real(pReal), dimension(:,:), intent(in) :: node_elem + real(pReal), dimension(:,:), intent(out) :: node_cell !< cell node coordinates + type(tCellNodeDefinition), dimension(:), intent(in) :: definition !< cell node definition (weights and parents) + real(pReal), dimension(:,:), intent(in) :: node_elem !< element nodes integer :: i, j, k, n n = size(node_elem,2) - node_cell(:,1:n) = node_elem + node_cell(:,1:n) = node_elem !< initial nodes coincide with element nodes do i = 1, size(cellNodeDefinition,1) do j = 1, size(cellNodeDefinition(i)%parents,1) @@ -866,20 +870,23 @@ subroutine buildCellNodes(node_cell,& end subroutine buildCellNodes -subroutine buildIPcoordinates(IPcoordinates,& +!-------------------------------------------------------------------------------------------------- +!> @brief Calculates IP coordinates as center of cell +!-------------------------------------------------------------------------------------------------- +subroutine buildIPcoordinates(IPcoordinates, & connectivity_cell,node_cell) - real(pReal), dimension(:,:), intent(out) :: IPcoordinates - integer,dimension(:,:), intent(in) :: connectivity_cell !< cell connectivity for each element,ip/cell - real(pReal), dimension(:,:), intent(in) :: node_cell + real(pReal), dimension(:,:), intent(out):: IPcoordinates !< cell-center/IP coordinates + integer, dimension(:,:), intent(in) :: connectivity_cell !< connectivity for each cell + real(pReal), dimension(:,:), intent(in) :: node_cell !< cell node coordinates - integer :: i, j + integer :: i, n do i = 1, size(connectivity_cell,2) IPcoordinates(:,i) = 0.0_pReal - do j = 1, size(connectivity_cell,1) + do n = 1, size(connectivity_cell,1) IPcoordinates(:,i) = IPcoordinates(:,i) & - + node_cell(:,connectivity_cell(j,i)) + + node_cell(:,connectivity_cell(n,i)) enddo IPcoordinates(:,i) = IPcoordinates(:,i)/real(size(connectivity_cell,1),pReal) enddo From 57fef8fa574d900a459ac7e5f947a2507a8dece0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 14 Oct 2019 10:08:35 +0200 Subject: [PATCH 43/48] consistent naming --- src/mesh_marc.f90 | 82 +++++++++++++++++++++++------------------------ 1 file changed, 40 insertions(+), 42 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 653d81fea..70da8688a 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -102,27 +102,27 @@ subroutine mesh_init(ip,el) inputFile = IO_read_ASCII(trim(modelName)//trim(InputFileExtension)) ! parsing Marc input file - fileFormatVersion = mesh_marc_get_fileFormat(inputFile) - call mesh_marc_get_tableStyles(initialcondTableStyle,hypoelasticTableStyle,inputFile) + fileFormatVersion = inputRead_fileFormat(inputFile) + call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle,inputFile) if (fileFormatVersion > 12) & - marc_matNumber = mesh_marc_get_matNumber(hypoelasticTableStyle,inputFile) - call mesh_marc_count_nodesAndElements(mesh_nNodes, mesh_nElems, inputFile) + marc_matNumber = inputRead_matNumber(hypoelasticTableStyle,inputFile) + call inputRead_NnodesAndElements(mesh_nNodes, mesh_nElems, inputFile) call IO_open_inputFile(FILEUNIT,modelName) - call mesh_marc_count_elementSets(mesh_NelemSets,mesh_maxNelemInSet,FILEUNIT) + 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 mesh_marc_map_elementSets(mesh_nameElemSet,mesh_mapElemSet,FILEUNIT) + call inputRead_mapElemSets(mesh_nameElemSet,mesh_mapElemSet,FILEUNIT) allocate (mesh_mapFEtoCPelem(2,mesh_nElems), source = 0) - call mesh_marc_map_elements(hypoelasticTableStyle,mesh_nameElemSet,mesh_mapElemSet,& + call inputRead_mapElems(hypoelasticTableStyle,mesh_nameElemSet,mesh_mapElemSet,& mesh_nElems,fileFormatVersion,marc_matNumber,FILEUNIT) allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes),source=0) - call mesh_marc_map_nodes(mesh_Nnodes,inputFile) !ToDo: don't work on global variables + call inputRead_mapNodes(mesh_Nnodes,inputFile) !ToDo: don't work on global variables - node0_elem = mesh_marc_build_nodes(mesh_Nnodes,inputFile) + node0_elem = inputRead_elemNodes(mesh_Nnodes,inputFile) - elemType = mesh_marc_getElemType(mesh_nElems,FILEUNIT) + elemType = inputRead_elemType(mesh_nElems,FILEUNIT) call theMesh%init('mesh',elemType,node0_elem) call theMesh%setNelems(mesh_nElems) @@ -137,13 +137,6 @@ subroutine mesh_init(ip,el) close (FILEUNIT) -#if defined(DAMASK_HDF5) - call results_openJobFile - call HDF5_closeGroup(results_addGroup('geometry')) - call results_writeDataset('geometry',connectivity_elem,'C',& - 'connectivity of the elements','-') - call results_closeJobFile -#endif allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal) allocate(cellNodeDefinition(theMesh%elem%nNodes-1)) @@ -178,6 +171,9 @@ subroutine mesh_init(ip,el) node0_elem) #if defined(DAMASK_HDF5) call results_openJobFile + call HDF5_closeGroup(results_addGroup('geometry')) + call results_writeDataset('geometry',connectivity_elem,'C',& + 'connectivity of the elements','-') call results_writeDataset('geometry',ip_reshaped,'x_c', & 'cell center coordinates','m') call results_writeDataset('geometry',node0_elem,'x_n', & @@ -191,7 +187,7 @@ end subroutine mesh_init !-------------------------------------------------------------------------------------------------- !> @brief Figures out version of Marc input file format !-------------------------------------------------------------------------------------------------- -integer function mesh_marc_get_fileFormat(fileContent) +integer function inputRead_fileFormat(fileContent) character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines @@ -201,18 +197,18 @@ integer function mesh_marc_get_fileFormat(fileContent) do l = 1, size(fileContent) chunkPos = IO_stringPos(fileContent(l)) if ( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'version') then - mesh_marc_get_fileFormat = IO_intValue(fileContent(l),chunkPos,2) + inputRead_fileFormat = IO_intValue(fileContent(l),chunkPos,2) exit endif enddo -end function mesh_marc_get_fileFormat +end function inputRead_fileFormat !-------------------------------------------------------------------------------------------------- !> @brief Figures out table styles for initial cond and hypoelastic !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_get_tableStyles(initialcond,hypoelastic,fileContent) +subroutine inputRead_tableStyles(initialcond,hypoelastic,fileContent) integer, intent(out) :: initialcond, hypoelastic character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines @@ -232,18 +228,18 @@ subroutine mesh_marc_get_tableStyles(initialcond,hypoelastic,fileContent) endif enddo -end subroutine mesh_marc_get_tableStyles +end subroutine inputRead_tableStyles !-------------------------------------------------------------------------------------------------- !> @brief Figures out material number of hypoelastic material !-------------------------------------------------------------------------------------------------- -function mesh_marc_get_matNumber(tableStyle,fileContent) +function inputRead_matNumber(tableStyle,fileContent) integer, intent(in) :: tableStyle character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines - integer, dimension(:), allocatable :: mesh_marc_get_matNumber + integer, dimension(:), allocatable :: inputRead_matNumber integer, allocatable, dimension(:) :: chunkPos integer :: i, j, data_blocks, l @@ -257,23 +253,24 @@ function mesh_marc_get_matNumber(tableStyle,fileContent) else data_blocks = 1 endif - allocate(mesh_marc_get_matNumber(data_blocks), source = 0) + allocate(inputRead_matNumber(data_blocks), source = 0) do i = 0, data_blocks - 1 j = i*(2+tableStyle) + 1 chunkPos = IO_stringPos(fileContent(l+1+j)) - mesh_marc_get_matNumber(i+1) = IO_intValue(fileContent(l+1+j),chunkPos,1) + inputRead_matNumber(i+1) = IO_intValue(fileContent(l+1+j),chunkPos,1) enddo exit endif enddo -end function mesh_marc_get_matNumber +end function inputRead_matNumber !-------------------------------------------------------------------------------------------------- !> @brief Count overall number of nodes and elements !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_count_nodesAndElements(nNodes, nElems, fileContent) +subroutine inputRead_NnodesAndElements(nNodes,nElems,& + fileContent) integer, intent(out) :: nNodes, nElems character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines @@ -294,13 +291,14 @@ subroutine mesh_marc_count_nodesAndElements(nNodes, nElems, fileContent) endif enddo -end subroutine mesh_marc_count_nodesAndElements +end subroutine inputRead_NnodesAndElements !-------------------------------------------------------------------------------------------------- !> @brief Count overall number of element sets in mesh. !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_count_elementSets(nElemSets,maxNelemInSet,fileUnit) +subroutine inputRead_NelemSets(nElemSets,maxNelemInSet,& + fileUnit) integer, intent(out) :: nElemSets, maxNelemInSet integer, intent(in) :: fileUnit @@ -323,13 +321,13 @@ subroutine mesh_marc_count_elementSets(nElemSets,maxNelemInSet,fileUnit) endif enddo -620 end subroutine mesh_marc_count_elementSets +620 end subroutine inputRead_NelemSets !-------------------------------------------------------------------------------------------------- !> @brief map element sets !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_map_elementSets(nameElemSet,mapElemSet,fileUnit) +subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,fileUnit) character(len=64), dimension(:), intent(out) :: nameElemSet integer, dimension(:,:), intent(out) :: mapElemSet @@ -353,14 +351,14 @@ subroutine mesh_marc_map_elementSets(nameElemSet,mapElemSet,fileUnit) endif enddo -620 end subroutine mesh_marc_map_elementSets +620 end subroutine inputRead_mapElemSets !-------------------------------------------------------------------------------------------------- !> @brief Maps elements from FE ID to internal (consecutive) representation. !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_map_elements(tableStyle,nameElemSet,mapElemSet,nElems,fileFormatVersion,matNumber,fileUnit) +subroutine inputRead_mapElems(tableStyle,nameElemSet,mapElemSet,nElems,fileFormatVersion,matNumber,fileUnit) integer, intent(in) :: fileUnit,tableStyle,nElems,fileFormatVersion integer, dimension(:), intent(in) :: matNumber @@ -418,13 +416,13 @@ subroutine mesh_marc_map_elements(tableStyle,nameElemSet,mapElemSet,nElems,fileF call math_sort(mesh_mapFEtoCPelem) -end subroutine mesh_marc_map_elements +end subroutine inputRead_mapElems !-------------------------------------------------------------------------------------------------- !> @brief Maps node from FE ID to internal (consecutive) representation. !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_map_nodes(nNodes,fileContent) +subroutine inputRead_mapNodes(nNodes,fileContent) integer, intent(in) :: nNodes character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines @@ -445,13 +443,13 @@ subroutine mesh_marc_map_nodes(nNodes,fileContent) call math_sort(mesh_mapFEtoCPnode) -end subroutine mesh_marc_map_nodes +end subroutine inputRead_mapNodes !-------------------------------------------------------------------------------------------------- !> @brief store x,y,z coordinates of all nodes in mesh. !-------------------------------------------------------------------------------------------------- -function mesh_marc_build_nodes(nNode,fileContent) result(nodes) +function inputRead_elemNodes(nNode,fileContent) result(nodes) integer, intent(in) :: nNode character(len=pStringLen), dimension(:), intent(in) :: fileContent !< file content, separated per lines @@ -474,13 +472,13 @@ function mesh_marc_build_nodes(nNode,fileContent) result(nodes) endif enddo -end function mesh_marc_build_nodes +end function inputRead_elemNodes !-------------------------------------------------------------------------------------------------- !> @brief Gets element type (and checks if the whole mesh comprises of only one type) !-------------------------------------------------------------------------------------------------- -integer function mesh_marc_getElemType(nElem,fileUnit) +integer function inputRead_elemType(nElem,fileUnit) integer, intent(in) :: & nElem, & @@ -505,7 +503,7 @@ integer function mesh_marc_getElemType(nElem,fileUnit) if (t == -1) then t = mapElemtype(IO_stringValue(line,chunkPos,2)) call tempEl%init(t) - mesh_marc_getElemType = t + inputRead_elemType = t else if (t /= mapElemtype(IO_stringValue(line,chunkPos,2))) call IO_error(191,el=t,ip=i) endif @@ -561,7 +559,7 @@ integer function mesh_marc_getElemType(nElem,fileUnit) end function mapElemtype -620 end function mesh_marc_getElemType +620 end function inputRead_elemType !-------------------------------------------------------------------------------------------------- From 57af822396c6de440b08e98fc446feaf25cbf778 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 14 Oct 2019 10:36:59 +0200 Subject: [PATCH 44/48] better readable --- src/mesh_marc.f90 | 66 +++++++++++++++++++++++++++++++++++------------ 1 file changed, 50 insertions(+), 16 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 70da8688a..0e5607028 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -130,8 +130,8 @@ subroutine mesh_init(ip,el) allocate(microstructureAt(theMesh%nElems), source=0) allocate(homogenizationAt(theMesh%nElems), source=0) - connectivity_elem = mesh_marc_buildElements(mesh_nElems,theMesh%elem%nNodes,FILEUNIT) - call mesh_marc_buildElements2(microstructureAt,homogenizationAt, & + connectivity_elem = inputRead_connectivityElem(mesh_nElems,theMesh%elem%nNodes,FILEUNIT) + call inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, & mesh_nElems,theMesh%elem%nNodes,mesh_nameElemSet,mesh_mapElemSet,& initialcondTableStyle,FILEUNIT) close (FILEUNIT) @@ -169,19 +169,53 @@ subroutine mesh_init(ip,el) call discretization_init(microstructureAt,homogenizationAt,& ip_reshaped,& node0_elem) + + call writeGeometry(0,connectivity_elem,& + reshape(connectivity_cell,[theMesh%elem%NcellNodesPerCell,theMesh%elem%nIPs*theMesh%nElems]),& + node0_cell,ip_reshaped) + +end subroutine mesh_init + + +subroutine writeGeometry(elemType,connectivity_elem,connectivity_cell,coordinates_nodes,coordinates_points) + + integer, intent(in) :: elemType + integer, dimension(:,:), intent(in) :: & + connectivity_elem, & + connectivity_cell + real(pReal), dimension(:,:), intent(in) :: & + coordinates_nodes, & + coordinates_points + + integer, dimension(:,:), allocatable :: & + connectivity_temp + real(pReal), dimension(:,:), allocatable :: & + coordinates_temp + #if defined(DAMASK_HDF5) call results_openJobFile call HDF5_closeGroup(results_addGroup('geometry')) - call results_writeDataset('geometry',connectivity_elem,'C',& + + connectivity_temp = connectivity_elem + call results_writeDataset('geometry',connectivity_temp,'T_e',& 'connectivity of the elements','-') - call results_writeDataset('geometry',ip_reshaped,'x_c', & - 'cell center coordinates','m') - call results_writeDataset('geometry',node0_elem,'x_n', & - 'nodal coordinates','m') + + connectivity_temp = connectivity_cell + call results_writeDataset('geometry',connectivity_temp,'T_c', & + 'connectivity of the cells','-') + + coordinates_temp = coordinates_nodes + call results_writeDataset('geometry',coordinates_temp,'x_n', & + 'coordinates of the nodes','m') + + coordinates_temp = coordinates_points + call results_writeDataset('geometry',coordinates_temp,'x_p', & + 'coordinates of the material points','m') + call results_closeJobFile() -#endif +#endif -end subroutine mesh_init +end subroutine writeGeometry !-------------------------------------------------------------------------------------------------- @@ -565,7 +599,7 @@ end function mapElemtype !-------------------------------------------------------------------------------------------------- !> @brief Stores node IDs !-------------------------------------------------------------------------------------------------- -function mesh_marc_buildElements(nElem,nNodes,fileUnit) +function inputRead_connectivityElem(nElem,nNodes,fileUnit) integer, intent(in) :: & nElem, & @@ -573,7 +607,7 @@ function mesh_marc_buildElements(nElem,nNodes,fileUnit) fileUnit integer, dimension(nElem,nNodes) :: & - mesh_marc_buildElements + inputRead_connectivityElem integer, allocatable, dimension(:) :: chunkPos character(len=300) line @@ -594,7 +628,7 @@ function mesh_marc_buildElements(nElem,nNodes,fileUnit) if (e /= 0) then ! disregard non CP elems nNodesAlreadyRead = 0 do j = 1,chunkPos(1)-2 - mesh_marc_buildElements(j,e) = & + inputRead_connectivityElem(j,e) = & mesh_FEasCP('node',IO_IntValue(line,chunkPos,j+2)) enddo nNodesAlreadyRead = chunkPos(1) - 2 @@ -602,7 +636,7 @@ function mesh_marc_buildElements(nElem,nNodes,fileUnit) read (fileUnit,'(A300)',END=620) line chunkPos = IO_stringPos(line) do j = 1,chunkPos(1) - mesh_marc_buildElements(nNodesAlreadyRead+j,e) = & + inputRead_connectivityElem(nNodesAlreadyRead+j,e) = & mesh_FEasCP('node',IO_IntValue(line,chunkPos,j)) enddo nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1) @@ -613,13 +647,13 @@ function mesh_marc_buildElements(nElem,nNodes,fileUnit) endif enddo -620 end function mesh_marc_buildElements +620 end function inputRead_connectivityElem !-------------------------------------------------------------------------------------------------- !> @brief Stores homogenization and microstructure ID !-------------------------------------------------------------------------------------------------- -subroutine mesh_marc_buildElements2(microstructureAt,homogenizationAt, & +subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, & nElem,nNodes,nameElemSet,mapElemSet,initialcondTableStyle,fileUnit) integer, dimension(:), intent(out) :: & @@ -677,7 +711,7 @@ subroutine mesh_marc_buildElements2(microstructureAt,homogenizationAt, & endif enddo -630 end subroutine mesh_marc_buildElements2 +630 end subroutine inputRead_microstructureAndHomogenization From dcbd7624dd9d9f6c7108399184f96282d9f74018 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 14 Oct 2019 11:25:48 +0200 Subject: [PATCH 45/48] array order was wrong --- src/mesh_marc.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 0e5607028..8109806ba 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -606,7 +606,7 @@ function inputRead_connectivityElem(nElem,nNodes,fileUnit) nNodes, & !< number of nodes per element fileUnit - integer, dimension(nElem,nNodes) :: & + integer, dimension(nNodes,nElem) :: & inputRead_connectivityElem integer, allocatable, dimension(:) :: chunkPos @@ -732,7 +732,7 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, & integer :: e, n, c, p, s,i,m,j,nParentNodes,nCellNode,Nelem,candidateID - Nelem = size(connectivity_elem,1) + Nelem = size(connectivity_elem,2) !--------------------------------------------------------------------------------------------------- ! initialize global connectivity to negative local connectivity From d73e653f501b4c9784deffade566e404331373cc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 14 Oct 2019 14:13:03 +0200 Subject: [PATCH 46/48] support for unstructured grid (FEM) --- processing/post/DADF5_vtk_cells.py | 35 +++++++++++++++++++----------- 1 file changed, 22 insertions(+), 13 deletions(-) diff --git a/processing/post/DADF5_vtk_cells.py b/processing/post/DADF5_vtk_cells.py index 473e887f8..915479e26 100755 --- a/processing/post/DADF5_vtk_cells.py +++ b/processing/post/DADF5_vtk_cells.py @@ -3,6 +3,7 @@ import os import argparse +import h5py import numpy as np import vtk from vtk.util import numpy_support @@ -41,20 +42,29 @@ for filename in options.filenames: results = damask.DADF5(filename) if results.structured: # for grid solvers use rectilinear grid - rGrid = vtk.vtkRectilinearGrid() + grid = vtk.vtkRectilineagrid() coordArray = [vtk.vtkDoubleArray(), vtk.vtkDoubleArray(), vtk.vtkDoubleArray(), ] - rGrid.SetDimensions(*(results.grid+1)) + grid.SetDimensions(*(results.grid+1)) for dim in [0,1,2]: for c in np.linspace(0,results.size[dim],1+results.grid[dim]): coordArray[dim].InsertNextValue(c) - rGrid.SetXCoordinates(coordArray[0]) - rGrid.SetYCoordinates(coordArray[1]) - rGrid.SetZCoordinates(coordArray[2]) + grid.SetXCoordinates(coordArray[0]) + grid.SetYCoordinates(coordArray[1]) + grid.SetZCoordinates(coordArray[2]) + else: + nodes = vtk.vtkPoints() + with h5py.File(filename) as f: + nodes.SetData(numpy_support.numpy_to_vtk(f['/geometry/x_n'][()],deep=True)) + grid = vtk.vtkUnstructuredGrid() + grid.SetPoints(nodes) + grid.Allocate(f['/geometry/T_c'].shape[0]) + for i in f['/geometry/T_c']: + grid.InsertNextCell(vtk.VTK_HEXAHEDRON,8,i-1) for i,inc in enumerate(results.iter_visible('increments')): @@ -75,7 +85,7 @@ for filename in options.filenames: shape = [array.shape[0],np.product(array.shape[1:])] vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) - rGrid.GetCellData().AddArray(vtk_data[-1]) + grid.GetCellData().AddArray(vtk_data[-1]) else: x = results.get_dataset_location(label) if len(x) == 0: @@ -84,7 +94,7 @@ for filename in options.filenames: shape = [array.shape[0],np.product(array.shape[1:])] vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) - rGrid.GetCellData().AddArray(vtk_data[-1]) + grid.GetCellData().AddArray(vtk_data[-1]) results.set_visible('constituents', False) results.set_visible('materialpoints',True) @@ -99,7 +109,7 @@ for filename in options.filenames: shape = [array.shape[0],np.product(array.shape[1:])] vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) - rGrid.GetCellData().AddArray(vtk_data[-1]) + grid.GetCellData().AddArray(vtk_data[-1]) else: x = results.get_dataset_location(label) if len(x) == 0: @@ -108,10 +118,10 @@ for filename in options.filenames: shape = [array.shape[0],np.product(array.shape[1:])] vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) - rGrid.GetCellData().AddArray(vtk_data[-1]) + grid.GetCellData().AddArray(vtk_data[-1]) - if results.structured: - writer = vtk.vtkXMLRectilinearGridWriter() + writer = vtk.vtkXMLRectilineagridWriter() if results.structured else \ + vtk.vtkXMLUnstructuredGridWriter() dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir)) @@ -122,7 +132,6 @@ for filename in options.filenames: writer.SetCompressorTypeToZLib() writer.SetDataModeToBinary() writer.SetFileName(os.path.join(dirname,file_out)) - if results.structured: - writer.SetInputData(rGrid) + writer.SetInputData(grid) writer.Write() From 1726d95a56f951a0296bf10709d88cf8b51f3109 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 14 Oct 2019 14:39:53 +0200 Subject: [PATCH 47/48] wrong macro was defined --- src/plastic_isotropic.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index 2e9a6b57e..18b6f8e86 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -434,7 +434,7 @@ end function plastic_isotropic_postResults !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- subroutine plastic_isotropic_results(instance,group) -#if defined(PETSc) || defined(DAMASKHDF5) +#if defined(PETSc) || defined(DAMASK_HDF5) integer, intent(in) :: instance character(len=*), intent(in) :: group From 9ee709d2147c6ce4a4981f9fcc0a3185fdd05a4e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 14 Oct 2019 15:09:16 +0200 Subject: [PATCH 48/48] polishing --- src/mesh_marc.f90 | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 8109806ba..482cde16c 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -176,8 +176,12 @@ subroutine mesh_init(ip,el) end subroutine mesh_init - -subroutine writeGeometry(elemType,connectivity_elem,connectivity_cell,coordinates_nodes,coordinates_points) +!-------------------------------------------------------------------------------------------------- +!> @brief Writes all information needed for the DADF5 geometry +!-------------------------------------------------------------------------------------------------- +subroutine writeGeometry(elemType, & + connectivity_elem,connectivity_cell, & + coordinates_nodes,coordinates_points) integer, intent(in) :: elemType integer, dimension(:,:), intent(in) :: &