Télécharger mesuvo.eso

Retour à la liste

Numérotation des lignes :

  1. C MESUVO SOURCE BP208322 16/11/18 21:19:14 9177
  2. C mesure du volume contenu dans une surface fermee
  3. C
  4. subroutine mesuvo(meleme,xmesu)
  5. IMPLICIT INTEGER(I-N)
  6. IMPLICIT REAL*8(A-H,O-Z)
  7. -INC CCOPTIO
  8. -INC SMELEME
  9. -INC CCGEOME
  10. -INC SMCOORD
  11. C Segments locaux
  12. C icpr tableau de compression des noeuds concernes
  13. C kon (numero point compresse par icpr,rang d'element,1) => pointeur su meleme
  14. C kon (numero point compresse par icpr,rang d'element,2) => numero element
  15. C ifut (rang d'element,1) => pointeur sur meleme
  16. C ifut (rang d'element,2) => element traite
  17. C ifut (rang d'element,3) => orientation a utiliser
  18. C
  19. segment kon(ite,nbk,3)
  20. segment ikon(ite)
  21. segment icpr(xcoor(/1)/(idim+1))
  22. segment ifut(nbel,3)
  23. xmesu=0.d0
  24. nbel=0
  25. segini icpr
  26. ite=0
  27. segact meleme*mod
  28. ipt1=meleme
  29. do 10 isous=1,max(1,lisous(/1))
  30. if (lisous(/1).ne.0) then
  31. ipt1=lisous(isous)
  32. segact ipt1*mod
  33. if (ipt1.itypel.ne.ksurf(ipt1.itypel)) call erreur(25)
  34. if (ierr.ne.0) return
  35. endif
  36. * cas des faces quaf tri7/qua9 : pas de prise en compte du noeud central
  37. if (ipt1.itypel.eq.7.or.ipt1.itypel.eq.11) then
  38. idec1=1
  39. else
  40. idec1=0
  41. endif
  42. nbel=nbel+ipt1.num(/2)
  43. do 15 ipt=1,ipt1.num(/1)-idec1
  44. do 151 iel=1,ipt1.num(/2)
  45. if (icpr(ipt1.num(ipt,iel)).ne.0) goto 151
  46. ite=ite+1
  47. icpr(ipt1.num(ipt,iel))=ite
  48. 151 continue
  49. 15 continue
  50. 10 continue
  51. nbk=8
  52. C on remplit kon qui donne les elements touchant un noeud (compacte)
  53. segini kon,ikon
  54. do 20 isous=1,max(1,lisous(/1))
  55. if (lisous(/1).ne.0) ipt1=lisous(isous)
  56. * cas des faces quaf tri7/qua9 : pas de prise en compte du noeud central
  57. if (ipt1.itypel.eq.7.or.ipt1.itypel.eq.11) then
  58. idec1=1
  59. else
  60. idec1=0
  61. endif
  62. do 25 ipt=1,ipt1.num(/1)-idec1
  63. do 251 iel=1,ipt1.num(/2)
  64. ip=icpr(ipt1.num(ipt,iel))
  65. ikon(ip)=ikon(ip)+1
  66. if (ikon(ip).gt.nbk) then
  67. nbk=nbk+1
  68. segadj kon
  69. endif
  70. kon(ip,ikon(ip),1)=ipt1
  71. kon(ip,ikon(ip),2)=iel
  72. 251 continue
  73. 25 continue
  74. 20 continue
  75. C on commence a travailler
  76. segini ifut
  77. do 30 i=1,nbel
  78. ifut(i,3)=1
  79. 30 continue
  80. ilcou=0
  81. iltop=1
  82. ifut(1,1)=ipt1
  83. ifut(1,2)=1
  84. ifut(1,3)=1
  85. ipt1.num(1,1)=-ipt1.num(1,1)
  86. 100 continue
  87. ilcou=ilcou+1
  88. if (ilcou.gt.iltop) goto 500
  89. ipt2=ifut(ilcou,1)
  90. iel=ifut(ilcou,2)
  91. isens=ifut(ilcou,3)
  92. * cas des faces quaf tri7/qua9 : pas de prise en compte du noeud central
  93. if (ipt2.itypel.eq.7.or.ipt2.itypel.eq.11) then
  94. idec2=1
  95. else
  96. idec2=0
  97. endif
  98.  
  99. C on cherche les elements contactant iel et on les oriente correctement
  100. C on note qu'ils sont traites en mettant en neg le premier noeud
  101. do 110 i=1,ipt2.num(/1)-idec2
  102. ip1=abs(ipt2.num(i,iel))
  103. i1=i+1
  104. if (i1.gt.ipt2.num(/1)-idec2) i1=1
  105. i0=i-1
  106. if (i0.lt.1) i0=ipt2.num(/1)-idec2
  107. if (isens.eq.1) ip2=abs(ipt2.num(i1,iel))
  108. if (isens.eq.-1) ip2=abs(ipt2.num(i0,iel))
  109. do 120 j=1,ikon(icpr(ip1))
  110. ipt3=kon(icpr(ip1),j,1)
  111. jel=kon(icpr(ip1),j,2)
  112. jsens=isign(1,ipt3.num(2,jel))
  113. * write (6,*) ' isens,jsens,iel,jel ',isens,jsens,iel,jel,ip1,ip2
  114. if (ipt3.eq.ipt2.and.jel.eq.iel) goto 120
  115. * cas des faces quaf tri7/qua9 : pas de prise en compte du noeud central
  116. if (ipt3.itypel.eq.7.or.ipt3.itypel.eq.11) then
  117. idec3=1
  118. else
  119. idec3=0
  120. endif
  121. do 130 in=1,ipt3.num(/1)-idec3
  122. if (abs(ipt3.num(in,jel)).ne.ip1) goto 130
  123. in1=in+1
  124. if (in1.gt.ipt3.num(/1)-idec3) in1=1
  125. in0=in-1
  126. if (in0.lt.1) in0=ipt3.num(/1)-idec3
  127. if ((abs(ipt3.num(in1,jel)).eq.ip2.and.jsens.eq.1).or.
  128. * (abs(ipt3.num(in0,jel)).eq.ip2.and.jsens.eq.-1))
  129. $ then
  130. C Il faut reorienter l'element
  131. if (iimpi.eq.5) write(6,*) ' element ',ipt3,jel
  132. $ ,' change de sens'
  133. if (ipt3.num(1,jel).lt.0) goto 999
  134. if (ierr.ne.0) return
  135. ipt3.num(1,jel)=-abs(ipt3.num(1,jel))
  136. iltop=iltop+1
  137. ifut(iltop,1)=ipt3
  138. ifut(iltop,2)=jel
  139. ifut(iltop,3)=-jsens
  140. ipt3.num(2,jel)=-abs(ipt3.num(2,jel))
  141. goto 110
  142. elseif ((abs(ipt3.num(in1,jel)).eq.ip2.and.jsens.eq.-1)
  143. * .or.(abs(ipt3.num(in0
  144. $ ,jel)).eq.ip2.and.jsens.eq.1)) then
  145. C Element OK
  146. if (iimpi.eq.5) write(6,*) ' element ',ipt3,jel,
  147. $ ' bon sens '
  148. if (ipt3.num(1,jel).lt.0) goto 110
  149. ipt3.num(1,jel)=-abs(ipt3.num(1,jel))
  150. iltop=iltop+1
  151. ifut(iltop,1)=ipt3
  152. ifut(iltop,2)=jel
  153. ifut(iltop,3)=jsens
  154. goto 110
  155. else
  156. endif
  157. 130 continue
  158. 120 continue
  159. goto 999
  160. 110 continue
  161. goto 100
  162. 500 continue
  163. C il ne reste plus qu'a calculer le volume et remettre les elements corrects
  164. if (iimpi.eq.5) write (6,*) ' iltop nbel ',iltop,nbel
  165. if (iltop.ne.nbel) goto 999
  166. if (ierr.ne.0) return
  167. do 600 i=1,nbel
  168. ipt1=ifut(i,1)
  169. iel=ifut(i,2)
  170. isens=ifut(i,3)
  171. IAD=LTEL(2,IPT1.ITYPEL)-1
  172. kdeg=kdegre(IPT1.ITYPEL)-1
  173. ITYP=LDEL(1,IAD+1)
  174. NPFAC=KDFAC(1,ITYP)
  175. JAD=LDEL(2,IAD+1)-1
  176. IDEP=KDFAC(2,ITYP)
  177. IFEP=IDEP+3*(KDFAC(3,ITYP)-1)
  178. dxmesu=0.d0
  179. DO 610 ITRIAN=IDEP,IFEP,3
  180. IAFA=LFAC(JAD+KFAC(ITRIAN))
  181. IBFA=LFAC(JAD+KFAC(ITRIAN+1))
  182. ICFA=LFAC(JAD+KFAC(ITRIAN+2))
  183. IA=abs(IPT1.NUM(IAFA,IEL))
  184. IB=abs(IPT1.NUM(IBFA,IEL))
  185. IC=abs(IPT1.NUM(ICFA,IEL))
  186. * if (iimpi.eq.5) write (6,*) ' triangle ',ia,ib,ic
  187. dxmesu=dxmesu+isens*voltri(0.d0,0.d0,0.d0,ia,ib,ic)/12
  188. 610 continue
  189. C SG 20160712 : comme dans mesu.eso, pour les faces TRI7, QUA9 (ityp
  190. C =7,8), il n'y a pas besoin de dedoubler car le decoupage de la
  191. C face est symetrique (barycentrique)
  192. IF (ITYP.NE.7.AND.ITYP.NE.8) THEN
  193. DO 611 ITRIAN=IDEP,IFEP,3
  194. ifa=KFAC(ITRIAN)+kdeg
  195. if (ifa.gt.npfac) ifa=ifa-npfac
  196. ifb=KFAC(ITRIAN+1)+kdeg
  197. if (ifb.gt.npfac) ifb=ifb-npfac
  198. ifc=KFAC(ITRIAN+2)+kdeg
  199. if (ifc.gt.npfac) ifc=ifc-npfac
  200. IAFA=LFAC(JAD+ifa)
  201. IBFA=LFAC(JAD+ifb)
  202. ICFA=LFAC(JAD+ifc)
  203. IA=abs(IPT1.NUM(IAFA,IEL))
  204. IB=abs(IPT1.NUM(IBFA,IEL))
  205. IC=abs(IPT1.NUM(ICFA,IEL))
  206. dxmesu=dxmesu+isens*voltri(0.d0,0.d0,0.d0,ia,ib,ic)/12
  207. 611 continue
  208. ELSE
  209. dxmesu=dxmesu*2
  210. ENDIF
  211. xmesu=xmesu+dxmesu
  212. 600 continue
  213. goto 998
  214. 999 continue
  215. call erreur(127)
  216. 998 continue
  217. C maintenant on desactive les segments et on remet les elements correctement
  218. ipt1=meleme
  219. do 620 isous=1,max(1,lisous(/1))
  220. if (lisous(/1).ne.0) ipt1=lisous(isous)
  221. do 630 iel=1,ipt1.num(/2)
  222. ipt1.num(1,iel)=abs(ipt1.num(1,iel))
  223. ipt1.num(2,iel)=abs(ipt1.num(2,iel))
  224. 630 continue
  225. segdes ipt1
  226. 620 continue
  227. segdes meleme
  228. segsup icpr,kon,ikon,ifut
  229. xmesu=abs(xmesu)
  230. end
  231.  
  232.  
  233.  
  234.  
  235.  
  236.  
  237.  
  238.  
  239.  
  240.  
  241.  
  242.  
  243.  
  244.  
  245.  
  246.  
  247.  

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