1 : C MESUVO SOURCE CHAT 07/10/29 21:24:21 5940 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. 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 32 : if (ipt1.itypel.ne.ksurf(ipt1.itypel)) call erreur(25) 33 : if (ierr.ne.0) return 34 : endif 35 : nbel=nbel+ipt1.num(/2) 36 : do 15 ipt=1,ipt1.num(/1) 37 : do 15 iel=1,ipt1.num(/2) 38 : if (icpr(ipt1.num(ipt,iel)).ne.0) goto 15 39 : ite=ite+1 40 : icpr(ipt1.num(ipt,iel))=ite 41 : 15 continue 42 : 10 continue 43 : nbk=8 44 : C on remplit kon qui donne les elements touchant un noeud (compacte) 45 : segini kon,ikon 46 : do 20 isous=1,max(1,lisous(/1)) 47 : if (lisous(/1).ne.0) ipt1=lisous(isous) 48 : do 25 ipt=1,ipt1.num(/1) 49 : do 25 iel=1,ipt1.num(/2) 50 : ip=icpr(ipt1.num(ipt,iel)) 51 : ikon(ip)=ikon(ip)+1 52 : if (ikon(ip).gt.nbk) then 53 : nbk=nbk+1 54 : segadj kon 55 : endif 56 : kon(ip,ikon(ip),1)=ipt1 57 : kon(ip,ikon(ip),2)=iel 58 : 25 continue 59 : 20 continue 60 : C on commence a travailler 61 : segini ifut 62 : do 30 i=1,nbel 63 : ifut(i,3)=1 64 : 30 continue 65 : ilcou=0 66 : iltop=1 67 : ifut(1,1)=ipt1 68 : ifut(1,2)=1 69 : ifut(1,3)=1 70 : ipt1.num(1,1)=-ipt1.num(1,1) 71 : 100 continue 72 : ilcou=ilcou+1 73 : if (ilcou.gt.iltop) goto 500 74 : ipt2=ifut(ilcou,1) 75 : iel=ifut(ilcou,2) 76 : isens=ifut(ilcou,3) 77 : C on cherche les elements contactant iel et on les oriente correctement 78 : C on note qu'ils sont traites en mettant en neg le premier noeud 79 : do 110 i=1,ipt2.num(/1) 80 : ip1=abs(ipt2.num(i,iel)) 81 : i1=i+1 82 : if (i1.gt.ipt2.num(/1)) i1=1 83 : i0=i-1 84 : if (i0.lt.1) i0=ipt2.num(/1) 85 : if (isens.eq.1) ip2=abs(ipt2.num(i1,iel)) 86 : if (isens.eq.-1) ip2=abs(ipt2.num(i0,iel)) 87 : do 120 j=1,ikon(icpr(ip1)) 88 : ipt3=kon(icpr(ip1),j,1) 89 : jel=kon(icpr(ip1),j,2) 90 : jsens=isign(1,ipt3.num(2,jel)) 91 : * write (6,*) ' isens,jsens,iel,jel ',isens,jsens,iel,jel,ip1,ip2 92 : if (ipt3.eq.ipt2.and.jel.eq.iel) goto 120 93 : do 130 in=1,ipt3.num(/1) 94 : if (abs(ipt3.num(in,jel)).ne.ip1) goto 130 95 : in1=in+1 96 : if (in1.gt.ipt3.num(/1)) in1=1 97 : in0=in-1 98 : if (in0.lt.1) in0=ipt3.num(/1) 99 : if ((abs(ipt3.num(in1,jel)).eq.ip2.and.jsens.eq.1).or. 100 : * (abs(ipt3.num(in0,jel)).eq.ip2.and.jsens.eq.-1)) then 101 : C Il faut reorienter l'element 102 : if (iimpi.eq.5) write(6,*) ' element ',ipt3,jel,' change de sens' 103 : if (ipt3.num(1,jel).lt.0) goto 999 104 : if (ierr.ne.0) return 105 : ipt3.num(1,jel)=-abs(ipt3.num(1,jel)) 106 : iltop=iltop+1 107 : ifut(iltop,1)=ipt3 108 : ifut(iltop,2)=jel 109 : ifut(iltop,3)=-jsens 110 : ipt3.num(2,jel)=-abs(ipt3.num(2,jel)) 111 : goto 110 112 : elseif ((abs(ipt3.num(in1,jel)).eq.ip2.and.jsens.eq.-1) 113 : * .or.(abs(ipt3.num(in0,jel)).eq.ip2.and.jsens.eq.1)) then 114 : C Element OK 115 : if (iimpi.eq.5) write(6,*) ' element ',ipt3,jel, ' bon sens ' 116 : if (ipt3.num(1,jel).lt.0) goto 110 117 : ipt3.num(1,jel)=-abs(ipt3.num(1,jel)) 118 : iltop=iltop+1 119 : ifut(iltop,1)=ipt3 120 : ifut(iltop,2)=jel 121 : ifut(iltop,3)=jsens 122 : goto 110 123 : else 124 : endif 125 : 130 continue 126 : 120 continue 127 : goto 999 128 : 110 continue 129 : goto 100 130 : 500 continue 131 : C il ne reste plus qu'a calculer le volume et remettre les elements corrects 132 : if (iimpi.eq.5) write (6,*) ' iltop nbel ',iltop,nbel 133 : if (iltop.ne.nbel) goto 999 134 : if (ierr.ne.0) return 135 : do 600 i=1,nbel 136 : ipt1=ifut(i,1) 137 : iel=ifut(i,2) 138 : isens=ifut(i,3) 139 : IAD=LTEL(2,IPT1.ITYPEL)-1 140 : kdeg=kdegre(IPT1.ITYPEL)-1 141 : ITYP=LDEL(1,IAD+1) 142 : NPFAC=KDFAC(1,ITYP) 143 : JAD=LDEL(2,IAD+1)-1 144 : IDEP=KDFAC(2,ITYP) 145 : IFEP=IDEP+3*(KDFAC(3,ITYP)-1) 146 : DO 610 ITRIAN=IDEP,IFEP,3 147 : IAFA=LFAC(JAD+KFAC(ITRIAN)) 148 : IBFA=LFAC(JAD+KFAC(ITRIAN+1)) 149 : ICFA=LFAC(JAD+KFAC(ITRIAN+2)) 150 : IA=abs(IPT1.NUM(IAFA,IEL)) 151 : IB=abs(IPT1.NUM(IBFA,IEL)) 152 : IC=abs(IPT1.NUM(ICFA,IEL)) 153 : * if (iimpi.eq.5) write (6,*) ' triangle ',ia,ib,ic 154 : xmesu=xmesu+isens*voltri(0.d0,0.d0,0.d0,ia,ib,ic)/12 155 : 610 continue 156 : DO 611 ITRIAN=IDEP,IFEP,3 157 : ifa=KFAC(ITRIAN)+kdeg 158 : if (ifa.gt.npfac) ifa=ifa-npfac 159 : ifb=KFAC(ITRIAN+1)+kdeg 160 : if (ifb.gt.npfac) ifb=ifb-npfac 161 : ifc=KFAC(ITRIAN+2)+kdeg 162 : if (ifc.gt.npfac) ifc=ifc-npfac 163 : IAFA=LFAC(JAD+ifa) 164 : IBFA=LFAC(JAD+ifb) 165 : ICFA=LFAC(JAD+ifc) 166 : IA=abs(IPT1.NUM(IAFA,IEL)) 167 : IB=abs(IPT1.NUM(IBFA,IEL)) 168 : IC=abs(IPT1.NUM(ICFA,IEL)) 169 : xmesu=xmesu+isens*voltri(0.d0,0.d0,0.d0,ia,ib,ic)/12 170 : 611 continue 171 : 600 continue 172 : goto 998 173 : 999 continue 173 : call erreur(127) 174 : 998 continue 175 : C maintenant on desactive les segments et on remet les elements correctement 176 : ipt1=meleme 177 : do 620 isous=1,max(1,lisous(/1)) 178 : if (lisous(/1).ne.0) ipt1=lisous(isous) 179 : do 630 iel=1,ipt1.num(/2) 180 : ipt1.num(1,iel)=abs(ipt1.num(1,iel)) 181 : ipt1.num(2,iel)=abs(ipt1.num(2,iel)) 182 : 630 continue 183 : segdes ipt1 184 : 620 continue 185 : segdes meleme 186 : segsup icpr,kon,ikon,ifut 187 : xmesu=abs(xmesu) 188 : end 189 : 190 : 191 : 192 : 193 : 194 : 195 : 196 : 197 : 198 : 199 : 200 : 201 :
© Cast3M 2003 - All rights reserved.
Disclaimer