1 : C REDCON SOURCE CHAT 12/01/17 21:15:53 7248 2 : subroutine redcon(mmodel,mmode2) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : -INC CCOPTIO 6 : -INC SMCOORD 7 : -INC SMELEME 8 : -INC CCREEL 9 : -INC SMMODEL 10 : character*4 modepl(3),moforc(3) 11 : data modepl /'UX ','UY ','UZ '/ 12 : data moforc /'FX ','FY ','FZ '/ 13 : * 14 : * lecture du maillage support des conditions 15 : * 16 : segact mmodel 17 : segini,mmode2=mmodel 18 : imoin=0 19 : do 1 imo=1,kmodel(/1) 20 : imodel=kmodel(imo) 21 : segact imodel 22 : do ifo=1,formod(/2) 23 : if( formod(ifo).eq.'CONTACT' ) go to 2 24 : enddo 25 : mmode2.kmodel(imo-imoin)=imodel 26 : go to 1 27 : 2 continue 28 : * on va reduire le maillage 29 : ipt1=imamod 30 : segact ipt1 31 : segini,ipt2=ipt1 32 : igard=0 33 : nbelem=ipt1.num(/2) 34 : if(idim.eq.3) then 35 : * 36 : * boucle sur le maillage support en 3D 37 : * 38 : 39 : do 10 iel=1,nbelem 40 : ip1=ipt1.num(2,iel) 41 : ip2=ipt1.num(3,iel) 42 : ip3=ipt1.num(4,iel) 43 : xp1=xcoor((ip1-1)*(idim+1)+1) 44 : yp1=xcoor((ip1-1)*(idim+1)+2) 45 : zp1=xcoor((ip1-1)*(idim+1)+3) 46 : xp2=xcoor((ip2-1)*(idim+1)+1) 47 : yp2=xcoor((ip2-1)*(idim+1)+2) 48 : zp2=xcoor((ip2-1)*(idim+1)+3) 49 : xp3=xcoor((ip3-1)*(idim+1)+1) 50 : yp3=xcoor((ip3-1)*(idim+1)+2) 51 : zp3=xcoor((ip3-1)*(idim+1)+3) 52 : * 53 : * verification si autre point dans la zone de selection 54 : jp=ipt1.num(5,iel) 55 : * write (6,*) 'mail ',ip1,ip2,ip3,jp 56 : xp =xcoor((jp-1)*(idim+1)+1) 57 : yp =xcoor((jp-1)*(idim+1)+2) 58 : zp =xcoor((jp-1)*(idim+1)+3) 59 : * verification que pas relation du point sur lui meme 60 : if (jp.eq.ip1) goto 10 61 : if (jp.eq.ip2) goto 10 62 : if (jp.eq.ip3) goto 10 63 : * verif distance au centre de gravite 64 : xg=(xp1+xp2+xp3)/3d0 65 : yg=(yp1+yp2+yp3)/3d0 66 : zg=(zp1+zp2+zp3)/3d0 67 : d1=sqrt((xg-xp1)**2+(yg-yp1)**2+(zg-zp1)**2) 68 : d2=sqrt((xg-xp2)**2+(yg-yp2)**2+(zg-zp2)**2) 69 : d3=sqrt((xg-xp3)**2+(yg-yp3)**2+(zg-zp3)**2) 70 : dp=sqrt((xg-xp)**2+(yg-yp)**2+(zg-zp)**2)/2d0 71 : if (dp.gt.d1.and.dp.gt.d2.and.dp.gt.d3) goto 10 72 : * write (6,*) 'contact test distance ok',dp,d1,d2,d3 73 : * verif position par rapport aux cotés 74 : * cote 1 2 75 : xv=xp2-xp1 76 : yv=yp2-yp1 77 : zv=zp2-zp1 78 : dv=sqrt((xv**2)+(yv**2)+(zv**2)) 79 : if(.not.(dv.lt.1d-10).and. 80 : > .not.(dv.gt.1d-10)) 81 : > write (6,*) 'prob 1 dans impo32' 82 : if (abs(dv).le.xpetit) write (6,*) 'prob 1 imp32' 83 : xv=xv/dv 84 : yv=yv/dv 85 : zv=zv/dv 86 : scal=(xp-xp1)*xv+(yp-yp1)*yv+(zp-zp1)*zv 87 : xm=xp1+scal*xv 88 : ym=yp1+scal*yv 89 : zm=zp1+scal*zv 90 : scal=(xg-xm)*xv+(yg-ym)*yv+(zg-zm)*zv 91 : xn=(xg-xm)-scal*xv 92 : yn=(yg-ym)-scal*yv 93 : zn=(zg-zm)-scal*zv 94 : dn=sqrt(xn**2+yn**2+zn**2) 95 : scal=(xp-xm)*xn+(yp-ym)*yn+(zp-zm)*zn 96 : dpm=sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2) 97 : if (dpm.gt.1d-3*dv.and.scal .lt.-0.2d0*dn*dpm) goto 10 98 : * write (6,*) 'contact test position 1 ok',scal/(dn*dpm),scal, 99 : * > dn,dpm 100 : * cote 2 3 101 : xv=xp3-xp2 102 : yv=yp3-yp2 103 : zv=zp3-zp2 104 : dv=sqrt((xv**2)+(yv**2)+(zv**2)) 105 : if(.not.(dv.lt.1d-10).and. 106 : > .not.(dv.gt.1d-10)) 107 : > write (6,*) 'prob 2 dans impo32' 108 : if (abs(dv).le.xpetit) write (6,*) 'prob 2 imp32' 109 : xv=xv/dv 110 : yv=yv/dv 111 : zv=zv/dv 112 : scal=(xp-xp2)*xv+(yp-yp2)*yv+(zp-zp2)*zv 113 : xm=xp2+scal*xv 114 : ym=yp2+scal*yv 115 : zm=zp2+scal*zv 116 : scal=(xg-xm)*xv+(yg-ym)*yv+(zg-zm)*zv 117 : xn=(xg-xm)-scal*xv 118 : yn=(yg-ym)-scal*yv 119 : zn=(zg-zm)-scal*zv 120 : dn=sqrt(xn**2+yn**2+zn**2) 121 : scal=(xp-xm)*xn+(yp-ym)*yn+(zp-zm)*zn 122 : dpm=sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2) 123 : if (dpm.gt.1d-3*dv.and.scal .lt.-0.2d0*dn*dpm) goto 10 124 : * write (6,*) 'contact test position 2 ok',scal/(dn*dpm),scal, 125 : * > dn,dpm 126 : * cote 3 1 127 : xv=xp1-xp3 128 : yv=yp1-yp3 129 : zv=zp1-zp3 130 : dv=sqrt((xv**2)+(yv**2)+(zv**2)) 131 : if(.not.(dv.lt.1d-10).and. 132 : > .not.(dv.gt.1d-10)) 133 : > write (6,*) 'prob 3 dans impo32' 134 : if (abs(dv).le.xpetit) write (6,*) 'prob 3 imp32' 135 : xv=xv/dv 136 : yv=yv/dv 137 : zv=zv/dv 138 : scal=(xp-xp3)*xv+(yp-yp3)*yv+(zp-zp3)*zv 139 : xm=xp3+scal*xv 140 : ym=yp3+scal*yv 141 : zm=zp3+scal*zv 142 : scal=(xg-xm)*xv+(yg-ym)*yv+(zg-zm)*zv 143 : xn=(xg-xm)-scal*xv 144 : yn=(yg-ym)-scal*yv 145 : zn=(zg-zm)-scal*zv 146 : dn=sqrt(xn**2+yn**2+zn**2) 147 : scal=(xp-xm)*xn+(yp-ym)*yn+(zp-zm)*zn 148 : dpm=sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2) 149 : if (dpm.gt.1d-3*dv.and.scal .lt.-0.2d0*dn*dpm) goto 10 150 : * write (6,*) 'contact test position 3 ok',scal/(dn*dpm),scal, 151 : * > dn,dpm 152 : * on a un point ou imposer la relation 153 : * Quelle est la relation ??? 154 : * write (6,*) 'contact potentiel ' 155 : * 156 : * direction de la relation 157 : * 158 : 159 : xn=(yp1-yp2)*(zp2-zp3)-(zp1-zp2)*(yp2-yp3) 160 : yn=(zp1-zp2)*(xp2-xp3)-(xp1-xp2)*(zp2-zp3) 161 : zn=(xp1-xp2)*(yp2-yp3)-(yp1-yp2)*(xp2-xp3) 162 : dn=sqrt(xn**2+yn**2+zn**2) 163 : if(.not.(dn.lt.1d-10).and. 164 : > .not.(dn.gt.1d-10)) 165 : > write (6,*) 'prob 4 dans impo32' 166 : if (abs(dn).le.xpetit) write (6,*) 'prob 4 imp32' 167 : xn=xn/dn 168 : yn=yn/dn 169 : zn=zn/dn 170 : if (abs(xn).lt.1d-10) xn=0.d0 171 : if (abs(yn).lt.1d-10) yn=0.d0 172 : if (abs(zn).lt.1d-10) zn=0.d0 173 : * 174 : * recherche du point a mettre en relation (coordonnees barycentriques) 175 : * 176 : scal=(xp-xg)*xn+(yp-yg)*yn+(zp-zg)*zn 177 : * 178 : * verif bon cote (a la tolerance pres) 179 : * 180 : if (scal.lt.-0.2d0*d1.and.scal.lt.-0.2d0*d2.and. 181 : > scal.lt.-0.2d0*d3) goto 10 182 : 183 : * 184 : * on range dans le meleme 185 : * 186 : 187 : igard=igard+1 188 : ipt2.num(1,igard)=ipt1.num(1,iel) 189 : ipt2.num(2,igard)=ip1 190 : ipt2.num(3,igard)=ip2 191 : ipt2.num(4,igard)=ip3 192 : ipt2.num(5,igard)=jp 193 : 10 continue 194 : segdes imodel,ipt1 195 : if(igard.eq.nbelem) then 196 : mmode2.kmodel(imo-imoin)=imodel 197 : segsup ipt2 198 : go to 1 199 : endif 200 : if( igard.eq.0) then 201 : imoin=imoin+1 202 : segsup ipt2 203 : go to 1 204 : endif 205 : nbelem=igard 206 : nbnn=ipt2.num(/1) 207 : nbref=0 208 : nbsous=0 209 : segadj ipt2 210 : segini,imode2=imodel 211 : imode2.imamod=ipt2 212 : kmodel(imo-imoin)=imode2 213 : segdes imode2,ipt2 214 : go to 1 215 : elseif(idim.eq.2) then 216 : if( ipt1.num(/1).eq.5) then 217 : * formulation faible 218 : do 210 iel=1,ipt1.num(/2) 219 : ip1=ipt1.num(2,iel) 220 : ip2=ipt1.num(3,iel) 221 : ip3=ipt1.num(4,iel) 222 : ip4=ipt1.num(5,iel) 223 : xp1=xcoor((ip1-1)*(idim+1)+1) 224 : yp1=xcoor((ip1-1)*(idim+1)+2) 225 : xp2=xcoor((ip2-1)*(idim+1)+1) 226 : yp2=xcoor((ip2-1)*(idim+1)+2) 227 : xp3=xcoor((ip3-1)*(idim+1)+1) 228 : yp3=xcoor((ip3-1)*(idim+1)+2) 229 : xp4=xcoor((ip4-1)*(idim+1)+1) 230 : yp4=xcoor((ip4-1)*(idim+1)+2) 231 : * 232 : * critere d'acceptation de l'éléments: 233 : * une extrémité d'un segment dans le cercle de sélection de l'autre segment 234 : * orientations respectives correctes des 2 segments 235 : * 236 : * verification (provisoire) que pas de noeuds doubles. 237 : if (ip1.eq.ip3) goto 210 238 : if (ip1.eq.ip4) goto 210 239 : if (ip2.eq.ip3) goto 210 240 : if (ip2.eq.ip4) goto 210 241 : * 242 : xl13=sqrt((xp3-xp1)**2+(yp3-yp1)**2) 243 : xl14=sqrt((xp4-xp1)**2+(yp4-yp1)**2) 244 : xl23=sqrt((xp3-xp2)**2+(yp3-yp2)**2) 245 : xl24=sqrt((xp4-xp2)**2+(yp4-yp2)**2) 246 : sca312=(xp3-xp1)*(xp3-xp2)+(yp3-yp1)*(yp3-yp2) 247 : sca412=(xp4-xp1)*(xp4-xp2)+(yp4-yp1)*(yp4-yp2) 248 : sca134=(xp1-xp3)*(xp1-xp4)+(yp1-yp3)*(yp1-yp4) 249 : sca234=(xp2-xp3)*(xp2-xp4)+(yp2-yp3)*(yp2-yp4) 250 : if (sca312/(xl13*xl23).gt.0.7.and. 251 : > sca412/(xl14*xl24).gt.0.7.and. 252 : > sca134/(xl13*xl14).gt.0.7.and. 253 : > sca234/(xl23*xl24).gt.0.7) goto 210 254 : * 255 : scal=(xp1-xp2)*(xp3-xp4)+(yp1-yp2)*(yp3-yp4) 256 : if (scal.gt.0.) goto 210 257 : 258 : * on a un element ou imposer la relation 259 : * Quelle est la relation ??? 260 : * 261 : * direction du contact 262 : * 263 : xr=(xp2-xp1+xp3-xp4) 264 : yr=(yp2-yp1+yp3-yp4) 265 : sr2 = xr**2 + yr**2 266 : sr = sqrt (sr2) 267 : xr=xr/sr 268 : yr=yr/sr 269 : * 270 : * projection des points sur la direction du contact 271 : * 272 : xl1=xp1*xr+yp1*yr 273 : xl2=xp2*xr+yp2*yr 274 : xl3=xp3*xr+yp3*yr 275 : xl4=xp4*xr+yp4*yr 276 : * 277 : * distance des points à leurs projections 278 : * 279 : d1=(xp1-xl1*xr)*yr-(yp1-xl1*yr)*xr 280 : d2=(xp2-xl2*xr)*yr-(yp2-xl2*yr)*xr 281 : d3=(xp3-xl3*xr)*yr-(yp3-xl3*yr)*xr 282 : d4=(xp4-xl4*xr)*yr-(yp4-xl4*yr)*xr 283 : * 284 : * calcul de l'intersection des projections. 285 : * 286 : xm1=min(xl1,xl2) 287 : xm2=max(xl1,xl2) 288 : xm3=min(xl3,xl4) 289 : xm4=max(xl3,xl4) 290 : * 291 : * critere 292 : * 293 : xcr=max(xm2-xm1,xm4-xm3)*1d-4 294 : xi=max(xm1,xm3) 295 : xj=min(xm2,xm4) 296 : * 297 : if (xi.ge.xj-xcr) goto 210 298 : 299 : * on range dans le meleme 300 : * 301 : igard=igard+1 302 : ipt2.num(1,igard)=ipt1.num(1,iel) 303 : ipt2.num(2,igard)=ip1 304 : ipt2.num(3,igard)=ip2 305 : ipt2.num(4,igard)=ip3 306 : ipt2.num(5,igard)=ip4 307 : * 308 : 210 continue 309 : segdes imodel,ipt1 310 : if(igard.eq.nbelem) then 311 : mmode2.kmodel(imo-imoin)=imodel 312 : segsup ipt2 313 : go to 1 314 : endif 315 : if( igard.eq.0) then 316 : imoin=imoin+1 317 : segsup ipt2 318 : go to 1 319 : endif 320 : nbelem=igard 321 : nbnn=ipt2.num(/1) 322 : nbref=0 323 : nbsous=0 324 : segadj ipt2 325 : segini,imode2=imodel 326 : imode2.imamod=ipt2 327 : kmodel(imo-imoin)=imode2 328 : segdes imode2,ipt2 329 : go to 1 330 : elseif(ipt1.num(/1).eq.4) then 331 : * formulation 1 point - 1 segment 332 : do 110 iel=1,ipt1.num(/2) 333 : ip1=ipt1.num(2,iel) 334 : ip2=ipt1.num(3,iel) 335 : xp1=xcoor((ip1-1)*(idim+1)+1) 336 : yp1=xcoor((ip1-1)*(idim+1)+2) 337 : xp2=xcoor((ip2-1)*(idim+1)+1) 338 : yp2=xcoor((ip2-1)*(idim+1)+2) 339 : * 340 : * verification si autre point dans le cercle de selection 341 : jp=ipt1.num(4,iel) 342 : * verification que pas relation du point sur lui meme 343 : if (jp.eq.ip1) goto 110 344 : if (jp.eq.ip2) goto 110 345 : xp=xcoor((jp-1)*(idim+1)+1) 346 : yp=xcoor((jp-1)*(idim+1)+2) 347 : scal=(xp-xp1)*(xp-xp2)+(yp-yp1)*(yp-yp2) 348 : xl1=sqrt((xp-xp1)**2+(yp-yp1)**2) 349 : xl2=sqrt((xp-xp2)**2+(yp-yp2)**2) 350 : scal=scal/(xl1*xl2) 351 : if (scal.gt.0.7) goto 110 352 : igard=igard+1 353 : ipt2.num(1,igard)=ipt1.num(1,iel) 354 : ipt2.num(2,igard)=ip1 355 : ipt2.num(3,igard)=ip2 356 : ipt2.num(4,igard)=jp 357 : 110 continue 358 : segdes imodel,ipt1 359 : if(igard.eq.nbelem) then 360 : mmode2.kmodel(imo-imoin)=imodel 361 : segsup ipt2 362 : go to 1 363 : endif 364 : if( igard.eq.0) then 365 : imoin=imoin+1 366 : segsup ipt2 367 : go to 1 368 : endif 369 : nbelem=igard 370 : nbnn=ipt2.num(/1) 371 : nbref=0 372 : nbsous=0 373 : segadj ipt2 374 : segini,imode2=imodel 375 : imode2.imamod=ipt2 376 : mmode2.kmodel(imo-imoin)=imode2 377 : segdes imode2,ipt2 378 : go to 1 379 : endif 380 : endif 381 : 1 continue 382 : segdes mmodel 383 : if(imoin.ne.0) then 384 : n1=kmodel(/1)-imoin 385 : if(n1.eq.0) then 385 : call erreur (21) 386 : return 387 : endif 388 : segadj mmode2 389 : endif 390 : segdes mmode2 391 : return 392 : end 393 :
© Cast3M 2003 - All rights reserved.
Disclaimer