Télécharger psiphi.eso

Retour à la liste

Numérotation des lignes :

  1. C PSIPHI SOURCE CB215821 19/07/31 21:17:16 10277
  2. SUBROUTINE PSIPHI
  3. ***********************************************************************
  4. c OPERATEUR : PSIP
  5. c
  6. c APPEL :
  7. c CHP1 (CHP2) (CHP3) = PSIP MAIL1 MAIL2 (CRIT1) (|'DEUX' | | P1 (P2) | );
  8. c |'TROI' | | MAIL3 |
  9. c
  10. c FONCTION : calcule les fonctions distances signées (level set)
  11. c relatives aux maillages MAIL2 (surface de fissure)
  12. c et MAIL3 en 3D (front de fissure)
  13. c ou P1 en 2D (pointe de fissure)
  14. c aux noeuds de MAIL1 (pas trop éloignés).
  15. c
  16. c CREATION : chat (16/07/2007)
  17. c MODIFS : bp (2009 -> 2012) : diverses corrections/ameliorations
  18. c repertoriees par les fiches d'anomalies/developpement
  19. c bp 14/03/2012 : on split psiphi en psiphi psip2d psip3d et
  20. c zonag2 (evolué seulement le 18/12/2013)
  21. c
  22. c TO DO : cas 3d avec plusieurs fronts (ou 1 front discontinu)
  23. c
  24. ***********************************************************************
  25. IMPLICIT REAL*8(A-H,O-Z)
  26. IMPLICIT INTEGER(I-N)
  27. -INC CCOPTIO
  28. -INC SMCOORD
  29. -INC SMELEME
  30. -INC SMCHPOI
  31. -INC CCREEL
  32. PARAMETER(NMOCLE=2)
  33. CHARACTER*4 MOCLE(NMOCLE)
  34. DATA MOCLE/'DEUX','TROI'/
  35.  
  36. if(iimpi.ge.1) write(ioimp,*)'==== ENTREE DANS PSIPHI ====='
  37.  
  38. ***********************************************************************
  39. * INITIALISATIONs et LECTUREs
  40. ***********************************************************************
  41. *
  42. * calcule t'on phi? phi ET psi? phi, psi et tau?
  43. ideux=0
  44. call lirmot (mocle,nmocle,ideux,0)
  45. if (ideux.eq.2.and.idim.lt.3) then
  46. write(*,*) ' ON NE PEUT CALCULER TROIS LEVEL SET QU EN 3D !'
  47. write(*,*) ' ON CONTINUE AVEC LE CALCUL DE 2 LEVEL SET ...'
  48. ideux=1
  49. endif
  50.  
  51. * lecture des points pour lesquels on veut calculer phi et psi
  52. call LIROBJ('MAILLAGE',ipt1,1,iretou)
  53. melmai=ipt1
  54. if(ierr.ne.0) return
  55. call change(ipt1,1)
  56. c segact ipt1
  57. c rem: inutile car change laisse ipt1 actif
  58.  
  59. * lecture du maillage de la fissure
  60. call LIROBJ('MAILLAGE',meleme,1,iretou)
  61. melfis=meleme
  62. if(ierr.ne.0) return
  63.  
  64. * lecture de la pointe de fissure (objet de type point) ou du front
  65. ip1=0
  66. ip2=0
  67. melfro=0
  68. 1 continue
  69. if (ideux.ge.1) then
  70. if (idim.eq.2) then
  71. call LIROBJ('POINT',ipt,0,iretou)
  72. if (iretou.ne.0) then
  73. if (ip1.eq.0) then
  74. ip1 = ipt
  75. go to 1
  76. else
  77. ip2 = ipt
  78. endif
  79. endif
  80. else
  81. call LIROBJ('MAILLAGE', melfro ,1,iretou)
  82. if(iretou.eq.0) write(ioimp,*) 'Il manque le maillage du front'
  83. if(IERR.NE.0) return
  84. endif
  85. endif
  86.  
  87. * lecture facultative d'une longueur max (=plus grande taille des elements
  88. * concernés par la fissure)
  89. xcrit=0.D0
  90. call lirree(xcrit,0,ircrit)
  91.  
  92.  
  93.  
  94. ***********************************************************************
  95. * CRÉATION DES MCHPOI DE SORTIE (1 POUR PSI , 2 POUR PHI, 3 POUR TAU)
  96. ***********************************************************************
  97. mpova1=0
  98. mpova2=0
  99. mpova3=0
  100.  
  101. *-----PHI-----
  102. c if(idebug.eq.1) write(6,*) '----Creation de PHI----'
  103. nat=1
  104. nsoupo=1
  105. nc=1
  106. n=ipt1.num(/2)
  107. segini,mchpo2
  108. segini,msoup2
  109. mchpo2.jattri(1)=1
  110. mchpo2.ipchp(1)=msoup2
  111. mchpo2.ifopoi=ifour
  112. segdes mchpo2
  113. msoup2.igeoc=ipt1
  114. msoup2.nocomp(1)='PHI'
  115. msoup2.noharm(1)=nifour
  116. segini,mpova2
  117. msoup2.ipoval=mpova2
  118. c segdes,msoup2
  119.  
  120. *-----PSI-----
  121. if (ideux.ge.1) then
  122. c if(idebug.eq.1) write(6,*) '----Creation de PSI----'
  123. segini,mchpo1
  124. segini,msoup1
  125. mchpo1.jattri(1)=1
  126. mchpo1.ipchp(1)=msoup1
  127. mchpo1.ifopoi=ifour
  128. segdes,mchpo1
  129. msoup1.igeoc=ipt1
  130. msoup1.nocomp(1)='PSI'
  131. msoup1.noharm(1)=nifour
  132. segini,mpova1
  133. msoup1.ipoval=mpova1
  134. c segdes,msoup1
  135. else
  136. msoup1 = 0
  137. endif
  138.  
  139. *-----TAU-----
  140. if (ideux.ge.2) then
  141. c if(idebug.eq.1) write(6,*) '----Creation de TAU----'
  142. segini,mchpo3
  143. segini,msoup3
  144. mchpo3.jattri(1)=1
  145. mchpo3.ipchp(1)=msoup3
  146. mchpo3.ifopoi=ifour
  147. segdes,mchpo3
  148. msoup3.igeoc=ipt1
  149. msoup3.nocomp(1)='TAU'
  150. msoup3.noharm(1)=nifour
  151. segini,mpova3
  152. msoup3.ipoval=mpova3
  153. c segdes,msoup3
  154. else
  155. msoup3 = 0
  156. endif
  157.  
  158.  
  159. ***********************************************************************
  160. * PSIPHI 3D/2D
  161. ***********************************************************************
  162. if (idim.eq.3) then
  163. c write(ioimp,*)'appel PSIP3D (',ideux,ipt1,melfis,melfro,xcrit
  164. call PSIP3D(ideux,ipt1,melfis,melfro,xcrit,
  165. & msoup1,msoup2,msoup3)
  166. else
  167. c write(ioimp,*)'appel PSIP3D (',ideux,ipt1,melfis,ip1,ip2,xcrit
  168. call PSIP2D(ideux,ipt1,melfis,ip1,ip2,xcrit,
  169. & msoup1,msoup2)
  170. endif
  171.  
  172.  
  173. ***********************************************************************
  174. * ECRITURE ET FIN DU PROGRAMME
  175. ***********************************************************************
  176.  
  177. c-----ecriture de (PSI) PHI ((TAU)) --------------
  178. if(ideux.ge.2) then
  179. call ACTOBJ('CHPOINT ',mchpo3,1)
  180. call ECROBJ('CHPOINT ',mchpo3)
  181. endif
  182. call ACTOBJ('CHPOINT ',mchpo2,1)
  183. call ECROBJ('CHPOINT ',mchpo2)
  184. if(ideux.ge.1) then
  185. call ACTOBJ('CHPOINT',mchpo1,1)
  186. call ECROBJ('CHPOINT',mchpo1)
  187. endif
  188.  
  189. end
  190.  
  191.  
  192.  

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