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