223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

cleaning

This commit is contained in:
Martin Diehl 2020-06-16 13:27:21 +02:00
parent 512e54a7ee
commit 753fbb70fd
2 changed files with 67 additions and 81 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/discretization.f90
View File

@ -49,7 +49,7 @@ subroutine discretization_init(homogenizationAt,microstructureAt,&
IPcoords0, & IPcoords0, &
NodeCoords0 NodeCoords0
integer, optional, intent(in) :: & integer, optional, intent(in) :: &
sharedNodesBegin sharedNodesBegin ! index of first node shared among different processes (MPI)
write(6,'(/,a)') ' <<<+- discretization init -+>>>'; flush(6) write(6,'(/,a)') ' <<<+- discretization init -+>>>'; flush(6)

146
src/mesh/discretization_mesh.f90
View File

@ -8,53 +8,50 @@ module discretization_mesh
#include <petsc/finclude/petscdmplex.h> #include <petsc/finclude/petscdmplex.h>
#include <petsc/finclude/petscis.h> #include <petsc/finclude/petscis.h>
#include <petsc/finclude/petscdmda.h> #include <petsc/finclude/petscdmda.h>
use PETScdmplex use PETScdmplex
use PETScdmda use PETScdmda
use PETScis use PETScis
use DAMASK_interface use DAMASK_interface
use IO use IO
use debug use debug
use discretization use discretization
use numerics use numerics
use FEsolving use FEsolving
use FEM_quadrature use FEM_quadrature
use prec use prec
implicit none
private
integer, public, protected :: &
mesh_Nboundaries, &
mesh_NcpElemsGlobal
integer :: & implicit none
mesh_NcpElems !< total number of CP elements in mesh private
integer, public, protected :: &
mesh_Nboundaries, &
mesh_NcpElemsGlobal
integer :: &
mesh_NcpElems !< total number of CP elements in mesh
!!!! BEGIN DEPRECATED !!!!! !!!! BEGIN DEPRECATED !!!!!
integer, public, protected :: & integer, public, protected :: &
mesh_maxNips !< max number of IPs in any CP element mesh_maxNips !< max number of IPs in any CP element
!!!! BEGIN DEPRECATED !!!!! !!!! BEGIN DEPRECATED !!!!!
integer, dimension(:,:), allocatable :: & real(pReal), dimension(:,:), allocatable :: &
mesh_element !DEPRECATED mesh_ipVolume, & !< volume associated with IP (initially!)
mesh_node0 !< node x,y,z coordinates (initially!)
real(pReal), dimension(:,:), allocatable :: & real(pReal), dimension(:,:,:), allocatable :: &
mesh_ipVolume, & !< volume associated with IP (initially!) mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!)
mesh_node0 !< node x,y,z coordinates (initially!)
real(pReal), dimension(:,:,:), allocatable :: &
mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!)
DM, public :: geomMesh DM, public :: geomMesh
PetscInt, dimension(:), allocatable, public, protected :: &
mesh_boundaries
public :: & PetscInt, dimension(:), allocatable, public, protected :: &
discretization_mesh_init, & mesh_boundaries
mesh_FEM_build_ipVolumes, &
mesh_FEM_build_ipCoordinates public :: &
discretization_mesh_init, &
mesh_FEM_build_ipVolumes, &
mesh_FEM_build_ipCoordinates
contains contains
@ -67,39 +64,32 @@ subroutine discretization_mesh_init(restart)
logical, intent(in) :: restart logical, intent(in) :: restart
integer, dimension(1), parameter:: FE_geomtype = [1] !< geometry type of particular element type
integer, dimension(1) :: FE_Nips !< number of IPs in a specific type of element
integer, allocatable, dimension(:) :: chunkPos integer, allocatable, dimension(:) :: chunkPos
integer :: dimPlex, & integer :: dimPlex, &
mesh_Nnodes, & !< total number of nodes in mesh mesh_Nnodes, & !< total number of nodes in mesh
j, l j, l
integer, parameter :: &
mesh_ElemType=1 !< Element type of the mesh (only support homogeneous meshes)
PetscSF :: sf PetscSF :: sf
DM :: globalMesh DM :: globalMesh
PetscInt :: nFaceSets PetscInt :: nFaceSets
PetscInt, pointer, dimension(:) :: pFaceSets PetscInt, pointer, dimension(:) :: pFaceSets
character(len=pStringLen), dimension(:), allocatable :: fileContent character(len=pStringLen), dimension(:), allocatable :: fileContent
IS :: faceSetIS IS :: faceSetIS
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer, dimension(:,:), allocatable :: &
mesh_element
write(6,'(/,a)') ' <<<+- mesh init -+>>>' write(6,'(/,a)') ' <<<+- mesh init -+>>>'
! read in file
call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr) call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
! get spatial dimension (2 or 3?)
call DMGetDimension(globalMesh,dimPlex,ierr) call DMGetDimension(globalMesh,dimPlex,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
write(6,*) 'dimension',dimPlex;flush(6)
call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr) call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
! get number of IDs in face sets (for boundary conditions?) ! get number of IDs in face sets (for boundary conditions?)
call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr) call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
write(6,*) 'number of "Face Sets"',mesh_Nboundaries;flush(6)
call MPI_Bcast(mesh_Nboundaries,1,MPI_INTEGER,0,PETSC_COMM_WORLD,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(mesh_NcpElemsGlobal,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
call MPI_Bcast(dimPlex,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) call MPI_Bcast(dimPlex,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
@ -142,48 +132,44 @@ subroutine discretization_mesh_init(restart)
enddo enddo
call DMClone(globalMesh,geomMesh,ierr) call DMClone(globalMesh,geomMesh,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
else else
call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr) call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
endif endif
call DMDestroy(globalMesh,ierr); CHKERRQ(ierr) call DMDestroy(globalMesh,ierr); CHKERRQ(ierr)
call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_NcpElems,ierr) call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_NcpElems,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr) call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
FE_Nips(FE_geomtype(1)) = FEM_nQuadrature(dimPlex,integrationOrder) mesh_maxNips = FEM_nQuadrature(dimPlex,integrationOrder)
mesh_maxNips = FE_Nips(1)
write(6,*) 'mesh_maxNips',mesh_maxNips
call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,integrationOrder)%p) call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,integrationOrder)%p)
call mesh_FEM_build_ipVolumes(dimPlex) call mesh_FEM_build_ipVolumes(dimPlex)
allocate (mesh_element (4,mesh_NcpElems)); mesh_element = 0 allocate (mesh_element (4,mesh_NcpElems)); mesh_element = 0
do j = 1, mesh_NcpElems do j = 1, mesh_NcpElems
mesh_element( 1,j) = -1 ! DEPRECATED mesh_element( 1,j) = -1
mesh_element( 2,j) = mesh_elemType ! elem type mesh_element( 2,j) = -1
mesh_element( 3,j) = 1 ! homogenization mesh_element( 3,j) = 1 ! homogenization
call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr) call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
end do end do
if (debug_e < 1 .or. debug_e > mesh_NcpElems) call IO_error(602,ext_msg='element')
if (debug_i < 1 .or. debug_i > mesh_maxNips) call IO_error(602,ext_msg='IP')
FEsolving_execElem = [1,mesh_NcpElems] ! parallel loop bounds set to comprise all DAMASK elements
FEsolving_execIP = [1,mesh_maxNips]
if (debug_e < 1 .or. debug_e > mesh_NcpElems) &
call IO_error(602,ext_msg='element') ! selected element does not exist
if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2,debug_e)))) &
call IO_error(602,ext_msg='IP') ! selected element does not have requested IP
FEsolving_execElem = [ 1,mesh_NcpElems ] ! parallel loop bounds set to comprise all DAMASK elements
FEsolving_execIP = [1,FE_Nips(FE_geomtype(mesh_element(2,1)))]
allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal) allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal)
call discretization_init(mesh_element(3,:),mesh_element(4,:),& call discretization_init(mesh_element(3,:),mesh_element(4,:),&
reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), & reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), &
mesh_node0) mesh_node0)
end subroutine discretization_mesh_init end subroutine discretization_mesh_init
@ -191,23 +177,23 @@ end subroutine discretization_mesh_init
!> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume' !> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_FEM_build_ipVolumes(dimPlex) subroutine mesh_FEM_build_ipVolumes(dimPlex)
PetscInt,intent(in):: dimPlex PetscInt,intent(in):: dimPlex
PetscReal :: vol PetscReal :: vol
PetscReal, pointer,dimension(:) :: pCent, pNorm PetscReal, pointer,dimension(:) :: pCent, pNorm
PetscInt :: cellStart, cellEnd, cell PetscInt :: cellStart, cellEnd, cell
PetscErrorCode :: ierr PetscErrorCode :: ierr
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems),source=0.0_pReal) allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
allocate(pCent(dimPlex)) allocate(pCent(dimPlex))
allocate(pNorm(dimPlex)) allocate(pNorm(dimPlex))
do cell = cellStart, cellEnd-1 do cell = cellStart, cellEnd-1
call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr) call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
mesh_ipVolume(:,cell+1) = vol/real(mesh_maxNips,pReal) mesh_ipVolume(:,cell+1) = vol/real(mesh_maxNips,pReal)
enddo enddo
end subroutine mesh_FEM_build_ipVolumes end subroutine mesh_FEM_build_ipVolumes
@ -219,20 +205,20 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
PetscInt, intent(in) :: dimPlex PetscInt, intent(in) :: dimPlex
PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex) PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex)
PetscReal, pointer,dimension(:) :: pV0, pCellJ, pInvcellJ PetscReal, pointer,dimension(:) :: pV0, pCellJ, pInvcellJ
PetscReal :: detJ PetscReal :: detJ
PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset
PetscErrorCode :: ierr PetscErrorCode :: ierr
allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal) allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
allocate(pV0(dimPlex)) allocate(pV0(dimPlex))
allocatE(pCellJ(dimPlex**2)) allocatE(pCellJ(dimPlex**2))
allocatE(pinvCellJ(dimPlex**2)) allocatE(pinvCellJ(dimPlex**2))
call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
qOffset = 0 qOffset = 0
@ -246,7 +232,7 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
enddo enddo
qOffset = qOffset + dimPlex qOffset = qOffset + dimPlex
enddo enddo
enddo enddo
end subroutine mesh_FEM_build_ipCoordinates end subroutine mesh_FEM_build_ipCoordinates