subroutine aeq4(coorr,coefr,prmt,estif,emass,edamp,eload,num)
implicit real*8 (a-h,o-z)
dimension estif(4,4),elump(4),emass(4),
&eload(4)
dimension prmt(*),coef(2),coefr(4,2),
& coorr(2,4),coor(2)
common /raeq4/ru(4,12),
& cu(4,3)
common /vaeq4/rctr(2,2),crtr(2,2),coefd(2,5),coefc(2,5)
common /daeq4/ refc(2,4),gaus(4),
& nnode,ngaus,ndisp,nrefc,ncoor,nvar,
& nvard(1),kdord(1),kvord(4,1)
ea=prmt(1)
eb=prmt(2)
gama=prmt(3)
time=prmt(4)
dt=prmt(5)
imate=prmt(6)+0.5
ielem=prmt(7)+0.5
if (num.eq.1) call aeq4i
do 10 i=1,nvar
emass(i)=0.0
eload(i)=0.0
do 10 j=1,nvar
estif(i,j)=0.0
10 continue
do 999 igaus=1,ngaus
call aeq4t(nnode,nrefc,ncoor,refc(1,igaus),coor,coorr,
& rctr,crtr,det,coefr)
x=coor(1)
y=coor(2)
rx=refc(1,igaus)
ry=refc(2,igaus)
call eaeq4(refc(1,igaus),coef,coorr,coefr,coefd)
iu=(igaus-1)*3+1
if (num.gt.1) goto 2
call aeq41(refc(1,igaus),ru(1,iu),rctr,crtr)
2 continue
call shapn(nrefc,ncoor,4,ru(1,iu),cu,crtr,1,3,3)
call shapc(nrefc,ncoor,2,coefd,coefc,crtr,2,5,5)
un=coef(1)
vn=coef(2)
weigh=det*gaus(igaus)
do 202 i=1,4
iv=kvord(i,1)
do 201 j=1,4
jv=kvord(j,1)
stif=+cu(i,2)*cu(j,2)
& +cu(i,3)*cu(j,3) 
estif(jv,iv)=estif(jv,iv)+stif*weigh
201 continue
202 continue
stif= 1.
elump(1)=stif*weigh
stif= 1.
elump(2)=stif*weigh
stif= 1.
elump(3)=stif*weigh
stif= 1.
elump(4)=stif*weigh
do 301 i=1,nnode
if (nvard(1).lt.i) goto 301
iv = kvord(i,1)
emass(iv)=emass(iv)+elump(iv)*cu(i,1)
301 continue
do 501 i=1,4
iv=kvord(i,1)
stif=+cu(i,1)*rc(un,vn,gama)*ea
eload(iv)=eload(iv)+stif*weigh
501 continue
999 continue
return
end

real*8 function rc(un,vn,gama)
implicit real*8 (a-h,o-z)
rc=(1-vn)*exp(gama*(1-1/(un+273.0)))
return
end

subroutine aeq4i
implicit real*8 (a-h,o-z)
common /daeq4/ refc(2,4),gaus(4),
& nnode,ngaus,ndisp,nrefc,ncoor,nvar,
& nvard(1),kdord(1),kvord(4,1)
ngaus= 4
ndisp= 1
nrefc= 2
ncoor= 2
nvar = 4
nnode= 4
kdord(1)=1
nvard(1)=4
kvord(1,1)=1
kvord(2,1)=2
kvord(3,1)=3
kvord(4,1)=4
refc(1,1)=-1.
refc(2,1)=-1.
gaus(1)=1.
refc(1,2)=1.
refc(2,2)=-1.
gaus(2)=1.
refc(1,3)=1.
refc(2,3)=1.
gaus(3)=1.
refc(1,4)=-1.
refc(2,4)=1.
gaus(4)=1.
end


