diff --git a/src/mesh_abaqus.f90 b/src/mesh_abaqus.f90 index 98bdda4ef..62ece4c93 100644 --- a/src/mesh_abaqus.f90 +++ b/src/mesh_abaqus.f90 @@ -389,7 +389,7 @@ integer(pInt), dimension(:,:), allocatable, private :: & mesh_abaqus_build_elements - type, public, extends(tMesh) :: tMesh_Abaqus + type, public, extends(tMesh) :: tMesh_abaqus integer(pInt):: & mesh_Nelems, & !< total number of elements in mesh (including non-DAMASK elements) @@ -406,16 +406,22 @@ integer(pInt), dimension(:,:), allocatable, private :: & logical:: noPart !< for cases where the ABAQUS input file does not use part/assembly information contains - procedure :: init=>tMesh_abaqus_init - end type tMesh_Abaqus + procedure, pass(self) :: tMesh_abaqus_init + generic, public :: init => tMesh_abaqus_init + end type tMesh_abaqus - type(tMesh_Abaqus), public, protected :: theMesh + type(tMesh_abaqus), public, protected :: theMesh contains -subroutine tMesh_abaqus_init(self) +subroutine tMesh_abaqus_init(self,elemType,nodes) + implicit none class(tMesh_abaqus) :: self + real(pReal), dimension(:,:), intent(in) :: nodes + integer(pInt), intent(in) :: elemType + + call self%tMesh%init('mesh',elemType,nodes) end subroutine tMesh_abaqus_init @@ -537,7 +543,7 @@ subroutine mesh_init(ip,el) mesh_microstructureAt = mesh_element(4,:) mesh_CPnodeID = mesh_element(5:4+mesh_NipsPerElem,:) !!!!!!!!!!!!!!!!!!!!!!!! - + call theMesh%init(mesh_element(2,1),mesh_node0) contains !-------------------------------------------------------------------------------------------------- !> @brief check if the input file for Abaqus contains part info @@ -859,97 +865,6 @@ pure function mesh_cellCenterCoordinates(ip,el) end function mesh_cellCenterCoordinates - -!-------------------------------------------------------------------------------------------------- -!> @brief builds mesh of (distorted) cubes for given coordinates (= center of the cubes) -!-------------------------------------------------------------------------------------------------- -function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes) - use debug, only: & - debug_mesh, & - debug_level, & - debug_levelBasic - use math, only: & - math_mul33x3 - - implicit none - real(pReal), intent(in), dimension(:,:,:,:) :: & - centres - real(pReal), dimension(3,size(centres,2)+1,size(centres,3)+1,size(centres,4)+1) :: & - nodes - real(pReal), intent(in), dimension(3) :: & - gDim - real(pReal), intent(in), dimension(3,3) :: & - Favg - real(pReal), dimension(3,size(centres,2)+2,size(centres,3)+2,size(centres,4)+2) :: & - wrappedCentres - - integer(pInt) :: & - i,j,k,n - integer(pInt), dimension(3), parameter :: & - diag = 1_pInt - integer(pInt), dimension(3) :: & - shift = 0_pInt, & - lookup = 0_pInt, & - me = 0_pInt, & - iRes = 0_pInt - integer(pInt), dimension(3,8) :: & - neighbor = reshape([ & - 0_pInt, 0_pInt, 0_pInt, & - 1_pInt, 0_pInt, 0_pInt, & - 1_pInt, 1_pInt, 0_pInt, & - 0_pInt, 1_pInt, 0_pInt, & - 0_pInt, 0_pInt, 1_pInt, & - 1_pInt, 0_pInt, 1_pInt, & - 1_pInt, 1_pInt, 1_pInt, & - 0_pInt, 1_pInt, 1_pInt ], [3,8]) - -!-------------------------------------------------------------------------------------------------- -! initializing variables - iRes = [size(centres,2),size(centres,3),size(centres,4)] - nodes = 0.0_pReal - wrappedCentres = 0.0_pReal - -!-------------------------------------------------------------------------------------------------- -! report - if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then - write(6,'(a)') ' Meshing cubes around centroids' - write(6,'(a,3(e12.5))') ' Dimension: ', gDim - write(6,'(a,3(i5))') ' Resolution:', iRes - endif - -!-------------------------------------------------------------------------------------------------- -! building wrappedCentres = centroids + ghosts - wrappedCentres(1:3,2_pInt:iRes(1)+1_pInt,2_pInt:iRes(2)+1_pInt,2_pInt:iRes(3)+1_pInt) = centres - do k = 0_pInt,iRes(3)+1_pInt - do j = 0_pInt,iRes(2)+1_pInt - do i = 0_pInt,iRes(1)+1_pInt - if (k==0_pInt .or. k==iRes(3)+1_pInt .or. & ! z skin - j==0_pInt .or. j==iRes(2)+1_pInt .or. & ! y skin - i==0_pInt .or. i==iRes(1)+1_pInt ) then ! x skin - me = [i,j,k] ! me on skin - shift = sign(abs(iRes+diag-2_pInt*me)/(iRes+diag),iRes+diag-2_pInt*me) - lookup = me-diag+shift*iRes - wrappedCentres(1:3,i+1_pInt, j+1_pInt, k+1_pInt) = & - centres(1:3,lookup(1)+1_pInt,lookup(2)+1_pInt,lookup(3)+1_pInt) & - - math_mul33x3(Favg, real(shift,pReal)*gDim) - endif - enddo; enddo; enddo - -!-------------------------------------------------------------------------------------------------- -! averaging - do k = 0_pInt,iRes(3); do j = 0_pInt,iRes(2); do i = 0_pInt,iRes(1) - do n = 1_pInt,8_pInt - nodes(1:3,i+1_pInt,j+1_pInt,k+1_pInt) = & - nodes(1:3,i+1_pInt,j+1_pInt,k+1_pInt) + wrappedCentres(1:3,i+1_pInt+neighbor(1,n), & - j+1_pInt+neighbor(2,n), & - k+1_pInt+neighbor(3,n) ) - enddo - enddo; enddo; enddo - nodes = nodes/8.0_pReal - -end function mesh_nodesAroundCentres - - !-------------------------------------------------------------------------------------------------- !> @brief Count overall number of nodes and elements in mesh and stores them in !! 'mesh_Nelems' and 'mesh_Nnodes' diff --git a/src/mesh_base.f90 b/src/mesh_base.f90 index 477fc3aed..f9a076f03 100644 --- a/src/mesh_base.f90 +++ b/src/mesh_base.f90 @@ -29,7 +29,7 @@ module mesh_base node !< node x,y,z coordinates (deformed) integer(pInt), dimension(:,:), allocatable, public :: & cellnodeParent !< cellnode's parent element ID, cellnode's intra-element ID - character(pStringLen) :: solver = "undefined" + character(pStringLen) :: type = "n/a" integer(pInt) :: & Nnodes, & !< total number of nodes in mesh Nelems = -1_pInt, & @@ -43,6 +43,24 @@ module mesh_base microstructureAt integer(pInt), dimension(:,:), allocatable, public :: & connectivity + contains + procedure, pass(self) :: tMesh_base_init + generic, public :: init => tMesh_base_init end type tMesh +contains +subroutine tMesh_base_init(self,meshType,elemType,nodes) + + implicit none + class(tMesh) :: self + character(len=*), intent(in) :: meshType + integer(pInt), intent(in) :: elemType + real(pReal), dimension(:,:), intent(in) :: nodes + + self%type = meshType + call self%elem%init(elemType) + self%node0 = nodes + +end subroutine tMesh_base_init + end module mesh_base diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index a2a041955..cff0dbc21 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -270,14 +270,11 @@ integer(pInt), dimension(:,:), allocatable, private :: & mesh_init, & mesh_cellCenterCoordinates - private :: & mesh_build_cellconnectivity, & mesh_build_ipAreas, & mesh_build_FEdata, & mesh_spectral_getHomogenization, & - mesh_spectral_count, & - mesh_spectral_count_cpSizes, & mesh_spectral_build_nodes, & mesh_spectral_build_elements, & mesh_spectral_build_ipNeighborhood, & @@ -302,19 +299,21 @@ integer(pInt), dimension(:,:), allocatable, private :: & size3offset contains - procedure :: init => tMesh_grid_init + procedure, pass(self) :: tMesh_grid_init + generic, public :: init => tMesh_grid_init end type tMesh_grid type(tMesh_grid), public, protected :: theMesh contains -subroutine tMesh_grid_init(self) +subroutine tMesh_grid_init(self,nodes) implicit none class(tMesh_grid) :: self + real(pReal), dimension(:,:), intent(in) :: nodes - call self%elem%init(10_pInt) + call self%tMesh%init('grid',10_pInt,nodes) end subroutine tMesh_grid_init @@ -364,7 +363,8 @@ subroutine mesh_init(ip,el) write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - call theMesh%init + + call mesh_build_FEdata ! get properties of the different types of elements mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh @@ -389,13 +389,23 @@ subroutine mesh_init(ip,el) grid3Offset = int(local_K_offset,pInt) size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal) size3Offset = geomSize(3)*real(grid3Offset,pReal)/real(grid(3),pReal) - if (myDebug) write(6,'(a)') ' Grid partitioned'; flush(6) - call mesh_spectral_count() - if (myDebug) write(6,'(a)') ' Counted nodes/elements'; flush(6) - call mesh_spectral_count_cpSizes - if (myDebug) write(6,'(a)') ' Built CP statistics'; flush(6) + mesh_NcpElems= product(grid(1:2))*grid3 + mesh_NcpElemsGlobal = product(grid) + + mesh_Nnodes = product(grid(1:2) + 1_pInt)*(grid3 + 1_pInt) + call mesh_spectral_build_nodes() + if (myDebug) write(6,'(a)') ' Built nodes'; flush(6) + call theMesh%init(mesh_node) + ! For compatibility + + mesh_maxNips = theMesh%elem%nIPs + mesh_maxNipNeighbors = theMesh%elem%nIPneighbors + mesh_maxNcellnodes = theMesh%elem%Ncellnodes + + + call mesh_spectral_build_elements(FILEUNIT) if (myDebug) write(6,'(a)') ' Built elements'; flush(6) call mesh_build_cellconnectivity @@ -434,8 +444,6 @@ subroutine mesh_init(ip,el) mesh_microstructureAt = mesh_element(4,:) mesh_CPnodeID = mesh_element(5:4+mesh_NipsPerElem,:) !!!!!!!!!!!!!!!!!!!!!!!! - - end subroutine mesh_init @@ -563,10 +571,9 @@ subroutine mesh_build_ipVolumes integer(pInt) :: e,t,g,c,i,m,f,n real(pReal), dimension(FE_maxNcellnodesPerCellface,FE_maxNcellfaces) :: subvolume - if (.not. allocated(mesh_ipVolume)) then - allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems)) - mesh_ipVolume = 0.0_pReal - endif + + allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems),source=0.0_pReal) + !$OMP PARALLEL DO PRIVATE(t,g,c,m,subvolume) do e = 1_pInt,mesh_NcpElems ! loop over cpElems @@ -894,43 +901,6 @@ integer(pInt) function mesh_spectral_getHomogenization(fileUnit) end function mesh_spectral_getHomogenization -!-------------------------------------------------------------------------------------------------- -!> @brief Count overall number of nodes and elements in mesh and stores them in -!! 'mesh_Nelems', 'mesh_Nnodes' and 'mesh_NcpElems' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_count() - - implicit none - - mesh_NcpElems= product(grid(1:2))*grid3 - mesh_Nnodes = product(grid(1:2) + 1_pInt)*(grid3 + 1_pInt) - - mesh_NcpElemsGlobal = product(grid) - -end subroutine mesh_spectral_count - - -!-------------------------------------------------------------------------------------------------- -!> @brief Gets maximum count of nodes, IPs, IP neighbors, and subNodes among cpElements. -!! Sets global values 'mesh_maxNips', 'mesh_maxNipNeighbors', -!! and 'mesh_maxNcellnodes' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_count_cpSizes - - implicit none - integer(pInt) :: t,g,c - - t = 10_pInt - g = FE_geomtype(t) - c = FE_celltype(g) - - mesh_maxNips = FE_Nips(g) - mesh_maxNipNeighbors = FE_NipNeighbors(c) - mesh_maxNcellnodes = FE_Ncellnodes(g) - -end subroutine mesh_spectral_count_cpSizes - - !-------------------------------------------------------------------------------------------------- !> @brief Store x,y,z coordinates of all nodes in mesh. !! Allocates global arrays 'mesh_node0' and 'mesh_node' @@ -941,7 +911,6 @@ subroutine mesh_spectral_build_nodes() integer(pInt) :: n allocate (mesh_node0 (3,mesh_Nnodes), source = 0.0_pReal) - allocate (mesh_node (3,mesh_Nnodes), source = 0.0_pReal) forall (n = 0_pInt:mesh_Nnodes-1_pInt) mesh_node0(1,n+1_pInt) = mesh_unitlength * & @@ -986,7 +955,6 @@ subroutine mesh_spectral_build_elements(fileUnit) headerLength = 0_pInt, & maxDataPerLine, & homog, & - elemType, & elemOffset integer(pInt), dimension(:), allocatable :: & microstructures, & @@ -1047,13 +1015,13 @@ subroutine mesh_spectral_build_elements(fileUnit) enddo enddo - elemType = 10_pInt + elemOffset = product(grid(1:2))*grid3Offset e = 0_pInt do while (e < mesh_NcpElems) ! fill expected number of elements, stop at end of data (or blank line!) e = e+1_pInt ! valid element entry mesh_element( 1,e) = -1_pInt ! DEPRECATED - mesh_element( 2,e) = elemType ! elem type + mesh_element( 2,e) = 10_pInt mesh_element( 3,e) = homog ! homogenization mesh_element( 4,e) = microGlobal(e+elemOffset) ! microstructure mesh_element( 5,e) = e + (e-1_pInt)/grid(1) + & diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 2cca47239..5607791fb 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -417,8 +417,9 @@ type, public, extends(tMesh) :: tMesh_marc integer(pInt), dimension(:,:), allocatable :: & mesh_mapElemSet !< list of elements in elementSet - contains - procedure :: init => tMesh_marc_init + contains + procedure, pass(self) :: tMesh_marc_init + generic, public :: init => tMesh_marc_init end type tMesh_marc type(tMesh_marc), public, protected :: theMesh @@ -426,10 +427,15 @@ end type tMesh_marc contains -subroutine tMesh_marc_init(self) +subroutine tMesh_marc_init(self,elemType,nodes) + implicit none class(tMesh_marc) :: self - + real(pReal), dimension(:,:), intent(in) :: nodes + integer(pInt), intent(in) :: elemType + + call self%tMesh%init('mesh',elemType,nodes) + end subroutine tMesh_marc_init !-------------------------------------------------------------------------------------------------- @@ -553,7 +559,8 @@ subroutine mesh_init(ip,el) mesh_microstructureAt = mesh_element(4,:) mesh_CPnodeID = mesh_element(5:4+mesh_NipsPerElem,:) !!!!!!!!!!!!!!!!!!!!!!!! - + call theMesh%init(mesh_element(2,1),mesh_node0) + end subroutine mesh_init