1 : C DEDU SOURCE CHAT 11/03/16 21:18:50 6902 2 : 3 : C----------------------------------------------------------------------- 4 : C 1ere possibilite 5 : C OBJ1 = DEDU OBJ2 MAIT_ANC MAIT_NOU (REGU) 6 : C Deduction d'un maillage a partir d'un autre et du deplacement de 7 : C noeuds maitres 8 : C AUTEUR : PV 9 : C----------------------------------------------------------------------- 10 : C 2eme possibilite 11 : C OBJ2 = OBJ1 DEDU 'TRAN' GEO1 GEO2 12 : C 3eme possibilite 13 : C OBJ2 = OBJ1 DEDU FLOT1 POIN1 (POIN2) 'ROTA' GEO1 GEO2 14 : C OBJ1 : type POINT, MAILLAGE, CHPOINT, MCHAML, MMODEL 15 : C Passage du support de OBJ1 a celui de OBJ2 comme celui de GEO1 a GEO2 16 : C Les points du support de OBJ1 doivent appartenir a GEO1 17 : C On realise eventuellement la rotation des composantes (envoi a PROPER) 18 : C AUTEUR : KICH 19 : C----------------------------------------------------------------------- 20 : C 10/2003 : Prise en compte du cas IDIM=1 et des modes associes. 21 : C----------------------------------------------------------------------- 22 : C Remarques : 23 : C ----------- 24 : C DEDU est appele par l'operateur DEPL (deplac.eso) avec ICAS=1. 25 : C Ce cas correspond a une transformation affine de l'espace (noeuds 26 : C maitres = un maillage a (IDIM+1) points). 27 : C----------------------------------------------------------------------- 28 : 29 : SUBROUTINE DEDU(ICAS) 30 : 31 : IMPLICIT INTEGER(I-N) 32 : IMPLICIT REAL*8 (A-H,O-Z) 33 : 34 : -INC CCOPTIO 35 : -INC SMCOORD 36 : -INC SMELEME 37 : -INC CCGEOME 38 : -INC SMRIGID 39 : -INC SMCHPOI 40 : -INC SMMODEL 41 : -INC SMCHAML 42 : -INC SMTEXTE 43 : DIMENSION iuniq(10) 44 : CHARACTER*4 mcle(5),strcomp 45 : 46 : SEGMENT kon(nkon,numnp,2) 47 : SEGMENT xkon(nkon+1,numnp) 48 : SEGMENT ipxkon(nkon+1,numnp) 49 : SEGMENT ikon(numnp) 50 : SEGMENT icpr(xcoor(/1)/(IDIM+1)) 51 : SEGMENT icp2(xcoor(/1)/(IDIM+1)) 52 : SEGMENT idcp(numnp) 53 : SEGMENT icorr(numnp) 54 : SEGMENT xxtra 55 : REAL*8 xdep(3,4),xarr(3,4),ved1(3),ved2(3),ved3(3),vea1(3), 56 : . vea2(3),vea3(3),pa(3,3),pa2(2,2),pa3(3,3),vb(3) 57 : END SEGMENT 58 : 59 : DATA mcle / 'REGU','ELIM','TRAN','ROTA','ADAP' / 60 : 61 : idimp1=IDIM+1 62 : 63 : C Mot cle optionnel REGU : coordonnees barycentriques identiques 63 : CALL LIRMOT(mcle,5,iregu,0) 64 : 65 : C 2eme possibilite : option TRAN 66 : C -------------------- 67 : IF (iregu.EQ.3) THEN 67 : CALL PROPER(4) 68 : RETURN 69 : ENDIF 70 : C 3eme possibilite : option ROTA 71 : C -------------------- 72 : C Option non disponible en dimension 1 (pas de signification) 73 : IF (iregu.EQ.4) THEN 74 : IF (IDIM.EQ.1) THEN 75 : MOTERR(1:4)=mcle(iregu) 76 : INTERR(1)=IDIM 76 : CALL ERREUR(971) 77 : ELSE 77 : CALL PROPER(5) 78 : ENDIF 79 : RETURN 80 : ENDIF 81 : C 5eme possibilite : option ADAP 82 : C On passe la main à la procédure DEDUADAP 83 : C -------------------- 84 : IF (iregu.EQ.5) THEN 85 : SEGINI MTEXTE 86 : LTT=8 87 : MTEXT(1:LTT) ='DEDUADAP' 88 : NCART=8 89 : SEGDES MTEXTE 89 : CALL ECROBJ('TEXTE',MTEXTE) 90 : RETURN 91 : ENDIF 92 : 93 : C Par defaut : option ELIM (ignore les parties concaves) 94 : C Correspond au cas ICAS=1 (operateur DEPL 'DEDU') 95 : IF (iregu.EQ.0) iregu=2 96 : 97 : C 1ere possibilite : options REGU et ELIM / appel par DEPLAC 98 : C -------------------- 99 : C Lecture du maillage initial 99 : CALL LIROBJ('MAILLAGE',ipt1,1,iretou) 100 : IF (ierr.NE.0) RETURN 101 : C Lecture des points maitres initiaux 101 : CALL LIROBJ('MAILLAGE',ipt2,0,iretou) 102 : IF (iretou.NE.0) THEN 103 : ichpo=0 104 : C lecture ddes points maitres finaux 104 : CALL LIROBJ('MAILLAGE',ipt3,1,iretou) 105 : IF (ierr.NE.0) RETURN 106 : ELSE 107 : C Lecture d'un chpoint de deplacement 108 : ichpo=1 108 : CALL LIROBJ('CHPOINT ',mchpoi,1,iretou) 109 : IF (ierr.NE.0) RETURN 110 : C Création de ipt2 et ipt3 pour se ramener au cas précédent 111 : nbpts=xcoor(/1)/idimp1 112 : nbpts0=nbpts 113 : C Comptage des points du chpoint concernes par UX UY UZ UR 114 : SEGACT mchpoi 115 : jatan=jattri(1) 116 : nbpcon=0 117 : DO nso=1,ipchp(/1) 118 : msoupo=ipchp(nso) 119 : SEGACT msoupo 120 : ipt6=igeoc 121 : SEGACT ipt6 122 : ipc=0 123 : DO i=1,nocomp(/2) 124 : strcomp=nocomp(i) 125 : IF ((strcomp.EQ.'UX ').OR.(strcomp.EQ.'UR ').OR. 126 : . ((strcomp.EQ.'UY ').AND.(IDIM.NE.1)).OR. 127 : . ((strcomp.EQ.'UZ ').AND.(IDIM.NE.1))) ipc=ipc+i 128 : ENDDO 129 : IF (ipc.NE.0) nbpcon=nbpcon+ipt6.num(/2) 130 : ENDDO 131 : C Test sur nbpcon=0 a ecrire 132 : nbref=0 133 : nbsous=0 134 : nbelem=nbpcon 135 : nbnn=1 136 : SEGINI ipt2 137 : SEGINI ipt3 138 : ipt2.itypel=1 139 : ipt3.itypel=1 140 : C Remplissage de ipt2 et ipt3 - Creation des points de ipt3 141 : idej=0 142 : DO 9 nso=1,ipchp(/1) 143 : msoupo=ipchp(nso) 144 : ipt6=igeoc 145 : ipox=0 146 : ipoy=0 147 : ipoz=0 148 : DO i=1,nocomp(/2) 149 : strcomp=nocomp(i) 150 : IF ((strcomp.EQ.'UX ').OR.(strcomp.EQ.'UR ')) ipox=i 151 : IF ((strcomp.EQ.'UY ').AND.(IDIM.NE.1)) ipoy=i 152 : IF (strcomp.EQ.'UZ ') THEN 153 : IF (IDIM.EQ.2) THEN 154 : ipoy=i 155 : ELSE IF (IDIM.GE.3) THEN 156 : ipoz=i 157 : ENDIF 158 : ENDIF 159 : ENDDO 160 : IF ((ipox+ipoy+ipoz).EQ.0) GOTO 9 161 : mpoval=ipoval 162 : SEGACT mpoval 163 : nbpan=xcoor(/1)/idimp1 164 : nbpts=nbpts+ipt6.num(/2) 165 : SEGADJ mcoord 166 : DO nel=1,ipt6.num(/2) 167 : iep=ipt6.num(1,nel) 168 : ipt2.num(1,idej+nel)=iep 169 : ipt3.num(1,idej+nel)=nbpan+nel 170 : irefep=(iep-1)*idimp1 171 : iref=(nbpan+nel-1)*idimp1 172 : IF (ipox.NE.0) THEN 173 : xcoor(iref+1)=xcoor(irefep+1)+vpocha(nel,ipox) 174 : ELSE 175 : xcoor(iref+1)=xcoor(irefep+1) 176 : ENDIF 177 : IF (ipoy.NE.0) THEN 178 : xcoor(iref+2)=xcoor(irefep+2)+vpocha(nel,ipoy) 179 : ELSE 180 : IF (IDIM.GE.2) xcoor(iref+2)=xcoor(irefep+2) 181 : ENDIF 182 : IF (ipoz.NE.0) THEN 183 : xcoor(iref+3)=xcoor(irefep+3)+vpocha(nel,ipoz) 184 : ELSE 185 : IF (IDIM.GE.3) xcoor(iref+3)=xcoor(irefep+3) 186 : ENDIF 187 : xcoor(iref+idimp1)=xcoor(irefep+idimp1) 188 : ENDDO 189 : idej=idej+ipt6.num(/2) 190 : SEGDES ipt6,mpoval,msoupo 191 : 9 CONTINUE 192 : SEGDES mchpoi 193 : ENDIF 194 : 195 : C Verification compatibilite points maitres initiaux et finaux 196 : C Cas ichpo = 0 : cette verification est a faire lorsque la 197 : C transformation est definie par deux maillages. 198 : C Cas ichpo = 1 : les points maitres initiaux et finaux sont crees 199 : C en parallele de la meme maniere a partir du CHPO donne. La veri- 200 : C fication est donc inutile. 201 : IF (ichpo.EQ.1) GOTO 100 202 : SEGACT ipt2,ipt3 203 : IF (ipt2.lisous(/1).NE.ipt3.lisous(/1)) THEN 203 : CALL ERREUR(16) 204 : SEGDES ipt2,ipt3 205 : RETURN 206 : ENDIF 207 : ipt7=ipt2 208 : ipt8=ipt3 209 : DO i=1,MAX(1,ipt2.lisous(/1)) 210 : IF (ipt2.lisous(/1).NE.0) ipt7=ipt2.lisous(i) 211 : IF (ipt3.lisous(/1).NE.0) ipt8=ipt3.lisous(i) 212 : SEGACT ipt7,ipt8 213 : moterr=' ' 214 : IF (ipt7.num(/1).NE.ipt8.num(/1)) moterr(1:8)='noeuds' 215 : IF (ipt7.num(/2).NE.ipt8.num(/2)) moterr(1:8)='elements' 216 : SEGDES ipt7,ipt8 217 : IF (moterr.NE.' ') THEN 217 : CALL ERREUR(550) 218 : RETURN 219 : ENDIF 220 : ENDDO 221 : 222 : 100 SEGACT MCOORD*MOD 223 : C Creation du tableau de correspondance icpr donnant pour chaque 224 : C noeud de ipt2 son noeud image de ipt3 225 : SEGACT ipt2,ipt3 226 : SEGINI icpr 227 : numnp=0 228 : ipt7=ipt2 229 : ipt8=ipt3 230 : DO i=1,MAX(1,ipt2.lisous(/1)) 231 : IF (ipt2.lisous(/1).NE.0) ipt7=ipt2.lisous(i) 232 : IF (ipt3.lisous(/1).NE.0) ipt8=ipt3.lisous(i) 233 : SEGACT ipt7,ipt8 234 : DO k=1,ipt7.num(/2) 235 : DO l=1,ipt7.num(/1) 236 : ip=ipt7.num(l,k) 237 : IF (icpr(ip).EQ.0) THEN 238 : numnp=numnp+1 239 : icpr(ip)=ipt8.num(l,k) 240 : ENDIF 241 : ENDDO 242 : ENDDO 243 : SEGDES ipt7,ipt8 244 : ENDDO 245 : SEGDES ipt2,ipt3 246 : 247 : C Est-ce une transformation affine ? 248 : IF (numnp.LE.idimp1) THEN 249 : SEGINI xxtra 250 : numnp=0 251 : DO i=1,icpr(/1) 252 : IF (icpr(i).NE.0) THEN 253 : numnp=numnp+1 254 : iref=(i-1)*idimp1 255 : irefep=(icpr(i)-1)*idimp1 256 : DO j=1,idim 257 : xdep(j,numnp)=xcoor(iref+j) 258 : xarr(j,numnp)=xcoor(irefep+j) 259 : ENDDO 260 : ENDIF 261 : ENDDO 262 : SEGSUP icpr 263 : C Faut-il creer de nouveaux points ? 264 : IF (numnp.EQ.1) THEN 265 : C On en cree un par translation en x 266 : numnp=numnp+1 267 : DO j=1,IDIM 268 : xdep(j,numnp)=xdep(j,numnp-1) 269 : xarr(j,numnp)=xarr(j,numnp-1) 270 : ENDDO 271 : xdep(1,numnp)=xdep(1,numnp)+1. 272 : xarr(1,numnp)=xarr(1,numnp)+1. 273 : ENDIF 274 : DO j=1,IDIM 275 : ved1(j)=xdep(j,2)-xdep(j,1) 276 : vea1(j)=xarr(j,2)-xarr(j,1) 277 : ENDDO 278 : IF (numnp.EQ.2) THEN 279 : C En 1D, on ne fait rien. 280 : C En 2D, on cherche un point sur la perpendiculaire a la droite 281 : C definie par les deux premiers points, passant par le premier 282 : C point et a une distance unite (formant un repere orthonorme direct). 283 : C En 3D, on prend un point dans la direction y. 284 : numnp=numnp+1 285 : IF (IDIM.EQ.2) THEN 286 : xdep(1,3)=ved1(2) 287 : xdep(2,3)=-ved1(1) 288 : aa=1./SQRT(xdep(1,3)*xdep(1,3)+xdep(2,3)*xdep(2,3)) 289 : xdep(1,3)=xdep(1,1)+xdep(1,3)*aa 290 : xdep(2,3)=xdep(2,1)+xdep(2,3)*aa 291 : xarr(1,3)=vea1(2) 292 : xarr(2,3)=-vea1(2) 293 : aa=1./SQRT(xarr(1,3)*xarr(1,3)+xarr(2,3)*xarr(2,3)) 294 : xarr(1,3)=xarr(1,1)+xarr(1,3)*aa 295 : xarr(2,3)=xarr(2,1)+xarr(2,3)*aa 296 : ELSE IF (IDIM.EQ.3) THEN 297 : xdep(1,3)=xdep(1,1) 298 : xdep(2,3)=xdep(2,1)+1. 299 : xdep(3,3)=xdep(3,1) 300 : xarr(1,3)=xarr(1,1) 301 : xarr(2,3)=xarr(2,1)+1. 302 : xarr(3,3)=xarr(3,1) 303 : ENDIF 304 : ENDIF 305 : IF (IDIM.NE.1) THEN 306 : DO j=1,IDIM 307 : ved2(j)=xdep(j,3)-xdep(j,1) 308 : vea2(j)=xarr(j,3)-xarr(j,1) 309 : ENDDO 310 : ENDIF 311 : C En 3D, il faut encore un point hors du plan forme par les 3 points 312 : IF (numnp.LT.idimp1) THEN 313 : ved3(1)=ved1(2)*ved2(3)-ved1(3)*ved2(2) 314 : ved3(2)=ved1(3)*ved2(1)-ved1(1)*ved2(3) 315 : ved3(3)=ved1(1)*ved2(2)-ved1(2)*ved2(1) 316 : aa=1./SQRT(ved3(1)*ved3(1)+ved3(2)*ved3(2)+ved3(3)*ved3(3)) 317 : vea3(1)=vea1(2)*vea2(3)-vea1(3)*vea2(2) 318 : vea3(2)=vea1(3)*vea2(1)-vea1(1)*vea2(3) 319 : vea3(3)=vea1(1)*vea2(2)-vea1(2)*vea2(1) 320 : bb=1./SQRT(vea3(1)*vea3(1)+vea3(2)*vea3(2)+vea3(3)*vea3(3)) 321 : DO j=1,3 322 : ved3(j)=ved3(j)*aa 323 : vea3(j)=vea3(j)*bb 324 : xdep(j,4)=xdep(j,1)+ved3(j) 325 : xarr(j,4)=xarr(j,1)+vea3(j) 326 : ENDDO 327 : ENDIF 328 : IF (IDIM.EQ.3) THEN 329 : DO j=1,3 330 : ved3(j)=xdep(j,4)-xdep(j,1) 331 : vea3(j)=xarr(j,4)-xarr(j,1) 332 : ENDDO 333 : ENDIF 334 : C Calcul de la matrice de passage pa du repere ved au repere vea 335 : C On traite le cas IDIM=1 a part (sinon reso33 serait a modifier). 336 : IF (IDIM.EQ.1) THEN 337 : IF ((ved1(1).EQ.0.).OR.(vea1(1).EQ.0.)) THEN 337 : CALL ERREUR(21) 338 : SEGSUP xxtra 339 : IF (ichpo.EQ.1) THEN 340 : nbpts=nbpts0 341 : SEGADJ mcoord 342 : SEGSUP ipt2,ipt3 343 : ENDIF 344 : RETURN 345 : ENDIF 346 : pa(1,1)=vea1(1)/ved1(1) 347 : ELSE 348 : DO j=1,IDIM 349 : vb(1)=vea1(j) 350 : vb(2)=vea2(j) 351 : IF (IDIM.EQ.3) vb(3)=vea3(j) 352 : pa3(1,1)=ved1(1) 353 : pa3(1,2)=ved1(2) 354 : pa3(2,1)=ved2(1) 355 : pa3(2,2)=ved2(2) 356 : IF (IDIM.EQ.3) THEN 357 : pa3(1,3)=ved1(3) 358 : pa3(2,3)=ved2(3) 359 : pa3(3,1)=ved3(1) 360 : pa3(3,2)=ved3(2) 361 : pa3(3,3)=ved3(3) 362 : ENDIF 362 : CALL reso33(pa3,vb,3,IDIM,kerr) 363 : C Erreur : en 2D -> 3 points alignes, en 3D -> 4 points coplanaires 364 : IF (kerr.NE.0) THEN 364 : CALL ERREUR(21) 365 : SEGSUP xxtra 366 : IF (ichpo.EQ.1) THEN 367 : nbpts=nbpts0 368 : SEGADJ mcoord 369 : SEGSUP ipt2,ipt3 370 : ENDIF 371 : RETURN 372 : ENDIF 373 : pa(j,1)=vb(1) 374 : pa(j,2)=vb(2) 375 : IF (IDIM.EQ.3) pa(j,3)=vb(3) 376 : ENDDO 377 : ENDIF 378 : C Fabrication de la liste des points a transformer et du maillage 379 : C resultat (segment MELEME) 380 : SEGINI icpr 381 : nbpan=xcoor(/1)/idimp1 382 : SEGACT ipt1 383 : SEGINI,ipt4=ipt1 384 : nbnn=ipt4.num(/1) 385 : nbelem=ipt4.num(/2) 386 : nbsous=ipt4.lisous(/1) 387 : nbref=ipt4.lisref(/1) 388 : IF (ICAS.NE.1) THEN 389 : nbref=0 390 : SEGADJ ipt4 391 : ENDIF 392 : ipt8=ipt1 393 : ipt5=ipt4 394 : numnp=0 395 : DO i=1,MAX(1,ipt1.lisous(/1)) 396 : IF (ipt1.lisous(/1).NE.0) THEN 397 : ipt8=ipt1.lisous(i) 398 : SEGACT ipt8 399 : SEGINI,ipt5=ipt8 400 : nbnn=ipt5.num(/1) 401 : nbelem=ipt5.num(/2) 402 : nbsous=ipt5.lisous(/1) 403 : nbref=ipt5.lisref(/1) 404 : IF (ICAS.NE.1) THEN 405 : nbref=0 406 : SEGADJ ipt5 407 : ENDIF 408 : ipt4.lisous(i)=ipt5 409 : ENDIF 410 : DO iel=1,ipt8.num(/2) 411 : DO inn=1,ipt8.num(/1) 412 : ipt=ipt8.num(inn,iel) 413 : IF (icpr(ipt).EQ.0) THEN 414 : IF (ICAS.EQ.1) THEN 415 : icpr(ipt)=ipt 416 : ELSE 417 : numnp=numnp+1 418 : icpr(ipt)=numnp+nbpan 419 : ENDIF 420 : ENDIF 421 : ipt5.num(inn,iel)=icpr(ipt) 422 : ENDDO 423 : ENDDO 424 : SEGDES ipt8,ipt5 425 : ENDDO 426 : SEGDES ipt4,ipt1 427 : IF (ICAS.EQ.0) THEN 428 : nbpts=numnp+nbpan 429 : SEGADJ mcoord 430 : ENDIF 431 : DO ipt=1,nbpan 432 : IF (icpr(ipt).NE.0) THEN 433 : jno=(ipt-1)*idimp1 434 : ino=(icpr(ipt)-1)*idimp1 435 : DO i=1,IDIM 436 : vb(i)=xcoor(jno+i)-xdep(i,1) 437 : ENDDO 438 : DO i=1,IDIM 439 : aa=xarr(i,1) 440 : DO j=1,IDIM 441 : aa=aa+pa(i,j)*vb(j) 442 : ENDDO 443 : xcoor(ino+i)=aa 444 : ENDDO 445 : xcoor(ino+idimp1)=xcoor(jno+idimp1) 446 : ENDIF 447 : ENDDO 448 : SEGSUP icpr,xxtra 448 : CALL ECROBJ('MAILLAGE',ipt4) 449 : IF (ichpo.EQ.1) THEN 450 : SEGSUP ipt2,ipt3 451 : ENDIF 452 : RETURN 453 : ENDIF 454 : 455 : C Cas d'une transformation non affine 456 : C Interdit dans le cas ICAS=1 (appel par deplac, operateur DEPL 'DEDU') 457 : C On recupere le maillage initial si ichpo=1. 458 : IF (ICAS.EQ.1) THEN 458 : CALL ERREUR(21) 459 : IF (ichpo.EQ.1) THEN 460 : nbpts=nbpts0 461 : SEGADJ mcoord 462 : SEGSUP ipt2,ipt3 463 : ENDIF 464 : RETURN 465 : ENDIF 466 : C CHOIX : Pour l'instant, cas non implemente en DIMEnsion 1. 467 : IF (IDIM.EQ.1) THEN 468 : MOTERR(1:4)=mcle(iregu) 469 : INTERR(1)=IDIM 469 : CALL ERREUR(971) 470 : IF (ichpo.EQ.1) THEN 471 : nbpts=nbpts0 472 : SEGADJ mcoord 473 : SEGSUP ipt2,ipt3 474 : ENDIF 475 : RETURN 476 : ENDIF 477 : C Creation d'un tableau de stockage contenant pour chaque point les 478 : C aretes des elements le touchant sauf celles qui le contiennent 479 : SEGACT ipt1 480 : SEGINI icpr 481 : numnp=0 482 : ipt8=ipt1 483 : DO i=1,MAX(1,ipt1.lisous(/1)) 484 : IF (ipt1.lisous(/1).NE.0) ipt8=ipt1.lisous(i) 485 : SEGACT ipt8 486 : DO iel=1,ipt8.num(/2) 487 : DO inn=1,ipt8.num(/1) 488 : ipt=ipt8.num(inn,iel) 489 : IF (icpr(ipt).EQ.0) THEN 490 : numnp=numnp+1 491 : icpr(ipt)=numnp 492 : ENDIF 493 : ENDDO 494 : ENDDO 495 : ENDDO 496 : SEGINI idcp 497 : DO i=1,icpr(/1) 498 : IF (icpr(i).NE.0) idcp(icpr(i))=i 499 : ENDDO 500 : nkon=5 501 : SEGINI kon,ikon 502 : DO i=1,MAX(1,ipt1.lisous(/1)) 503 : IF (ipt1.lisous(/1).NE.0) ipt8=ipt1.lisous(i) 504 : DO iel=1,ipt8.num(/2) 505 : DO inn=1,ipt8.num(/1) 506 : ip=ipt8.num(inn,iel) 507 : ipc=icpr(ip) 508 : DO 41 inp=1,ipt8.num(/1) 509 : kpt1=ipt8.num(inp,iel) 510 : IF (kpt1.EQ.ip) GOTO 41 511 : inps=inp+1 512 : IF (inps.GT.ipt8.num(/1)) inps=1 513 : kpt2=ipt8.num(inps,iel) 514 : IF (kpt2.EQ.ip) GOTO 41 515 : DO idj=1,ikon(ipc) 516 : IF ((kon(idj,ipc,1).EQ.kpt1).AND. 517 : . (kon(idj,ipc,2).EQ.kpt2)) GOTO 41 518 : IF ((kon(idj,ipc,1).EQ.kpt2).AND. 519 : . (kon(idj,ipc,2).EQ.kpt1)) GOTO 41 520 : ENDDO 521 : ikon(ipc)=ikon(ipc)+1 522 : IF (ikon(ipc).GT.nkon) THEN 523 : nkon=ikon(ipc) 524 : SEGADJ kon 525 : ENDIF 526 : kon(ikon(ipc),ipc,1)=kpt1 527 : kon(ikon(ipc),ipc,2)=kpt2 528 : 41 CONTINUE 529 : ENDDO 530 : ENDDO 531 : ENDDO 532 : C Fermeture du contour si necessaire (avec 1 segment supplementaire) 533 : DO 80 ipc=1,kon(/2) 534 : iuni=0 535 : DO iko=1,ikon(ipc) 536 : kpt1=kon(iko,ipc,1) 537 : kpt2=kon(iko,ipc,2) 538 : DO ikq=1,ikon(ipc) 539 : IF (ikq.NE.iko) THEN 540 : IF (kpt1.EQ.kon(ikq,ipc,1)) GOTO 86 541 : IF (kpt1.EQ.kon(ikq,ipc,2)) GOTO 86 542 : ENDIF 543 : ENDDO 544 : C kpt1 unique 545 : iuni=iuni+1 546 : iuniq(iuni)=kpt1 547 : 86 DO 88 ikq=1,ikon(ipc) 548 : IF (ikq.EQ.iko) GOTO 88 549 : IF (kpt2.EQ.kon(ikq,ipc,1)) GOTO 90 550 : IF (kpt2.EQ.kon(ikq,ipc,2)) GOTO 90 551 : 88 CONTINUE 552 : C kpt1 unique 553 : iuni=iuni+1 554 : iuniq(iuni)=kpt2 555 : 90 CONTINUE 556 : ENDDO 557 : IF (iuni.EQ.0) GOTO 80 558 : ikon(ipc)=ikon(ipc)+1 559 : IF (ikon(ipc).GT.nkon) THEN 560 : nkon=ikon(ipc) 561 : SEGADJ kon 562 : ENDIF 563 : kon(ikon(ipc),ipc,1)=iuniq(1) 564 : kon(ikon(ipc),ipc,2)=iuniq(2) 565 : 80 CONTINUE 566 : C Calcul des coordonnees barycentriques 567 : SEGINI xkon,ipxkon 568 : DO i=1,kon(/2) 569 : ipcour=0 570 : jp=idcp(i) 571 : ipcini=ipcour 572 : 58 CONTINUE 573 : DO ipc=ipcini+1,ipcour 574 : ipxkon(ipc,i)=0 575 : xkon(ipc,i)=0. 576 : ENDDO 577 : ipcour=ipcini 578 : DO 59 j=1,ikon(i) 579 : C Creation d'un angle 580 : IF (kon(j,i,1).EQ.0) GOTO 59 581 : ip1=kon(j,i,1) 582 : ip2=kon(j,i,2) 583 : DO 60 k=j+1,ikon(i) 584 : IF (kon(k,i,1).EQ.0) GOTO 60 585 : jp1=0 586 : IF (ip1.EQ.kon(k,i,1)) THEN 587 : jp1=ip1 588 : jp2=ip2 589 : jp3=kon(k,i,2) 590 : ELSE IF (ip1.EQ.kon(k,i,2)) THEN 591 : jp1=ip1 592 : jp2=ip2 593 : jp3=kon(k,i,1) 594 : ELSE IF (ip2.EQ.kon(k,i,1)) THEN 595 : jp1=ip2 596 : jp2=ip1 597 : jp3=kon(k,i,2) 598 : ELSE IF (ip2.EQ.kon(k,i,2)) THEN 599 : jp1=ip2 600 : jp2=ip1 601 : jp3=kon(k,i,1) 602 : ELSE 603 : GOTO 60 604 : ENDIF 605 : * jp le point en cours de travail 606 : * jp1 le sommet jp2 et jp3 les deux autres points 607 : xp=xcoor((jp-1)*idimp1+1) 608 : yp=xcoor((jp-1)*idimp1+2) 609 : zp=0. 610 : IF (IDIM.EQ.3) zp=xcoor((jp-1)*idimp1+3) 611 : xp1=xcoor((jp1-1)*idimp1+1) 612 : yp1=xcoor((jp1-1)*idimp1+2) 613 : zp1=0. 614 : IF (IDIM.EQ.3) zp1=xcoor((jp1-1)*idimp1+3) 615 : xp2=xcoor((jp2-1)*idimp1+1) 616 : yp2=xcoor((jp2-1)*idimp1+2) 617 : zp2=0. 618 : IF (IDIM.EQ.3) zp2=xcoor((jp2-1)*idimp1+3) 619 : xp3=xcoor((jp3-1)*idimp1+1) 620 : yp3=xcoor((jp3-1)*idimp1+2) 621 : zp3=0. 622 : IF (IDIM.EQ.3) zp3=xcoor((jp3-1)*idimp1+3) 623 : * le sinus (trafique) 624 : xv1=xp2-xp1 625 : yv1=yp2-yp1 626 : zv1=zp2-zp1 627 : aa=SQRT(xv1*xv1+yv1*yv1+zv1*zv1) 628 : xv1=xv1/aa 629 : yv1=yv1/aa 630 : zv1=zv1/aa 631 : xv2=xp3-xp1 632 : yv2=yp3-yp1 633 : zv2=zp3-zp1 634 : aa=SQRT(xv2*xv2+yv2*yv2+zv2*zv2) 635 : xv2=xv2/aa 636 : yv2=yv2/aa 637 : zv2=zv2/aa 638 : xv=yv1*zv2-zv1*yv2 639 : yv=zv1*xv2-xv1*zv2 640 : zv=xv1*yv2-yv1*xv2 641 : * provisoirement on prend le signe positif 642 : sinp=MAX(1.D-15,SQRT(xv*xv+yv*yv+zv*zv)) 643 : * recherche du signe comparaison avec jp jp1 vect jp2 jp3 644 : xw1=xp-xp1 645 : yw1=yp-yp1 646 : zw1=zp-zp1 647 : aa=SQRT(xw1*xw1+yw1*yw1+zw1*zw1) 648 : xw1=xw1/aa 649 : yw1=yw1/aa 650 : zw1=zw1/aa 651 : xw2=xp3-xp2 652 : yw2=yp3-yp2 653 : zw2=zp3-zp2 654 : aa=SQRT(xw2*xw2+yw2*yw2+zw2*zw2) 655 : xw2=xw2/aa 656 : yw2=yw2/aa 657 : zw2=zw2/aa 658 : xw=yw1*zw2-zw1*yw2 659 : yw=zw1*xw2-xw1*zw2 660 : zw=xw1*yw2-yw1*xw2 661 : IF (iregu.NE.2) THEN 662 : IF ((xv*xw+yv*yw+zv*zw).LT.0.) sinp=-sinp 663 : ELSE 664 : IF ((xv*xw+yv*yw+zv*zw).LT.-1.D-15) then 665 : * Essai si le signe est negatif on supprime le coin pour ne travailler 666 : * que dans un convexe 667 : kon(k,i,1)=0 668 : kon(k,i,2)=0 669 : kon(j,i,1)=jp2 670 : kon(j,i,2)=jp3 671 : GOTO 58 672 : ENDIF 673 : ENDIF 674 : ipcour=ipcour+1 675 : ipxkon(ipcour,i)=jp1 676 : * les distances aux cotes 677 : xv=xp-xp1 678 : yv=yp-yp1 679 : zv=zp-zp1 680 : scal=xv*xv1+yv*yv1+zv*zv1 681 : xv=xv-scal*xv1 682 : yv=yv-scal*yv1 683 : zv=zv-scal*zv1 684 : dis1=MAX(1.D-20,SQRT(xv*xv+yv*yv+zv*zv)) 685 : xv=xp-xp1 686 : yv=yp-yp1 687 : zv=zp-zp1 688 : scal=xv*xv2+yv*yv2+zv*zv2 689 : xv=xv-scal*xv2 690 : yv=yv-scal*yv2 691 : zv=zv-scal*zv2 692 : dis2=MAX(1.D-20,SQRT(xv*xv+yv*yv+zv*zv)) 693 : C la coordonnee barycentrique : 694 : xkon(ipcour,i)=sinp/(dis1*dis2) 695 : IF (iregu.EQ.1) xkon(ipcour,i)=1. 696 : 60 CONTINUE 697 : 59 CONTINUE 698 : * on normalise 699 : skon=0. 700 : DO j=1,ikon(i)+1 701 : skon=skon+xkon(j,i) 702 : ENDDO 703 : DO j=1,ikon(i)+1 704 : xkon(j,i)=xkon(j,i)/skon 705 : ENDDO 706 : ENDDO 707 : SEGSUP kon,ikon 708 : * OK on a les coefficients barycentriques. 709 : * On duplique le maillage puis on cree les nouveaux noeuds 710 : * Mais d'abord voir les points imposes 711 : SEGACT ipt2,ipt3 712 : SEGINI icp2 713 : ipt7=ipt2 714 : ipt8=ipt3 715 : DO isous=1,MAX(1,ipt2.lisous(/1)) 716 : IF (ipt2.lisous(/1).NE.0) ipt7=ipt2.lisous(isous) 717 : IF (ipt3.lisous(/1).NE.0) ipt8=ipt3.lisous(isous) 718 : SEGACT ipt7,ipt8 719 : DO i=1,ipt7.num(/1) 720 : DO j=1,ipt7.num(/2) 721 : ip=ipt7.num(i,j) 722 : IF (icp2(ip).EQ.0) icp2(ip)=ipt8.num(i,j) 723 : ENDDO 724 : ENDDO 725 : SEGDES ipt7,ipt8 726 : ENDDO 727 : SEGDES ipt2,ipt3 728 : SEGINI icorr 729 : DO 210 i=1,icpr(/1) 730 : IF (icpr(i).EQ.0) GOTO 210 731 : IF (icp2(i).EQ.0) THEN 732 : nbpts=xcoor(/1)/idimp1+1 733 : SEGADJ mcoord 734 : icorr(icpr(i))=nbpts 735 : DO j=1,idimp1 736 : xcoor((nbpts-1)*idimp1+j)=xcoor((i-1)*idimp1+j) 737 : ENDDO 738 : ELSE 739 : icorr(icpr(i))=icp2(i) 740 : ENDIF 741 : 210 CONTINUE 742 : * On duplique les maillages 743 : SEGINI,ipt4=ipt1 744 : nbref=0 745 : nbsous=ipt4.lisous(/1) 746 : nbnn=ipt4.num(/1) 747 : nbelem=ipt4.num(/2) 748 : SEGADJ ipt4 749 : ipt7=ipt1 750 : ipt8=ipt4 751 : DO isous=1,MAX(1,ipt1.lisous(/1)) 752 : IF (ipt1.lisous(/1).NE.0) THEN 753 : ipt7=ipt1.lisous(isous) 754 : SEGINI,ipt8=ipt7 755 : ipt4.lisous(isous)=ipt8 756 : ENDIF 757 : DO i=1,ipt8.num(/1) 758 : DO j=1,ipt8.num(/2) 759 : ipt8.num(i,j)=icorr(icpr(ipt8.num(i,j))) 760 : ENDDO 761 : ENDDO 762 : SEGDES ipt7,ipt8 763 : ENDDO 764 : SEGDES ipt1,ipt4 765 : C Maintenant on peut deplacer les points 766 : nfois=50 766 : CALL LIRENT(nfois,0,iretou) 767 : DO ifois=1,nfois 768 : DO i=1,numnp 769 : IF (icp2(idcp(i)).EQ.0) THEN 770 : DO j=1,IDIM 771 : xaux=0. 772 : DO k=1,nkon+1 773 : IF (ipxkon(k,i).NE.0) xaux=xaux+xkon(k,i) 774 : . *xcoor((icorr(icpr(ipxkon(k,i)))-1)*idimp1+j) 775 : ENDDO 776 : xcoor((icorr(i)-1)*idimp1+j)=xaux 777 : ENDDO 778 : ENDIF 779 : ENDDO 780 : ENDDO 781 : 782 : IF (ichpo.EQ.0) THEN 782 : CALL ECROBJ('MAILLAGE',ipt4) 783 : ELSE 784 : C Sortie d'un CHPOINT, les points concernes sont reperes dans icpr 785 : SEGSUP ipt2,ipt3 786 : nbelem=icorr(/1) 787 : nbnn=1 788 : nbsous=0 789 : nbref=0 790 : SEGINI ipt2 791 : ipt2.itypel=1 792 : nat=1 793 : nsoupo=1 794 : SEGINI mchpo1 795 : mchpo1.jattri(1)=jatan 796 : IF (IDIM.EQ.1) nc=1 797 : IF (IDIM.EQ.2) nc=2 798 : IF (IDIM.EQ.3) nc=3 799 : SEGINI msoupo 800 : mchpo1.ipchp(1)=msoupo 801 : igeoc=ipt2 802 : n=nbelem 803 : SEGINI mpoval 804 : ipoval=mpoval 805 : IF (IDIM.EQ.1) THEN 806 : IF (NIFOUR.LE.9) THEN 807 : nocomp(1)='UX' 808 : ELSE 809 : nocomp(1)='UR' 810 : ENDIF 811 : ELSE IF (IDIM.EQ.2) THEN 812 : IF (IFOUR.GE.0) THEN 813 : nocomp(1)='UR' 814 : nocomp(2)='UZ' 815 : ELSE 816 : nocomp(1)='UX' 817 : nocomp(2)='UY' 818 : ENDIF 819 : ELSE IF (IDIM.EQ.3) THEN 820 : nocomp(1)='UX' 821 : nocomp(2)='UY' 822 : nocomp(3)='UZ' 823 : ENDIF 824 : SEGDES msoupo 825 : DO i=1,idcp(/1) 826 : ipt2.num(1,i)=idcp(i) 827 : ifi=(icorr(i)-1)*idimp1 828 : ide=(idcp(i)-1)*idimp1 829 : DO j=1,IDIM 830 : vpocha(i,j)=xcoor(ifi+j)-xcoor(ide+j) 831 : ENDDO 832 : ENDDO 833 : SEGDES mpoval,ipt2,mchpo1 833 : CALL ECROBJ('CHPOINT ',mchpo1) 834 : ENDIF 835 : SEGSUP icpr,idcp,xkon,ipxkon,icp2,icorr 836 : 837 : RETURN 838 : END 839 : 840 : 841 :
© Cast3M 2003 - All rights reserved.
Disclaimer