1 : C RAVC SOURCE PV 13/12/06 21:15:01 7877 2 : subroutine ravc 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8(A-H,O-Z) 5 : -INC CCOPTIO 6 : -INC SMCHPOI 7 : -INC SMELEME 8 : -INC SMCOORD 9 : segment icpr(xcoor(/1)/(idim+1)) 9 : call lirobj ( 'CHPOINT',MCHPO1, 1,iretou) 10 : if(ierr . ne.0) return 10 : call lirobj ( 'CHPOINT',MCHPO2, 1,iretou) 11 : if(ierr . ne.0) return 11 : call lirobj ( 'CHPOINT',MCHPO3, 1,iretou) 12 : if(ierr . ne.0) return 13 : segact mchpo1, mchpo2,mchpo3 14 : if(mchpo1.ipchp(/1). ne . 1) then 15 : * write(6,*) ' erreur 1 ' 15 : call erreur(19) 16 : return 17 : endif 18 : if(mchpo2.ipchp(/1). ne . 1) then 19 : * write(6,*) ' erreur 2 ' 19 : call erreur(19) 20 : return 21 : endif 22 : if(mchpo3.ipchp(/1). ne . 1) then 23 : * write(6,*) ' erreur 3 ' 23 : call erreur(19) 24 : return 25 : endif 26 : msoup1= mchpo1.ipchp(1) 27 : msoup2= mchpo2.ipchp(1) 28 : msoup3= mchpo3.ipchp(1) 29 : segact msoup1,msoup2,msoup3 30 : ipt1=msoup1.igeoc 31 : ipt2=msoup2.igeoc 32 : ipt3=msoup3.igeoc 33 : segact ipt1,ipt2,ipt3 34 : segini icpr 35 : n1 = ipt1.num(/2) 36 : n2 = ipt2.num(/2) 37 : n3 = ipt3.num(/2) 38 : ipo=0 39 : do 100 i=1,n1 40 : ia=ipt1.num(1,i) 41 : if(icpr(ia).eq.0) then 42 : ipo=ipo+1 43 : icpr(ia)=ipo 44 : endif 45 : 100 continue 46 : do 101 i=1,n2 47 : ia=ipt2.num(1,i) 48 : if(icpr(ia).eq.0) then 49 : ipo=ipo+1 50 : icpr(ia)=ipo 51 : endif 52 : 101 continue 53 : do 102 i=1,n3 54 : ia=ipt3.num(1,i) 55 : if(icpr(ia).eq.0) then 56 : ipo=ipo+1 57 : icpr(ia)=ipo 58 : endif 59 : 102 continue 60 : nbelem=ipo 61 : nbnn=1 62 : nbsous=0 63 : nbref=0 64 : segini meleme 65 : itypel=1 66 : do 103 i=1,icpr(/1) 67 : if(icpr(I).ne.0) then 68 : num(1,icpr(i))=i 69 : endif 70 : 103 continue 71 : 72 : * write(6,*) ' ipt1',(ipt1.num(1,iou),iou=1,ipt1.num(/2)) 73 : * write(6,*) ' ipt2',(ipt2.num(1,iou),iou=1,ipt2.num(/2)) 74 : * write(6,*) ' ipt3',(ipt3.num(1,iou),iou=1,ipt3.num(/2)) 75 : * if( (n1.ne.n2). or .(n1.ne.n3))then 76 : * write(6,*) ' erreur 4 ' 76 : * call erreur(19) 77 : * return 78 : * endif 79 : if( msoup1.noharm(/1) . ne . 1) then 80 : * write(6,*) ' erreur 5 ' 80 : call erreur(19) 81 : return 82 : endif 83 : if( msoup2.noharm(/1) . ne . 1) then 84 : * write(6,*) ' erreur 6 ' 84 : call erreur(19) 85 : return 86 : endif 87 : if( msoup3.noharm(/1) . ne . 1) then 88 : * write(6,*) ' erreur 7 ' 88 : call erreur(19) 89 : return 90 : endif 91 : segini,mchpoi=mchpo1 92 : segini,msoupo=msoup1 93 : ipchp(1)=msoupo 94 : mpova1=msoup1.ipoval 95 : mpova2=msoup2.ipoval 96 : mpova3=msoup3.ipoval 97 : segact mpova1,mpova2,mpova3 98 : n=nbelem 99 : nc=1 100 : segini mpoval 101 : ipoval=mpoval 102 : igeoc=meleme 103 : * write(6,*) ' meleme' , ( num(1,iou),iou=1,nbelem) 104 : * 105 : * on cherche les max de chaque champs 106 : * 107 : * amma=0. 108 : * bmma=0. 109 : * cmma=0. 110 : * do 10 i=1,nbelem 111 : * if(amma.le.(abs(mpova1.vpocha(i,1))))amma=abs(mpova1.vpocha(i,1)) 112 : * if(bmma.le.(abs(mpova2.vpocha(i,1))))bmma=abs(mpova2.vpocha(i,1)) 113 : * if(cmma.le.(abs(mpova3.vpocha(i,1))))cmma=abs(mpova3.vpocha(i,1)) 114 : * 10 continue 115 : * 116 : * calcul des racines 117 : * 118 : do 1 j=1,nbelem 119 : ia = num(1,j) 120 : * write(6,*) ' ia ' , ia 121 : a = 0.d0 122 : do 8 i=1,n1 123 : if( ipt1.num(1,i).eq.ia)then 124 : a = mpova1.vpocha(i,1) 125 : go to 9 126 : endif 127 : 8 continue 128 : 9 continue 129 : b=0.d0 130 : do 2 k=1,n2 131 : if(ipt2.num(1,k).eq.ia) then 132 : b = mpova2.vpocha(k,1) 133 : go to 3 134 : endif 135 : 2 continue 136 : 3 continue 137 : c=0.d0 138 : do 4 l=1,n3 139 : if(ipt3.num(1,l).eq.ia)then 140 : c = mpova3.vpocha(l,1) 141 : go to 5 142 : endif 143 : 4 continue 144 : 5 continue 145 : 146 : * test sur les grandeurs 147 : * write (6,*) ' a b c ',a,b,c 148 : if( a . eq. 0.d0. or . abs(a).lt. (abs(b)*1.e-13))then 149 : * on resoud une seule equation 150 : if( b . ne . 0.d0) then 151 : xx = -c / b 152 : if( xx .ge. 0.d0 ) xx = 0.d0 153 : else 154 : xx = 0.d0 155 : endif 156 : else 157 : xx=0.d0 158 : delta2 = ( b*b - 4*a*c) 159 : if (delta2.lt. b*b*1d-13) then 160 : xx = -b/(2.d0*a) 161 : xx = min(xx,0.d0) 162 : xx = 0.d0 163 : ** write (6,*) ' ravc discriminant negatif ',delta2,b 164 : else 165 : delta2=max(0.d0,delta2) 166 : delta = sqrt (delta2) 167 : s1= (-b + delta )/ (2.*a) 168 : s2= (-b - delta) / (2.*a) 169 : s1=min(s1,0.d0) 170 : s2=min(s2,0.d0) 171 : * write (6,*) ' s1 s2 ' , s1,s2 172 : if( s1. gt . s2) then 173 : if( s1 . lt. 0.d0) then 174 : xx = s1 175 : else 176 : if(s2.lt.0.d0) xx = s2 177 : endif 178 : else 179 : if( s2 . lt. 0.d0) then 180 : xx = s2 181 : else 182 : if(s1.lt.0.d0) xx = s1 183 : endif 184 : endif 185 : endif 186 : endif 187 : ** write (6,*) ' xx ',xx 188 : vpocha(j,1)=xx 189 : 1 continue 190 : segdes mpova1,mpova2,mpova3,mpoval 191 : segdes ipt1,ipt2,ipt3,meleme 192 : segdes msoupo,msoup1,msoup2,msoup3 193 : segdes mchpoi,mchpo1,mchpo2,mchpo3 194 : segsup icpr 194 : call ecrobj('CHPOINT',mchpoi) 195 : return 196 : end 197 : 198 : 199 : 200 : 201 : 202 : 203 : 204 : 205 : 206 :
© Cast3M 2003 - All rights reserved.
Disclaimer