subroutine aeq4t(nnode,nrefc,ncoor,refc,coor,coorr,rc,cr,det,
& coefr)
implicit real*8 (a-h,o-z)
dimension refc(nrefc),rc(ncoor,nrefc),cr(nrefc,ncoor),a(5,10),
& coorr(ncoor,nnode),coor(ncoor),coefr(nnode,*)
call taeq4(refc,coor,coorr,coefr,rc)
n=nrefc
m=n*2
det = 1.0
do 10 i=1,n
do 10 j=1,n
if (i.le.ncoor) a(i,j) = rc(i,j)
if (i.gt.ncoor) a(i,j)=1.0
a(i,n+j)=0.0
if (i.eq.j) a(i,n+i) = 1.0
10 continue
c write(*,*) 'a ='
c do 21 i=1,n
c21 write(*,8) (a(i,j),j=1,m)
do 400 i=1,n
amax = 0.0
l = 0
do 50 j=i,n
c = a(j,i)
if (c.lt.0.0) c = -c
if (c.le.amax) goto 50
amax = c
l = j
50 continue
do 60 k=1,m
c = a(l,k)
a(l,k) = a(i,k)
a(i,k) = c
60 continue
c = a(i,i)
det = c*det
do 100 k=i+1,m
100 a(i,k) = a(i,k)/c
do 300 j=1,n
if (i.eq.j) goto 300
do 200 k=i+1,m
200 a(j,k) = a(j,k)-a(i,k)*a(j,i)
c write(*,*) 'i =',i,' j =',j,' a ='
c do 11 ii=1,n
c11 write(*,8) (a(ii,jj),jj=1,m)
300 continue
400 continue
do 500 i=1,nrefc
do 500 j=1,ncoor
500 cr(i,j) = a(i,n+j)
c write(*,*) 'a ='
c do 22 i=1,n
c22 write(*,8) (a(i,j),j=1,m)
c write(*,*) 'rc ='
c do 24 i=1,ncoor
c24 write(*,8) (rc(i,j),j=1,nrefc)
c write(*,*) 'cr ='
c do 23 i=1,nrefc
c23 write(*,8) (cr(i,j),j=1,ncoor)
c write(*,*) 'det =',det
if (det.lt.0.0) det=-det
c write(*,*) 'det =',det
8 format(1x,6f12.3)
end

subroutine aeq41(refc,shpr,rctr,crtr)
implicit real*8 (a-h,o-z)
dimension refc(2),shpr(4,3),rctr(2,2),crtr(2,2)
external faeq41
rx=refc(1)
ry=refc(2)
call dshap(faeq41,refc,shpr,2,4,1)
return
end

real*8 function faeq41(refc,n)
implicit real*8 (a-h,o-z)
common /ccaeq4/ xa(4),ya(4),una(4),vna(4)
common /vaeq4/ rctr(2,2),crtr(2,2),coefd(2,5),coefc(2,5)
dimension refc(2)
common /coord/ coor(3),coora(27,3)
x=coor(1)
y=coor(2)
rx=refc(1)
ry=refc(2)
goto (1,2,3,4) n
1 faeq41=+(+1.-rx)/2.*(+1.-ry)/2. 
goto 1000
2 faeq41=+(+1.+rx)/2.*(+1.-ry)/2. 
goto 1000
3 faeq41=+(+1.+rx)/2.*(+1.+ry)/2. 
goto 1000
4 faeq41=+(+1.-rx)/2.*(+1.+ry)/2. 
goto 1000
1000 return
end

subroutine taeq4(refc,coor,coorr,coefr,rc)
implicit real*8 (a-h,o-z)
dimension refc(2),coor(2),coorr(2,4),coefr(4,2),rc(2,2)
common /ccaeq4/ x(4),y(4),un(4),vn(4)
external ftaeq4
do 100 n=1,4
x(n)=coorr(1,n)
y(n)=coorr(2,n)
100 continue
do 200 n=1,4
un(n)=coefr(n,1)
vn(n)=coefr(n,2)
200 continue
rx=refc(1)
ry=refc(2)
call dcoor(ftaeq4,refc,coor,rc,2,2,1)
return
end

real*8 function ftaeq4(refc,n)
implicit real*8 (a-h,o-z)
dimension refc(2)
common /ccaeq4/ x(4),y(4),un(4),vn(4)
common /vaeq4/ rctr(2,2),crtr(2,2),coefd(2,5),coefc(2,5)
rx=refc(1)
ry=refc(2)
goto (1,2) n
1 ftaeq4=+(+(+1.-rx)/2.*(+1.-ry)/2.)*x(1)+(+(+1.+rx)/
& 2.*(+1.-ry)/2.)*x(2)+(+(+1.+rx)/2.*(+1.+ry)/2.)*x(3)+(+
& (+1.-rx)/2.*(+1.+ry)/2.)*x(4)
goto 1000
2 ftaeq4=+(+(+1.-rx)/2.*(+1.-ry)/2.)*y(1)+(+(+1.+rx)/
& 2.*(+1.-ry)/2.)*y(2)+(+(+1.+rx)/2.*(+1.+ry)/2.)*y(3)+(+
& (+1.-rx)/2.*(+1.+ry)/2.)*y(4)
goto 1000
1000 return
end

