You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
246 lines
7.8 KiB
Fortran
246 lines
7.8 KiB
Fortran
c----------------------------------------------------------------blocks
|
|
subroutine blocks(nx,ny,nz,np,ne,nop,cord,nxyz,egrp,nbeg,nend,
|
|
* nelg,nbb,nee,nelb,nblk,egps,itype,ebox,m,num)
|
|
c---------------------------------------------------------------------c
|
|
c Purpose: c
|
|
c To divide FE mesh into blocks. c
|
|
c---------------------------------------------------------------------c
|
|
common /etype/inode(4)
|
|
integer*2 nop
|
|
dimension nop(20,1),cord(3,1),egrp(6,1),nbeg(1),nend(1),
|
|
* nelg(1),nee(1),nbb(1),nelb(1),nblk(2,1),egps(6,1),
|
|
* itype(1),ebox(6,1)
|
|
c
|
|
xmax=cord(1,1)
|
|
xmin=cord(1,1)
|
|
ymax=cord(2,1)
|
|
ymin=cord(2,1)
|
|
zmax=cord(3,1)
|
|
zmin=cord(3,1)
|
|
disd=0.0
|
|
do i=1,ne
|
|
if(nop(9,i).eq.0) goto 10
|
|
c
|
|
c----------store element type
|
|
c
|
|
if(nop(20,i).ne.0) then
|
|
itype(i)=3
|
|
else if(nop(15,i).ne.0) then
|
|
itype(i)=2
|
|
else if(nop(13,i).ne.0) then
|
|
itype(i)=4
|
|
else
|
|
itype(i)=1
|
|
endif
|
|
c
|
|
c----------find B.L. and T.R. conners of each element
|
|
c
|
|
ebox(1,i)=cord(1,nop(1,i))
|
|
ebox(4,i)=cord(1,nop(1,i))
|
|
ebox(2,i)=cord(2,nop(1,i))
|
|
ebox(5,i)=cord(2,nop(1,i))
|
|
ebox(3,i)=cord(3,nop(1,i))
|
|
ebox(6,i)=cord(3,nop(1,i))
|
|
do k=2,inode(itype(i))
|
|
ebox(1,i)=min(ebox(1,i),cord(1,nop(k,i)))
|
|
ebox(4,i)=max(ebox(4,i),cord(1,nop(k,i)))
|
|
ebox(2,i)=min(ebox(2,i),cord(2,nop(k,i)))
|
|
ebox(5,i)=max(ebox(5,i),cord(2,nop(k,i)))
|
|
ebox(3,i)=min(ebox(3,i),cord(3,nop(k,i)))
|
|
ebox(6,i)=max(ebox(6,i),cord(3,nop(k,i)))
|
|
enddo
|
|
ebox(1,i)=ebox(1,i)-1.0e-5*abs(ebox(1,i))
|
|
ebox(4,i)=ebox(4,i)+1.0e-5*abs(ebox(4,i))
|
|
ebox(2,i)=ebox(2,i)-1.0e-5*abs(ebox(2,i))
|
|
ebox(5,i)=ebox(5,i)+1.0e-5*abs(ebox(5,i))
|
|
ebox(3,i)=ebox(3,i)-1.0e-5*abs(ebox(3,i))
|
|
ebox(6,i)=ebox(6,i)+1.0e-5*abs(ebox(6,i))
|
|
c
|
|
c----------calculate the maximum diagonal distance
|
|
c
|
|
disd0=sqrt((ebox(4,i)-ebox(1,i))*(ebox(4,i)-ebox(1,i))+
|
|
* (ebox(5,i)-ebox(2,i))*(ebox(5,i)-ebox(2,i)))
|
|
disd=max(disd,disd0)
|
|
c
|
|
c-------find B.L. and T.R. conners of the FE mesh
|
|
c
|
|
xmax=max(xmax,ebox(4,i))
|
|
xmin=min(xmin,ebox(1,i))
|
|
ymax=max(ymax,ebox(5,i))
|
|
ymin=min(ymin,ebox(2,i))
|
|
zmax=max(zmax,ebox(6,i))
|
|
zmin=min(zmin,ebox(3,i))
|
|
10 enddo
|
|
|
|
m=1
|
|
dx=(xmax-xmin)/float(nx)
|
|
if(dx.lt.disd) then
|
|
nx=(xmax-xmin)/disd
|
|
if(nx.lt.1) nx=1
|
|
dx=(xmax-xmin)/float(nx)
|
|
endif
|
|
dy=(ymax-ymin)/float(ny)
|
|
if(dy.lt.disd) then
|
|
ny=(ymax-ymin)/disd
|
|
if(ny.lt.1) ny=1
|
|
dy=(ymax-ymin)/float(ny)
|
|
endif
|
|
dz=(zmax-zmin)/float(nz)
|
|
x=xmin
|
|
c
|
|
c-------divide into blocks
|
|
c
|
|
do i=1,nx
|
|
y=ymin
|
|
do j=1,ny
|
|
z=zmin
|
|
do k=1,nz
|
|
egrp(1,m)=x-dx*0.01
|
|
egrp(2,m)=y-dy*0.01
|
|
egrp(3,m)=z-dz*0.01
|
|
egrp(4,m)=x+dx*1.01
|
|
egrp(5,m)=y+dy*1.01
|
|
egrp(6,m)=z+dz*1.01
|
|
z=z+dz
|
|
m=m+1
|
|
enddo
|
|
y=y+dy
|
|
enddo
|
|
x=x+dx
|
|
enddo
|
|
nxyz=nx*ny*nz
|
|
nxyz0=nxyz
|
|
num=1
|
|
do i=1,nxyz
|
|
nbeg(i)=num
|
|
do j=1,ne
|
|
if(nop(9,j).ne.0) then
|
|
do k=1,inode(itype(j))
|
|
if(cord(1,nop(k,j)).gt.egrp(4,i).or.
|
|
* cord(1,nop(k,j)).lt.egrp(1,i)) goto 30
|
|
if(cord(2,nop(k,j)).gt.egrp(5,i).or.
|
|
* cord(2,nop(k,j)).lt.egrp(2,i)) goto 30
|
|
if(cord(3,nop(k,j)).gt.egrp(6,i).or.
|
|
* cord(3,nop(k,j)).lt.egrp(3,i)) goto 30
|
|
nelg(num)=j
|
|
num=num+1
|
|
goto 40
|
|
30 enddo
|
|
endif
|
|
40 enddo
|
|
nend(i)=num-1
|
|
enddo
|
|
i=1
|
|
45 if(nend(i)-nbeg(i).lt.0) then
|
|
nxyz=nxyz-1
|
|
do j=i,nxyz
|
|
nbeg(j)=nbeg(j+1)
|
|
nend(j)=nend(j+1)
|
|
do k=1,6
|
|
egrp(k,j)=egrp(k,j+1)
|
|
enddo
|
|
enddo
|
|
else
|
|
i=i+1
|
|
endif
|
|
if(i.le.nxyz) goto 45
|
|
c
|
|
c-------divide into sub-blocks
|
|
c
|
|
m=1
|
|
num=1
|
|
do n=1,nxyz
|
|
num_elem=nend(n)-nbeg(n)+1
|
|
if(num_elem.gt.nxyz+nxyz) then
|
|
a_multiplier=sqrt(float(num_elem))/float(nx*ny)
|
|
nxs=int(nx*a_multiplier)
|
|
if(nxs.lt.1) nxs=1
|
|
nys=int(ny*a_multiplier)
|
|
if(nys.lt.1) nys=1
|
|
dx=(egrp(4,n)-egrp(1,n))/float(nxs)
|
|
dy=(egrp(5,n)-egrp(2,n))/float(nys)
|
|
if(dx.lt.disd) then
|
|
nxs=(egrp(4,n)-egrp(1,n))/disd
|
|
if(nxs.lt.1) nxs=1
|
|
dx=(egrp(4,n)-egrp(1,n))/float(nxs)
|
|
endif
|
|
if(dy.lt.disd) then
|
|
dy=disd
|
|
nys=(egrp(5,n)-egrp(2,n))/disd
|
|
if(nys.lt.1) nys=1
|
|
dy=(egrp(5,n)-egrp(2,n))/float(nys)
|
|
endif
|
|
nxyzs=nxs*nys
|
|
if(nxyzs.le.1) goto 70
|
|
dz=(egrp(6,n)-egrp(3,n))
|
|
|
|
nblk(1,n)=nxyzs
|
|
nblk(2,n)=m
|
|
x=egrp(1,n)
|
|
do i=1,nxs
|
|
y=egrp(2,n)
|
|
do j=1,nys
|
|
z=egrp(3,n)
|
|
egps(1,m)=x-0.01*dx
|
|
egps(2,m)=y-0.01*dy
|
|
egps(3,m)=z-0.01*dz
|
|
egps(4,m)=x+dx*1.01
|
|
egps(5,m)=y+dy*1.01
|
|
egps(6,m)=z+dz*1.01
|
|
y=y+dy
|
|
m=m+1
|
|
enddo
|
|
x=x+dx
|
|
enddo
|
|
do i=nblk(2,n),nblk(2,n)+nblk(1,n)-1
|
|
nbb(i)=num
|
|
do j=nbeg(n),nend(n)
|
|
if(nop(9,nelg(j)).ne.0) then
|
|
do k=1,inode(itype(nelg(j)))
|
|
if(cord(1,nop(k,nelg(j))).gt.egps(4,i).or.
|
|
* cord(1,nop(k,nelg(j))).lt.egps(1,i)) goto 50
|
|
if(cord(2,nop(k,nelg(j))).gt.egps(5,i).or.
|
|
* cord(2,nop(k,nelg(j))).lt.egps(2,i)) goto 50
|
|
if(cord(3,nop(k,nelg(j))).gt.egps(6,i).or.
|
|
* cord(3,nop(k,nelg(j))).lt.egps(3,i)) goto 50
|
|
nelb(num)=nelg(j)
|
|
num=num+1
|
|
goto 60
|
|
50 enddo
|
|
endif
|
|
60 enddo
|
|
nee(i)=num-1
|
|
enddo
|
|
i=nblk(2,n)
|
|
65 if(nee(i)-nbb(i).lt.0) then
|
|
nblk(1,n)=nblk(1,n)-1
|
|
m=m-1
|
|
do j=i,nblk(2,n)+nblk(1,n)-1
|
|
nbb(j)=nbb(j+1)
|
|
nee(j)=nee(j+1)
|
|
do k=1,6
|
|
egps(k,j)=egps(k,j+1)
|
|
enddo
|
|
enddo
|
|
else
|
|
i=i+1
|
|
endif
|
|
if(i.le.nblk(2,n)+nblk(1,n)-1) goto 65
|
|
else
|
|
70 nblk(1,n)=1
|
|
nblk(2,n)=m
|
|
nbb(m)=num
|
|
do i=nbeg(n),nend(n)
|
|
nelb(num)=nelg(i)
|
|
num=num+1
|
|
enddo
|
|
do j=1,6
|
|
egps(j,m)=egrp(j,n)
|
|
enddo
|
|
nee(m)=num-1
|
|
m=m+1
|
|
endif
|
|
enddo
|
|
return
|
|
end
|