1 : C IMPO31 SOURCE PV 13/12/13 21:15:01 7883 2 : 3 : C Operateur 'IMPO' 'MAIL' ( | 'SYME' | ) Surf1 Surf2 ( Crit) ; 4 : C 3D | 'MESC' | 5 : C | 'FAIB' | 6 : 7 : SUBROUTINE IMPO31 8 : * 9 : * modif octobre 2012 on ne cree pas les noeuds support des mult ici 10 : * on le fera dynamiquement dans impo32 11 : * 12 : 13 : 14 : IMPLICIT INTEGER(I-N) 15 : IMPLICIT REAL*8 (A-H,O-Z) 16 : 17 : -INC CCOPTIO 18 : 19 : -INC SMCOORD 20 : -INC SMELEME 21 : -INC CCGEOME 22 : 23 : SEGMENT IAUXTB(NBAUX) 24 : PARAMETER ( X1s4 = 0.25000000000000000000000000000000000000000D0 , 25 : & X1s6 = 0.16666666666666666666666666666666666666667D0 ) 26 : * 27 : LOGICAL faible 28 : PARAMETER (ncle = 3) 29 : CHARACTER*4 mcle(ncle) 30 : DATA mcle / 'MESC','FAIB','SYME' / 31 : * 32 : idimp1 = IDIM + 1 33 : * 34 : * Lecture d'un mot-cle eventuel : SYME(trique) ou MESC ou FAIB(le) 35 : * 36 : mesc = 0 36 : call lirmot(mcle,ncle,mesc,0) 37 : if (ierr.ne.0) return 38 : * 39 : * Lecture des deux surfaces en contact 40 : * 40 : call lirobj('MAILLAGE',ipt1,1,iretou) 41 : if (ierr.ne.0) return 41 : call lirobj('MAILLAGE',ipt2,1,iretou) 42 : if (ierr.ne.0) return 43 : * 44 : * Lecture eventuelle de la distance limite max(crit) 45 : * 46 : crit = 1.e+35 46 : call lirree(crit,0,iretou) 47 : if (ierr.ne.0) return 48 : cri2 = crit*crit 49 : * crit = abs(crit) 50 : *D write(ioimp,*) ' critere ',crit,cri2 51 : * 52 : * Lecture eventuelle de la distance limite min(critmin) 53 : * 54 : critmin = 0.d0 54 : call lirree(critmin,0,iretou) 55 : if (ierr.ne.0) return 56 : crimi2 = critmin*critmin 57 : * crit = abs(crit) 58 : *D write(ioimp,*) ' critere ',crit,cri2 59 : * 60 : * Activation et quelques verifications sur les maillages des surfaces 61 : * 62 : segact ipt1,ipt2 62 : if (ipt1.itypel.ne.4) call erreur(16) 62 : if (ipt2.itypel.ne.4) call erreur(16) 63 : if (ierr.ne.0) goto 900 64 : * 65 : nbpts0 = xcoor(/1) / idimp1 66 : * 67 : nel1 = ipt1.num(/2) 68 : nel2 = ipt2.num(/2) 69 : * noe1 = ipt1.num(/1) 70 : noe2 = ipt2.num(/1) 71 : * En fait : noe1 = noe2 = 3 car elements de type TRI3 ! 72 : * 73 : * Increment sur le nombre d'elements de contact retenus pour agrandir 74 : * les longueurs des segments (maillage, mcoord) de maniere adequate 75 : incnel = 1000000 76 : * 77 : faible = mesc.eq.2 78 : if (mesc .ne. 1) mesc = 0 79 : * creation si il y a lieu du point support des multiplicateurs de lagrange non instanciées 80 : if (ilgni.eq.0) then 81 : nbpts = nbpts0 +1 82 : nbpts0 = nbpts 83 : segadj,mcoord 84 : ilgni = nbpts 85 : endif 86 : * 87 : *----------------------------------------------------------------------- 88 : *- FORMULATION "STANDARD" DU CONTACT - 89 : *----------------------------------------------------------------------- 90 : IF (faible) GOTO 500 91 : * 92 : * Nombre de points supports (multiplicateurs de Lagrange) associes 93 : * au contact (1 par element de contact) 94 : * Ajustement du segment mcoord 95 : * 96 : nbpts = nbpts0 + (incnel * 1) 97 : segadj,mcoord 98 : * 99 : * tableau auxiliaire pour accelerer le test element deja cree 100 : * 101 : nbaux = nbpts0 102 : * 103 : * Creation du maillage de contact (meleme) 104 : * 105 : nbsous = 0 106 : nbref = 0 107 : nbnn = 5 108 : nbelem = incnel 109 : segini,meleme 110 : itypel = 22 111 : * 112 : nblag = nbpts0 113 : nbcoo = (nblag-1)*idimp1 114 : nbelc = 0 115 : * 116 : ipas = 1 117 : * 118 : 5 continue 119 : segini iauxtb 120 : * 121 : * Boucle sur les elements de la 1ere surface de contact 122 : * 123 : do 10 iel=1,nel1 124 : nbeini = nbelc 125 : ip1 = ipt1.num(1,iel) 126 : ip2 = ipt1.num(2,iel) 127 : ip3 = ipt1.num(3,iel) 128 : ipv = (ip1-1)*idimp1 129 : xp1 = xcoor(ipv+1) 130 : yp1 = xcoor(ipv+2) 131 : zp1 = xcoor(ipv+3) 132 : ipv = (ip2-1)*idimp1 133 : xp2 = xcoor(ipv+1) 134 : yp2 = xcoor(ipv+2) 135 : zp2 = xcoor(ipv+3) 136 : ipv = (ip3-1)*idimp1 137 : xp3 = xcoor(ipv+1) 138 : yp3 = xcoor(ipv+2) 139 : zp3 = xcoor(ipv+3) 140 : xg1 = xp1+xp2+xp3 141 : yg1 = yp1+yp2+yp3 142 : zg1 = zp1+zp2+zp3 143 : * 144 : * Boucle sur les noeuds de la 2eme surface de contact 145 : * 146 : do 20 jel=1,nel2 147 : do 20 jnn=1,noe2 148 : jp=ipt2.num(jnn,jel) 149 : * verification que pas relation du point sur lui meme 150 : if (jp.eq.ip1) goto 20 151 : if (jp.eq.ip2) goto 20 152 : if (jp.eq.ip3) goto 20 153 : * verification que pas deja la relation 154 : if (iauxtb(jp).eq.iel) goto 20 155 : iauxtb(jp) = iel 156 : ipv = (jp-1)*idimp1 157 : xp=xcoor(ipv+1) 158 : yp=xcoor(ipv+2) 159 : zp=xcoor(ipv+3) 160 : * test distance 161 : d1=(xp1-xp)**2+(yp1-yp)**2+(zp1-zp)**2 162 : if (d1.le.cri2.and.d1.ge.crimi2) goto 22 163 : d1=(xp2-xp)**2+(yp2-yp)**2+(zp2-zp)**2 164 : if (d1.le.cri2.and.d1.ge.crimi2) goto 22 165 : d1=(xp3-xp)**2+(yp3-yp)**2+(zp3-zp)**2 166 : if (d1.le.cri2.and.d1.ge.crimi2) goto 22 167 : goto 20 168 : * 169 : 22 continue 170 : * 171 : * Mise a jour de la dimension des segments si necessaire 172 : * 173 : if (nbelc.eq.nbelem) then 174 : nbelem = nbelem + incnel 175 : write (6,*) ' ajustement impo31 ',nbelem 176 : segadj,meleme 177 : nbpts = nbpts + incnel 178 : segadj,mcoord 179 : endif 180 : * 181 : * on incremente mcoord avec le pt support du multiplicateur 182 : * 183 : nblag = nblag + 1 184 : nbcoo = nbcoo + idimp1 185 : nbelc = nbelc + 1 186 : * 187 : xcoor(nbcoo+1) = (xg1+xp) * X1s4 188 : xcoor(nbcoo+2) = (yg1+yp) * X1s4 189 : xcoor(nbcoo+3) = (zg1+zp) * X1s4 190 : xcoor(nbcoo+4) = 0. 191 : * 192 : * on range dans le meleme 193 : * 194 : * num(1,nbelc) = ilgni 195 : num(1,nbelc) = nblag 196 : num(2,nbelc) = ip1 197 : num(3,nbelc) = ip2 198 : num(4,nbelc) = ip3 199 : num(5,nbelc) = jp 200 : * 201 : 20 continue 202 : * 203 : 10 continue 204 : segsup iauxtb 205 : * 206 : * On recommence en inversant les deux surfaces en cas de traitement 207 : * symetrique du contact (option MESC non activee) 208 : * 209 : if (ipas.eq.1.and.mesc.ne.1) then 210 : ipas=ipas+1 211 : iaux=ipt1 212 : ipt1=ipt2 213 : ipt2=iaux 214 : nel1 = ipt1.num(/2) 215 : nel2 = ipt2.num(/2) 216 : noe2 = ipt2.num(/1) 217 : goto 5 218 : endif 219 : * 220 : GOTO 1000 221 : * 222 : *----------------------------------------------------------------------- 223 : *- FORMULATION "FAIBLE" DU CONTACT - 224 : *----------------------------------------------------------------------- 225 : 500 CONTINUE 226 : * 227 : * Le nb d'elements maximal a creer pour le maillage de contact FAIBLE : 228 : * 1 entre chaque element de chacune des surfaces 229 : * 230 : * Creation du maillage de contact (meleme) 231 : * 232 : nbnn = 7 233 : nbsous = 0 234 : nbref = 0 235 : nbelem = incnel 236 : segini,meleme 237 : itypel = 22 238 : * 239 : * Nombre de points supports (multiplicateurs de Lagrange) associes 240 : * au contact : 1 par element de contact 241 : * Ajustement du segment mcoord 242 : * 243 : nbpts = nbpts0 + (incnel * 1) 244 : segadj,mcoord 245 : * 246 : nblag = nbpts0 247 : nbcoo = (nblag-1)*idimp1 248 : nbelc = 0 249 : * 250 : * Boucle sur les elements du maillage de la 1ere surface de contact 251 : * 252 : do 510 iel = 1, nel1 253 : ip1 = ipt1.num(1,iel) 254 : ip2 = ipt1.num(2,iel) 255 : ip3 = ipt1.num(3,iel) 256 : ipv = (ip1-1)*idimp1 257 : xp1 = xcoor(ipv+1) 258 : yp1 = xcoor(ipv+2) 259 : zp1 = xcoor(ipv+3) 260 : ipv = (ip2-1)*idimp1 261 : xp2 = xcoor(ipv+1) 262 : yp2 = xcoor(ipv+2) 263 : zp2 = xcoor(ipv+3) 264 : ipv = (ip3-1)*idimp1 265 : xp3 = xcoor(ipv+1) 266 : yp3 = xcoor(ipv+2) 267 : zp3 = xcoor(ipv+3) 268 : xg1 = xp1 + xp2 + xp3 269 : yg1 = yp1 + yp2 + yp3 270 : zg1 = zp1 + zp2 + zp3 271 : * 272 : * Boucle sur les elements du maillage de la 2e surface de contact 273 : * 274 : do 520 jel = 1, nel2 275 : jp1 = ipt2.num(1,jel) 276 : jp2 = ipt2.num(2,jel) 277 : jp3 = ipt2.num(3,jel) 278 : ipv = (jp1-1)*idimp1 279 : xp4 = xcoor(ipv+1) 280 : yp4 = xcoor(ipv+2) 281 : zp4 = xcoor(ipv+3) 282 : ipv = (jp2-1)*idimp1 283 : xp5 = xcoor(ipv+1) 284 : yp5 = xcoor(ipv+2) 285 : zp5 = xcoor(ipv+3) 286 : ipv = (jp3-1)*idimp1 287 : xp6 = xcoor(ipv+1) 288 : yp6 = xcoor(ipv+2) 289 : zp6 = xcoor(ipv+3) 290 : 291 : * test distance 292 : d1 = (xp1-xp4)**2+(yp1-yp4)**2+(zp1-zp4)**2 293 : d2 = (xp2-xp4)**2+(yp2-yp4)**2+(zp2-zp4)**2 294 : d3 = (xp3-xp4)**2+(yp3-yp4)**2+(zp3-zp4)**2 295 : if ((d1.le.cri2.or.d2.le.cri2.or.d3.le.cri2).and. 296 : > (d1.ge.crimi2.or.d2.ge.crimi2.or.d3.ge.crimi2)) goto 550 297 : 298 : d1 = (xp1-xp5)**2+(yp1-yp5)**2+(zp1-zp5)**2 299 : d2 = (xp2-xp5)**2+(yp2-yp5)**2+(zp2-zp5)**2 300 : d3 = (xp3-xp5)**2+(yp3-yp5)**2+(zp3-zp5)**2 301 : if ((d1.le.cri2.or.d2.le.cri2.or.d3.le.cri2).and. 302 : > (d1.ge.crimi2.or.d2.ge.crimi2.or.d3.ge.crimi2)) goto 550 303 : 304 : d1 = (xp1-xp6)**2+(yp1-yp6)**2+(zp1-zp6)**2 305 : d2 = (xp2-xp6)**2+(yp2-yp6)**2+(zp2-zp6)**2 306 : d3 = (xp3-xp6)**2+(yp3-yp6)**2+(zp3-zp6)**2 307 : if ((d1.le.cri2.or.d2.le.cri2.or.d3.le.cri2).and. 308 : > (d1.ge.crimi2.or.d2.ge.crimi2.or.d3.ge.crimi2)) goto 550 309 : 310 : goto 520 311 : * 312 : 550 continue 313 : * 314 : * Mise a jour de la dimension des segments si necessaire 315 : * 316 : if (nbelc.eq.nbelem) then 317 : nbelem = nbelem + incnel 318 : segadj,meleme 319 : nbpts = nbpts + incnel 320 : segadj,mcoord 321 : endif 322 : * 323 : * on incremente mcoord avec les pts support des multiplicateurs 324 : * 325 : nblag = nblag + 1 326 : nbcoo = nbcoo + idimp1 327 : nbelc = nbelc + 1 328 : * 329 : xg2 = xp4 + xp5 + xp6 330 : yg2 = yp4 + yp5 + yp6 331 : zg2 = zp4 + zp5 + zp6 332 : * 333 : xcoor(nbcoo+1) = (xg1+xg2) * X1s6 334 : xcoor(nbcoo+2) = (yg1+yg2) * X1s6 335 : xcoor(nbcoo+3) = (zg1+zg2) * X1s6 336 : xcoor(nbcoo+4) = 0. 337 : * 338 : * on range dans le meleme 339 : * 340 : * num(1,nbelc) = ilgni 341 : num(1,nbelc) = nblag 342 : num(2,nbelc) = ip1 343 : num(3,nbelc) = ip2 344 : num(4,nbelc) = ip3 345 : num(5,nbelc) = jp1 346 : num(6,nbelc) = jp3 347 : num(7,nbelc) = jp2 348 : * 349 : 520 continue 350 : * 351 : 510 continue 352 : * 353 : * GOTO 1000 354 : * 355 : *----------------------------------------------------------------------- 356 : * FIN : 357 : *----------------------------------------------------------------------- 358 : 1000 CONTINUE 359 : *D write(ioimp,*) ' nbpts nbpts0 nblag',nbpts,nbpts0,nblag 360 : *D write(ioimp,*) ' nbelem nbelc',nbelem,nbelc 361 : * 362 : * On renvoie le resultat mais au prealable on ajuste les longueurs 363 : * 364 : if (nbelem.lt.nbelc) then 364 : call erreur(5) 365 : segsup,meleme 366 : goto 900 367 : else if (nbelem.gt.nbelc) then 368 : nbelem = nbelc 369 : segadj,meleme 370 : nbpts = nblag 371 : segadj,mcoord 372 : endif 373 : segdes,meleme 373 : call ecrobj('MAILLAGE',meleme) 374 : 375 : 900 CONTINUE 376 : segdes,ipt1,ipt2 377 : 378 : return 379 : end 380 : 381 : 382 : 383 :
© Cast3M 2003 - All rights reserved.
Disclaimer