polishing
This commit is contained in:
Martin Diehl
parent
74927d9622
commit
56ab4f723d
224
src/mesh_FEM.f90
224
src/mesh_FEM.f90
|
|
@ -66,128 +66,124 @@ contains
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine mesh_init
|
subroutine mesh_init
|
||||||
|
|
||||||
integer, dimension(1), parameter:: FE_geomtype = [1] !< geometry type of particular element type
|
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, parameter :: FILEUNIT = 222
|
||||||
|
integer :: j
|
||||||
|
integer, allocatable, dimension(:) :: chunkPos
|
||||||
|
integer :: dimPlex, &
|
||||||
|
mesh_Nnodes !< total number of nodes in mesh
|
||||||
|
integer, parameter :: &
|
||||||
|
mesh_ElemType=1 !< Element type of the mesh (only support homogeneous meshes)
|
||||||
|
character(len=512) :: &
|
||||||
|
line
|
||||||
|
logical :: flag
|
||||||
|
PetscSF :: sf
|
||||||
|
DM :: globalMesh
|
||||||
|
PetscInt :: face, nFaceSets
|
||||||
|
PetscInt, pointer :: pFaceSets(:)
|
||||||
|
IS :: faceSetIS
|
||||||
|
PetscErrorCode :: ierr
|
||||||
|
|
||||||
integer, dimension(1) :: FE_Nips !< number of IPs in a specific type of element
|
|
||||||
|
write(6,'(/,a)') ' <<<+- mesh init -+>>>'
|
||||||
|
|
||||||
|
! read in file
|
||||||
integer, parameter :: FILEUNIT = 222
|
call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr)
|
||||||
integer :: j
|
CHKERRQ(ierr)
|
||||||
integer, allocatable, dimension(:) :: chunkPos
|
! get spatial dimension (2 or 3?)
|
||||||
integer :: dimPlex, &
|
call DMGetDimension(globalMesh,dimPlex,ierr)
|
||||||
mesh_Nnodes !< total number of nodes in mesh
|
CHKERRQ(ierr)
|
||||||
integer, parameter :: &
|
write(6,*) 'dimension',dimPlex;flush(6)
|
||||||
mesh_ElemType=1 !< Element type of the mesh (only support homogeneous meshes)
|
call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr)
|
||||||
character(len=512) :: &
|
CHKERRQ(ierr)
|
||||||
line
|
! get number of IDs in face sets (for boundary conditions?)
|
||||||
logical :: flag
|
call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr)
|
||||||
PetscSF :: sf
|
CHKERRQ(ierr)
|
||||||
DM :: globalMesh
|
write(6,*) 'number of "Face Sets"',mesh_Nboundaries;flush(6)
|
||||||
PetscInt :: face, nFaceSets
|
call MPI_Bcast(mesh_Nboundaries,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
|
||||||
PetscInt, pointer :: pFaceSets(:)
|
call MPI_Bcast(mesh_NcpElemsGlobal,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
|
||||||
IS :: faceSetIS
|
call MPI_Bcast(dimPlex,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
|
||||||
PetscErrorCode :: ierr
|
|
||||||
|
|
||||||
|
allocate(mesh_boundaries(mesh_Nboundaries), source = 0)
|
||||||
write(6,'(/,a)') ' <<<+- mesh init -+>>>'
|
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)
|
||||||
|
|
||||||
! read in file
|
! this read in function should ignore C and C++ style comments
|
||||||
call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr)
|
! it is used for BC only?
|
||||||
CHKERRQ(ierr)
|
if (worldrank == 0) then
|
||||||
! get spatial dimension (2 or 3?)
|
j = 0
|
||||||
call DMGetDimension(globalMesh,dimPlex,ierr)
|
flag = .false.
|
||||||
CHKERRQ(ierr)
|
call IO_open_file(FILEUNIT,trim(geometryFile))
|
||||||
write(6,*) 'dimension',dimPlex;flush(6)
|
do
|
||||||
call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr)
|
read(FILEUNIT,'(A)') line
|
||||||
CHKERRQ(ierr)
|
if (trim(line) == IO_EOF) exit ! skip empty lines
|
||||||
! get number of IDs in face sets (for boundary conditions?)
|
if (trim(line) == '$Elements') then
|
||||||
call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr)
|
read(FILEUNIT,'(A)') line ! number of elements (ignore)
|
||||||
CHKERRQ(ierr)
|
read(FILEUNIT,'(A)') line
|
||||||
write(6,*) 'number of "Face Sets"',mesh_Nboundaries;flush(6)
|
flag = .true.
|
||||||
call MPI_Bcast(mesh_Nboundaries,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
|
endif
|
||||||
call MPI_Bcast(mesh_NcpElemsGlobal,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
|
if (trim(line) == '$EndElements') exit
|
||||||
call MPI_Bcast(dimPlex,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
|
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)
|
||||||
|
call DMClone(globalMesh,geomMesh,ierr)
|
||||||
|
CHKERRQ(ierr)
|
||||||
|
else
|
||||||
|
call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr)
|
||||||
|
CHKERRQ(ierr)
|
||||||
|
endif
|
||||||
|
|
||||||
allocate(mesh_boundaries(mesh_Nboundaries), source = 0)
|
call DMDestroy(globalMesh,ierr); CHKERRQ(ierr)
|
||||||
call DMGetLabelSize(globalMesh,'Face Sets',nFaceSets,ierr)
|
|
||||||
CHKERRQ(ierr)
|
call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_NcpElems,ierr)
|
||||||
call DMGetLabelIdIS(globalMesh,'Face Sets',faceSetIS,ierr)
|
CHKERRQ(ierr)
|
||||||
CHKERRQ(ierr)
|
call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr)
|
||||||
if (nFaceSets > 0) call ISGetIndicesF90(faceSetIS,pFaceSets,ierr)
|
CHKERRQ(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)
|
|
||||||
|
|
||||||
! this read in function should ignore C and C++ style comments
|
FE_Nips(FE_geomtype(1)) = FEM_Zoo_nQuadrature(dimPlex,integrationOrder)
|
||||||
! it is used for BC only?
|
mesh_maxNips = FE_Nips(1)
|
||||||
if (worldrank == 0) then
|
|
||||||
j = 0
|
write(6,*) 'mesh_maxNips',mesh_maxNips
|
||||||
flag = .false.
|
call mesh_FEM_build_ipCoordinates(dimPlex,FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p)
|
||||||
call IO_open_file(FILEUNIT,trim(geometryFile))
|
call mesh_FEM_build_ipVolumes(dimPlex)
|
||||||
do
|
|
||||||
read(FILEUNIT,'(a512)') line
|
allocate (mesh_element (4,mesh_NcpElems)); mesh_element = 0
|
||||||
if (trim(line) == IO_EOF) exit ! skip empty lines
|
do j = 1, mesh_NcpElems
|
||||||
if (trim(line) == '$Elements') then
|
mesh_element( 1,j) = -1 ! DEPRECATED
|
||||||
read(FILEUNIT,'(a512)') line ! number of elements (ignore)
|
mesh_element( 2,j) = mesh_elemType ! elem type
|
||||||
read(FILEUNIT,'(a512)') line
|
mesh_element( 3,j) = 1 ! homogenization
|
||||||
flag = .true.
|
call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr)
|
||||||
endif
|
CHKERRQ(ierr)
|
||||||
if (trim(line) == '$EndElements') exit
|
end do
|
||||||
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)
|
|
||||||
call DMClone(globalMesh,geomMesh,ierr)
|
|
||||||
CHKERRQ(ierr)
|
|
||||||
else
|
|
||||||
call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr)
|
|
||||||
CHKERRQ(ierr)
|
|
||||||
endif
|
|
||||||
|
|
||||||
call DMDestroy(globalMesh,ierr); CHKERRQ(ierr)
|
if (debug_e < 1 .or. debug_e > mesh_NcpElems) &
|
||||||
|
call IO_error(602,ext_msg='element') ! selected element does not exist
|
||||||
call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_NcpElems,ierr)
|
if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2,debug_e)))) &
|
||||||
CHKERRQ(ierr)
|
call IO_error(602,ext_msg='IP') ! selected element does not have requested IP
|
||||||
call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr)
|
|
||||||
CHKERRQ(ierr)
|
FEsolving_execElem = [ 1,mesh_NcpElems ] ! parallel loop bounds set to comprise all DAMASK elements
|
||||||
|
FEsolving_execIP = spread([1,FE_Nips(FE_geomtype(mesh_element(2,1))),2,nElems)
|
||||||
|
|
||||||
|
allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal)
|
||||||
|
|
||||||
FE_Nips(FE_geomtype(1)) = FEM_Zoo_nQuadrature(dimPlex,integrationOrder)
|
call discretization_init(mesh_element(3,:),mesh_element(4,:),&
|
||||||
mesh_maxNips = FE_Nips(1)
|
|
||||||
|
|
||||||
write(6,*) 'mesh_maxNips',mesh_maxNips
|
|
||||||
call mesh_FEM_build_ipCoordinates(dimPlex,FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p)
|
|
||||||
call mesh_FEM_build_ipVolumes(dimPlex)
|
|
||||||
|
|
||||||
allocate (mesh_element (4,mesh_NcpElems)); mesh_element = 0
|
|
||||||
do j = 1, mesh_NcpElems
|
|
||||||
mesh_element( 1,j) = -1 ! DEPRECATED
|
|
||||||
mesh_element( 2,j) = mesh_elemType ! elem type
|
|
||||||
mesh_element( 3,j) = 1 ! 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,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
|
|
||||||
if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP)
|
|
||||||
allocate(FEsolving_execIP(2,mesh_NcpElems)); FEsolving_execIP = 1 ! parallel loop bounds set to comprise from first IP...
|
|
||||||
forall (j = 1:mesh_NcpElems) FEsolving_execIP(2,j) = FE_Nips(FE_geomtype(mesh_element(2,j))) ! ...up to own IP count for each element
|
|
||||||
|
|
||||||
allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal)
|
|
||||||
|
|
||||||
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)
|
||||||
|
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue