Télécharger ravc.eso

Retour à la liste

Numérotation des lignes :

  1. C RAVC SOURCE PV 13/12/06 21:15:01 7877
  2. subroutine ravc
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8(A-H,O-Z)
  5. -INC CCOPTIO
  6. -INC SMCHPOI
  7. -INC SMELEME
  8. -INC SMCOORD
  9. segment icpr(xcoor(/1)/(idim+1))
  10. call lirobj ( 'CHPOINT',MCHPO1, 1,iretou)
  11. if(ierr . ne.0) return
  12. call lirobj ( 'CHPOINT',MCHPO2, 1,iretou)
  13. if(ierr . ne.0) return
  14. call lirobj ( 'CHPOINT',MCHPO3, 1,iretou)
  15. if(ierr . ne.0) return
  16. segact mchpo1, mchpo2,mchpo3
  17. if(mchpo1.ipchp(/1). ne . 1) then
  18. * write(6,*) ' erreur 1 '
  19. call erreur(19)
  20. return
  21. endif
  22. if(mchpo2.ipchp(/1). ne . 1) then
  23. * write(6,*) ' erreur 2 '
  24. call erreur(19)
  25. return
  26. endif
  27. if(mchpo3.ipchp(/1). ne . 1) then
  28. * write(6,*) ' erreur 3 '
  29. call erreur(19)
  30. return
  31. endif
  32. msoup1= mchpo1.ipchp(1)
  33. msoup2= mchpo2.ipchp(1)
  34. msoup3= mchpo3.ipchp(1)
  35. segact msoup1,msoup2,msoup3
  36. ipt1=msoup1.igeoc
  37. ipt2=msoup2.igeoc
  38. ipt3=msoup3.igeoc
  39. segact ipt1,ipt2,ipt3
  40. segini icpr
  41. n1 = ipt1.num(/2)
  42. n2 = ipt2.num(/2)
  43. n3 = ipt3.num(/2)
  44. ipo=0
  45. do 100 i=1,n1
  46. ia=ipt1.num(1,i)
  47. if(icpr(ia).eq.0) then
  48. ipo=ipo+1
  49. icpr(ia)=ipo
  50. endif
  51. 100 continue
  52. do 101 i=1,n2
  53. ia=ipt2.num(1,i)
  54. if(icpr(ia).eq.0) then
  55. ipo=ipo+1
  56. icpr(ia)=ipo
  57. endif
  58. 101 continue
  59. do 102 i=1,n3
  60. ia=ipt3.num(1,i)
  61. if(icpr(ia).eq.0) then
  62. ipo=ipo+1
  63. icpr(ia)=ipo
  64. endif
  65. 102 continue
  66. nbelem=ipo
  67. nbnn=1
  68. nbsous=0
  69. nbref=0
  70. segini meleme
  71. itypel=1
  72. do 103 i=1,icpr(/1)
  73. if(icpr(I).ne.0) then
  74. num(1,icpr(i))=i
  75. endif
  76. 103 continue
  77.  
  78. * write(6,*) ' ipt1',(ipt1.num(1,iou),iou=1,ipt1.num(/2))
  79. * write(6,*) ' ipt2',(ipt2.num(1,iou),iou=1,ipt2.num(/2))
  80. * write(6,*) ' ipt3',(ipt3.num(1,iou),iou=1,ipt3.num(/2))
  81. * if( (n1.ne.n2). or .(n1.ne.n3))then
  82. * write(6,*) ' erreur 4 '
  83. * call erreur(19)
  84. * return
  85. * endif
  86. if( msoup1.noharm(/1) . ne . 1) then
  87. * write(6,*) ' erreur 5 '
  88. call erreur(19)
  89. return
  90. endif
  91. if( msoup2.noharm(/1) . ne . 1) then
  92. * write(6,*) ' erreur 6 '
  93. call erreur(19)
  94. return
  95. endif
  96. if( msoup3.noharm(/1) . ne . 1) then
  97. * write(6,*) ' erreur 7 '
  98. call erreur(19)
  99. return
  100. endif
  101. segini,mchpoi=mchpo1
  102. segini,msoupo=msoup1
  103. ipchp(1)=msoupo
  104. mpova1=msoup1.ipoval
  105. mpova2=msoup2.ipoval
  106. mpova3=msoup3.ipoval
  107. segact mpova1,mpova2,mpova3
  108. n=nbelem
  109. nc=1
  110. segini mpoval
  111. ipoval=mpoval
  112. igeoc=meleme
  113. * write(6,*) ' meleme' , ( num(1,iou),iou=1,nbelem)
  114. *
  115. * on cherche les max de chaque champs
  116. *
  117. * amma=0.
  118. * bmma=0.
  119. * cmma=0.
  120. * do 10 i=1,nbelem
  121. * if(amma.le.(abs(mpova1.vpocha(i,1))))amma=abs(mpova1.vpocha(i,1))
  122. * if(bmma.le.(abs(mpova2.vpocha(i,1))))bmma=abs(mpova2.vpocha(i,1))
  123. * if(cmma.le.(abs(mpova3.vpocha(i,1))))cmma=abs(mpova3.vpocha(i,1))
  124. * 10 continue
  125. *
  126. * calcul des racines
  127. *
  128. do 1 j=1,nbelem
  129. ia = num(1,j)
  130. * write(6,*) ' ia ' , ia
  131. a = 0.d0
  132. do 8 i=1,n1
  133. if( ipt1.num(1,i).eq.ia)then
  134. a = mpova1.vpocha(i,1)
  135. go to 9
  136. endif
  137. 8 continue
  138. 9 continue
  139. b=0.d0
  140. do 2 k=1,n2
  141. if(ipt2.num(1,k).eq.ia) then
  142. b = mpova2.vpocha(k,1)
  143. go to 3
  144. endif
  145. 2 continue
  146. 3 continue
  147. c=0.d0
  148. do 4 l=1,n3
  149. if(ipt3.num(1,l).eq.ia)then
  150. c = mpova3.vpocha(l,1)
  151. go to 5
  152. endif
  153. 4 continue
  154. 5 continue
  155.  
  156. * test sur les grandeurs
  157. * write (6,*) ' a b c ',a,b,c
  158. if( a . eq. 0.d0. or . abs(a).lt. (abs(b)*1.e-13))then
  159. * on resoud une seule equation
  160. if( b . ne . 0.d0) then
  161. xx = -c / b
  162. if( xx .ge. 0.d0 ) xx = 0.d0
  163. else
  164. xx = 0.d0
  165. endif
  166. else
  167. xx=0.d0
  168. delta2 = ( b*b - 4*a*c)
  169. if (delta2.lt. b*b*1d-13) then
  170. xx = -b/(2.d0*a)
  171. xx = min(xx,0.d0)
  172. xx = 0.d0
  173. ** write (6,*) ' ravc discriminant negatif ',delta2,b
  174. else
  175. delta2=max(0.d0,delta2)
  176. delta = sqrt (delta2)
  177. s1= (-b + delta )/ (2.*a)
  178. s2= (-b - delta) / (2.*a)
  179. s1=min(s1,0.d0)
  180. s2=min(s2,0.d0)
  181. * write (6,*) ' s1 s2 ' , s1,s2
  182. if( s1. gt . s2) then
  183. if( s1 . lt. 0.d0) then
  184. xx = s1
  185. else
  186. if(s2.lt.0.d0) xx = s2
  187. endif
  188. else
  189. if( s2 . lt. 0.d0) then
  190. xx = s2
  191. else
  192. if(s1.lt.0.d0) xx = s1
  193. endif
  194. endif
  195. endif
  196. endif
  197. ** write (6,*) ' xx ',xx
  198. vpocha(j,1)=xx
  199. 1 continue
  200. segdes mpova1,mpova2,mpova3,mpoval
  201. segdes ipt1,ipt2,ipt3,meleme
  202. segdes msoupo,msoup1,msoup2,msoup3
  203. segdes mchpoi,mchpo1,mchpo2,mchpo3
  204. segsup icpr
  205. call ecrobj('CHPOINT',mchpoi)
  206. return
  207. end
  208.  
  209.  
  210.  
  211.  
  212.  
  213.  
  214.  
  215.  
  216.  
  217.  
  218.  

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