From 18b8e71f69f19091a55835a773095d39ac618b46 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 5 Oct 2019 20:09:01 +0200 Subject: [PATCH] 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