1 : C VOLSHB SOURCE CHAT 07/10/22 21:18:15 5921 2 : subroutine volshb(ipt1,epai,mchpoi,ip1,meleme) 3 : implicit real*8(a-h,o-z) 4 : implicit integer(i-n) 5 : -INC CCOPTIO 6 : -INC SMCHPOI 7 : -INC SMELEME 8 : -INC SMCOORD 9 : segment icpr(xcoor(/1)/(idim+1)) 10 : segment ipos(nbp) 11 : segment vec(3,nbp) 12 : segment icprp(xcoor(/1)/(idim+1)) 13 : segact ipt1 14 : if( ipt1.lisous(/1).ne.0.or.ipt1.itypel.ne.8) then 14 : call erreur (1002) 15 : return 16 : endif 17 : icle=1 18 : if(mchpoi.ne.0) then 19 : icle=2 20 : segact mchpoi 21 : if(ipchp(/1).ne.1) then 21 : call erreur(180) 22 : return 23 : endif 24 : msoupo=ipchp(1) 25 : segact msoupo 26 : if(nocomp(/2).ne.1) then 26 : call erreur(180) 27 : return 28 : endif 29 : mpoval=ipoval 30 : ipt2=igeoc 31 : segact ipt2,mpoval 32 : segini icprp 33 : do i=1,ipt2.num(/2) 34 : ia=ipt2.num(1,i) 35 : icprp(ia)=i 36 : enddo 37 : segdes ipt2 38 : endif 39 : * 40 : * pour chaque point on regarde quels sont les éléments qui le touchent 41 : * 42 : xp1=xcoor((ip1-1)*(idim+1)+1) 43 : yp1=xcoor((ip1-1)*(idim+1)+2) 44 : zp1=xcoor((ip1-1)*(idim+1)+3) 45 : nbp=0 46 : segini icpr 47 : do i=1,ipt1.num(/2) 48 : do j=1,4 49 : ia= ipt1.num(j,i) 50 : if(icpr(ia).eq.0) then 51 : nbp=nbp+1 52 : icpr(ia)=nbp 53 : endif 54 : enddo 55 : enddo 56 : segini ipos,vec 57 : * 58 : * on calcule pour chaque point le vecteur normale 59 : * 60 : isens=0 61 : do i=1,ipt1.num(/2) 62 : ia=(ipt1.num(1,i)-1)*4 63 : ib=(ipt1.num(2,i)-1)*4 64 : ic=(ipt1.num(3,i)-1)*4 65 : xa=xcoor(ia+1) 66 : ya=xcoor(ia+2) 67 : za=xcoor(ia+3) 68 : xb=xcoor(ib+1) 69 : yb=xcoor(ib+2) 70 : zb=xcoor(ib+3) 71 : xc=xcoor(ic+1) 72 : yc=xcoor(ic+2) 73 : zc=xcoor(ic+3) 74 : xn=(yb-ya)*(zc-za)-(yc-ya)*(zb-za) 75 : yn=(zb-za)*(xc-xa)-(zc-za)*(xb-xa) 76 : zn=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya) 77 : xyz=sqrt(xn*xn+yn*yn+zn*zn) 78 : xn=xn/xyz 79 : yn=yn/xyz 80 : zn=zn/xyz 81 : do j=1,4 82 : ia=ipt1.num(j,i) 83 : ib=icpr(ia) 84 : ipos(ib)=ipos(ib)+1 85 : vec(1,ib)=vec(1,ib)+xn 86 : vec(2,ib)=vec(2,ib)+yn 87 : vec(3,ib)=vec(3,ib)+zn 88 : enddo 89 : enddo 90 : do i=1,nbp 91 : na= ipos(i) 92 : ipos(i)=0 93 : vec(1,i)=vec(1,i)/na 94 : vec(2,i)=vec(2,i)/na 95 : vec(3,i)=vec(3,i)/na 96 : xv= vec(1,i)*vec(1,i)+vec(2,i)*vec(2,i)+vec(3,i)*vec(3,i) 97 : xv=sqrt(xv) 98 : vec(1,i)=vec(1,i)/xv 99 : vec(2,i)=vec(2,i)/xv 100 : vec(3,i)=vec(3,i)/xv 101 : xps= vec(1,i)*xp1+vec(2,i)*yp1+vec(3,i)*zp1 102 : if(xps.ge.0.d0) then 103 : if( isens.lt.0) then 103 : call erreur (19) 104 : return 105 : endif 106 : isens=1 107 : elseif(xps.le.0.d0) then 108 : if( isens.gt.0) then 108 : call erreur (19) 109 : return 110 : endif 111 : isens=-1 112 : endif 113 : enddo 114 : * 115 : * il faut maintenant créer le maillage 116 : * 117 : nbnn=8 118 : nbelem=ipt1.num(/2) 119 : nbref=2 120 : nbsous=0 121 : segini meleme 122 : itypel=14 123 : ep=epai 124 : nbpts=xcoor(/1)/(idim+1) 125 : nbpini=nbpts 126 : nbpts=nbpts+2*nbp 127 : segadj mcoord 128 : nbnn=4 129 : nbref=0 130 : segini ipt5,ipt6 131 : lisref(1)=ipt5 132 : lisref(2)=ipt6 133 : ipt5.itypel=8 134 : ipt6.itypel=8 135 : do i=1,ipt1.num(/2) 136 : do j=1,4 137 : ia=ipt1.num(j,i) 138 : ib=icpr(ia) 139 : * faut-il créer les points ? 140 : if( ipos(ib).eq.0) then 141 : iadec=(ia-1)*(idim+1) 142 : if(icle.eq.2) then 143 : ic=icprp(ia) 144 : if(ic.eq.0) then 145 : meleme=0 145 : call erreur(19) 146 : return 147 : endif 148 : ep=vpocha(ic,1) 149 : endif 150 : nbpini=nbpini+1 151 : nbpdec=(nbpini-1)*(idim+1) 152 : ie=nbpini 153 : eps=ep*isens 154 : do k=1,3 155 : xcoor(nbpdec+k)= xcoor(iadec+k)-vec(k,ib)*eps 156 : enddo 157 : nbpini=nbpini+1 158 : nbpdec=nbpdec+4 159 : do k=1,3 160 : xcoor(nbpdec+k)= xcoor(iadec+k)+vec(k,ib)*eps 161 : enddo 162 : ipos(ib)=ie 163 : else 164 : ie=ipos(ib) 165 : endif 166 : num(J,i)=ie 167 : num(J+4,i)=ie+1 168 : ipt5.num(j,i)=ie 169 : ipt6.num(j,i)=ie+1 170 : enddo 171 : enddo 172 : if( icle.eq.2)then 173 : segdes mpoval,ipt2,msoupo,mchpoi 174 : segsup icprp 175 : endif 176 : segsup icpr,ipos,vec 177 : segdes ipt1 178 : segdes meleme,ipt5,ipt6 179 : return 180 : end 181 : 182 : 183 : 184 : 185 : 186 : 187 : 188 :
© Cast3M 2003 - All rights reserved.
Disclaimer