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

added function mesh_FEtoCPelement(FEid)

This commit is contained in:
Philip Eisenlohr 2007-03-21 16:18:33 +00:00
parent 2b05803684
commit b91ab5ae61
1 changed files with 217 additions and 175 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

392
trunk/mesh.f90
View File

@ -1,210 +1,252 @@
!############################################################## !##############################################################
MODULE mesh MODULE mesh
!############################################################## !##############################################################
use prec, only: pRe,pIn use prec, only: pReal,pInt
implicit none implicit none
! --------------------------- ! ---------------------------
! _Nelems : total number of elements in mesh ! _Nelems : total number of elements in mesh
! _Nnodes : total number of nodes in mesh ! _Nnodes : total number of nodes in mesh
! _maxNnodes : max number of nodes in any element ! _maxNnodes : max number of nodes in any element
! _maxNips : max number of IPs in any element ! _maxNips : max number of IPs in any element
! _element : type, material, node indices ! _mapFEtoCPelement : [sorted FEid, corresponding CPid]
! _node : x,y,z coordinates (initially!) ! _element : FEid, type, material, texture, node indices
! _nodeIndex : count of elements containing node, ! _node : x,y,z coordinates (initially!)
! [element_num, node_index], ... ! _nodeIndex : count of elements containing node,
! _envIP : 6 neighboring IPs as [element_num, IP_index] ! [element_num, node_index], ...
! order is +x, +y,+z, -x, -y, -z in local coord ! _envIP : 6 neighboring IPs as [element_num, IP_index]
! --------------------------- ! order is +x, +y,+z, -x, -y, -z in local coord
integer(pIn) mesh_Nelems, mesh_Nnodes, mesh_maxNnodes,mesh_maxNips ! ---------------------------
integer(pIn), allocatable :: mesh_element (:,:) integer(pInt) mesh_Nelems, mesh_Nnodes, mesh_maxNnodes,mesh_maxNips
integer(pIn), allocatable :: mesh_nodeIndex (:,:) integer(pInt), dimension(2,:), allocatable :: mesh_mapFEtoCPelement
integer(pIn), allocatable :: mesh_envIP (:,:) integer(pInt), dimension(:,:), allocatable :: mesh_element, mesh_nodeIndex, mesh_envIP
real(pRe), allocatable :: mesh_node (:,:) real(pReal), allocatable :: mesh_node (:,:)
CONTAINS CONTAINS
! --------------------------- ! ---------------------------
! subroutine mesh_init() ! subroutine mesh_init()
! subroutine mesh_parse_inputFile() ! function mesh_FEtoCPelement(FEid)
! --------------------------- ! function mesh_build_IPenvironment()
! subroutine mesh_parse_inputFile()
! ---------------------------
! *********************************************************** ! ***********************************************************
! initialization ! FE to CP id mapping by binary search thru lookup array
! *********************************************************** ! ***********************************************************
SUBROUTINE mesh_init () FUNCTION mesh_FEtoCPelement(FEid)
use prec, only: pInt
implicit none
integer(pInt), intent(in) :: FEid
integer(pInt) mesh_FEtoCPelement, lower,upper,center
mesh_FEtoCPelement = 0_pInt
lower = 1_pInt
upper = size(mesh_mapFEtoCPelement,2)
if (mesh_mapFEtoCPelement(lower,1) == FEid) then
mesh_FEtoCPelement = mesh_mapFEtoCPelement(lower,2)
return
elseif (mesh_mapFEtoCPelement(upper,1) == FEid) then
mesh_FEtoCPelement = mesh_mapFEtoCPelement(upper,2)
return
endif
do while (upper-lower > 0)
center = (lower+upper)/2
if (mesh_mapFEtoCPelement(center,1) < FEid) then
lower = center
elseif (mesh_mapFEtoCPelement(center,1) > FEid) then
upper = center
else
mesh_FEtoCPelement = mesh_mapFEtoCPelement(center,2)
exit
end if
end do
return
END FUNCTION
mesh_Nelems = 0_pIn ! ***********************************************************
mesh_Nnodes = 0_pIn ! initialization
mesh_maxNips = 0_pIn ! ***********************************************************
mesh_maxNnodes = 0_pIn SUBROUTINE mesh_init ()
call mesh_parse_inputFile ()
END SUBROUTINE mesh_Nelems = 0_pInt
mesh_Nnodes = 0_pInt
! *********************************************************** mesh_maxNips = 0_pInt
! parsing of input file mesh_maxNnodes = 0_pInt
! *********************************************************** call mesh_parse_inputFile ()
FUNCTION mesh_parse_inputFile ()
use prec, only: pRe,pIn END SUBROUTINE
use IO
implicit none ! ***********************************************************
! parsing of input file
! ***********************************************************
FUNCTION mesh_parse_inputFile ()
logical mesh_parse_inputFile use prec, only: pReal,pInt
integer(pIn) i,j,positions(10*2+1) use IO
integer(pIn) elem_num,elem_type,Nnodes,node_num,num_ip,mat,tp(70,2)
implicit none
! Set a format to read the entire line (max. len is 80 characters) logical mesh_parse_inputFile
integer(pInt) i,j,positions(10*2+1)
integer(pInt) elem_num,elem_type,Nnodes,node_num,num_ip,mat,tp(70,2)
! Set a format to read the entire line (max. len is 80 characters)
610 FORMAT(A80) 610 FORMAT(A80)
if (.not. IO_open_inputFile(600)) then if (.not. IO_open_inputFile(600)) then
mesh_parse_inputFile = .false. mesh_parse_inputFile = .false.
return return
endif endif
do while(.true.) do while(.true.)
read(600,610,end=620) line read(600,610,end=620) line
positions = IO_stringPos(line,3) positions = IO_stringPos(line,3)
select case (IO_stringValue(line,positions,1) select case (IO_stringValue(line,positions,1)
!----------------------------------- !-----------------------------------
case ('sizing') case ('sizing')
!----------------------------------- !-----------------------------------
mesh_Nelems = IO_intValue(line,positions,2) mesh_Nelems = IO_intValue(line,positions,2)
mesh_Nnodes = IO_intValue(line,positions,3) mesh_Nnodes = IO_intValue(line,positions,3)
!----------------------------------- !-----------------------------------
case ('elements') case ('elements')
!----------------------------------- !-----------------------------------
select case (IO_intValue(line,positions,2)) ! elem type select case (IO_intValue(line,positions,2)) ! elem type
case (3) ! 2D Triangle case (3) ! 2D Triangle
mesh_maxNips = max(3,mesh_maxNips) mesh_maxNips = max(3,mesh_maxNips)
mesh_maxNnodes = max(3,mesh_maxNnodes) mesh_maxNnodes = max(3,mesh_maxNnodes)
case (6) ! 2D Quad. case (6) ! 2D Quad.
mesh_maxNips = max(4,mesh_maxNips) mesh_maxNips = max(4,mesh_maxNips)
mesh_maxNnodes = max(4,mesh_maxNnodes) mesh_maxNnodes = max(4,mesh_maxNnodes)
case (7) ! 3D hexahedral case (7) ! 3D hexahedral
mesh_maxNips = max(8,mesh_maxNips) mesh_maxNips = max(8,mesh_maxNips)
mesh_maxNnodes = max(8,mesh_maxNnodes) mesh_maxNnodes = max(8,mesh_maxNnodes)
case default case default
mesh_maxNips = max(8,mesh_maxNips) mesh_maxNips = max(8,mesh_maxNips)
mesh_maxNnodes = max(8,mesh_maxNnodes) mesh_maxNnodes = max(8,mesh_maxNnodes)
end select end select
!----------------------------------- !-----------------------------------
case ('connectivity') case ('connectivity')
!----------------------------------- !-----------------------------------
allocate (mesh_element(mesh_Nelems,2+mesh_maxNips)) allocate (mesh_element(mesh_Nelems,2+mesh_maxNips))
allocate (mesh_nodeIndex (mesh_Nnodes,1+mesh_maxNnodes*2) allocate (mesh_nodeIndex (mesh_Nnodes,1+mesh_maxNnodes*2)
allocate (mesh_envIP(mesh_Nelems,mesh_maxNips,6,2)) allocate (mesh_envIP(mesh_Nelems,mesh_maxNips,6,2))
mesh_element = 0_pIn mesh_element = 0_pInt
mesh_nodeIndex = 0_pIn mesh_nodeIndex = 0_pInt
mesh_envIP = 0_pIn mesh_envIP = 0_pInt
! MISSING: setting up of envIP ! MISSING: setting up of envIP
read(600,610,end=620) line ! skip line ?? read(600,610,end=620) line ! skip line ??
do i=1,mesh_Nelems do i=1,mesh_Nelems
read(600,610,end=620) line read(600,610,end=620) line
positions = IO_stringPos(line,0) ! find all chunks positions = IO_stringPos(line,0) ! find all chunks
elem_num = IO_intValue(line,positions,1) elem_num = IO_intValue(line,positions,1)
elem_type = IO_intValue(line,positions,2) elem_type = IO_intValue(line,positions,2)
select case (elem_type) select case (elem_type)
case (3) ! 2D Triangle case (3) ! 2D Triangle
Nnodes = 3 Nnodes = 3
case (6) ! 2D Quad. case (6) ! 2D Quad.
Nnodes = 4 Nnodes = 4
case (7) ! 3D hexahedral case (7) ! 3D hexahedral
Nnodes = 8 Nnodes = 8
case default case default
Nnodes = 8 Nnodes = 8
end select end select
mesh_element(elem_num,1) = elem_type mesh_element(elem_num,1) = elem_type
do j=1,Nnodes ! store all node indices do j=1,Nnodes ! store all node indices
node_num = IO_intValue(line,positions,2+j) node_num = IO_intValue(line,positions,2+j)
mesh_element(elem_num,1+j) = node_num mesh_element(elem_num,1+j) = node_num
mesh_nodeIndex(node_num,1) = mesh_nodeIndex(node_num,1)+1 ! inc count mesh_nodeIndex(node_num,1) = mesh_nodeIndex(node_num,1)+1 ! inc count
mesh_nodeIndex(node_num,mesh_nodeIndex(node_num,1)*2 ) = elem_num mesh_nodeIndex(node_num,mesh_nodeIndex(node_num,1)*2 ) = elem_num
mesh_nodeIndex(node_num,mesh_nodeIndex(node_num,1)*2+1) = j mesh_nodeIndex(node_num,mesh_nodeIndex(node_num,1)*2+1) = j
end do end do
end do end do
!----------------------------------- !-----------------------------------
case ('coordinates') case ('coordinates')
!----------------------------------- !-----------------------------------
allocate (mesh_node(mesh_Nnodes,3)) ! x,y,z, per node allocate (mesh_node(mesh_Nnodes,3)) ! x,y,z, per node
read(600,610,end=620) line ! skip line ?? read(600,610,end=620) line ! skip line ??
do i=1,mesh_Nnodes do i=1,mesh_Nnodes
read(600,610,end=620) line read(600,610,end=620) line
positions = IO_stringPos(line,0) ! find all (4) chunks positions = IO_stringPos(line,0) ! find all (4) chunks
node_num = IO_intValue(line,positions,1) node_num = IO_intValue(line,positions,1)
do j=1,3 ! store x,y,z coordinates do j=1,3 ! store x,y,z coordinates
mesh_node(node_num,j) = IO_floatValue(line,positions,1+j) mesh_node(node_num,j) = IO_floatValue(line,positions,1+j)
end do end do
end do end do
!----------------------------------- !-----------------------------------
case ('hypoelastic') case ('hypoelastic')
!----------------------------------- !-----------------------------------
! *************************************************** ! ***************************************************
! Search for key word "hypoelastic". ! Search for key word "hypoelastic".
! This section contains the # of materials and ! This section contains the # of materials and
! the element range of each material ! the element range of each material
! *************************************************** ! ***************************************************
ELSE IF( line(1:11).eq.'hypoelastic' )THEN ELSE IF( line(1:11).eq.'hypoelastic' )THEN
mat=0 mat=0
flag=0 flag=0
DO WHILE( line(1:8).ne.'geometry' ) DO WHILE( line(1:8).ne.'geometry' )
READ(600,610,END=620) line READ(600,610,END=620) line
i=1 i=1
DO WHILE( i.le.len(line)-8 ) DO WHILE( i.le.len(line)-8 )
IF( line(i:i+2).eq.'mat' )THEN IF( line(i:i+2).eq.'mat' )THEN
mat=mat+1 mat=mat+1
flag=1 flag=1
END IF END IF
i=i+1 i=i+1
END DO END DO
IF( flag.eq.1 )THEN IF( flag.eq.1 )THEN
flag=0 flag=0
READ(600,610,END=620) line READ(600,610,END=620) line
READ(600,610,END=620) line READ(600,610,END=620) line
i=1 i=1
DO WHILE( line(i:i).eq.' ') DO WHILE( line(i:i).eq.' ')
i=i+1 i=i+1
END DO END DO
left=i left=i
DO WHILE( line(i:i).ne.' ') DO WHILE( line(i:i).ne.' ')
i=i+1 i=i+1
END DO END DO
right=i-1 right=i-1
READ(UNIT=line(left:right), FMT=' (I5) ') tp(mat,1) READ(UNIT=line(left:right), FMT=' (I5) ') tp(mat,1)
DO WHILE( (line(i:i).eq.' ').or. DO WHILE( (line(i:i).eq.' ').or.
& (line(i:i).eq.'t').or. & (line(i:i).eq.'t').or.
& (line(i:i).eq.'o') ) & (line(i:i).eq.'o') )
i=i+1 i=i+1
END DO END DO
left=i left=i
DO WHILE( line(i:i).ne.' ') DO WHILE( line(i:i).ne.' ')
i=i+1 i=i+1
END DO END DO
right=i-1 right=i-1
READ(UNIT=line(left:right), FMT=' (I5) ') tp(mat,2) READ(UNIT=line(left:right), FMT=' (I5) ') tp(mat,2)
END IF END IF
END DO END DO
WRITE(6,*) 'mat: ',mat,' ',tp(1,1),' ',tp(1,2) WRITE(6,*) 'mat: ',mat,' ',tp(1,1),' ',tp(1,2)
end select end select
END DO END DO
! Code jumps to 620 when it reaches the end of the file ! Code jumps to 620 when it reaches the end of the file
620 continue 620 continue
WRITE(6,*) 'Finished with .dat file' WRITE(6,*) 'Finished with .dat file'
CALL FLUSH(6) CALL FLUSH(6)
END FUNCTION END FUNCTION
END MODULE mesh END MODULE mesh