subroutine eaeq4(refc,coef,coorr,coefr,coefd)
implicit real*8 (a-h,o-z)
dimension refc(2),coef(2),coorr(2,4),coefr(4,2),coefd(2,2)
external feaeq4
rx=refc(1)
ry=refc(2)
call dcoef(feaeq4,refc,coef,coefd,2,2,2)
return
end

real*8 function feaeq4(refc,n)
implicit real*8 (a-h,o-z)
dimension refc(2)
common /ccaeq4/ xa(4),ya(4),un(4),vn(4)
common /vaeq4/ rctr(2,2),crtr(2,2),coefd(2,5),coefc(2,5)
common /coord/ coor(3),coora(27,3)
x=coor(1)
y=coor(2)
rx=refc(1)
ry=refc(2)
goto (1,2) n
2 feaeq4=+(+(+1.-rx)/2.*(+1.-ry)/2.)*vn(1)+(+(+1.+rx)/
& 2.*(+1.-ry)/2.)*vn(2)+(+(+1.+rx)/2.*(+1.+ry)/2.)*vn(3)
& +(+(+1.-rx)/2.*(+1.+ry)/2.)*vn(4)
goto 1000
1 feaeq4=+(+(+1.-rx)/2.*(+1.-ry)/2.)*un(1)+(+(+1.+rx)/
& 2.*(+1.-ry)/2.)*un(2)+(+(+1.+rx)/2.*(+1.+ry)/2.)*un(3)
& +(+(+1.-rx)/2.*(+1.+ry)/2.)*un(4)
goto 1000
1000 return
end

Concentration of chemical constituents

subroutine beq4(coorr,coefr,prmt,estif,emass,edamp,eload,num)
implicit real*8 (a-h,o-z)
dimension estif(4,4),elump(4),emass(4),
&eload(4)
dimension prmt(*),coef(2),coefr(4,2),
& coorr(2,4),coor(2)
common /rbeq4/rv(4,12),
& cv(4,3)
common /vbeq4/rctr(2,2),crtr(2,2),coefd(2,5),coefc(2,5)
common /dbeq4/ refc(2,4),gaus(4),
& nnode,ngaus,ndisp,nrefc,ncoor,nvar,
& nvard(1),kdord(1),kvord(4,1)
ea=prmt(1)
eb=prmt(2)
gama=prmt(3)
time=prmt(4)
dt=prmt(5)
imate=prmt(6)+0.5
ielem=prmt(7)+0.5
if (num.eq.1) call beq4i
do 10 i=1,nvar
emass(i)=0.0
eload(i)=0.0
do 10 j=1,nvar
estif(i,j)=0.0
10 continue
do 999 igaus=1,ngaus
call beq4t(nnode,nrefc,ncoor,refc(1,igaus),coor,coorr,
& rctr,crtr,det,coefr)
x=coor(1)
y=coor(2)
rx=refc(1,igaus)
ry=refc(2,igaus)
call ebeq4(refc(1,igaus),coef,coorr,coefr,coefd)
iv=(igaus-1)*3+1
if (num.gt.1) goto 2
call beq41(refc(1,igaus),rv(1,iv),rctr,crtr)
2 continue
call shapn(nrefc,ncoor,4,rv(1,iv),cv,crtr,1,3,3)
call shapc(nrefc,ncoor,2,coefd,coefc,crtr,2,5,5)
vn=coef(1)
un=coef(2)
weigh=det*gaus(igaus)
do 202 i=1,4
iv=kvord(i,1)
do 201 j=1,4
jv=kvord(j,1)
stif=+cv(i,2)*cv(j,2)
& +cv(i,3)*cv(j,3) 
estif(jv,iv)=estif(jv,iv)+stif*weigh
201 continue
202 continue
stif= 1.
elump(1)=stif*weigh
stif= 1.
elump(2)=stif*weigh
stif= 1.
elump(3)=stif*weigh
stif= 1.
elump(4)=stif*weigh
do 301 i=1,nnode
if (nvard(1).lt.i) goto 301
iv = kvord(i,1)
emass(iv)=emass(iv)+elump(iv)*cv(i,1)
301 continue
do 501 i=1,4
iv=kvord(i,1)
stif=+cv(i,1)*rc(un,vn,gama)*eb
eload(iv)=eload(iv)+stif*weigh
501 continue
999 continue
return
end

