1 : C ITGLDP SOURCE CHAT 05/01/13 00:44:44 5004 2 : C LINESEARCH for FEFP models 3 : C ********************************************************************** 4 : CCC 5 : cc Trabaja en direcciones principales y las incognitas son DEFO ELAS! 6 : 7 : SUBROUTINE INTEGRA_LS_DPRAL(XTRI,XINI,NDIMX,LAM,NMODEL,TOLREL, 8 : . jmax,nescri,ues,nnumer,deltax,kerre,j) 9 : 10 : IMPLICIT INTEGER(I-N) 11 : integer nmodel,nitmax,nescri,nnumer,ues,j,kerre, 12 : . ndimv,jmax,i,ii,k,n,ndimx,uesb 13 : real*8 xtri(ndimx),xini(ndimx),lam,tolrel,deltax, 14 : . x(5),dx(5),r(5),err,normainf,norm_l2_2,xnew(5),rnew(5), 15 : . cotaere,lam0,lam1,lam2,f0,fp0,f1,f2,f3,s1,s2,aa,bb,cc, 16 : . sig(4),vecm(4),err0,aux,beta 17 : logical escr 18 : c kmax=0 => no line search 19 : integer level,kmax 20 : common /linesearch/ level,kmax 21 : 22 : c muestra tensiones 23 : c if (nescri.eq.1) then 24 : c write(ues,*)' Variables iniciales' 25 : c write(ues,*)(xini(i),i=1,ndimx) 26 : c write(ues,*)' Variables trial' 27 : c write(ues,*)(xtri(i),i=1,ndimx) 28 : c write(ues,*)' Lambda inicial' 29 : c write(ues,*)lam 30 : c write(ues,*)' Variables varias' 31 : c write(ues,*)level,kmax,nmodel,nnumer,deltax 31 : c call der_enerelas_dpral(xini,x,nmodel) 32 : c write(ues,*)' Tensiones iniciales' 33 : c write(ues,*)(x(i),i=1,3) 33 : c call Invari_t(x,3,cc) 33 : c call Invari_J2(x,3,bb) 33 : c call Invari_I1(x,3,aa) 34 : c write(ues,999)0,aa,bb,cc*45.D0/atan(1.D0) 35 : c if (ndimx.eq.4) write(ues,*)' Variable interna inicial',x(4) 36 : c endif 37 : 38 : c muestra iteraciones de line-search 39 : uesb=199 40 : escr=.false. 41 : if (nescri.eq.1) escr=.true. 42 : 43 : c parametros generales 44 : kerre=0 45 : n=ndimx+1 46 : ndimv=1 46 : call equv(x,xini,ndimx) 47 : x(n)=lam 48 : beta=1.D-4 49 : j=0 50 : 51 : 10 continue 52 : j=j+1 52 : call residuo_dpral(x,r,n,ndimx,nmodel,xtri,sig,vecm) 53 : if (nescri.eq.2) then 54 : write(ues,*)' TRAS residuo_dpral' 55 : write(ues,*)(x(i),i=1,n) 56 : write(ues,*)(r(i),i=1,n) 57 : write(ues,*)(vecm(i),i=1,ndimx) 58 : endif 59 : 60 : c control constrain gamma positiva 61 : if (x(n).le.0.D0) then 62 : aux=0.D0 63 : do i=1,ndimx 64 : aux=aux+vecm(i)*(x(i)-xtri(i)) 65 : enddo 66 : if (aux.gt.0.D0) 67 : . write(ues,*)' ERROR: Active constrain gam>0 !',x(n),aux 68 : c if (x(n).lt.0.D0) 69 : c . write(ues,*)' GAMMA negative !',x(n),' => ',0.D0 70 : x(n)=0.D0 71 : endif 72 : 73 : c line-search cero 74 : f0=norm_l2_2(r,n) 75 : fp0=-2.D0*f0 75 : call direction_dpral(x,r,dx,n,ndimx,nmodel,sig,vecm,nnumer,deltax) 76 : c if (nescri.eq.2) then 77 : c write(ues,*)' TRAS direction_dpral' 78 : c write(ues,*)(x(i),i=1,n) 79 : c write(ues,*)(r(i),i=1,n) 80 : c write(ues,*)(dx(i),i=1,n) 81 : c endif 82 : err0=normainf(dx,ndimx+1) 83 : lam0=1.D0 83 : call addvec(xnew,x,lam0,dx,n) 83 : call residuo_dpral(xnew,rnew,n,ndimx,nmodel,xtri,sig,vecm) 84 : f1=norm_l2_2(rnew,n) 85 : cotaere=f0+beta*lam0*fp0 86 : lam1=lam0 87 : if (escr) write(uesb,*)' more ',0,lam0 88 : if (f1.lt.cotaere) goto 80 89 : if ((f0.lt.tolrel).and.(f1.lt.tolrel)) goto 80 90 : if (kmax.le.0) goto 70 91 : c line-search uno 92 : k=1 93 : lam1=(-fp0)/2.D0/(f1-f0-fp0) 94 : if (lam1.lt.0.1D0) lam1=0.1D0 94 : call addvec(xnew,x,lam1,dx,n) 94 : call residuo_dpral(xnew,rnew,n,ndimx,nmodel,xtri,sig,vecm) 95 : f2=norm_l2_2(rnew,n) 96 : cotaere=f0+beta*lam1*fp0 97 : if (escr) write(uesb,*)' more ',k,lam1 98 : if (f2.lt.cotaere) goto 80 99 : if (kmax.le.1) goto 70 100 : 20 continue 101 : c posteriores line-search 102 : k=k+1 103 : lam2=(-fp0*lam1**2)/2.D0/(f2-f0-fp0*lam1) 104 : if (lam2.lt.(0.1D0*lam1)) lam2=0.1D0*lam1 105 : if (lam2.gt.(0.5D0*lam1)) lam2=0.5D0*lam1 105 : call addvec(xnew,x,lam2,dx,n) 105 : call residuo_dpral(xnew,rnew,n,ndimx,nmodel,xtri,sig,vecm) 106 : f3=norm_l2_2(rnew,n) 107 : cotaere=f0+beta*lam2*fp0 108 : f1=f2 109 : f2=f3 110 : lam0=lam1 111 : lam1=lam2 112 : if (escr) write(uesb,*)' more ',k,lam2 113 : if (f3.lt.cotaere) goto 80 114 : if (k.lt.kmax) goto 20 115 : 70 continue 116 : c if (escr) 117 : if (escr) write(uesb,*)' more line search',kmax 118 : 80 continue 119 : c final line-search 120 : err=err0*lam1/max(normainf(xnew,ndimx),dabs(xnew(ndimx+1))) 120 : call equv(x,xnew,n) 121 : c muestra resultados iteraciones 122 : c if (nescri.eq.1) then 122 : c call der_enerelas_dpral(x,xnew,nmodel) 122 : c call Invari_t(xnew,3,cc) 122 : c call Invari_J2(xnew,3,bb) 122 : c call Invari_I1(xnew,3,aa) 123 : c write(ues,999)j,(x(i),i=1,n),err 124 : c write(ues,999)j,aa,bb,cc*45.D0/atan(1.D0),err,rnew(n) 125 : 999 format(I3,' ',5(E14.8,1x)) 126 : c endif 127 : c control de convergencia 128 : if (err.lt.tolrel) goto 99 129 : if (j.lt.jmax) goto 10 130 : kerre=1 131 : 99 continue 131 : call equv(xini,x,ndimx) 132 : lam=x(n) 133 : 134 : c muestra resultados finales 135 : c if (nescri.eq.2) then 136 : c write(ues,*)' Deform elas finales' 137 : c write(ues,*)(x(i),i=1,3) 138 : c write(ues,*)' Deform plas finales' 139 : c write(ues,*)(xtri(i)-x(i),i=1,3) 140 : c if (ndimx.eq.4) write(ues,*)' Variable interna final',x(4) 140 : c call der_enerelas_dpral(x,xini,nmodel) 141 : c write(ues,*)' Tensiones finales' 142 : c write(ues,*)(xini(i),i=1,3) 143 : c if (ndimx.eq.4) write(ues,*)' Variable interna inicial',x(4) 143 : c call equv(xini,x,ndimx) 144 : c endif 145 : end 146 : 147 : 148 : 149 :
© Cast3M 2003 - All rights reserved.
Disclaimer