diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index cb08a6482..17eb677d2 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -6,95 +6,104 @@ !> @brief Sets up the mesh for the solvers MSC.Marc, Abaqus and the spectral solver !-------------------------------------------------------------------------------------------------- module mesh - use, intrinsic :: iso_c_binding - use prec - use geometry_plastic_nonlocal - use mesh_base + use, intrinsic :: iso_c_binding + use prec + use geometry_plastic_nonlocal + use mesh_base +#include + use PETScsys + use DAMASK_interface + use IO + use debug + use numerics + use FEsolving - implicit none - private - - integer, public, protected :: & - mesh_Nnodes - - integer, dimension(:), allocatable, private :: & - microGlobal - integer, dimension(:), allocatable, private :: & - mesh_homogenizationAt - - integer, dimension(:,:), allocatable, public, protected :: & - mesh_element !< entryCount and list of elements containing node - - integer, dimension(:,:,:,:), allocatable, public, protected :: & - mesh_ipNeighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] - - real(pReal), public, protected :: & - 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!!!) - - - real(pReal), dimension(:,:), allocatable, public, protected :: & - mesh_ipVolume, & !< volume associated with IP (initially!) - mesh_node0 !< node x,y,z coordinates (initially!) - - real(pReal), dimension(:,:,:), allocatable, public, protected :: & - mesh_ipArea !< area of interface to neighboring IP (initially!) - - real(pReal), dimension(:,:,:), allocatable, public :: & - mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) - - real(pReal),dimension(:,:,:,:), allocatable, public, protected :: & - mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!) - - logical, dimension(3), public, parameter :: mesh_periodicSurface = .true. !< flag indicating periodic outer surfaces (used for fluxes) - - -! grid specific - integer, dimension(3), public, protected :: & - grid !< (global) grid - integer, public, protected :: & - mesh_NcpElemsGlobal, & !< total number of CP elements in global mesh - grid3, & !< (local) grid in 3rd direction - grid3Offset !< (local) grid offset in 3rd direction - real(pReal), dimension(3), public, protected :: & - geomSize - real(pReal), public, protected :: & - size3, & !< (local) size in 3rd direction - size3offset !< (local) size offset in 3rd direction - - public :: & - mesh_init - - private :: & - mesh_build_ipAreas, & - mesh_build_ipNormals, & - mesh_spectral_build_nodes, & - mesh_spectral_build_elements, & - mesh_spectral_build_ipNeighborhood, & - mesh_build_ipCoordinates - - type, public, extends(tMesh) :: tMesh_grid - integer, dimension(3), public :: & - grid !< (global) grid - integer, public :: & - mesh_NcpElemsGlobal, & !< total number of CP elements in global mesh - grid3, & !< (local) grid in 3rd direction - grid3Offset !< (local) grid offset in 3rd direction - real(pReal), dimension(3), public :: & - geomSize - real(pReal), public :: & - size3, & !< (local) size in 3rd direction - size3offset - - contains - procedure, pass(self) :: tMesh_grid_init - generic, public :: init => tMesh_grid_init - end type tMesh_grid + implicit none + private - type(tMesh_grid), public, protected :: theMesh + include 'fftw3-mpi.f03' + integer, public, protected :: & + mesh_Nnodes + + integer, dimension(:), allocatable, private :: & + microGlobal + integer, dimension(:), allocatable, private :: & + mesh_homogenizationAt + + integer, dimension(:,:), allocatable, public, protected :: & + mesh_element !< entryCount and list of elements containing node + + integer, dimension(:,:,:,:), allocatable, public, protected :: & + mesh_ipNeighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] + + real(pReal), public, protected :: & + 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!!!) + + + real(pReal), dimension(:,:), allocatable, public, protected :: & + mesh_ipVolume, & !< volume associated with IP (initially!) + mesh_node0 !< node x,y,z coordinates (initially!) + + real(pReal), dimension(:,:,:), allocatable, public, protected :: & + mesh_ipArea !< area of interface to neighboring IP (initially!) + + real(pReal), dimension(:,:,:), allocatable, public :: & + mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) + + real(pReal),dimension(:,:,:,:), allocatable, public, protected :: & + mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!) + + logical, dimension(3), public, parameter :: mesh_periodicSurface = .true. !< flag indicating periodic outer surfaces (used for fluxes) + + + ! grid specific + integer, dimension(3), public, protected :: & + grid !< (global) grid + integer, public, protected :: & + mesh_NcpElemsGlobal, & !< total number of CP elements in global mesh + grid3, & !< (local) grid in 3rd direction + grid3Offset !< (local) grid offset in 3rd direction + real(pReal), dimension(3), public, protected :: & + geomSize + real(pReal), public, protected :: & + size3, & !< (local) size in 3rd direction + size3offset !< (local) size offset in 3rd direction + + public :: & + mesh_init + + private :: & + mesh_build_ipAreas, & + mesh_build_ipNormals, & + mesh_spectral_build_nodes, & + mesh_spectral_build_elements, & + mesh_spectral_build_ipNeighborhood, & + mesh_build_ipCoordinates + + type, public, extends(tMesh) :: tMesh_grid + + integer, dimension(3), public :: & + grid !< (global) grid + integer, public :: & + mesh_NcpElemsGlobal, & !< total number of CP elements in global mesh + grid3, & !< (local) grid in 3rd direction + grid3Offset !< (local) grid offset in 3rd direction + real(pReal), dimension(3), public :: & + geomSize + real(pReal), public :: & + size3, & !< (local) size in 3rd direction + size3offset + + contains + procedure, pass(self) :: tMesh_grid_init + generic, public :: init => tMesh_grid_init + end type tMesh_grid + + type(tMesh_grid), public, protected :: theMesh contains @@ -103,7 +112,7 @@ subroutine tMesh_grid_init(self,nodes) class(tMesh_grid) :: self real(pReal), dimension(:,:), intent(in) :: nodes - call self%tMesh%init('grid',10_pInt,nodes) + call self%tMesh%init('grid',10,nodes) end subroutine tMesh_grid_init @@ -113,101 +122,82 @@ end subroutine tMesh_grid_init !-------------------------------------------------------------------------------------------------- subroutine mesh_init(ip,el) -#include - use PETScsys + integer(C_INTPTR_T) :: devNull, local_K, local_K_offset + integer :: ierr, worldsize, j + integer, intent(in), optional :: el, ip + logical :: myDebug - use DAMASK_interface - use IO, only: & - IO_error - use debug, only: & - debug_e, & - debug_i, & - debug_level, & - debug_mesh, & - debug_levelBasic - use numerics, only: & - numerics_unitlength - use FEsolving, only: & - FEsolving_execElem, & - FEsolving_execIP + write(6,'(/,a)') ' <<<+- mesh init -+>>>' - include 'fftw3-mpi.f03' - integer(C_INTPTR_T) :: devNull, local_K, local_K_offset - integer :: ierr, worldsize, j - integer, intent(in), optional :: el, ip - logical :: myDebug + mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh - write(6,'(/,a)') ' <<<+- mesh init -+>>>' + myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0) - mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh - - myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) - - call fftw_mpi_init() - call mesh_spectral_read_grid() + call fftw_mpi_init() + call mesh_spectral_read_grid() - call MPI_comm_size(PETSC_COMM_WORLD, worldsize, ierr) - if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_comm_size') - if(worldsize>grid(3)) call IO_error(894_pInt, ext_msg='number of processes exceeds grid(3)') + call MPI_comm_size(PETSC_COMM_WORLD, worldsize, ierr) + if(ierr /=0) call IO_error(894, ext_msg='MPI_comm_size') + if(worldsize>grid(3)) call IO_error(894, ext_msg='number of processes exceeds grid(3)') - devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T), & - int(grid(2),C_INTPTR_T), & - int(grid(1),C_INTPTR_T)/2+1, & - PETSC_COMM_WORLD, & - local_K, & ! domain grid size along z - local_K_offset) ! domain grid offset along z - grid3 = int(local_K,pInt) - 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) + devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T), & + int(grid(2),C_INTPTR_T), & + int(grid(1),C_INTPTR_T)/2+1, & + PETSC_COMM_WORLD, & + local_K, & ! domain grid size along z + local_K_offset) ! domain grid offset along z + grid3 = int(local_K,pInt) + 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) - mesh_NcpElemsGlobal = product(grid) + mesh_NcpElemsGlobal = product(grid) - mesh_Nnodes = product(grid(1:2) + 1_pInt)*(grid3 + 1_pInt) + mesh_Nnodes = product(grid(1:2) + 1)*(grid3 + 1) - mesh_node0 = mesh_spectral_build_nodes() - mesh_node = mesh_node0 - if (myDebug) write(6,'(a)') ' Built nodes'; flush(6) + mesh_node0 = mesh_spectral_build_nodes() + mesh_node = mesh_node0 + if (myDebug) write(6,'(a)') ' Built nodes'; flush(6) - call theMesh%init(mesh_node) - call theMesh%setNelems(product(grid(1:2))*grid3) - 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) - - + call theMesh%init(mesh_node) + call theMesh%setNelems(product(grid(1:2))*grid3) + 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) + + - 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) - mesh_ipArea = mesh_build_ipAreas() - mesh_ipAreaNormal = mesh_build_ipNormals() - if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) + 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) + 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 - call geometry_plastic_nonlocal_set_IPneighborhood(mesh_ipNeighborhood) + call mesh_spectral_build_ipNeighborhood + call geometry_plastic_nonlocal_set_IPneighborhood(mesh_ipNeighborhood) - if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6) + if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6) - if (debug_e < 1 .or. debug_e > theMesh%nElems) & - call IO_error(602_pInt,ext_msg='element') ! selected element does not exist - if (debug_i < 1 .or. debug_i > theMesh%elem%nIPs) & - call IO_error(602_pInt,ext_msg='IP') ! selected element does not have requested IP + 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 - FEsolving_execElem = [ 1_pInt,theMesh%nElems ] ! parallel loop bounds set to comprise all DAMASK elements - allocate(FEsolving_execIP(2_pInt,theMesh%nElems), source=1_pInt) ! parallel loop bounds set to comprise from first IP... - forall (j = 1_pInt:theMesh%nElems) FEsolving_execIP(2,j) = theMesh%elem%nIPs ! ...up to own IP count for each element + 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... + forall (j = 1:theMesh%nElems) FEsolving_execIP(2,j) = theMesh%elem%nIPs ! ...up to own IP count for each element !!!! COMPATIBILITY HACK !!!! - theMesh%homogenizationAt = mesh_element(3,:) - theMesh%microstructureAt = mesh_element(4,:) + theMesh%homogenizationAt = mesh_element(3,:) + theMesh%microstructureAt = mesh_element(4,:) !!!!!!!!!!!!!!!!!!!!!!!! end subroutine mesh_init @@ -219,22 +209,13 @@ end subroutine mesh_init ! supposed to be called only once! !-------------------------------------------------------------------------------------------------- subroutine mesh_spectral_read_grid() - use IO, only: & - IO_stringPos, & - IO_lc, & - IO_stringValue, & - IO_intValue, & - IO_floatValue, & - IO_error - use DAMASK_interface, only: & - geometryFile character(len=:), allocatable :: rawData character(len=65536) :: line integer, allocatable, dimension(:) :: chunkPos - integer :: h =- 1_pInt + integer :: h =- 1 integer :: & - headerLength = -1_pInt, & !< length of header (in lines) + headerLength = -1, & !< length of header (in lines) fileLength, & !< length of the geom file (in characters) fileUnit, & startPos, endPos, & @@ -245,7 +226,7 @@ subroutine mesh_spectral_read_grid() e, & !< "element", i.e. spectral collocation point i, j - grid = -1_pInt + grid = -1 geomSize = -1.0_pReal !-------------------------------------------------------------------------------------------------- @@ -253,7 +234,7 @@ subroutine mesh_spectral_read_grid() inquire(file = trim(geometryFile), size=fileLength) open(newunit=fileUnit, file=trim(geometryFile), access='stream',& status='old', position='rewind', action='read',iostat=myStat) - if(myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=trim(geometryFile)) + if(myStat /= 0) call IO_error(100,ext_msg=trim(geometryFile)) allocate(character(len=fileLength)::rawData) read(fileUnit) rawData close(fileUnit) @@ -263,106 +244,106 @@ subroutine mesh_spectral_read_grid() endPos = index(rawData,new_line('')) if(endPos <= index(rawData,'head')) then startPos = len(rawData) - call IO_error(error_ID=841_pInt, ext_msg='mesh_spectral_read_grid') + call IO_error(error_ID=841, ext_msg='mesh_spectral_read_grid') else chunkPos = IO_stringPos(rawData(1:endPos)) - if (chunkPos(1) < 2_pInt) call IO_error(error_ID=841_pInt, ext_msg='mesh_spectral_read_grid') - headerLength = IO_intValue(rawData(1:endPos),chunkPos,1_pInt) - startPos = endPos + 1_pInt + if (chunkPos(1) < 2) call IO_error(error_ID=841, ext_msg='mesh_spectral_read_grid') + headerLength = IO_intValue(rawData(1:endPos),chunkPos,1) + startPos = endPos + 1 endif !-------------------------------------------------------------------------------------------------- ! read and interprete header l = 0 do while (l < headerLength .and. startPos < len(rawData)) - endPos = startPos + index(rawData(startPos:),new_line('')) - 1_pInt + endPos = startPos + index(rawData(startPos:),new_line('')) - 1 if (endPos < startPos) endPos = len(rawData) ! end of file without new line line = rawData(startPos:endPos) - startPos = endPos + 1_pInt - l = l + 1_pInt + startPos = endPos + 1 + l = l + 1 chunkPos = IO_stringPos(trim(line)) if (chunkPos(1) < 2) cycle ! need at least one keyword value pair - select case ( IO_lc(IO_StringValue(trim(line),chunkPos,1_pInt,.true.)) ) + select case ( IO_lc(IO_StringValue(trim(line),chunkPos,1,.true.)) ) case ('grid') if (chunkPos(1) > 6) then - do j = 2_pInt,6_pInt,2_pInt + do j = 2,6,2 select case (IO_lc(IO_stringValue(line,chunkPos,j))) case('a') - grid(1) = IO_intValue(line,chunkPos,j+1_pInt) + grid(1) = IO_intValue(line,chunkPos,j+1) case('b') - grid(2) = IO_intValue(line,chunkPos,j+1_pInt) + grid(2) = IO_intValue(line,chunkPos,j+1) case('c') - grid(3) = IO_intValue(line,chunkPos,j+1_pInt) + grid(3) = IO_intValue(line,chunkPos,j+1) end select enddo endif case ('size') if (chunkPos(1) > 6) then - do j = 2_pInt,6_pInt,2_pInt + do j = 2,6,2 select case (IO_lc(IO_stringValue(line,chunkPos,j))) case('x') - geomSize(1) = IO_floatValue(line,chunkPos,j+1_pInt) + geomSize(1) = IO_floatValue(line,chunkPos,j+1) case('y') - geomSize(2) = IO_floatValue(line,chunkPos,j+1_pInt) + geomSize(2) = IO_floatValue(line,chunkPos,j+1) case('z') - geomSize(3) = IO_floatValue(line,chunkPos,j+1_pInt) + geomSize(3) = IO_floatValue(line,chunkPos,j+1) end select enddo endif case ('homogenization') - if (chunkPos(1) > 1) h = IO_intValue(line,chunkPos,2_pInt) + if (chunkPos(1) > 1) h = IO_intValue(line,chunkPos,2) end select enddo !-------------------------------------------------------------------------------------------------- ! sanity checks - if(h < 1_pInt) & - call IO_error(error_ID = 842_pInt, ext_msg='homogenization (mesh_spectral_read_grid)') - if(any(grid < 1_pInt)) & - call IO_error(error_ID = 842_pInt, ext_msg='grid (mesh_spectral_read_grid)') + if(h < 1) & + call IO_error(error_ID = 842, ext_msg='homogenization (mesh_spectral_read_grid)') + if(any(grid < 1)) & + call IO_error(error_ID = 842, ext_msg='grid (mesh_spectral_read_grid)') if(any(geomSize < 0.0_pReal)) & - call IO_error(error_ID = 842_pInt, ext_msg='size (mesh_spectral_read_grid)') + call IO_error(error_ID = 842, ext_msg='size (mesh_spectral_read_grid)') - allocate(microGlobal(product(grid)), source = -1_pInt) + allocate(microGlobal(product(grid)), source = -1) allocate(mesh_homogenizationAt(product(grid)), source = h) ! too large in case of MPI (shrink later, not very elegant) !-------------------------------------------------------------------------------------------------- ! read and interpret content - e = 1_pInt + e = 1 do while (startPos < len(rawData)) - endPos = startPos + index(rawData(startPos:),new_line('')) - 1_pInt + endPos = startPos + index(rawData(startPos:),new_line('')) - 1 if (endPos < startPos) endPos = len(rawData) ! end of file without new line line = rawData(startPos:endPos) - startPos = endPos + 1_pInt - l = l + 1_pInt + startPos = endPos + 1 + l = l + 1 chunkPos = IO_stringPos(trim(line)) noCompression: if (chunkPos(1) /= 3) then c = chunkPos(1) - microGlobal(e:e+c-1_pInt) = [(IO_intValue(line,chunkPos,i+1_pInt), i=0_pInt, c-1_pInt)] + microGlobal(e:e+c-1) = [(IO_intValue(line,chunkPos,i+1), i=0, c-1)] else noCompression compression: if (IO_lc(IO_stringValue(line,chunkPos,2)) == 'of') then c = IO_intValue(line,chunkPos,1) - microGlobal(e:e+c-1_pInt) = [(IO_intValue(line,chunkPos,3),i = 1_pInt,IO_intValue(line,chunkPos,1))] + microGlobal(e:e+c-1) = [(IO_intValue(line,chunkPos,3),i = 1,IO_intValue(line,chunkPos,1))] else if (IO_lc(IO_stringValue(line,chunkPos,2)) == 'to') then compression - c = abs(IO_intValue(line,chunkPos,3) - IO_intValue(line,chunkPos,1)) + 1_pInt - o = merge(+1_pInt, -1_pInt, IO_intValue(line,chunkPos,3) > IO_intValue(line,chunkPos,1)) - microGlobal(e:e+c-1_pInt) = [(i, i = IO_intValue(line,chunkPos,1),IO_intValue(line,chunkPos,3),o)] + c = abs(IO_intValue(line,chunkPos,3) - IO_intValue(line,chunkPos,1)) + 1 + o = merge(+1, -1, IO_intValue(line,chunkPos,3) > IO_intValue(line,chunkPos,1)) + microGlobal(e:e+c-1) = [(i, i = IO_intValue(line,chunkPos,1),IO_intValue(line,chunkPos,3),o)] else compression c = chunkPos(1) - microGlobal(e:e+c-1_pInt) = [(IO_intValue(line,chunkPos,i+1_pInt), i=0_pInt, c-1_pInt)] + microGlobal(e:e+c-1) = [(IO_intValue(line,chunkPos,i+1), i=0, c-1)] endif compression endif noCompression e = e+c end do - if (e-1 /= product(grid)) call IO_error(error_ID = 843_pInt, el=e) + if (e-1 /= product(grid)) call IO_error(error_ID = 843, el=e) end subroutine mesh_spectral_read_grid @@ -413,30 +394,30 @@ end function mesh_build_ipCoordinates !> @brief Store FEid, type, material, texture, and node list per element. !! Allocates global array 'mesh_element' !-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_build_elements() - integer :: & - e, & - elemOffset +subroutine mesh_spectral_build_elements + integer :: & + e, & + elemOffset - allocate(mesh_element (4_pInt+8_pInt,theMesh%nElems), source = 0_pInt) + allocate(mesh_element (4+8,theMesh%nElems), source = 0) - elemOffset = product(grid(1:2))*grid3Offset - do e=1, theMesh%nElems - mesh_element( 1,e) = -1_pInt ! DEPRECATED - 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) + & - ((e-1_pInt)/(grid(1)*grid(2)))*(grid(1)+1_pInt) ! base node - mesh_element( 6,e) = mesh_element(5,e) + 1_pInt - mesh_element( 7,e) = mesh_element(5,e) + grid(1) + 2_pInt - mesh_element( 8,e) = mesh_element(5,e) + grid(1) + 1_pInt - mesh_element( 9,e) = mesh_element(5,e) +(grid(1) + 1_pInt) * (grid(2) + 1_pInt) ! second floor base node - mesh_element(10,e) = mesh_element(9,e) + 1_pInt - mesh_element(11,e) = mesh_element(9,e) + grid(1) + 2_pInt - mesh_element(12,e) = mesh_element(9,e) + grid(1) + 1_pInt - enddo + elemOffset = product(grid(1:2))*grid3Offset + do e=1, theMesh%nElems + mesh_element( 1,e) = -1 ! DEPRECATED + mesh_element( 2,e) = -1 ! DEPRECATED + mesh_element( 3,e) = mesh_homogenizationAt(e) + mesh_element( 4,e) = microGlobal(e+elemOffset) ! microstructure + mesh_element( 5,e) = e + (e-1)/grid(1) + & + ((e-1)/(grid(1)*grid(2)))*(grid(1)+1) ! base node + mesh_element( 6,e) = mesh_element(5,e) + 1 + mesh_element( 7,e) = mesh_element(5,e) + grid(1) + 2 + mesh_element( 8,e) = mesh_element(5,e) + grid(1) + 1 + mesh_element( 9,e) = mesh_element(5,e) +(grid(1) + 1) * (grid(2) + 1) ! second floor base node + mesh_element(10,e) = mesh_element(9,e) + 1 + mesh_element(11,e) = mesh_element(9,e) + grid(1) + 2 + mesh_element(12,e) = mesh_element(9,e) + grid(1) + 1 + enddo end subroutine mesh_spectral_build_elements @@ -447,50 +428,50 @@ end subroutine mesh_spectral_build_elements !-------------------------------------------------------------------------------------------------- subroutine mesh_spectral_build_ipNeighborhood - integer :: & - x,y,z, & - e - allocate(mesh_ipNeighborhood(3,6,1,theMesh%nElems),source=0_pInt) + integer :: & + x,y,z, & + e + allocate(mesh_ipNeighborhood(3,6,1,theMesh%nElems),source=0) - e = 0_pInt - do z = 0_pInt,grid3-1_pInt - do y = 0_pInt,grid(2)-1_pInt - do x = 0_pInt,grid(1)-1_pInt - e = e + 1_pInt - mesh_ipNeighborhood(1,1,1,e) = z * grid(1) * grid(2) & - + y * grid(1) & - + modulo(x+1_pInt,grid(1)) & - + 1_pInt - mesh_ipNeighborhood(1,2,1,e) = z * grid(1) * grid(2) & - + y * grid(1) & - + modulo(x-1_pInt,grid(1)) & - + 1_pInt - mesh_ipNeighborhood(1,3,1,e) = z * grid(1) * grid(2) & - + modulo(y+1_pInt,grid(2)) * grid(1) & - + x & - + 1_pInt - mesh_ipNeighborhood(1,4,1,e) = z * grid(1) * grid(2) & - + modulo(y-1_pInt,grid(2)) * grid(1) & - + x & - + 1_pInt - mesh_ipNeighborhood(1,5,1,e) = modulo(z+1_pInt,grid3) * grid(1) * grid(2) & - + y * grid(1) & - + x & - + 1_pInt - mesh_ipNeighborhood(1,6,1,e) = modulo(z-1_pInt,grid3) * grid(1) * grid(2) & - + y * grid(1) & - + x & - + 1_pInt - mesh_ipNeighborhood(2,1:6,1,e) = 1_pInt - mesh_ipNeighborhood(3,1,1,e) = 2_pInt - mesh_ipNeighborhood(3,2,1,e) = 1_pInt - mesh_ipNeighborhood(3,3,1,e) = 4_pInt - mesh_ipNeighborhood(3,4,1,e) = 3_pInt - mesh_ipNeighborhood(3,5,1,e) = 6_pInt - mesh_ipNeighborhood(3,6,1,e) = 5_pInt - enddo - enddo - enddo + e = 0 + do z = 0,grid3-1 + do y = 0,grid(2)-1 + do x = 0,grid(1)-1 + e = e + 1 + mesh_ipNeighborhood(1,1,1,e) = z * grid(1) * grid(2) & + + y * grid(1) & + + modulo(x+1,grid(1)) & + + 1 + mesh_ipNeighborhood(1,2,1,e) = z * grid(1) * grid(2) & + + y * grid(1) & + + modulo(x-1,grid(1)) & + + 1 + mesh_ipNeighborhood(1,3,1,e) = z * grid(1) * grid(2) & + + modulo(y+1,grid(2)) * grid(1) & + + x & + + 1 + mesh_ipNeighborhood(1,4,1,e) = z * grid(1) * grid(2) & + + modulo(y-1,grid(2)) * grid(1) & + + x & + + 1 + mesh_ipNeighborhood(1,5,1,e) = modulo(z+1,grid3) * grid(1) * grid(2) & + + y * grid(1) & + + x & + + 1 + mesh_ipNeighborhood(1,6,1,e) = modulo(z-1,grid3) * grid(1) * grid(2) & + + y * grid(1) & + + x & + + 1 + mesh_ipNeighborhood(2,1:6,1,e) = 1 + mesh_ipNeighborhood(3,1,1,e) = 2 + mesh_ipNeighborhood(3,2,1,e) = 1 + mesh_ipNeighborhood(3,3,1,e) = 4 + mesh_ipNeighborhood(3,4,1,e) = 3 + mesh_ipNeighborhood(3,5,1,e) = 6 + mesh_ipNeighborhood(3,6,1,e) = 5 + enddo + enddo + enddo end subroutine mesh_spectral_build_ipNeighborhood @@ -499,41 +480,37 @@ end subroutine mesh_spectral_build_ipNeighborhood !> @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 - - 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 :: & - i,j,k,n - integer, dimension(3), parameter :: & - diag = 1_pInt - integer, dimension(3) :: & - shift = 0_pInt, & - lookup = 0_pInt, & - me = 0_pInt, & - iRes = 0_pInt - integer, 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]) + + 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 :: & + i,j,k,n + integer, dimension(3), parameter :: & + diag = 1 + integer, dimension(3) :: & + shift = 0, & + lookup = 0, & + me = 0, & + iRes = 0 + integer, dimension(3,8) :: & + neighbor = reshape([ & + 0, 0, 0, & + 1, 0, 0, & + 1, 1, 0, & + 0, 1, 0, & + 0, 0, 1, & + 1, 0, 1, & + 1, 1, 1, & + 0, 1, 1 ], [3,8]) !-------------------------------------------------------------------------------------------------- ! initializing variables @@ -541,43 +518,35 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes) 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) & - - matmul(Favg, real(shift,pReal)*gDim) - endif - enddo; enddo; enddo + wrappedCentres(1:3,2:iRes(1)+1,2:iRes(2)+1,2:iRes(3)+1) = centres + do k = 0,iRes(3)+1 + do j = 0,iRes(2)+1 + do i = 0,iRes(1)+1 + if (k==0 .or. k==iRes(3)+1 .or. & ! z skin + j==0 .or. j==iRes(2)+1 .or. & ! y skin + i==0 .or. i==iRes(1)+1 ) then ! x skin + me = [i,j,k] ! me on skin + shift = sign(abs(iRes+diag-2*me)/(iRes+diag),iRes+diag-2*me) + lookup = me-diag+shift*iRes + wrappedCentres(1:3,i+1, j+1, k+1) = & + centres(1:3,lookup(1)+1,lookup(2)+1,lookup(3)+1) & + - matmul(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 + do k = 0,iRes(3); do j = 0,iRes(2); do i = 0,iRes(1) + do n = 1,8 + nodes(1:3,i+1,j+1,k+1) = & + nodes(1:3,i+1,j+1,k+1) + wrappedCentres(1:3,i+1+neighbor(1,n), & + j+1+neighbor(2,n), & + k+1+neighbor(3,n) ) + enddo + enddo; enddo; enddo + nodes = nodes/8.0_pReal end function mesh_nodesAroundCentres