real*8 function rc(un,vn,gama)
implicit real*8 (a-h,o-z)
rc=(1-vn)*exp(gama*(1-1/(un+273.0)))
return
end

subroutine beq4i
implicit real*8 (a-h,o-z)
common /dbeq4/ refc(2,4),gaus(4),
& nnode,ngaus,ndisp,nrefc,ncoor,nvar,
& nvard(1),kdord(1),kvord(4,1)
ngaus= 4
ndisp= 1
nrefc= 2
ncoor= 2
nvar = 4
nnode= 4
kdord(1)=1
nvard(1)=4
kvord(1,1)=1
kvord(2,1)=2
kvord(3,1)=3
kvord(4,1)=4
refc(1,1)=-1.
refc(2,1)=-1.
gaus(1)=1.
refc(1,2)=1.
refc(2,2)=-1.
gaus(2)=1.
refc(1,3)=1.
refc(2,3)=1.
gaus(3)=1.
refc(1,4)=-1.
refc(2,4)=1.
gaus(4)=1.
end


subroutine beq4t(nnode,nrefc,ncoor,refc,coor,coorr,rc,cr,det,
& coefr)
implicit real*8 (a-h,o-z)
dimension refc(nrefc),rc(ncoor,nrefc),cr(nrefc,ncoor),a(5,10),
& coorr(ncoor,nnode),coor(ncoor),coefr(nnode,*)
call tbeq4(refc,coor,coorr,coefr,rc)
n=nrefc
m=n*2
det = 1.0
do 10 i=1,n
do 10 j=1,n
if (i.le.ncoor) a(i,j) = rc(i,j)
if (i.gt.ncoor) a(i,j)=1.0
a(i,n+j)=0.0
if (i.eq.j) a(i,n+i) = 1.0
10 continue
c write(*,*) 'a ='
c do 21 i=1,n
c21 write(*,8) (a(i,j),j=1,m)
do 400 i=1,n
amax = 0.0
l = 0
do 50 j=i,n
c = a(j,i)
if (c.lt.0.0) c = -c
if (c.le.amax) goto 50
amax = c
l = j
50 continue
do 60 k=1,m
c = a(l,k)
a(l,k) = a(i,k)
a(i,k) = c
60 continue
c = a(i,i)
det = c*det
do 100 k=i+1,m
100 a(i,k) = a(i,k)/c
do 300 j=1,n
if (i.eq.j) goto 300
do 200 k=i+1,m
200 a(j,k) = a(j,k)-a(i,k)*a(j,i)
c write(*,*) 'i =',i,' j =',j,' a ='
c do 11 ii=1,n
c11 write(*,8) (a(ii,jj),jj=1,m)
300 continue
400 continue
do 500 i=1,nrefc
do 500 j=1,ncoor
500 cr(i,j) = a(i,n+j)
c write(*,*) 'a ='
c do 22 i=1,n
c22 write(*,8) (a(i,j),j=1,m)
c write(*,*) 'rc ='
c do 24 i=1,ncoor
c24 write(*,8) (rc(i,j),j=1,nrefc)
c write(*,*) 'cr ='
c do 23 i=1,nrefc
c23 write(*,8) (cr(i,j),j=1,ncoor)
c write(*,*) 'det =',det
if (det.lt.0.0) det=-det
c write(*,*) 'det =',det
8 format(1x,6f12.3)
end

subroutine beq41(refc,shpr,rctr,crtr)
implicit real*8 (a-h,o-z)
dimension refc(2),shpr(4,3),rctr(2,2),crtr(2,2)
external fbeq41
rx=refc(1)
ry=refc(2)
call dshap(fbeq41,refc,shpr,2,4,1)
return
end

