From 344e6ca51a32ffd0e2ab3f93901fee8520e3237b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 12 May 2019 23:56:44 +0200 Subject: [PATCH 01/13] IP volume is trivial for spectral solver --- src/mesh_grid.f90 | 79 +++-------------------------------------------- 1 file changed, 4 insertions(+), 75 deletions(-) diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index 7274582a2..f8a874c55 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -127,7 +127,6 @@ integer(pInt), dimension(:,:), allocatable, private :: & mesh_spectral_build_elements, & mesh_spectral_build_ipNeighborhood, & mesh_build_cellnodes, & - mesh_build_ipVolumes, & mesh_build_ipCoordinates type, public, extends(tMesh) :: tMesh_grid @@ -190,7 +189,7 @@ subroutine mesh_init(ip,el) implicit none include 'fftw3-mpi.f03' integer(C_INTPTR_T) :: devNull, local_K, local_K_offset - integer :: ierr, worldsize + integer :: ierr, worldsize, i integer(pInt), intent(in), optional :: el, ip integer(pInt) :: j logical :: myDebug @@ -244,7 +243,7 @@ subroutine mesh_init(ip,el) if (myDebug) write(6,'(a)') ' Built cell nodes'; flush(6) call mesh_build_ipCoordinates if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6) - call mesh_build_ipVolumes + allocate(mesh_ipVolume(1,theMesh%nElems),source=product([geomSize(1:2),size3]/real([grid(1:2),grid3]))) if (myDebug) write(6,'(a)') ' Built IP volumes'; flush(6) call mesh_build_ipAreas if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) @@ -508,7 +507,7 @@ subroutine mesh_spectral_build_ipNeighborhood integer(pInt) :: & x,y,z, & e - allocate(mesh_ipNeighborhood(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems),source=0_pInt) + allocate(mesh_ipNeighborhood(3,6,1,theMesh%nElems),source=0_pInt) e = 0_pInt do z = 0_pInt,grid3-1_pInt @@ -768,75 +767,6 @@ function mesh_build_cellnodes(nodes,Ncellnodes) end function mesh_build_cellnodes -!-------------------------------------------------------------------------------------------------- -!> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume' -!> @details The IP volume is calculated differently depending on the cell type. -!> 2D cells assume an element depth of one in order to calculate the volume. -!> For the hexahedral cell we subdivide the cell into subvolumes of pyramidal -!> shape with a cell face as basis and the central ip at the tip. This subvolume is -!> calculated as an average of four tetrahedals with three corners on the cell face -!> and one corner at the central ip. -!-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipVolumes - use math, only: & - math_volTetrahedron, & - math_areaTriangle - - implicit none - integer(pInt) :: e,t,g,c,i,m,f,n - real(pReal), dimension(FE_maxNcellnodesPerCellface,FE_maxNcellfaces) :: subvolume - - - allocate(mesh_ipVolume(theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal) - - - !$OMP PARALLEL DO PRIVATE(t,g,c,m,subvolume) - do e = 1_pInt,theMesh%nElems ! loop over cpElems - select case (theMesh%elem%cellType) - - case (1_pInt) ! 2D 3node - forall (i = 1_pInt:theMesh%elem%nIPs) & ! loop over ips=cells in this element - mesh_ipVolume(i,e) = math_areaTriangle(mesh_cellnode(1:3,mesh_cell(1,i,e)), & - mesh_cellnode(1:3,mesh_cell(2,i,e)), & - mesh_cellnode(1:3,mesh_cell(3,i,e))) - - case (2_pInt) ! 2D 4node - forall (i = 1_pInt:theMesh%elem%nIPs) & ! loop over ips=cells in this element - mesh_ipVolume(i,e) = math_areaTriangle(mesh_cellnode(1:3,mesh_cell(1,i,e)), & ! here we assume a planar shape, so division in two triangles suffices - mesh_cellnode(1:3,mesh_cell(2,i,e)), & - mesh_cellnode(1:3,mesh_cell(3,i,e))) & - + math_areaTriangle(mesh_cellnode(1:3,mesh_cell(3,i,e)), & - mesh_cellnode(1:3,mesh_cell(4,i,e)), & - mesh_cellnode(1:3,mesh_cell(1,i,e))) - - case (3_pInt) ! 3D 4node - forall (i = 1_pInt:theMesh%elem%nIPs) & ! loop over ips=cells in this element - mesh_ipVolume(i,e) = math_volTetrahedron(mesh_cellnode(1:3,mesh_cell(1,i,e)), & - mesh_cellnode(1:3,mesh_cell(2,i,e)), & - mesh_cellnode(1:3,mesh_cell(3,i,e)), & - mesh_cellnode(1:3,mesh_cell(4,i,e))) - - case (4_pInt) - c = theMesh%elem%cellType ! 3D 8node - m = FE_NcellnodesPerCellface(c) - do i = 1_pInt,theMesh%elem%nIPs ! loop over ips=cells in this element - subvolume = 0.0_pReal - forall(f = 1_pInt:FE_NipNeighbors(c), n = 1_pInt:FE_NcellnodesPerCellface(c)) & - subvolume(n,f) = math_volTetrahedron(& - mesh_cellnode(1:3,mesh_cell(FE_cellface( n ,f,c),i,e)), & - mesh_cellnode(1:3,mesh_cell(FE_cellface(1+mod(n ,m),f,c),i,e)), & - mesh_cellnode(1:3,mesh_cell(FE_cellface(1+mod(n+1,m),f,c),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 - - end select - enddo - !$OMP END PARALLEL DO - -end subroutine mesh_build_ipVolumes - - !-------------------------------------------------------------------------------------------------- !> @brief Calculates IP Coordinates. Allocates global array 'mesh_ipCoordinates' ! Called by all solvers in mesh_init in order to initialize the ip coordinates. @@ -855,8 +785,7 @@ subroutine mesh_build_ipCoordinates integer(pInt) :: e,c,i,n real(pReal), dimension(3) :: myCoords - if (.not. allocated(mesh_ipCoordinates)) & - allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal) + allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal) !$OMP PARALLEL DO PRIVATE(c,myCoords) do e = 1_pInt,theMesh%nElems ! loop over cpElems From df01a3ff36ac227b41bef9ce13263835d0fcecda Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 13 May 2019 00:00:14 +0200 Subject: [PATCH 02/13] only for internal use --- src/mesh_grid.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index f8a874c55..094f7da4e 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -34,7 +34,7 @@ module mesh real(pReal), public, protected :: & mesh_unitlength !< physical length of one unit in mesh - real(pReal), dimension(:,:), allocatable, public :: & + real(pReal), dimension(:,:), allocatable, private :: & mesh_node, & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) mesh_cellnode !< cell node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) From 4c148f1a07f8cbb5e8d3b546c09ea06e4f91f448 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 13 May 2019 08:47:53 +0200 Subject: [PATCH 03/13] not needed --- src/mesh_grid.f90 | 38 +++----------------------------------- 1 file changed, 3 insertions(+), 35 deletions(-) diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index 094f7da4e..11afe0175 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -69,8 +69,6 @@ integer(pInt), dimension(:,:), allocatable, private :: & integer(pInt), parameter, private :: & FE_Ngeomtypes = 10_pInt, & FE_Ncelltypes = 4_pInt, & - FE_maxNmatchingNodesPerFace = 4_pInt, & - FE_maxNfaces = 6_pInt, & FE_maxNcellnodesPerCell = 8_pInt, & FE_maxNcellfaces = 6_pInt, & FE_maxNcellnodesPerCellface = 4_pInt @@ -116,8 +114,7 @@ integer(pInt), dimension(:,:), allocatable, private :: & size3offset !< (local) size offset in 3rd direction public :: & - mesh_init, & - mesh_cellCenterCoordinates + mesh_init private :: & mesh_build_cellconnectivity, & @@ -460,12 +457,8 @@ end subroutine mesh_spectral_build_nodes !-------------------------------------------------------------------------------------------------- !> @brief Store FEid, type, material, texture, and node list per element. !! Allocates global array 'mesh_element' -!> @todo does the IO_error makes sense? !-------------------------------------------------------------------------------------------------- subroutine mesh_spectral_build_elements() - use IO, only: & - IO_error - implicit none integer(pInt) :: & e, & elemOffset @@ -474,11 +467,9 @@ subroutine mesh_spectral_build_elements() allocate(mesh_element (4_pInt+8_pInt,theMesh%nElems), source = 0_pInt) elemOffset = product(grid(1:2))*grid3Offset - e = 0_pInt - do while (e < theMesh%nElems) ! fill expected number of elements, stop at end of data - e = e+1_pInt ! valid element entry + do e=1, theMesh%nElems mesh_element( 1,e) = -1_pInt ! DEPRECATED - mesh_element( 2,e) = 10_pInt + mesh_element( 2,e) = -1_pInt ! DEPRECATED mesh_element( 3,e) = mesh_homogenizationAt(e) mesh_element( 4,e) = microGlobal(e+elemOffset) ! microstructure mesh_element( 5,e) = e + (e-1_pInt)/grid(1) + & @@ -492,8 +483,6 @@ subroutine mesh_spectral_build_elements() mesh_element(12,e) = mesh_element(9,e) + grid(1) + 1_pInt enddo - if (e /= theMesh%nElems) call IO_error(880_pInt,e) - end subroutine mesh_spectral_build_elements @@ -803,27 +792,6 @@ subroutine mesh_build_ipCoordinates end subroutine mesh_build_ipCoordinates -!-------------------------------------------------------------------------------------------------- -!> @brief Calculates cell center coordinates. -!-------------------------------------------------------------------------------------------------- -pure function mesh_cellCenterCoordinates(ip,el) - - implicit none - integer(pInt), intent(in) :: el, & !< element number - ip !< integration point number - real(pReal), dimension(3) :: mesh_cellCenterCoordinates !< x,y,z coordinates of the cell center of the requested IP cell - integer(pInt) :: c,n - - c = theMesh%elem%cellType - mesh_cellCenterCoordinates = 0.0_pReal - do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell - mesh_cellCenterCoordinates = mesh_cellCenterCoordinates + mesh_cellnode(1:3,mesh_cell(n,ip,el)) - enddo - mesh_cellCenterCoordinates = mesh_cellCenterCoordinates / real(FE_NcellnodesPerCell(c),pReal) - -end function mesh_cellCenterCoordinates - - !-------------------------------------------------------------------------------------------------- !> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal' !-------------------------------------------------------------------------------------------------- From fb49acdb976c48fb088a0a057fcc3c7b46eda5de Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 13 May 2019 22:54:18 +0200 Subject: [PATCH 04/13] more explicit --- src/mesh_grid.f90 | 36 +++++++++++++++--------------------- 1 file changed, 15 insertions(+), 21 deletions(-) diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index 11afe0175..39e0111c1 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -221,7 +221,8 @@ subroutine mesh_init(ip,el) mesh_Nnodes = product(grid(1:2) + 1_pInt)*(grid3 + 1_pInt) - call mesh_spectral_build_nodes() + mesh_node0 = mesh_spectral_build_nodes() + mesh_node = mesh_node0 if (myDebug) write(6,'(a)') ' Built nodes'; flush(6) call theMesh%init(mesh_node) @@ -429,29 +430,22 @@ end subroutine mesh_spectral_read_grid !> @brief Store x,y,z coordinates of all nodes in mesh. !! Allocates global arrays 'mesh_node0' and 'mesh_node' !-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_build_nodes() +pure function mesh_spectral_build_nodes() - implicit none - integer(pInt) :: n + real(pReal), dimension(3,mesh_Nnodes) :: mesh_spectral_build_nodes + integer :: n,a,b,c - allocate (mesh_node0 (3,mesh_Nnodes), source = 0.0_pReal) + n = 0 + do c = 0, grid3 + do b = 0, grid(2) + do a = 0, grid(1) + n = n + 1 + mesh_spectral_build_nodes(1:3,n) = geomSize/real(grid,pReal) * real([a,b,grid3Offset+c],pReal) + enddo + enddo + enddo - forall (n = 0_pInt:mesh_Nnodes-1_pInt) - mesh_node0(1,n+1_pInt) = mesh_unitlength * & - geomSize(1)*real(mod(n,(grid(1)+1_pInt) ),pReal) & - / real(grid(1),pReal) - mesh_node0(2,n+1_pInt) = mesh_unitlength * & - geomSize(2)*real(mod(n/(grid(1)+1_pInt),(grid(2)+1_pInt)),pReal) & - / real(grid(2),pReal) - mesh_node0(3,n+1_pInt) = mesh_unitlength * & - size3*real(mod(n/(grid(1)+1_pInt)/(grid(2)+1_pInt),(grid3+1_pInt)),pReal) & - / real(grid3,pReal) + & - size3offset - end forall - - mesh_node = mesh_node0 - -end subroutine mesh_spectral_build_nodes +end function mesh_spectral_build_nodes !-------------------------------------------------------------------------------------------------- From c8794af3bbbf5a771ee703fc6e82e505881055dd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 13 May 2019 23:18:02 +0200 Subject: [PATCH 05/13] better readable --- src/mesh_grid.f90 | 92 ++++++++++++++--------------------------------- 1 file changed, 27 insertions(+), 65 deletions(-) diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index 39e0111c1..e752f5d41 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -239,7 +239,7 @@ subroutine mesh_init(ip,el) if (myDebug) write(6,'(a)') ' Built cell connectivity'; flush(6) mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes) if (myDebug) write(6,'(a)') ' Built cell nodes'; flush(6) - call mesh_build_ipCoordinates + mesh_ipCoordinates = mesh_build_ipCoordinates() if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6) allocate(mesh_ipVolume(1,theMesh%nElems),source=product([geomSize(1:2),size3]/real([grid(1:2),grid3]))) if (myDebug) write(6,'(a)') ' Built IP volumes'; flush(6) @@ -391,7 +391,7 @@ subroutine mesh_spectral_read_grid() allocate(mesh_homogenizationAt(product(grid)), source = h) ! too large in case of MPI (shrink later, not very elegant) !-------------------------------------------------------------------------------------------------- -! read and interprete content +! read and interpret content e = 1_pInt do while (startPos < len(rawData)) endPos = startPos + index(rawData(startPos:),new_line('')) - 1_pInt @@ -426,10 +426,9 @@ subroutine mesh_spectral_read_grid() end subroutine mesh_spectral_read_grid -!-------------------------------------------------------------------------------------------------- -!> @brief Store x,y,z coordinates of all nodes in mesh. -!! Allocates global arrays 'mesh_node0' and 'mesh_node' -!-------------------------------------------------------------------------------------------------- +!--------------------------------------------------------------------------------------------------- +!> @brief Calculates position of nodes (pretend to be an element) +!--------------------------------------------------------------------------------------------------- pure function mesh_spectral_build_nodes() real(pReal), dimension(3,mesh_Nnodes) :: mesh_spectral_build_nodes @@ -448,6 +447,27 @@ pure function mesh_spectral_build_nodes() end function mesh_spectral_build_nodes +!--------------------------------------------------------------------------------------------------- +!> @brief Calculates position of IPs/cell centres (pretend to be an element) +!--------------------------------------------------------------------------------------------------- +function mesh_build_ipCoordinates() + + real(pReal), dimension(3,1,theMesh%nElems) :: mesh_build_ipCoordinates + integer :: n,a,b,c + + n = 0 + do c = 1, grid3 + do b = 1, grid(2) + do a = 1, grid(1) + n = n + 1 + mesh_build_ipCoordinates(1:3,1,n) = geomSize/real(grid,pReal) * (real([a,b,grid3Offset+c],pReal) -0.5_pReal) + enddo + enddo + enddo + +end function mesh_build_ipCoordinates + + !-------------------------------------------------------------------------------------------------- !> @brief Store FEid, type, material, texture, and node list per element. !! Allocates global array 'mesh_element' @@ -544,7 +564,6 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes) debug_level, & debug_levelBasic - implicit none real(pReal), intent(in), dimension(:,:,:,:) :: & centres real(pReal), dimension(3,size(centres,2)+1,size(centres,3)+1,size(centres,4)+1) :: & @@ -623,16 +642,6 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes) end function mesh_nodesAroundCentres -!################################################################################################################# -!################################################################################################################# -!################################################################################################################# -! The following routines are not solver specific and should be included in mesh_base (most likely in modified form) -!################################################################################################################# -!################################################################################################################# -!################################################################################################################# - - - !-------------------------------------------------------------------------------------------------- !> @brief Split CP elements into cells. !> @details Build a mapping between cells and the corresponding cell nodes ('mesh_cell'). @@ -762,28 +771,6 @@ end function mesh_build_cellnodes ! HAS TO BE CHANGED IN A LATER VERSION. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipCoordinates - - implicit none - integer(pInt) :: e,c,i,n - real(pReal), dimension(3) :: myCoords - - allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal) - - !$OMP PARALLEL DO PRIVATE(c,myCoords) - do e = 1_pInt,theMesh%nElems ! loop over cpElems - c = theMesh%elem%cellType - do i = 1_pInt,theMesh%elem%nIPs ! loop over ips=cells in this element - myCoords = 0.0_pReal - do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell - myCoords = myCoords + mesh_cellnode(1:3,mesh_cell(n,i,e)) - enddo - mesh_ipCoordinates(1:3,i,e) = myCoords / real(FE_NcellnodesPerCell(c),pReal) - enddo - enddo - !$OMP END PARALLEL DO - -end subroutine mesh_build_ipCoordinates !-------------------------------------------------------------------------------------------------- @@ -806,32 +793,7 @@ subroutine mesh_build_ipAreas c = theMesh%elem%cellType select case (c) - case (1_pInt,2_pInt) ! 2D 3 or 4 node - do i = 1_pInt,theMesh%elem%nIPs ! loop over ips=cells in this element - do f = 1_pInt,FE_NipNeighbors(c) ! loop over cell faces - forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) & - nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),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 - 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 - - case (3_pInt) ! 3D 4node - do i = 1_pInt,theMesh%elem%nIPs ! loop over ips=cells in this element - do f = 1_pInt,FE_NipNeighbors(c) ! loop over cell faces - forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) & - nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e)) - normal = math_crossproduct(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 - enddo - enddo - - case (4_pInt) ! 3D 8node + case (4_pInt) ! 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 From cac472c506b135f5bb2bae8ddd8f792e5ef5ab63 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 14 May 2019 00:13:43 +0200 Subject: [PATCH 06/13] cleaning more complex initialization will be used by FEM solvers only --- src/mesh_grid.f90 | 295 ++++------------------------------------------ 1 file changed, 23 insertions(+), 272 deletions(-) diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index e752f5d41..a7da77bf0 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -35,8 +35,8 @@ module mesh mesh_unitlength !< physical length of one unit in mesh real(pReal), dimension(:,:), allocatable, private :: & - mesh_node, & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) - mesh_cellnode !< cell node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) + mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) + real(pReal), dimension(:,:), allocatable, public, protected :: & mesh_ipVolume, & !< volume associated with IP (initially!) @@ -53,53 +53,14 @@ module mesh logical, dimension(3), public, parameter :: mesh_periodicSurface = .true. !< flag indicating periodic outer surfaces (used for fluxes) -integer(pInt), dimension(:,:), allocatable, private :: & - mesh_cellnodeParent !< cellnode's parent element ID, cellnode's intra-element ID integer(pInt),dimension(:,:,:), allocatable, private :: & mesh_cell !< cell connectivity for each element,ip/cell - integer(pInt), dimension(:,:,:), allocatable, private :: & - FE_cellface !< list of intra-cell cell node IDs that constitute the cell faces of a specific type of cell - - -! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS) -! Hence, I suggest to prefix with "FE_" - integer(pInt), parameter, private :: & FE_Ngeomtypes = 10_pInt, & FE_Ncelltypes = 4_pInt, & - FE_maxNcellnodesPerCell = 8_pInt, & - FE_maxNcellfaces = 6_pInt, & - FE_maxNcellnodesPerCellface = 4_pInt - - - - integer(pInt), dimension(FE_Ncelltypes), parameter, private :: FE_NcellnodesPerCell = & !< number of cell nodes in a specific cell type - int([ & - 3, & ! (2D 3node) - 4, & ! (2D 4node) - 4, & ! (3D 4node) - 8 & ! (3D 8node) - ],pInt) - - integer(pInt), dimension(FE_Ncelltypes), parameter, private :: 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(pInt), dimension(FE_Ncelltypes), parameter, private :: 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) - + FE_maxNcellnodesPerCell = 8_pInt integer(pInt), dimension(3), public, protected :: & grid !< (global) grid @@ -117,13 +78,10 @@ integer(pInt), dimension(:,:), allocatable, private :: & mesh_init private :: & - mesh_build_cellconnectivity, & mesh_build_ipAreas, & - mesh_build_FEdata, & mesh_spectral_build_nodes, & mesh_spectral_build_elements, & mesh_spectral_build_ipNeighborhood, & - mesh_build_cellnodes, & mesh_build_ipCoordinates type, public, extends(tMesh) :: tMesh_grid @@ -234,16 +192,14 @@ subroutine mesh_init(ip,el) if (myDebug) write(6,'(a)') ' Built elements'; flush(6) - call mesh_build_FEdata ! get properties of the different types of elements - call mesh_build_cellconnectivity - if (myDebug) write(6,'(a)') ' Built cell connectivity'; flush(6) - mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes) + if (myDebug) write(6,'(a)') ' Built cell nodes'; flush(6) mesh_ipCoordinates = mesh_build_ipCoordinates() if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6) allocate(mesh_ipVolume(1,theMesh%nElems),source=product([geomSize(1:2),size3]/real([grid(1:2),grid3]))) if (myDebug) write(6,'(a)') ' Built IP volumes'; flush(6) - call mesh_build_ipAreas + mesh_ipArea = mesh_build_ipAreas() + call mesh_build_ipAreas2 if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) call mesh_spectral_build_ipNeighborhood @@ -261,9 +217,6 @@ subroutine mesh_init(ip,el) !!!! COMPATIBILITY HACK !!!! -! for a homogeneous mesh, all elements have the same number of IPs and and cell nodes. -! hence, xxPerElem instead of maxXX -! better name theMesh%homogenizationAt = mesh_element(3,:) theMesh%microstructureAt = mesh_element(4,:) !!!!!!!!!!!!!!!!!!!!!!!! @@ -643,236 +596,34 @@ end function mesh_nodesAroundCentres !-------------------------------------------------------------------------------------------------- -!> @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. -!> Also allocates the 'mesh_node' array. +!> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal' !-------------------------------------------------------------------------------------------------- -subroutine mesh_build_cellconnectivity +pure function mesh_build_ipAreas() - implicit none - integer(pInt), dimension(:), allocatable :: & - matchingNode2cellnode - integer(pInt), dimension(:,:), allocatable :: & - cellnodeParent - integer(pInt), dimension(theMesh%elem%Ncellnodes) :: & - localCellnode2globalCellnode - integer(pInt) :: & - e,n,i, & - matchingNodeID, & - localCellnodeID - - integer(pInt), 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) + real(pReal), dimension(6,1,theMesh%nElems) :: mesh_build_ipAreas - allocate(mesh_cell(FE_maxNcellnodesPerCell,theMesh%elem%nIPs,theMesh%nElems), source=0_pInt) - allocate(matchingNode2cellnode(theMesh%nNodes), source=0_pInt) - allocate(cellnodeParent(2_pInt,theMesh%elem%Ncellnodes*theMesh%nElems), source=0_pInt) - - mesh_Ncells = theMesh%nElems*theMesh%elem%nIPs -!-------------------------------------------------------------------------------------------------- -! Count cell nodes (including duplicates) and generate cell connectivity list - mesh_Ncellnodes = 0_pInt - - do e = 1_pInt,theMesh%nElems - localCellnode2globalCellnode = 0_pInt - do i = 1_pInt,theMesh%elem%nIPs - do n = 1_pInt,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_pInt+localCellnodeID,e) - if (matchingNode2cellnode(matchingNodeID) == 0_pInt) then ! if this matching node does not yet exist in the glbal cell node list ... - mesh_Ncellnodes = mesh_Ncellnodes + 1_pInt ! ... count it as cell node ... - matchingNode2cellnode(matchingNodeID) = mesh_Ncellnodes ! ... and remember its global ID - cellnodeParent(1_pInt,mesh_Ncellnodes) = e ! ... and where it belongs to - cellnodeParent(2_pInt,mesh_Ncellnodes) = localCellnodeID - endif - mesh_cell(n,i,e) = matchingNode2cellnode(matchingNodeID) - else ! this cell node is no matching node - if (localCellnode2globalCellnode(localCellnodeID) == 0_pInt) then ! if this local cell node does not yet exist in the global cell node list ... - mesh_Ncellnodes = mesh_Ncellnodes + 1_pInt ! ... count it as cell node ... - localCellnode2globalCellnode(localCellnodeID) = mesh_Ncellnodes ! ... and remember its global ID ... - cellnodeParent(1_pInt,mesh_Ncellnodes) = e ! ... and it belongs to - cellnodeParent(2_pInt,mesh_Ncellnodes) = localCellnodeID - endif - mesh_cell(n,i,e) = localCellnode2globalCellnode(localCellnodeID) - endif - enddo - enddo - enddo - - allocate(mesh_cellnodeParent(2_pInt,mesh_Ncellnodes)) - allocate(mesh_cellnode(3_pInt,mesh_Ncellnodes)) - - forall(n = 1_pInt: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(nodes,Ncellnodes) - - implicit none - integer(pInt), intent(in) :: Ncellnodes !< requested number of cellnodes - real(pReal), dimension(3,mesh_Nnodes), intent(in) :: nodes - real(pReal), dimension(3,Ncellnodes) :: mesh_build_cellnodes - - integer(pInt) :: & - e,n,m, & - localCellnodeID - real(pReal), dimension(3) :: & - myCoords - - mesh_build_cellnodes = 0.0_pReal -!$OMP PARALLEL DO PRIVATE(e,localCellnodeID,myCoords) - do n = 1_pInt,Ncellnodes ! loop over cell nodes - e = mesh_cellnodeParent(1,n) - localCellnodeID = mesh_cellnodeParent(2,n) - myCoords = 0.0_pReal - do m = 1_pInt,theMesh%elem%nNodes - myCoords = myCoords + nodes(1:3,mesh_element(4_pInt+m,e)) & - * theMesh%elem%cellNodeParentNodeWeights(m,localCellnodeID) - enddo - mesh_build_cellnodes(1:3,n) = myCoords / sum(theMesh%elem%cellNodeParentNodeWeights(:,localCellnodeID)) - enddo -!$OMP END PARALLEL DO - -end function mesh_build_cellnodes - - -!-------------------------------------------------------------------------------------------------- -!> @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, -! 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. -! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!-------------------------------------------------------------------------------------------------- + mesh_build_ipAreas(1:2,1,:) = geomSize(2)/real(grid(2)) * geomSize(3)/real(grid(3)) + mesh_build_ipAreas(3:4,1,:) = geomSize(3)/real(grid(3)) * geomSize(1)/real(grid(1)) + mesh_build_ipAreas(5:6,1,:) = geomSize(1)/real(grid(1)) * geomSize(2)/real(grid(2)) + +end function mesh_build_ipAreas !-------------------------------------------------------------------------------------------------- !> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal' !-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipAreas - use math, only: & - math_crossproduct +subroutine mesh_build_ipAreas2 - implicit none - integer(pInt) :: e,t,g,c,i,f,n,m - real(pReal), dimension (3,FE_maxNcellnodesPerCellface) :: nodePos, normals - real(pReal), dimension(3) :: normal + allocate(mesh_ipAreaNormal(3,6,1,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_pInt,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) - - !$OMP PARALLEL DO PRIVATE(t,g,c,nodePos,normal,normals) - do e = 1_pInt,theMesh%nElems ! loop over cpElems - c = theMesh%elem%cellType - select case (c) - - case (4_pInt) - ! 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 = FE_NcellnodesPerCellface(c) - do i = 1_pInt,theMesh%elem%nIPs ! loop over ips=cells in this element - do f = 1_pInt,FE_NipNeighbors(c) ! loop over cell faces - forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) & - nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e)) - forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) & - normals(1:3,n) = 0.5_pReal & - * math_crossproduct(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) - enddo - enddo - - end select - enddo - !$OMP END PARALLEL DO - -end subroutine mesh_build_ipAreas - - -!-------------------------------------------------------------------------------------------------- -!> @brief get properties of different types of finite elements -!> @details assign globals: FE_nodesAtIP, FE_ipNeighbor, FE_subNodeOnIPFace -!-------------------------------------------------------------------------------------------------- -subroutine mesh_build_FEdata - - implicit none - integer(pInt) :: me - allocate(FE_cellface(FE_maxNcellnodesPerCellface,FE_maxNcellfaces,FE_Ncelltypes), source=0_pInt) - - - ! *** FE_cellface *** - me = 0_pInt - - me = me + 1_pInt - FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 2D 3node, VTK_TRIANGLE (5) - reshape(int([& - 2,3, & - 3,1, & - 1,2 & - ],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)]) - - me = me + 1_pInt - FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 2D 4node, VTK_QUAD (9) - reshape(int([& - 2,3, & - 4,1, & - 3,4, & - 1,2 & - ],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)]) - - me = me + 1_pInt - FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 3D 4node, VTK_TETRA (10) - reshape(int([& - 1,3,2, & - 1,2,4, & - 2,3,4, & - 1,4,3 & - ],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)]) - - me = me + 1_pInt - FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 3D 8node, VTK_HEXAHEDRON (12) - reshape(int([& - 2,3,7,6, & - 4,1,5,8, & - 3,4,8,7, & - 1,2,6,5, & - 5,6,7,8, & - 1,4,3,2 & - ],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)]) + mesh_ipAreaNormal(1:3,1,1,:) = spread([+1.0_pReal, 0.0_pReal, 0.0_pReal],2,theMesh%nElems) + mesh_ipAreaNormal(1:3,2,1,:) = spread([-1.0_pReal, 0.0_pReal, 0.0_pReal],2,theMesh%nElems) + mesh_ipAreaNormal(1:3,3,1,:) = spread([ 0.0_pReal,+1.0_pReal, 0.0_pReal],2,theMesh%nElems) + mesh_ipAreaNormal(1:3,4,1,:) = spread([ 0.0_pReal,-1.0_pReal, 0.0_pReal],2,theMesh%nElems) + mesh_ipAreaNormal(1:3,5,1,:) = spread([ 0.0_pReal, 0.0_pReal,+1.0_pReal],2,theMesh%nElems) + mesh_ipAreaNormal(1:3,6,1,:) = spread([ 0.0_pReal, 0.0_pReal,-1.0_pReal],2,theMesh%nElems) - -end subroutine mesh_build_FEdata +end subroutine mesh_build_ipAreas2 end module mesh From c9f1f8756d48a632df05cb9bac2c35b500edd171 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 14 May 2019 06:06:01 +0200 Subject: [PATCH 07/13] better name/not needed --- src/mesh_grid.f90 | 36 +++++++++++++++--------------------- 1 file changed, 15 insertions(+), 21 deletions(-) diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index a7da77bf0..320024dc3 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -35,7 +35,7 @@ module mesh mesh_unitlength !< physical length of one unit in mesh real(pReal), dimension(:,:), allocatable, private :: & - mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) + mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) real(pReal), dimension(:,:), allocatable, public, protected :: & @@ -54,14 +54,7 @@ module mesh logical, dimension(3), public, parameter :: mesh_periodicSurface = .true. !< flag indicating periodic outer surfaces (used for fluxes) - integer(pInt),dimension(:,:,:), allocatable, private :: & - mesh_cell !< cell connectivity for each element,ip/cell - - integer(pInt), parameter, private :: & - FE_Ngeomtypes = 10_pInt, & - FE_Ncelltypes = 4_pInt, & - FE_maxNcellnodesPerCell = 8_pInt - +! grid specific integer(pInt), dimension(3), public, protected :: & grid !< (global) grid integer(pInt), public, protected :: & @@ -79,6 +72,7 @@ module mesh private :: & mesh_build_ipAreas, & + mesh_build_ipNormals, & mesh_spectral_build_nodes, & mesh_spectral_build_elements, & mesh_spectral_build_ipNeighborhood, & @@ -198,8 +192,8 @@ subroutine mesh_init(ip,el) if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6) allocate(mesh_ipVolume(1,theMesh%nElems),source=product([geomSize(1:2),size3]/real([grid(1:2),grid3]))) if (myDebug) write(6,'(a)') ' Built IP volumes'; flush(6) - mesh_ipArea = mesh_build_ipAreas() - call mesh_build_ipAreas2 + mesh_ipArea = mesh_build_ipAreas() + mesh_ipAreaNormal = mesh_build_ipNormals() if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) call mesh_spectral_build_ipNeighborhood @@ -220,7 +214,7 @@ subroutine mesh_init(ip,el) theMesh%homogenizationAt = mesh_element(3,:) theMesh%microstructureAt = mesh_element(4,:) !!!!!!!!!!!!!!!!!!!!!!!! - deallocate(mesh_cell) + end subroutine mesh_init @@ -612,18 +606,18 @@ end function mesh_build_ipAreas !-------------------------------------------------------------------------------------------------- !> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal' !-------------------------------------------------------------------------------------------------- -subroutine mesh_build_ipAreas2 +pure function mesh_build_ipNormals() - allocate(mesh_ipAreaNormal(3,6,1,theMesh%nElems), source=0.0_pReal) + real, dimension(3,6,1,theMesh%nElems) :: mesh_build_ipNormals - mesh_ipAreaNormal(1:3,1,1,:) = spread([+1.0_pReal, 0.0_pReal, 0.0_pReal],2,theMesh%nElems) - mesh_ipAreaNormal(1:3,2,1,:) = spread([-1.0_pReal, 0.0_pReal, 0.0_pReal],2,theMesh%nElems) - mesh_ipAreaNormal(1:3,3,1,:) = spread([ 0.0_pReal,+1.0_pReal, 0.0_pReal],2,theMesh%nElems) - mesh_ipAreaNormal(1:3,4,1,:) = spread([ 0.0_pReal,-1.0_pReal, 0.0_pReal],2,theMesh%nElems) - mesh_ipAreaNormal(1:3,5,1,:) = spread([ 0.0_pReal, 0.0_pReal,+1.0_pReal],2,theMesh%nElems) - mesh_ipAreaNormal(1:3,6,1,:) = spread([ 0.0_pReal, 0.0_pReal,-1.0_pReal],2,theMesh%nElems) + mesh_build_ipNormals(1:3,1,1,:) = spread([+1.0_pReal, 0.0_pReal, 0.0_pReal],2,theMesh%nElems) + mesh_build_ipNormals(1:3,2,1,:) = spread([-1.0_pReal, 0.0_pReal, 0.0_pReal],2,theMesh%nElems) + mesh_build_ipNormals(1:3,3,1,:) = spread([ 0.0_pReal,+1.0_pReal, 0.0_pReal],2,theMesh%nElems) + mesh_build_ipNormals(1:3,4,1,:) = spread([ 0.0_pReal,-1.0_pReal, 0.0_pReal],2,theMesh%nElems) + mesh_build_ipNormals(1:3,5,1,:) = spread([ 0.0_pReal, 0.0_pReal,+1.0_pReal],2,theMesh%nElems) + mesh_build_ipNormals(1:3,6,1,:) = spread([ 0.0_pReal, 0.0_pReal,-1.0_pReal],2,theMesh%nElems) -end subroutine mesh_build_ipAreas2 +end function mesh_build_ipNormals end module mesh From 8dea95879c22fadfde569e8b6704246f8ddcd45a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 14 May 2019 07:22:29 +0200 Subject: [PATCH 08/13] specific for nonlocal, can be calculated during post processing --- src/crystallite.f90 | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 02aec8b7a..41d1cb1c3 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -87,7 +87,6 @@ module crystallite enumerator :: undefined_ID, & phase_ID, & texture_ID, & - volume_ID, & orientation_ID, & grainrotation_ID, & defgrad_ID, & @@ -273,8 +272,6 @@ subroutine crystallite_init crystallite_outputID(o,c) = phase_ID case ('texture') outputName crystallite_outputID(o,c) = texture_ID - case ('volume') outputName - crystallite_outputID(o,c) = volume_ID case ('orientation') outputName crystallite_outputID(o,c) = orientation_ID case ('grainrotation') outputName @@ -323,7 +320,7 @@ subroutine crystallite_init do r = 1,size(config_crystallite) do o = 1,crystallite_Noutput(r) select case(crystallite_outputID(o,r)) - case(phase_ID,texture_ID,volume_ID) + case(phase_ID,texture_ID) mySize = 1 case(orientation_ID,grainrotation_ID) mySize = 4 @@ -948,7 +945,6 @@ function crystallite_postResults(ipc, ip, el) use mesh, only: & theMesh, & mesh_element, & - mesh_ipVolume, & mesh_ipNeighborhood use material, only: & plasticState, & @@ -999,11 +995,6 @@ function crystallite_postResults(ipc, ip, el) case (texture_ID) mySize = 1 crystallite_postResults(c+1) = real(material_texture(ipc,ip,el),pReal) ! textureID of grain - case (volume_ID) - mySize = 1 - detF = math_det33(crystallite_partionedF(1:3,1:3,ipc,ip,el)) ! V_current = det(F) * V_reference - crystallite_postResults(c+1) = detF * mesh_ipVolume(ip,el) & - / real(homogenization_Ngrains(mesh_element(3,el)),pReal) ! grain volume (not fraction but absolute) case (orientation_ID) mySize = 4 crystallite_postResults(c+1:c+mySize) = crystallite_orientation(ipc,ip,el)%asQuaternion() From a8b9b5d1c901ce6f5750a9c7807054d73794e3d8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 14 May 2019 08:11:23 +0200 Subject: [PATCH 09/13] homogenizationAt array had wrong shape for MPI (too large) --- src/crystallite.f90 | 2 +- src/mesh_grid.f90 | 16 +++++----------- 2 files changed, 6 insertions(+), 12 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 41d1cb1c3..e972ff162 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -299,7 +299,7 @@ subroutine crystallite_init case ('neighboringelement') outputName ! ToDo: this is not a result, it is static. Should be written out by mesh crystallite_outputID(o,c) = neighboringelement_ID case default outputName - call IO_error(105,ext_msg=trim(str(o))//' (Crystallite)') + !call IO_error(105,ext_msg=trim(str(o))//' (Crystallite)') end select outputName enddo enddo diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index 320024dc3..685a53ba8 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -13,12 +13,7 @@ module mesh implicit none private integer(pInt), public, protected :: & - mesh_Nnodes, & !< total number of nodes in mesh - mesh_Ncellnodes, & !< total number of cell nodes in mesh (including duplicates) - mesh_Ncells, & !< total number of cells in mesh - mesh_maxNipNeighbors, & !< max number of IP neighbors in any CP element - mesh_maxNsharedElems !< max number of CP elements sharing a node - + mesh_Nnodes integer(pInt), dimension(:), allocatable, private :: & microGlobal @@ -138,9 +133,8 @@ subroutine mesh_init(ip,el) implicit none include 'fftw3-mpi.f03' integer(C_INTPTR_T) :: devNull, local_K, local_K_offset - integer :: ierr, worldsize, i + integer :: ierr, worldsize, j integer(pInt), intent(in), optional :: el, ip - integer(pInt) :: j logical :: myDebug write(6,'(/,a)') ' <<<+- mesh init -+>>>' @@ -179,10 +173,10 @@ subroutine mesh_init(ip,el) call theMesh%init(mesh_node) call theMesh%setNelems(product(grid(1:2))*grid3) - mesh_homogenizationAt = mesh_homogenizationAt(product(grid(1:2))*grid3) ! reallocate/shrink in case of MPI - mesh_maxNipNeighbors = theMesh%elem%nIPneighbors - call mesh_spectral_build_elements() + mesh_homogenizationAt = mesh_homogenizationAt(product(grid(1:2))*grid3Offset+1: & + product(grid(1:2))*(grid3Offset+grid3)) ! reallocate/shrink in case of MPI + if (myDebug) write(6,'(a)') ' Built elements'; flush(6) From 4bce087d3d44350aa591b9bfefd103d2e3aa446d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 14 May 2019 11:32:25 +0200 Subject: [PATCH 10/13] separating functionality (stub only) - mesh (initialization depends on type of solver) only provides information about number of elements, IPs and writes out displacements - new module (will get setter functions for solver specific initialization) will provide information that is only used by the nonlocal model --- src/commercialFEM_fileList.f90 | 1 + src/geometry_plastic_nonlocal.f90 | 28 ++++++++++++++++++++++++++++ src/plastic_nonlocal.f90 | 11 ++++++++--- 3 files changed, 37 insertions(+), 3 deletions(-) create mode 100644 src/geometry_plastic_nonlocal.f90 diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 0ae1323f0..ef3870ece 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -28,6 +28,7 @@ #endif #include "material.f90" #include "lattice.f90" +#include "geometry_plastic_nonlocal.f90" #include "source_thermal_dissipation.f90" #include "source_thermal_externalheat.f90" #include "source_damage_isoBrittle.f90" diff --git a/src/geometry_plastic_nonlocal.f90 b/src/geometry_plastic_nonlocal.f90 new file mode 100644 index 000000000..579004b2f --- /dev/null +++ b/src/geometry_plastic_nonlocal.f90 @@ -0,0 +1,28 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!> @author Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @brief Geometric information about the IP cells needed for the nonlocal +! plasticity model +!-------------------------------------------------------------------------------------------------- +module geometry_plastic_nonlocal + use prec + + implicit none + private + logical, dimension(3), public, parameter :: & + geometry_plastic_nonlocal_periodicSurface = .true. !< flag indicating periodic outer surfaces (used for fluxes) NEEDED? + + integer, dimension(:,:,:,:), allocatable, public, protected :: & + geometry_plastic_nonlocal_IPneighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] + + real(pReal), dimension(:,:), allocatable, public, protected :: & + geometry_plastic_nonlocal_IPvolume !< volume associated with IP (initially!) + + real(pReal), dimension(:,:,:), allocatable, public, protected :: & + geometry_plastic_nonlocal_IParea !< area of interface to neighboring IP (initially!) + + real(pReal),dimension(:,:,:,:), allocatable, public, protected :: & + geometry_plastic_nonlocal_IPareaNormal !< area normal of interface to neighboring IP (initially!) + +end module geometry_plastic_nonlocal diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90 index 6097bbbc8..66e8f8980 100644 --- a/src/plastic_nonlocal.f90 +++ b/src/plastic_nonlocal.f90 @@ -5,10 +5,15 @@ !> @brief material subroutine for plasticity including dislocation flux !-------------------------------------------------------------------------------------------------- module plastic_nonlocal - use prec, only: & - pReal + use prec use future - + use geometry_plastic_nonlocal, only: & + periodicSurface => geometry_plastic_nonlocal_periodicSurface, & + IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, & + IPvolume => geometry_plastic_nonlocal_IPvolume, & + IParea => geometry_plastic_nonlocal_IParea, & + IPareaNormal => geometry_plastic_nonlocal_IPareaNormal + implicit none private real(pReal), parameter, private :: & From 0dbc6fb435f0dc503410672894c08deb89de516e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 14 May 2019 11:52:28 +0200 Subject: [PATCH 11/13] usage example --- src/commercialFEM_fileList.f90 | 2 +- src/geometry_plastic_nonlocal.f90 | 24 ++++++++++++++++++++++++ src/mesh_grid.f90 | 4 +++- 3 files changed, 28 insertions(+), 2 deletions(-) diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index ef3870ece..1c7287813 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -14,6 +14,7 @@ #include "Lambert.f90" #include "rotations.f90" #include "FEsolving.f90" +#include "geometry_plastic_nonlocal.f90" #include "element.f90" #include "mesh_base.f90" #ifdef Abaqus @@ -28,7 +29,6 @@ #endif #include "material.f90" #include "lattice.f90" -#include "geometry_plastic_nonlocal.f90" #include "source_thermal_dissipation.f90" #include "source_thermal_externalheat.f90" #include "source_damage_isoBrittle.f90" diff --git a/src/geometry_plastic_nonlocal.f90 b/src/geometry_plastic_nonlocal.f90 index 579004b2f..0b63b7f9c 100644 --- a/src/geometry_plastic_nonlocal.f90 +++ b/src/geometry_plastic_nonlocal.f90 @@ -24,5 +24,29 @@ module geometry_plastic_nonlocal real(pReal),dimension(:,:,:,:), allocatable, public, protected :: & geometry_plastic_nonlocal_IPareaNormal !< area normal of interface to neighboring IP (initially!) + + public :: & + geometry_plastic_nonlocal_set_IPneighborhood, & + geometry_plastic_nonlocal_set_IPvolume + + contains + +subroutine geometry_plastic_nonlocal_set_IPneighborhood(IPneighborhood) + + integer, dimension(:,:,:,:), intent(in) :: IPneighborhood + + geometry_plastic_nonlocal_IPneighborhood = IPneighborhood + +end subroutine geometry_plastic_nonlocal_set_IPneighborhood + + +subroutine geometry_plastic_nonlocal_set_IPvolume(IPvolume) + + real(pReal), dimension(:,:), intent(in) :: IPvolume + + geometry_plastic_nonlocal_IPvolume = IPvolume + +end subroutine geometry_plastic_nonlocal_set_IPvolume + end module geometry_plastic_nonlocal diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index 685a53ba8..d873e3542 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -7,7 +7,8 @@ !-------------------------------------------------------------------------------------------------- module mesh use, intrinsic :: iso_c_binding - use prec, only: pReal, pInt + use prec + use geometry_plastic_nonlocal use mesh_base implicit none @@ -191,6 +192,7 @@ subroutine mesh_init(ip,el) if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) call mesh_spectral_build_ipNeighborhood + call geometry_plastic_nonlocal_set_IPneighborhood(mesh_ipNeighborhood) if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6) From d4392dc1bdb875cdd50b77af8bcc7d000841af36 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 30 May 2019 23:58:59 +0200 Subject: [PATCH 12/13] one implicit none is enough --- src/grid/spectral_utilities.f90 | 22 ---------------------- 1 file changed, 22 deletions(-) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 1d5e42070..f545eab4e 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -196,7 +196,6 @@ subroutine utilities_init grid3Offset, & geomSize - implicit none PetscErrorCode :: ierr integer :: i, j, k, & FFTW_planner_flag @@ -425,7 +424,6 @@ subroutine utilities_updateGamma(C,saveReference) math_det33, & math_invert2 - implicit none real(pReal), intent(in), dimension(3,3,3,3) :: C !< input stiffness to store as reference stiffness logical , intent(in) :: saveReference !< save reference stiffness to file for restart complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx @@ -473,7 +471,6 @@ end subroutine utilities_updateGamma !> @details Does an unweighted filtered FFT transform from real to complex !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTtensorForward - implicit none call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) @@ -485,7 +482,6 @@ end subroutine utilities_FFTtensorForward !> @details Does an weighted inverse FFT transform from complex to real !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTtensorBackward - implicit none call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real) tensorField_real = tensorField_real * wgt ! normalize the result by number of elements @@ -497,7 +493,6 @@ end subroutine utilities_FFTtensorBackward !> @details Does an unweighted filtered FFT transform from real to complex !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTscalarForward - implicit none call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier) @@ -509,7 +504,6 @@ end subroutine utilities_FFTscalarForward !> @details Does an weighted inverse FFT transform from complex to real !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTscalarBackward - implicit none call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real) scalarField_real = scalarField_real * wgt ! normalize the result by number of elements @@ -522,7 +516,6 @@ end subroutine utilities_FFTscalarBackward !> @details Does an unweighted filtered FFT transform from real to complex. !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTvectorForward - implicit none call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier) @@ -534,7 +527,6 @@ end subroutine utilities_FFTvectorForward !> @details Does an weighted inverse FFT transform from complex to real !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTvectorBackward - implicit none call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real) vectorField_real = vectorField_real * wgt ! normalize the result by number of elements @@ -554,7 +546,6 @@ subroutine utilities_fourierGammaConvolution(fieldAim) grid, & grid3Offset - implicit none real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx real(pReal), dimension(6,6) :: A, A_inv @@ -615,7 +606,6 @@ subroutine utilities_fourierGreenConvolution(D_ref, mobility_ref, deltaT) grid, & grid3 - implicit none real(pReal), dimension(3,3), intent(in) :: D_ref real(pReal), intent(in) :: mobility_ref, deltaT complex(pReal) :: GreenOp_hat @@ -644,7 +634,6 @@ real(pReal) function utilities_divergenceRMS() grid, & grid3 - implicit none integer :: i, j, k, ierr complex(pReal), dimension(3) :: rescaledGeom @@ -694,7 +683,6 @@ real(pReal) function utilities_curlRMS() grid, & grid3 - implicit none integer :: i, j, k, l, ierr complex(pReal), dimension(3,3) :: curl_fourier complex(pReal), dimension(3) :: rescaledGeom @@ -766,7 +754,6 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) math_rotate_forward33, & math_invert2 - implicit none real(pReal), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance real(pReal), intent(in) , dimension(3,3,3,3) :: C !< current average stiffness real(pReal), intent(in) , dimension(3,3) :: rot_BC !< rotation of load frame @@ -861,7 +848,6 @@ subroutine utilities_fourierScalarGradient() grid3, & grid - implicit none integer :: i, j, k vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) @@ -879,7 +865,6 @@ subroutine utilities_fourierVectorDivergence() grid3, & grid - implicit none integer :: i, j, k scalarField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) @@ -898,7 +883,6 @@ subroutine utilities_fourierVectorGradient() grid3, & grid - implicit none integer :: i, j, k, m, n tensorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) @@ -919,7 +903,6 @@ subroutine utilities_fourierTensorDivergence() grid3, & grid - implicit none integer :: i, j, k, m, n vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) @@ -954,7 +937,6 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& materialpoint_dPdF, & materialpoint_stressAndItsTangent - implicit none real(pReal),intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress @@ -1032,7 +1014,6 @@ pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate) grid3, & grid - implicit none real(pReal), intent(in), dimension(3,3) :: & avRate !< homogeneous addon real(pReal), intent(in) :: & @@ -1063,7 +1044,6 @@ function utilities_forwardField(timeinc,field_lastInc,rate,aim) grid3, & grid - implicit none real(pReal), intent(in) :: & timeinc !< timeinc of current step real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & @@ -1100,7 +1080,6 @@ pure function utilities_getFreqDerivative(k_s) geomSize, & grid - implicit none integer, intent(in), dimension(3) :: k_s !< indices of frequency complex(pReal), dimension(3) :: utilities_getFreqDerivative @@ -1158,7 +1137,6 @@ subroutine utilities_updateIPcoords(F) grid3Offset, & geomSize, & mesh_ipCoordinates - implicit none real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F integer :: i, j, k, m, ierr From 30826c5c8666d5c50d1e42d1304a51c64b1ef258 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 31 May 2019 09:05:58 +0200 Subject: [PATCH 13/13] not supported anymore is part of the geometry --- examples/SpectralMethod/Polycrystal/material.config | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/SpectralMethod/Polycrystal/material.config b/examples/SpectralMethod/Polycrystal/material.config index 71d7e07d7..8103e7128 100644 --- a/examples/SpectralMethod/Polycrystal/material.config +++ b/examples/SpectralMethod/Polycrystal/material.config @@ -11,7 +11,6 @@ mech none [almostAll] (output) phase (output) texture -(output) volume (output) orientation # quaternion (output) grainrotation # deviation from initial orientation as axis (1-3) and angle in degree (4) (output) f # deformation gradient tensor; synonyms: "defgrad"