FEPG






Displacement field subroutine aeq4(coorr,coefr,prmt,estif,emass,edamp,eload,num) implicit real8 (a-h,o-z) dimension estif(8,8),elump(8),emass(8), &eload(8) dimension prmt(), & eex(8),eey(8),eexy(8), & coorr(2,4),coor(2) common /raeq4/ru(4,12),rv(4,12), & cu(4,3),cv(4,3) common /vaeq4/rctr(2,2),crtr(2,2) common /daeq4/ refc(2,4),gaus(4), & nnode,ngaus,ndisp,nrefc,ncoor,nvar, & nvard(2),kdord(2),kvord(8,2) pe=prmt(1) pv=prmt(2) fx=prmt(3) fy=prmt(4) time=prmt(5) dt=prmt(6) imate=prmt(7)+0.5 ielem=prmt(8)+0.5 fact = pe/(1.+pv)/(1.-pv) if (num.eq.1) call aeq4i do 10 i=1,nvar 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) x=coor(1) y=coor(2) rx=refc(1,igaus) ry=refc(2,igaus) iu=(igaus-1)3+1 iv=(igaus-1)3+1 if (num.gt.1) goto 2 call aeq41(refc(1,igaus),ru(1,iu),rctr,crtr) call aeq42(refc(1,igaus),rv(1,iv),rctr,crtr) 2 continue call shapn(nrefc,ncoor,4,ru(1,iu),cu,crtr,1,3,3) call shapn(nrefc,ncoor,4,rv(1,iv),cv,crtr,1,3,3) weigh=detgaus(igaus) do 100 i=1, 8 eex(i) = 0.0 eey(i) = 0.0 eexy(i) = 0.0 100 continue do 101 i=1,4 iv=kvord(i,1) stif=+cu(i,2) eex(iv)=eex(iv)+stif 101 continue do 102 i=1,4 iv=kvord(i,2) stif=+cv(i,3) eey(iv)=eey(iv)+stif 102 continue do 103 i=1,4 iv=kvord(i,1) stif=+cu(i,3) eexy(iv)=eexy(iv)+stif 103 continue do 104 i=1,4 iv=kvord(i,2) stif=+cv(i,2) eexy(iv)=eexy(iv)+stif 104 continue do 202 iv=1,8 do 201 jv=1,8 stif=+eex(iv)eex(jv)1fact & +eex(iv)eey(jv)pvfact & +eey(iv)eex(jv)pvfact & +eey(iv)eey(jv)1fact & +eexy(iv)eexy(jv)fact((1.-pv)/2) estif(iv,jv)=estif(iv,jv)+stifweigh 201 continue 202 continue do 501 i=1,4 iv=kvord(i,1) stif=+cu(i,1)fx eload(iv)=eload(iv)+stifweigh 501 continue do 502 i=1,4 iv=kvord(i,2) stif=+cv(i,1)fy eload(iv)=eload(iv)+stifweigh 502 continue 999 continue return end subroutine aeq4i implicit real8 (a-h,o-z) common /daeq4/ refc(2,4),gaus(4), & nnode,ngaus,ndisp,nrefc,ncoor,nvar, & nvard(2),kdord(2),kvord(8,2) ngaus= 4 ndisp= 2 nrefc= 2 ncoor= 2 nvar = 8 nnode= 4 kdord(1)=1 nvard(1)=4 kvord(1,1)=1 kvord(2,1)=3 kvord(3,1)=5 kvord(4,1)=7 kdord(2)=1 nvard(2)=4 kvord(1,2)=2 kvord(2,2)=4 kvord(3,2)=6 kvord(4,2)=8 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) implicit real8 (a-h,o-z) dimension refc(nrefc),rc(ncoor,nrefc),cr(nrefc,ncoor),a(5,10), coorr(ncoor,nnode),coor(ncoor) call taeq4(refc,coor,coorr,rc) n=nrefc m=n2 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 = cdet 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 real8 (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 real8 function faeq41(refc,n) implicit real8 (a-h,o-z) common /ccaeq4/ xa(4),ya(4) common /vaeq4/ rctr(2,2),crtr(2,2) 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 aeq42(refc,shpr,rctr,crtr) implicit real8 (a-h,o-z) dimension refc(2),shpr(4,3),rctr(2,2),crtr(2,2) external faeq42 rx=refc(1) ry=refc(2) call dshap(faeq42,refc,shpr,2,4,1) return end real8 function faeq42(refc,n) implicit real8 (a-h,o-z) common /ccaeq4/ xa(4),ya(4) common /vaeq4/ rctr(2,2),crtr(2,2) 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 faeq42=+(+1.-rx)/2.(+1.-ry)/2. goto 1000 2 faeq42=+(+1.+rx)/2.(+1.-ry)/2. goto 1000 3 faeq42=+(+1.+rx)/2.(+1.+ry)/2. goto 1000 4 faeq42=+(+1.-rx)/2.(+1.+ry)/2. goto 1000 1000 return end subroutine taeq4(refc,coor,coorr,rc) implicit real8 (a-h,o-z) dimension refc(2),coor(2),coorr(2,4),rc(2,2) common /ccaeq4/ x(4),y(4) external ftaeq4 do 100 n=1,4 x(n)=coorr(1,n) y(n)=coorr(2,n) 100 continue rx=refc(1) ry=refc(2) call dcoor(ftaeq4,refc,coor,rc,2,2,1) return end real8 function ftaeq4(refc,n) implicit real8 (a-h,o-z) dimension refc(2) common /ccaeq4/ x(4),y(4) common /vaeq4/ rctr(2,2),crtr(2,2) 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 Stress field subroutine beq4(coorr,coefr,prmt,estif,emass,edamp,eload,num) implicit real8 (a-h,o-z) dimension estif(12,12),elump(12),emass(12), &eload(12) dimension prmt(),coef(2),coefr(4,2), & coorr(2,4),coor(2) common /rbeq4/rsx(4,12),rsy(4,12),rsxy(4,12), & csx(4,3),csy(4,3),csxy(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(3),kdord(3),kvord(12,3) pe=prmt(1) pv=prmt(2) fx=prmt(3) fy=prmt(4) time=prmt(5) dt=prmt(6) imate=prmt(7)+0.5 ielem=prmt(8)+0.5 fact = pe/(1.+pv)/(1.-pv) 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) isx=(igaus-1)3+1 isy=(igaus-1)3+1 isxy=(igaus-1)3+1 if (num.gt.1) goto 2 call beq41(refc(1,igaus),rsx(1,isx),rctr,crtr) call beq42(refc(1,igaus),rsy(1,isy),rctr,crtr) call beq43(refc(1,igaus),rsxy(1,isxy),rctr,crtr) 2 continue call shapn(nrefc,ncoor,4,rsx(1,isx),csx,crtr,1,3,3) call shapn(nrefc,ncoor,4,rsy(1,isy),csy,crtr,1,3,3) call shapn(nrefc,ncoor,4,rsxy(1,isxy),csxy,crtr,1,3,3) call shapc(nrefc,ncoor,2,coefd,coefc,crtr,2,5,5) u=coef(1) v=coef(2) weigh=detgaus(igaus) ex=+coefc(1,1) ey=+coefc(2,2) fsx = +1exfact+pveyfact fsy = +pvexfact+1eyfact exy=+coefc(1,2)+coefc(2,1) fsxy=exy((1.-pv)/2)fact do 202 i=1,4 iv=kvord(i,1) do 201 j=1,4 jv=kvord(j,1) stif=+csx(i,1)csx(j,1)0.0 estif(jv,iv)=estif(jv,iv)+stifweigh 201 continue 202 continue stif= 1. elump(1)=stifweigh stif= 1. elump(4)=stifweigh stif= 1. elump(7)=stifweigh stif= 1. elump(10)=stifweigh stif= 1. elump(2)=stifweigh stif= 1. elump(5)=stifweigh stif= 1. elump(8)=stifweigh stif= 1. elump(11)=stifweigh stif= 1. elump(3)=stifweigh stif= 1. elump(6)=stifweigh stif= 1. elump(9)=stifweigh stif= 1. elump(12)=stifweigh do 301 i=1,nnode if (nvard(1).lt.i) goto 301 iv = kvord(i,1) emass(iv)=emass(iv)+elump(iv)csx(i,1) if (nvard(2).lt.i) goto 301 iv = kvord(i,2) emass(iv)=emass(iv)+elump(iv)csy(i,1) if (nvard(3).lt.i) goto 301 iv = kvord(i,3) emass(iv)=emass(iv)+elump(iv)csxy(i,1) 301 continue do 501 i=1,4 iv=kvord(i,1) stif=+csx(i,1)fsx eload(iv)=eload(iv)+stifweigh 501 continue do 502 i=1,4 iv=kvord(i,2) stif=+csy(i,1)fsy eload(iv)=eload(iv)+stifweigh 502 continue do 503 i=1,4 iv=kvord(i,3) stif=+csxy(i,1)fsxy eload(iv)=eload(iv)+stifweigh 503 continue 999 continue return end subroutine beq4i implicit real8 (a-h,o-z) common /dbeq4/ refc(2,4),gaus(4), & nnode,ngaus,ndisp,nrefc,ncoor,nvar, & nvard(3),kdord(3),kvord(12,3) ngaus= 4 ndisp= 3 nrefc= 2 ncoor= 2 nvar = 12 nnode= 4 kdord(1)=1 nvard(1)=4 kvord(1,1)=1 kvord(2,1)=4 kvord(3,1)=7 kvord(4,1)=10 kdord(2)=1 nvard(2)=4 kvord(1,2)=2 kvord(2,2)=5 kvord(3,2)=8 kvord(4,2)=11 kdord(3)=1 nvard(3)=4 kvord(1,3)=3 kvord(2,3)=6 kvord(3,3)=9 kvord(4,3)=12 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 real8 (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=n2 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 = cdet 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 real8 (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 real8 function fbeq41(refc,n) implicit real8 (a-h,o-z) common /ccbeq4/ xa(4),ya(4),ua(4),va(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 beq42(refc,shpr,rctr,crtr) implicit real8 (a-h,o-z) dimension refc(2),shpr(4,3),rctr(2,2),crtr(2,2) external fbeq42 rx=refc(1) ry=refc(2) call dshap(fbeq42,refc,shpr,2,4,1) return end real8 function fbeq42(refc,n) implicit real8 (a-h,o-z) common /ccbeq4/ xa(4),ya(4),ua(4),va(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 fbeq42=+(+1.-rx)/2.(+1.-ry)/2. goto 1000 2 fbeq42=+(+1.+rx)/2.(+1.-ry)/2. goto 1000 3 fbeq42=+(+1.+rx)/2.(+1.+ry)/2. goto 1000 4 fbeq42=+(+1.-rx)/2.(+1.+ry)/2. goto 1000 1000 return end subroutine beq43(refc,shpr,rctr,crtr) implicit real8 (a-h,o-z) dimension refc(2),shpr(4,3),rctr(2,2),crtr(2,2) external fbeq43 rx=refc(1) ry=refc(2) call dshap(fbeq43,refc,shpr,2,4,1) return end real8 function fbeq43(refc,n) implicit real8 (a-h,o-z) common /ccbeq4/ xa(4),ya(4),ua(4),va(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 fbeq43=+(+1.-rx)/2.(+1.-ry)/2. goto 1000 2 fbeq43=+(+1.+rx)/2.(+1.-ry)/2. goto 1000 3 fbeq43=+(+1.+rx)/2.(+1.+ry)/2. goto 1000 4 fbeq43=+(+1.-rx)/2.(+1.+ry)/2. goto 1000 1000 return end subroutine tbeq4(refc,coor,coorr,coefr,rc) implicit real8 (a-h,o-z) dimension refc(2),coor(2),coorr(2,4),coefr(4,2),rc(2,2) common /ccbeq4/ x(4),y(4),u(4),v(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 u(n)=coefr(n,1) v(n)=coefr(n,2) 200 continue rx=refc(1) ry=refc(2) call dcoor(ftbeq4,refc,coor,rc,2,2,1) return end real8 function ftbeq4(refc,n) implicit real8 (a-h,o-z) dimension refc(2) common /ccbeq4/ x(4),y(4),u(4),v(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 real8 (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 real8 function febeq4(refc,n) implicit real8 (a-h,o-z) dimension refc(2) common /ccbeq4/ xa(4),ya(4),u(4),v(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 1 febeq4=+(+(+1.-rx)/2.(+1.-ry)/2.)u(1)+(+(+1.+rx)/ & 2.(+1.-ry)/2.)u(2)+(+(+1.+rx)/2.(+1.+ry)/2.)u(3)+(+ & (+1.-rx)/2.(+1.+ry)/2.)u(4) goto 1000 2 febeq4=+(+(+1.-rx)/2.(+1.-ry)/2.)v(1)+(+(+1.+rx)/ & 2.(+1.-ry)/2.)v(2)+(+(+1.+rx)/2.(+1.+ry)/2.)v(3)+(+ & (+1.-rx)/2.(+1.+ry)/2.)*v(4) goto 1000 1000 return end 【Close】