Télécharger mesuvo.eso

Retour à la liste

Numérotation des lignes :

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

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