Télécharger mesu.eso

Retour à la liste

Numérotation des lignes :

  1. C MESU SOURCE GOUNAND 16/08/01 21:15:21 9043
  2.  
  3. C fournit la mesure d'un maillage
  4.  
  5. SUBROUTINE mesu
  6.  
  7. IMPLICIT INTEGER(I-N)
  8. IMPLICIT REAL*8 (A-H,O-Z)
  9.  
  10. -INC CCOPTIO
  11. -INC CCREEL
  12. -INC CCGEOME
  13. -INC SMCOORD
  14. -INC SMELEME
  15.  
  16. SEGMENT /FER/ (NFI(ITT),MAI(IPP),ITOUR)
  17. CHARACTER*4 mcle(3)
  18. DATA mcle / 'LONG','SURF','VOLU' /
  19.  
  20. *dbg write(ioimp,*) 'coucou mesu'
  21. ncle=IDIM
  22. icle=0
  23. CALL LIRMOT(mcle,ncle,icle,0)
  24. IF (IERR.NE.0) RETURN
  25. CALL LIROBJ('MAILLAGE',meleme,1,iOK)
  26. IF (IERR.NE.0) RETURN
  27.  
  28. idimp1=IDIM+1
  29. SEGACT,meleme
  30. ipt1=meleme
  31. itp=0
  32. xmesu=xzero
  33. nsous=lisous(/1)
  34. DO isou=1,MAX(nsous,1)
  35. IF (nsous.NE.0) THEN
  36. ipt1=lisous(isou)
  37. SEGACT,ipt1
  38. ENDIF
  39. jtyp=3
  40. IF (ksurf(ipt1.itypel).EQ.ipt1.itypel) jtyp=2
  41. IF (kdegre(ipt1.itypel).EQ.ipt1.itypel) jtyp=1
  42. IF (jtyp.EQ.1.AND.nsous.GT.1) CALL ERREUR(25)
  43. IF (itp.EQ.0) itp=jtyp
  44. IF (itp.NE.jtyp) CALL ERREUR(25)
  45. IF (icle.EQ.0) icle=itp
  46. IF (icle.EQ.1.AND.itp.NE.1) CALL ERREUR(21)
  47. IF (icle.EQ.2.AND.itp.EQ.3) CALL ERREUR(21)
  48. IF (icle.EQ.3.AND.itp.EQ.1) CALL ERREUR(21)
  49. IF (IERR.NE.0) RETURN
  50. GOTO (100,200,300),icle
  51. C calcul de la longueur d'une ligne
  52. 100 IF (IIMPI.NE.0) WRITE(6,*) ' Mesure d''une ligne'
  53. nbnn=ipt1.num(/1)
  54. IF (IDIM.EQ.1) THEN
  55. DO iel=1,ipt1.num(/2)
  56. ip1=idimp1*(ipt1.num(1,iel)-1)
  57. ip2=idimp1*(ipt1.num(nbnn,iel)-1)
  58. xp12=XCOOR(ip2+1)-XCOOR(ip1+1)
  59. xmesu=xmesu+ABS(xp12)
  60. ENDDO
  61. ELSE IF (IDIM.EQ.2) THEN
  62. DO iel=1,ipt1.num(/2)
  63. ip1=idimp1*(ipt1.num(1,iel)-1)
  64. ip2=idimp1*(ipt1.num(nbnn,iel)-1)
  65. xp12=XCOOR(ip2+1)-XCOOR(ip1+1)
  66. yp12=XCOOR(ip2+2)-XCOOR(ip1+2)
  67. xmesu=xmesu+SQRT(xp12*xp12+yp12*yp12)
  68. ENDDO
  69. ELSE IF (IDIM.EQ.3) THEN
  70. DO iel=1,ipt1.num(/2)
  71. ip1=idimp1*(ipt1.num(1,iel)-1)
  72. ip2=idimp1*(ipt1.num(nbnn,iel)-1)
  73. xp12=XCOOR(ip2+1)-XCOOR(ip1+1)
  74. yp12=XCOOR(ip2+2)-XCOOR(ip1+2)
  75. zp12=XCOOR(ip2+3)-XCOOR(ip1+3)
  76. xmesu=xmesu+SQRT(xp12*xp12+yp12*yp12+zp12*zp12)
  77. ENDDO
  78. ENDIF
  79. GOTO 500
  80. C calcul de la surface d'une aire
  81. 200 GOTO (220,240),jtyp
  82. C cas d'une ligne
  83. 220 IF (IIMPI.NE.0) WRITE(6,*) ' Mesure de surface dans un contour'
  84. CALL AVTRSF(ipt1,fer,ipt5)
  85. SEGSUP,ipt5
  86. ib=0
  87. xb=xzero
  88. yb=xzero
  89. zb=xzero
  90. DO it=1,itour
  91. imdeb=mai(it)+1
  92. imfin=mai(it+1)
  93. DO im=imdeb,imfin
  94. ib=ib+1
  95. ip=idimp1*(nfi(im)-1)
  96. xb=xb+XCOOR(ip+1)
  97. yb=yb+XCOOR(ip+2)
  98. IF (IDIM.EQ.3) zb=zb+XCOOR(ip+3)
  99. ENDDO
  100. ENDDO
  101. xb=xb/ib
  102. yb=yb/ib
  103. zb=zb/ib
  104. IF (IIMPI.NE.0) WRITE(6,*) ' Barycentre ',xb,yb,zb
  105. xs=xzero
  106. ys=xzero
  107. zs=xzero
  108. zv1=xzero
  109. zv2=xzero
  110. DO it=1,itour
  111. imdeb=mai(it)+1
  112. imfin=mai(it+1)
  113. DO im=imdeb,imfin
  114. im1=im+1
  115. IF (im1.GT.imfin) im1=imdeb
  116. ip1=idimp1*(nfi(im)-1)
  117. ip2=idimp1*(nfi(im1)-1)
  118. xv1=XCOOR(ip1+1)-xb
  119. yv1=XCOOR(ip1+2)-yb
  120. xv2=XCOOR(ip2+1)-xb
  121. yv2=XCOOR(ip2+2)-yb
  122. IF (IDIM.GE.3) THEN
  123. zv1=XCOOR(ip1+3)-zb
  124. zv2=XCOOR(ip2+3)-zb
  125. ENDIF
  126. xs=xs+(yv1*zv2-zv1*yv2)
  127. ys=ys+(zv1*xv2-xv1*zv2)
  128. zs=zs+(xv1*yv2-yv1*xv2)
  129. ENDDO
  130. ENDDO
  131. xmesu=xmesu+0.5D0*SQRT(xs*xs+ys*ys+zs*zs)
  132. SEGSUP,fer
  133. GOTO 500
  134. C cas d'une surface
  135. 240 IF (IIMPI.NE.0) WRITE(6,*) ' Mesure d''une surface'
  136. iad=LTEL(2,ipt1.itypel)-1
  137. jad=LDEL(2,iad+1)-1
  138. ityp=LDEL(1,iad+1)
  139. npfac=KDFAC(1,ityp)
  140. idep=KDFAC(2,ityp)
  141. ifep=idep+3*(KDFAC(3,ityp)-1)
  142. xs=xzero
  143. ys=xzero
  144. zs=xzero
  145. Cc=xzero
  146. DO iel=1,ipt1.num(/2)
  147. DO i=idep,ifep,3
  148. iafa=LFAC(jad+KFAC(i))
  149. ibfa=LFAC(jad+KFAC(i+1))
  150. icfa=LFAC(jad+KFAC(i+2))
  151. ia=ipt1.num(iafa,iel)
  152. ib=ipt1.num(ibfa,iel)
  153. ic=ipt1.num(icfa,iel)
  154. IF (IIMPI.EQ.5) WRITE(6,*) ' Triangle ',ia,ib,ic
  155. ia=idimp1*(ia-1)
  156. ib=idimp1*(ib-1)
  157. ic=idimp1*(ic-1)
  158. xb=XCOOR(ib+1)
  159. yb=XCOOR(ib+2)
  160. xv1=XCOOR(ia+1)-xb
  161. yv1=XCOOR(ia+2)-yb
  162. xv2=XCOOR(ic+1)-xb
  163. yv2=XCOOR(ic+2)-yb
  164. IF (IDIM.GE.3) THEN
  165. zb=XCOOR(ib+3)
  166. zv1=XCOOR(ia+3)-zb
  167. zv2=XCOOR(ic+3)-zb
  168. xs=yv1*zv2-zv1*yv2
  169. ys=zv1*xv2-xv1*zv2
  170. ENDIF
  171. zs=xv1*yv2-yv1*xv2
  172. Cc=Cc+SQRT(xs*xs+ys*ys+zs*zs)
  173. ENDDO
  174. ENDDO
  175. xmesu=xmesu+0.5D0*Cc
  176. GOTO 500
  177. C calcul d'un volume
  178. 300 GOTO (320,320,340),jtyp
  179. C cas du volume interieur a une surface
  180. C du fait du probleme d'orientation des elements on ne peut pas
  181. C traiter le probleme sous-maillage par sous maillage
  182. 320 CALL MESUVO(meleme,xmesu)
  183. GOTO 11
  184. C cas d'un volume. On va utiliser le decoupage en triangle de la
  185. C surface des elements mais en le dedoublant de facon a etre sur
  186. C de correspondre d'un element au voisin
  187. C SG 20160712 : pour les faces TRI7, QUA9 (ityp=7,8), il n'y a pas
  188. C besoin de dedoubler car le decoupage de la face est symetrique
  189. C (barycentrique)
  190. 340 nbfac=LTEL(1,ipt1.itypel)
  191. IF (nbfac.EQ.0) GOTO 500
  192. iad=LTEL(2,ipt1.itypel)-1
  193. kdeg=KDEGRE(ipt1.itypel)-1
  194. Cc=xzero
  195. DO iel=1,ipt1.num(/2)
  196. C centre de gravite du volume
  197. xb=xzero
  198. yb=xzero
  199. zb=xzero
  200. nbnn=ipt1.num(/1)
  201. DO i=1,nbnn
  202. ipt=idimp1*(ipt1.num(i,iel)-1)
  203. xb=xb+XCOOR(ipt+1)
  204. yb=yb+XCOOR(ipt+2)
  205. zb=zb+XCOOR(ipt+3)
  206. ENDDO
  207. xb=xb/nbnn
  208. yb=yb/nbnn
  209. zb=zb/nbnn
  210. DO ifac=1,nbfac
  211. v=xzero
  212. ityp=LDEL(1,iad+ifac)
  213. npfac=KDFAC(1,ityp)
  214. jad=LDEL(2,iad+ifac)-1
  215. idep=KDFAC(2,ityp)
  216. ifep=idep+3*(KDFAC(3,ityp)-1)
  217. DO i=idep,ifep,3
  218. iafa=LFAC(jad+KFAC(i))
  219. ibfa=LFAC(jad+KFAC(i+1))
  220. icfa=LFAC(jad+KFAC(i+2))
  221. ia=ipt1.num(iafa,iel)
  222. ib=ipt1.num(ibfa,iel)
  223. ic=ipt1.num(icfa,iel)
  224. v=v+VOLTRI(xb,yb,zb,ia,ib,ic)
  225. ENDDO
  226. * SG
  227. IF (ITYP.NE.7.AND.ITYP.NE.8) THEN
  228. DO i=idep,ifep,3
  229. ifa=KFAC(i)+kdeg
  230. IF (ifa.GT.npfac) ifa=ifa-npfac
  231. ifb=KFAC(i+1)+kdeg
  232. IF (ifb.GT.npfac) ifb=ifb-npfac
  233. ifc=KFAC(i+2)+kdeg
  234. IF (ifc.GT.npfac) ifc=ifc-npfac
  235. iafa=LFAC(jad+ifa)
  236. ibfa=LFAC(jad+ifb)
  237. icfa=LFAC(jad+ifc)
  238. ia=ipt1.num(iafa,iel)
  239. ib=ipt1.num(ibfa,iel)
  240. ic=ipt1.num(icfa,iel)
  241. v=v+VOLTRI(xb,yb,zb,ia,ib,ic)
  242. ENDDO
  243. ELSE
  244. v=v*2
  245. ENDIF
  246. IF (IIMPI.EQ.1) WRITE(6,*) 'v-2 ',v
  247. Cc=Cc+ABS(v)
  248. ENDDO
  249. ENDDO
  250. xmesu=xmesu+Cc/12.D0
  251. C * GOTO 500
  252. 500 SEGDES,ipt1
  253. ENDDO
  254.  
  255. 11 CONTINUE
  256. CALL ECRREE(xmesu)
  257. SEGDES,meleme
  258.  
  259. RETURN
  260. END
  261.  
  262.  
  263.  
  264.  
  265.  

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