real*8 function fbeq41(refc,n)
implicit real*8 (a-h,o-z)
common /ccbeq4/ xa(4),ya(4),vna(4),una(4)
common /vbeq4/ rctr(2,2),crtr(2,2),coefd(2,5),coefc(2,5)
dimension refc(2)
common /coord/ coor(3),coora(27,3)
x=coor(1)
y=coor(2)
rx=refc(1)
ry=refc(2)
goto (1,2,3,4) n
1 fbeq41=+(+1.-rx)/2.*(+1.-ry)/2. 
goto 1000
2 fbeq41=+(+1.+rx)/2.*(+1.-ry)/2. 
goto 1000
3 fbeq41=+(+1.+rx)/2.*(+1.+ry)/2. 
goto 1000
4 fbeq41=+(+1.-rx)/2.*(+1.+ry)/2. 
goto 1000
1000 return
end

subroutine tbeq4(refc,coor,coorr,coefr,rc)
implicit real*8 (a-h,o-z)
dimension refc(2),coor(2),coorr(2,4),coefr(4,2),rc(2,2)
common /ccbeq4/ x(4),y(4),vn(4),un(4)
external ftbeq4
do 100 n=1,4
x(n)=coorr(1,n)
y(n)=coorr(2,n)
100 continue
do 200 n=1,4
vn(n)=coefr(n,1)
un(n)=coefr(n,2)
200 continue
rx=refc(1)
ry=refc(2)
call dcoor(ftbeq4,refc,coor,rc,2,2,1)
return
end

real*8 function ftbeq4(refc,n)
implicit real*8 (a-h,o-z)
dimension refc(2)
common /ccbeq4/ x(4),y(4),vn(4),un(4)
common /vbeq4/ rctr(2,2),crtr(2,2),coefd(2,5),coefc(2,5)
rx=refc(1)
ry=refc(2)
goto (1,2) n
1 ftbeq4=+(+(+1.-rx)/2.*(+1.-ry)/2.)*x(1)+(+(+1.+rx)/
& 2.*(+1.-ry)/2.)*x(2)+(+(+1.+rx)/2.*(+1.+ry)/2.)*x(3)+(+
& (+1.-rx)/2.*(+1.+ry)/2.)*x(4)
goto 1000
2 ftbeq4=+(+(+1.-rx)/2.*(+1.-ry)/2.)*y(1)+(+(+1.+rx)/
& 2.*(+1.-ry)/2.)*y(2)+(+(+1.+rx)/2.*(+1.+ry)/2.)*y(3)+(+
& (+1.-rx)/2.*(+1.+ry)/2.)*y(4)
goto 1000
1000 return
end

subroutine ebeq4(refc,coef,coorr,coefr,coefd)
implicit real*8 (a-h,o-z)
dimension refc(2),coef(2),coorr(2,4),coefr(4,2),coefd(2,2)
external febeq4
rx=refc(1)
ry=refc(2)
call dcoef(febeq4,refc,coef,coefd,2,2,2)
return
end

real*8 function febeq4(refc,n)
implicit real*8 (a-h,o-z)
dimension refc(2)
common /ccbeq4/ xa(4),ya(4),vn(4),un(4)
common /vbeq4/ rctr(2,2),crtr(2,2),coefd(2,5),coefc(2,5)
common /coord/ coor(3),coora(27,3)
x=coor(1)
y=coor(2)
rx=refc(1)
ry=refc(2)
goto (1,2) n
2 febeq4=+(+(+1.-rx)/2.*(+1.-ry)/2.)*un(1)+(+(+1.+rx)/
& 2.*(+1.-ry)/2.)*un(2)+(+(+1.+rx)/2.*(+1.+ry)/2.)*un(3)
& +(+(+1.-rx)/2.*(+1.+ry)/2.)*un(4)
goto 1000
1 febeq4=+(+(+1.-rx)/2.*(+1.-ry)/2.)*vn(1)+(+(+1.+rx)/
& 2.*(+1.-ry)/2.)*vn(2)+(+(+1.+rx)/2.*(+1.+ry)/2.)*vn(3)
& +(+(+1.-rx)/2.*(+1.+ry)/2.)*vn(4)
goto 1000
1000 return
end

【Close】