!-------------------------------------------------------------------------------------------------- !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @brief Driver controlling inner and outer load case looping of the FEM solver !> @details doing cutbacking, forwarding in case of restart, reporting statistics, writing !> results !-------------------------------------------------------------------------------------------------- module mesh #include #include use prec, only: pReal, pInt use PETScdmda use PETScis implicit none private integer(pInt), public, parameter :: & mesh_ElemType=1_pInt !< Element type of the mesh (only support homogeneous meshes) integer(pInt), public, protected :: & mesh_Nboundaries, & mesh_NcpElems, & !< total number of CP elements in mesh mesh_NcpElemsGlobal, & mesh_Nnodes, & !< total number of nodes in mesh mesh_maxNipNeighbors, & !< max number of IP neighbors in any CP element mesh_maxNipNeighbors !!!! BEGIN DEPRECATED !!!!! integer(pInt), public, protected :: & mesh_maxNips, & !< max number of IPs in any CP element !!!! BEGIN DEPRECATED !!!!! real(pReal), public, protected :: charLength integer(pInt), dimension(:,:), allocatable, public, protected :: & mesh_element !< FEid, type(internal representation), material, texture, node indices as CP IDs real(pReal), dimension(:,:), allocatable, public :: & 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 :: & mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) real(pReal), dimension(:,:,:), allocatable, public, protected :: & mesh_ipArea !< area of interface to neighboring IP (initially!) real(pReal),dimension(:,:,:,:), allocatable, public, protected :: & mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!) integer(pInt), dimension(:,:,:,:), allocatable, public, protected :: & mesh_ipNeighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] logical, dimension(3), public, protected :: mesh_periodicSurface !< flag indicating periodic outer surfaces (used for fluxes) DM, public :: geomMesh integer(pInt), dimension(:), allocatable, public, protected :: & mesh_boundaries integer(pInt), dimension(1_pInt), parameter, public :: FE_geomtype = & !< geometry type of particular element type int([1],pInt) integer(pInt), dimension(1_pInt), parameter, public :: FE_celltype = & !< cell type that is used by each geometry type int([1],pInt) integer(pInt), dimension(1_pInt), public :: FE_Nips = & !< number of IPs in a specific type of element int([0],pInt) integer(pInt), dimension(1_pInt), parameter, public :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type int([6],pInt) public :: & mesh_init, & mesh_FEM_build_ipVolumes, & mesh_FEM_build_ipCoordinates, & mesh_cellCenterCoordinates external :: & DMPlexCreateFromFile, & DMPlexDistribute, & DMPlexCopyCoordinates, & DMGetStratumSize, & DMPlexGetHeightStratum, & DMPlexGetLabelValue, & DMPlexSetLabelValue contains !-------------------------------------------------------------------------------------------------- !> @brief initializes the mesh by calling all necessary private routines the mesh module !! Order and routines strongly depend on type of solver !-------------------------------------------------------------------------------------------------- subroutine mesh_init(ip,el) use DAMASK_interface use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use IO, only: & IO_timeStamp, & IO_error, & IO_open_file, & IO_stringPos, & IO_intValue, & IO_EOF, & IO_read, & IO_isBlank use debug, only: & debug_e, & debug_i use numerics, only: & usePingPong, & integrationOrder, & worldrank, & worldsize use FEsolving, only: & FEsolving_execElem, & FEsolving_execIP use FEM_Zoo, only: & FEM_Zoo_nQuadrature, & FEM_Zoo_QuadraturePoints implicit none integer(pInt), parameter :: FILEUNIT = 222_pInt integer(pInt), intent(in) :: el, ip integer(pInt) :: j integer(pInt), allocatable, dimension(:) :: chunkPos integer :: dimPlex character(len=512) :: & line logical :: flag PetscSF :: sf DM :: globalMesh PetscInt :: face, nFaceSets PetscInt, pointer :: pFaceSets(:) IS :: faceSetIS PetscErrorCode :: ierr write(6,'(/,a)') ' <<<+- mesh init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr) CHKERRQ(ierr) call DMGetDimension(globalMesh,dimPlex,ierr) CHKERRQ(ierr) call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr) CHKERRQ(ierr) call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr) CHKERRQ(ierr) call MPI_Bcast(mesh_Nboundaries,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) call MPI_Bcast(mesh_NcpElemsGlobal,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) call MPI_Bcast(dimPlex,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) allocate(mesh_boundaries(mesh_Nboundaries), source = 0_pInt) call DMGetLabelSize(globalMesh,'Face Sets',nFaceSets,ierr) CHKERRQ(ierr) call DMGetLabelIdIS(globalMesh,'Face Sets',faceSetIS,ierr) CHKERRQ(ierr) if (nFaceSets > 0) call ISGetIndicesF90(faceSetIS,pFaceSets,ierr) do face = 1, nFaceSets mesh_boundaries(face) = pFaceSets(face) enddo if (nFaceSets > 0) call ISRestoreIndicesF90(faceSetIS,pFaceSets,ierr) call MPI_Bcast(mesh_boundaries,mesh_Nboundaries,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) if (worldrank == 0) then j = 0 flag = .false. call IO_open_file(FILEUNIT,trim(geometryFile)) do read(FILEUNIT,'(a512)') line if (trim(line) == IO_EOF) exit ! skip empty lines if (trim(line) == '$Elements') then read(FILEUNIT,'(a512)') line read(FILEUNIT,'(a512)') line flag = .true. endif if (trim(line) == '$EndElements') exit if (flag) then chunkPos = IO_stringPos(line) if (chunkPos(1) == 3+IO_intValue(line,chunkPos,3)+dimPlex+1) then call DMSetLabelValue(globalMesh,'material',j,IO_intValue(line,chunkPos,4),ierr) CHKERRQ(ierr) j = j + 1 endif ! count all identifiers to allocate memory and do sanity check endif enddo close (FILEUNIT) endif if (worldsize > 1) then call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr) CHKERRQ(ierr) else call DMClone(globalMesh,geomMesh,ierr) CHKERRQ(ierr) endif call DMDestroy(globalMesh,ierr); CHKERRQ(ierr) call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_NcpElems,ierr) CHKERRQ(ierr) call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr) CHKERRQ(ierr) FE_Nips(FE_geomtype(1_pInt)) = FEM_Zoo_nQuadrature(dimPlex,integrationOrder) mesh_maxNips = FE_Nips(1_pInt) call mesh_FEM_build_ipCoordinates(dimPlex,FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p) call mesh_FEM_build_ipVolumes(dimPlex) allocate (mesh_element (4_pInt,mesh_NcpElems)); mesh_element = 0_pInt do j = 1, mesh_NcpElems mesh_element( 1,j) = j mesh_element( 2,j) = mesh_elemType ! elem type mesh_element( 3,j) = 1_pInt ! homogenization call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr) CHKERRQ(ierr) end do if (debug_e < 1 .or. debug_e > mesh_NcpElems) & call IO_error(602_pInt,ext_msg='element') ! selected element does not exist if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2_pInt,debug_e)))) & call IO_error(602_pInt,ext_msg='IP') ! selected element does not have requested IP FEsolving_execElem = [ 1_pInt,mesh_NcpElems ] ! parallel loop bounds set to comprise all DAMASK elements if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP) allocate(FEsolving_execIP(2_pInt,mesh_NcpElems)); FEsolving_execIP = 1_pInt ! parallel loop bounds set to comprise from first IP... forall (j = 1_pInt:mesh_NcpElems) FEsolving_execIP(2,j) = FE_Nips(FE_geomtype(mesh_element(2,j))) ! ...up to own IP count for each element !!!! COMPATIBILITY HACK !!!! ! for a homogeneous mesh, all elements have the same number of IPs and and cell nodes. ! hence, xxPerElem instead of maxXX mesh_NipsPerElem = mesh_maxNips !!!!!!!!!!!!!!!!!!!!!!!! end subroutine mesh_init !-------------------------------------------------------------------------------------------------- !> @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 end function mesh_cellCenterCoordinates !-------------------------------------------------------------------------------------------------- !> @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_FEM_build_ipVolumes(dimPlex) use math, only: & math_I3, & math_det33 implicit none PetscInt :: dimPlex PetscReal :: vol PetscReal, target :: cent(dimPlex), norm(dimPlex) PetscReal, pointer :: pCent(:), pNorm(:) PetscInt :: cellStart, cellEnd, cell PetscErrorCode :: ierr if (.not. allocated(mesh_ipVolume)) then allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems)) mesh_ipVolume = 0.0_pReal endif call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) pCent => cent pNorm => norm do cell = cellStart, cellEnd-1 call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr) CHKERRQ(ierr) mesh_ipVolume(:,cell+1) = vol/real(mesh_maxNips,pReal) enddo end subroutine mesh_FEM_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. ! 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. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !-------------------------------------------------------------------------------------------------- subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints) implicit none PetscInt, intent(in) :: dimPlex PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex) PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), invcellJ(dimPlex*dimPlex) PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:) PetscReal :: detJ PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset PetscErrorCode :: ierr allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal) pV0 => v0 pCellJ => cellJ pInvcellJ => invcellJ call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) do cell = cellStart, cellEnd-1 !< loop over all elements call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) CHKERRQ(ierr) qOffset = 0 do qPt = 1, mesh_maxNips do dirI = 1, dimPlex mesh_ipCoordinates(dirI,qPt,cell+1) = pV0(dirI) do dirJ = 1, dimPlex mesh_ipCoordinates(dirI,qPt,cell+1) = mesh_ipCoordinates(dirI,qPt,cell+1) + & pCellJ((dirI-1)*dimPlex+dirJ)*(qPoints(qOffset+dirJ) + 1.0) enddo enddo qOffset = qOffset + dimPlex enddo enddo end subroutine mesh_FEM_build_ipCoordinates end module mesh