diff --git a/src/FEM_mech.f90 b/src/FEM_mech.f90 index 50bb68edd..d05e3a184 100644 --- a/src/FEM_mech.f90 +++ b/src/FEM_mech.f90 @@ -63,7 +63,6 @@ use PETScDMplex FEM_mech_destroy external :: & - MPI_Allreduce, & MatZeroRowsColumnsLocalIS, & PetscQuadratureCreate, & PetscFECreateDefault, & @@ -189,7 +188,7 @@ subroutine FEM_mech_init(fieldBC) do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries if (fieldBC%componentBC(field)%Mask(faceSet)) then numBC = numBC + 1 - call ISCreateGeneral(PETSC_COMM_WORLD,1,field-1,PETSC_COPY_VALUES,bcComps(numBC),ierr) + call ISCreateGeneral(PETSC_COMM_WORLD,1,[field-1],PETSC_COPY_VALUES,bcComps(numBC),ierr) CHKERRQ(ierr) call DMGetStratumSize(mech_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,ierr) CHKERRQ(ierr) diff --git a/src/FEM_mesh.f90 b/src/FEM_mesh.f90 deleted file mode 100644 index 82b91ddc9..000000000 --- a/src/FEM_mesh.f90 +++ /dev/null @@ -1,446 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @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 - use, intrinsic :: iso_c_binding - use prec, only: pReal, pInt - - implicit none -#include - private - 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_maxNnodes, & !< max number of nodes in any CP element - mesh_maxNips, & !< max number of IPs in any CP element - mesh_maxNipNeighbors, & - mesh_Nelems !< total number of elements in mesh - - 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) - - integer(pInt), dimension(:,:), allocatable, target, private :: & - mesh_mapFEtoCPelem, & !< [sorted FEid, corresponding CPid] - mesh_mapFEtoCPnode !< [sorted FEid, corresponding CPid] - - DM, public :: geomMesh - - integer(pInt), dimension(:), allocatable, public, protected :: & - mesh_boundaries - -! 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, public :: & - FE_Nelemtypes = 1_pInt, & - FE_Ngeomtypes = 1_pInt, & - FE_Ncelltypes = 1_pInt, & - FE_maxNnodes = 1_pInt, & - FE_maxNips = 14_pInt - - integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_geomtype = & !< geometry type of particular element type - int([1],pInt) - - integer(pInt), dimension(FE_Ngeomtypes), parameter, public :: FE_celltype = & !< cell type that is used by each geometry type - int([1],pInt) - - integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_Nnodes = & !< number of nodes that constitute a specific type of element - int([0],pInt) - - integer(pInt), dimension(FE_Ngeomtypes), public :: FE_Nips = & !< number of IPs in a specific type of element - int([0],pInt) - - integer(pInt), dimension(FE_Ncelltypes), parameter, public :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type - int([6],pInt) - - - public :: & - mesh_init, & - mesh_FEasCP, & - mesh_FEM_build_ipVolumes, & - mesh_FEM_build_ipCoordinates, & - mesh_cellCenterCoordinates - - external :: & - MPI_abort, & - MPI_Bcast, & - DMClone, & - DMGetDimension, & - DMPlexCreateFromFile, & - DMPlexDistribute, & - DMPlexCopyCoordinates, & - DMGetStratumSize, & - DMPlexGetHeightStratum, & - DMPlexGetLabelValue, & - DMPlexSetLabelValue, & - DMDestroy - -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, & - calcMode - 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 - - - if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- mesh init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - endif - - if (allocated(mesh_mapFEtoCPelem)) deallocate(mesh_mapFEtoCPelem) - if (allocated(mesh_mapFEtoCPnode)) deallocate(mesh_mapFEtoCPnode) - if (allocated(mesh_node0)) deallocate(mesh_node0) - if (allocated(mesh_node)) deallocate(mesh_node) - if (allocated(mesh_element)) deallocate(mesh_element) - if (allocated(mesh_ipCoordinates)) deallocate(mesh_ipCoordinates) - if (allocated(mesh_ipVolume)) deallocate(mesh_ipVolume) - - 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_Nelems,ierr) - CHKERRQ(ierr) - call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr) - CHKERRQ(ierr) - mesh_NcpElems = mesh_Nelems - call mesh_FEM_mapNodesAndElems - - FE_Nips(FE_geomtype(1_pInt)) = FEM_Zoo_nQuadrature(dimPlex,integrationOrder) - mesh_maxNnodes = FE_Nnodes(1_pInt) - 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_maxNnodes,mesh_NcpElems)); mesh_element = 0_pInt - do j = 1, mesh_NcpElems - mesh_element( 1,j) = j - mesh_element( 2,j) = 1_pInt ! 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 (usePingPong .and. (mesh_Nelems /= mesh_NcpElems)) & - call IO_error(600_pInt) ! ping-pong must be disabled when having non-DAMASK elements - 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 - - if (allocated(calcMode)) deallocate(calcMode) - allocate(calcMode(mesh_maxNips,mesh_NcpElems)) - 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" - -end subroutine mesh_init - -!-------------------------------------------------------------------------------------------------- -!> @brief Gives the FE to CP ID mapping by binary search through lookup array -!! valid questions (what) are 'elem', 'node' -!-------------------------------------------------------------------------------------------------- -integer(pInt) function mesh_FEasCP(what,myID) - use IO, only: & - IO_lc - - implicit none - character(len=*), intent(in) :: what - integer(pInt), intent(in) :: myID - - integer(pInt), dimension(:,:), pointer :: lookupMap - integer(pInt) :: lower,upper,center - - mesh_FEasCP = 0_pInt - select case(IO_lc(what(1:4))) - case('elem') - lookupMap => mesh_mapFEtoCPelem - case('node') - lookupMap => mesh_mapFEtoCPnode - case default - return - endselect - - lower = 1_pInt - upper = int(size(lookupMap,2_pInt),pInt) - - if (lookupMap(1_pInt,lower) == myID) then ! check at bounds QUESTION is it valid to extend bounds by 1 and just do binary search w/o init check at bounds? - mesh_FEasCP = lookupMap(2_pInt,lower) - return - elseif (lookupMap(1_pInt,upper) == myID) then - mesh_FEasCP = lookupMap(2_pInt,upper) - return - endif - - binarySearch: do while (upper-lower > 1_pInt) - center = (lower+upper)/2_pInt - if (lookupMap(1_pInt,center) < myID) then - lower = center - elseif (lookupMap(1_pInt,center) > myID) then - upper = center - else - mesh_FEasCP = lookupMap(2_pInt,center) - exit - endif - enddo binarySearch - -end function mesh_FEasCP - - -!-------------------------------------------------------------------------------------------------- -!> @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 - - if (.not. allocated(mesh_ipCoordinates)) then - allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems)) - mesh_ipCoordinates = 0.0_pReal - endif - - 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 - - -!-------------------------------------------------------------------------------------------------- -!> @brief fake map node from FE ID to internal (consecutive) representation for node and element -!! Allocates global array 'mesh_mapFEtoCPnode' and 'mesh_mapFEtoCPelem' -!-------------------------------------------------------------------------------------------------- -subroutine mesh_FEM_mapNodesAndElems - use math, only: & - math_range - - implicit none - allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes), source = 0_pInt) - allocate (mesh_mapFEtoCPelem(2_pInt,mesh_NcpElems), source = 0_pInt) - - mesh_mapFEtoCPnode = spread(math_range(mesh_Nnodes),1,2) - mesh_mapFEtoCPelem = spread(math_range(mesh_NcpElems),1,2) - -end subroutine mesh_FEM_mapNodesAndElems - - -end module mesh diff --git a/src/FEM_utilities.f90 b/src/FEM_utilities.f90 index 1b1c33b3a..4947fb0c7 100644 --- a/src/FEM_utilities.f90 +++ b/src/FEM_utilities.f90 @@ -141,7 +141,6 @@ use PETScis COMPONENT_MGTWIN_PHI_ID external :: & - MPI_Allreduce, & PetscOptionsInsertString, & PetscObjectSetName, & DMPlexGetHeightStratum, & diff --git a/src/meshFEM.f90 b/src/meshFEM.f90 index 7dc5c93af..ee11a37bd 100644 --- a/src/meshFEM.f90 +++ b/src/meshFEM.f90 @@ -97,7 +97,6 @@ use PETScis mesh_cellCenterCoordinates external :: & - MPI_Bcast, & DMPlexCreateFromFile, & DMPlexDistribute, & DMPlexCopyCoordinates, &