Télécharger itgldp.eso

Retour à la liste

Numérotation des lignes :

  1. C ITGLDP SOURCE CB215821 16/04/21 21:17:18 8920
  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
  32. c call der_enerelas_dpral(xini,x,nmodel)
  33. c write(ues,*)' Tensiones iniciales'
  34. c write(ues,*)(x(i),i=1,3)
  35. c call Invari_t(x,3,cc)
  36. c call Invari_J2(x,3,bb)
  37. c call Invari_I1(x,3,aa)
  38. c write(ues,999)0,aa,bb,cc*45.D0/atan(1.D0)
  39. c if (ndimx.eq.4) write(ues,*)' Variable interna inicial',x(4)
  40. c endif
  41.  
  42. c muestra iteraciones de line-search
  43. uesb=199
  44. escr=.false.
  45. if (nescri.eq.1) escr=.true.
  46.  
  47. c parametros generales
  48. kerre=0
  49. n=ndimx+1
  50. ndimv=1
  51. call equv(x,xini,ndimx)
  52. x(n)=lam
  53. beta=1.D-4
  54. j=0
  55.  
  56. 10 continue
  57. j=j+1
  58. call residuo_dpral(x,r,n,ndimx,nmodel,xtri,sig,vecm)
  59. if (nescri.eq.2) then
  60. write(ues,*)' TRAS residuo_dpral'
  61. write(ues,*)(x(i),i=1,n)
  62. write(ues,*)(r(i),i=1,n)
  63. write(ues,*)(vecm(i),i=1,ndimx)
  64. endif
  65.  
  66. c control constrain gamma positiva
  67. if (x(n).le.0.D0) then
  68. aux=0.D0
  69. do i=1,ndimx
  70. aux=aux+vecm(i)*(x(i)-xtri(i))
  71. enddo
  72. if (aux.gt.0.D0)
  73. . write(ues,*)' ERROR: Active constrain gam>0 !',x(n),aux
  74. c if (x(n).lt.0.D0)
  75. c . write(ues,*)' GAMMA negative !',x(n),' => ',0.D0
  76. x(n)=0.D0
  77. endif
  78.  
  79. c line-search cero
  80. f0=norm_l2_2(r,n)
  81. fp0=-2.D0*f0
  82. call direction_dpral(x,r,dx,n,ndimx,nmodel,sig,vecm,nnumer,deltax)
  83. c if (nescri.eq.2) then
  84. c write(ues,*)' TRAS direction_dpral'
  85. c write(ues,*)(x(i),i=1,n)
  86. c write(ues,*)(r(i),i=1,n)
  87. c write(ues,*)(dx(i),i=1,n)
  88. c endif
  89. err0=normainf(dx,ndimx+1)
  90. lam0=1.D0
  91. call addvec(xnew,x,lam0,dx,n)
  92. call residuo_dpral(xnew,rnew,n,ndimx,nmodel,xtri,sig,vecm)
  93. f1=norm_l2_2(rnew,n)
  94. cotaere=f0+beta*lam0*fp0
  95. lam1=lam0
  96. if (escr) write(uesb,*)' more ',0,lam0
  97. if (f1.lt.cotaere) goto 80
  98. if ((f0.lt.tolrel).and.(f1.lt.tolrel)) goto 80
  99. if (kmax.le.0) goto 70
  100. c line-search uno
  101. k=1
  102. lam1=(-fp0)/2.D0/(f1-f0-fp0)
  103. if (lam1.lt.0.1D0) lam1=0.1D0
  104. call addvec(xnew,x,lam1,dx,n)
  105. call residuo_dpral(xnew,rnew,n,ndimx,nmodel,xtri,sig,vecm)
  106. f2=norm_l2_2(rnew,n)
  107. cotaere=f0+beta*lam1*fp0
  108. if (escr) write(uesb,*)' more ',k,lam1
  109. if (f2.lt.cotaere) goto 80
  110. if (kmax.le.1) goto 70
  111. 20 continue
  112. c posteriores line-search
  113. k=k+1
  114. lam2=(-fp0*lam1**2)/2.D0/(f2-f0-fp0*lam1)
  115. if (lam2.lt.(0.1D0*lam1)) lam2=0.1D0*lam1
  116. if (lam2.gt.(0.5D0*lam1)) lam2=0.5D0*lam1
  117. call addvec(xnew,x,lam2,dx,n)
  118. call residuo_dpral(xnew,rnew,n,ndimx,nmodel,xtri,sig,vecm)
  119. f3=norm_l2_2(rnew,n)
  120. cotaere=f0+beta*lam2*fp0
  121. f1=f2
  122. f2=f3
  123. lam0=lam1
  124. lam1=lam2
  125. if (escr) write(uesb,*)' more ',k,lam2
  126. if (f3.lt.cotaere) goto 80
  127. if (k.lt.kmax) goto 20
  128. 70 continue
  129. c if (escr)
  130. if (escr) write(uesb,*)' more line search',kmax
  131. 80 continue
  132. c final line-search
  133. err=err0*lam1/max(normainf(xnew,ndimx),ABS(xnew(ndimx+1)))
  134. call equv(x,xnew,n)
  135. c muestra resultados iteraciones
  136. c if (nescri.eq.1) then
  137. c call der_enerelas_dpral(x,xnew,nmodel)
  138. c call Invari_t(xnew,3,cc)
  139. c call Invari_J2(xnew,3,bb)
  140. c call Invari_I1(xnew,3,aa)
  141. c write(ues,999)j,(x(i),i=1,n),err
  142. c write(ues,999)j,aa,bb,cc*45.D0/atan(1.D0),err,rnew(n)
  143. 999 format(I3,' ',5(E14.8,1x))
  144. c endif
  145. c control de convergencia
  146. if (err.lt.tolrel) goto 99
  147. if (j.lt.jmax) goto 10
  148. kerre=1
  149. 99 continue
  150. call equv(xini,x,ndimx)
  151. lam=x(n)
  152.  
  153. c muestra resultados finales
  154. c if (nescri.eq.2) then
  155. c write(ues,*)' Deform elas finales'
  156. c write(ues,*)(x(i),i=1,3)
  157. c write(ues,*)' Deform plas finales'
  158. c write(ues,*)(xtri(i)-x(i),i=1,3)
  159. c if (ndimx.eq.4) write(ues,*)' Variable interna final',x(4)
  160. c call der_enerelas_dpral(x,xini,nmodel)
  161. c write(ues,*)' Tensiones finales'
  162. c write(ues,*)(xini(i),i=1,3)
  163. c if (ndimx.eq.4) write(ues,*)' Variable interna inicial',x(4)
  164. c call equv(xini,x,ndimx)
  165. c endif
  166. end
  167.  
  168.  
  169.  
  170.  
  171.  
  172.  
  173.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales