Télécharger accro.eso

Retour à la liste

Numérotation des lignes :

  1. C ACCRO SOURCE AM 18/03/13 21:15:01 9776
  2.  
  3. C========================================================================
  4. C ACCROCHAGE D UN MAILLAGE SUR UN MASSIF PAR L INTERMEDIAIRE
  5. C DES FONCTIONS DE FORME
  6. C SUBROUTINE APPELLEE PAR RELA2
  7. C JMB OCTOBRE 1996
  8. C========================================================================
  9. SUBROUTINE ACCRO(igliss)
  10.  
  11. IMPLICIT INTEGER(I-N)
  12. IMPLICIT REAL*8 (A-H,O-Z)
  13.  
  14. -INC CCOPTIO
  15. -INC CCREEL
  16. -INC CCHAMP
  17. -INC CCGEOME
  18.  
  19. -INC SMCOORD
  20. -INC SMELEME
  21. -INC SMRIGID
  22. -INC SMLMOTS
  23. -INC SMMODEL
  24.  
  25. SEGMENT iexclu
  26. INTEGER IDEJVU(NBPTI)
  27. ENDSEGMENT
  28.  
  29. SEGMENT icount
  30. INTEGER NBPTOT
  31. INTEGER MAXPZ(nbsz)
  32. INTEGER IEINT(3,NODES)
  33. REAL*8 QQQ(3,NODES)
  34. ENDSEGMENT
  35.  
  36. SEGMENT ioubli
  37. INTEGER iptou(nboub)
  38. ENDSEGMENT
  39.  
  40. SEGMENT mtrav
  41. INTEGER inomil(NBNOE1)
  42. REAL*8 QSI(3)
  43. REAL*8 SHPP(6,NBNOE1),XEL(3,NBNOE1)
  44. ENDSEGMENT
  45.  
  46. SEGMENT itange
  47. INTEGER iextr(nbpti)
  48. REAL*8 xperp1(3,nbpti), xperp2(3,nbpti)
  49. ENDSEGMENT
  50.  
  51. SEGMENT icpr(nbpti)
  52. SEGMENT imul(idim,nbpt)
  53. SEGMENT icomu(nbpt)
  54.  
  55. DIMENSION ICOR(100),IMEL(100),imtt(100)
  56. DIMENSION xpo(3),CCOE(100)
  57.  
  58. LOGICAL ltelq
  59.  
  60. idimp1 = idim + 1
  61.  
  62. MAIL1 = 0
  63. MAIL2 = 0
  64. IFROCA = 0
  65. MLMOTX = 0
  66. CCRIT = 1.D-5
  67.  
  68. C**** LECTURE des MAILLAGE et MMODEL ***********************************
  69. CALL LIROBJ('MAILLAGE',MAIL1,1,iretou)
  70. IF (IERR.NE.0) RETURN
  71. CALL LIROBJ('MAILLAGE',MAIL2,0,iret2)
  72. IF (IERR.NE.0) RETURN
  73.  
  74. IF (iret2.NE.0) THEN
  75. MELEMI = MAIL2
  76. IPP2 = MAIL1
  77. ELSE
  78. CALL LIROBJ('MMODEL',IPMODL,1,IFROCA)
  79. IF (IERR.NE.0) RETURN
  80. c IPP2 sera extrait du modele IPMODL
  81. MELEMI = MAIL1
  82. IPP2 = 0
  83. ENDIF
  84.  
  85. C**** LECTURE DES NOMS D'INCONNUE ET DU CRITERE ************************
  86. CALL LIROBJ('LISTMOTS',MLMOTX,0,iretou)
  87. IF (IERR.NE.0) RETURN
  88.  
  89. CALL LIRREE(CCRIT,0,iretou)
  90. IF (IERR.NE.0) RETURN
  91. CCRIT = -ABS(CCRIT)
  92.  
  93. C LISTE DES NOMS D'INCONNUES
  94. C Options par defaut
  95. IF (MLMOTX.EQ.0) THEN
  96. NINC = IDIM
  97. IF (IFOUR.EQ.1) NINC = NINC + 1
  98. JGN = 4
  99. JGM = NINC
  100. SEGINI,mlmots
  101. C 2D axisymetrique
  102. IF (IFOUR.EQ.0) THEN
  103. mlmots.MOTS(1) = 'UR '
  104. mlmots.MOTS(2) = 'UZ '
  105. C 2D Fourier
  106. ELSE IF (IFOUR.EQ.1) THEN
  107. mlmots.MOTS(1) = 'UR '
  108. mlmots.MOTS(2) = 'UZ '
  109. mlmots.MOTS(3) = 'UT '
  110. ELSE
  111. mlmots.MOTS(1) = 'UX '
  112. mlmots.MOTS(2) = 'UY '
  113. IF (IDIM.EQ.3) mlmots.MOTS(3) = 'UZ '
  114. ENDIF
  115. C On a lu les noms d'inconnues
  116. ELSE
  117. IF (igliss.EQ.1) THEN
  118. CALL ERREUR(19)
  119. RETURN
  120. ENDIF
  121. mlmots = MLMOTX
  122. SEGACT,mlmots
  123. NINC = mlmots.MOTS(/2)
  124. IF (NINC.LT.1 .OR. NINC.GT.100) THEN
  125. SEGDES,mlmots
  126. CALL ERREUR(19)
  127. RETURN
  128. ENDIF
  129. ENDIF
  130. C Identification des duals des ddl en question
  131. c bp 18.01.2011 : mieux vaut utiliser LNOMDD renseigné dans cchamp+bdata
  132. DO i = 1, NINC
  133. k = 0
  134. CALL PLACE(NOMDD,LNOMDD,k,mlmots.MOTS(i))
  135. IF (k.EQ.0) THEN
  136. MOTERR(1:8) = mlmots.MOTS(i)
  137. CALL ERREUR(931)
  138. RETURN
  139. ENDIF
  140. ICOR(i) = k
  141. ENDDO
  142. SEGDES,mlmots
  143. IF (MLMOTX.EQ.0) SEGSUP,mlmots
  144.  
  145. C**** Nombre de noeuds du segment MCOORD *******************************
  146. ** SEGACT,mcoord*NOMOD
  147. NBPTI = mcoord.XCOOR(/1) / idimp1
  148.  
  149. C**** VERIFICATION DU MAILLAGE MELEMI **********************************
  150. C CE MAILLAGE DOIT ETRE MASSIF
  151. C TETRAEDRES ET PYRAMIDES SONT EXCLUS
  152. IPT1 = MELEMI
  153. SEGACT,IPT1
  154. NBSZ1 = IPT1.LISOUS(/1)
  155. nbsz = MAX(1,NBSZ1)
  156. NBNOE1 = 0
  157. meleme = IPT1
  158. DO iob = 1, nbsz
  159. IF (NBSZ1.NE.0) THEN
  160. meleme = ipt1.LISOUS(iob)
  161. SEGACT,meleme
  162. ENDIF
  163. ity = meleme.ITYPEL
  164. nbnn = meleme.NUM(/1)
  165. NBNOE1 = MAX(NBNOE1,nbnn)
  166. * SEGDES,meleme
  167. IF ((IDIM.EQ.2.AND.ity.LT.4) .OR.
  168. & (IDIM.EQ.3.AND.ity.LT.14)) THEN
  169. CALL ERREUR(16)
  170. RETURN
  171. ENDIF
  172. ENDDO
  173.  
  174. C**** CAS OU UN MODELE A ETE LU * EXTRATION DE IPP2 ********************
  175. IF (IFROCA.NE.0) THEN
  176. mmode1 = IPMODL
  177. segact,mmode1
  178. c* nbpti = mcoord.xcoor(/1) / idimp1
  179. segini,icpr
  180. nbpt = 0
  181. C creation des l'ensemble des points a accrocher et des multiplicateurs
  182. C correspondant. imul(i,j) est le ieme multiplicateur a accrocher au j eme
  183. C point local ( icpr(ia)=j)
  184. C --- puis Assemblage du maillage global (tous les câbles)
  185. ltelq=.false.
  186. do 500 isouf = 1, mmode1.kmodel(/1)
  187. imode1 = mmode1.kmodel(isouf)
  188. segact imode1
  189. if (imode1.tymode(1).ne.'MMODEL') then
  190. call erreur(975)
  191. return
  192. endif
  193. ipt6 = imode1.imamod
  194. segact ipt6
  195. if (ipt6.itypel.ne.22) then
  196. call erreur(19)
  197. return
  198. endif
  199. do iel = 1,ipt6.num(/2)
  200. ia = ipt6.num(2,iel)
  201. if (icpr(ia).eq.0) then
  202. nbpt = nbpt+1
  203. icpr(ia) = nbpt
  204. endif
  205. enddo
  206. if (imode1.tymode(1).ne.'MMODEL') then
  207. call erreur(975)
  208. return
  209. endif
  210. mmode2 = imode1.ivamod(1)
  211. segact mmode2
  212. imode2 = mmode2.kmodel(1)
  213. segact imode2
  214. ipt7 = imode2.imamod
  215. segdes,imode2,mmode2
  216. if (isouf.eq.1) then
  217. ipt5 = ipt7
  218. else
  219. call fuse(ipt5,ipt7,ipp4,ltelq)
  220. ipt5 = ipp4
  221. endif
  222. 500 continue
  223. segini,imul,icomu
  224. do 502 isouf = 1, mmode1.kmodel(/1)
  225. imode1 = mmode1.kmodel(isouf)
  226. ipt6 = imode1.imamod
  227. do iel = 1, ipt6.num(/2)
  228. ia = ipt6.num(2,iel)
  229. ib = icpr(ia)
  230. icomu(ib) = icomu(ib)+1
  231. imul(icomu(ib),ib) = ipt6.num(1,iel)
  232. enddo
  233. segdes,imode1,ipt6
  234. 502 continue
  235. segdes,mmode1
  236. segsup,icomu
  237. IPP2 = ipt5
  238. ENDIF
  239. C**** FIN DU PRE-TRAITEMENT DANS LE CAS D'UN MODELE ********************
  240.  
  241. C**** TRAVAIL **********************************************************
  242.  
  243. IPT2 = IPP2
  244. SEGACT,IPT2
  245. C--------------------------------------
  246. *---- cas de l'accrochage glissant ----
  247. IF (igliss.EQ.1) THEN
  248. * on va attacher a chaque point un vecteur tangent au cable
  249. * si on est arrive avec un modele, on prend le maillage du modele
  250. * en reference
  251. segini itange
  252. nbsz2 = ipt2.lisous(/1)
  253. meleme = ipt2
  254. do iob = 1, MAX(1,nbsz2)
  255. if (nbsz2.ne.0) then
  256. meleme = ipt2.lisous(i)
  257. segact,meleme
  258. endif
  259. if (meleme.itypel.ne.2) then
  260. call erreur(946)
  261. return
  262. endif
  263. do j = 1, meleme.num(/2)
  264. ia = num(1,j)
  265. ib = num(2,j)
  266. if (iextr(ia).gt.2.or.iextr(ib).gt.2) then
  267. call erreur(947)
  268. return
  269. endif
  270. ja = (ia-1)*idimp1
  271. jb = (ib-1)*idimp1
  272. xab = xcoor(jb+1) - xcoor(ja+1)
  273. yab = xcoor(jb+2) - xcoor(ja+2)
  274. zab = 0.D0
  275. if (idim.eq.3) zab = xcoor(jb+3) - xcoor(ja+3)
  276. xlo = sqrt(xab * xab + yab * yab + zab * zab)
  277. if (xlo.eq.0.D0) then
  278. call erreur(945)
  279. return
  280. endif
  281. xab = xab/xlo
  282. yab = yab/xlo
  283. zab = zab/xlo
  284. if (iextr(ia).eq.0) then
  285. xperp1(1,ia) = xab
  286. xperp1(2,ia) = yab
  287. xperp1(3,ia) = zab
  288. else
  289. xxx = 1.d0
  290. if (iextr(ia).eq.1) xxx = -1.d0
  291. zz1 = xab*xxx + xperp1(1,ia)
  292. zz2 = yab*xxx + xperp1(2,ia)
  293. zz3 = zab*xxx + xperp1(3,ia)
  294. xlo = sqrt(zz1*zz1 + zz2*zz2 + zz3*zz3)
  295. if (xlo.le.0.d0) then
  296. call erreur(945)
  297. return
  298. endif
  299. xperp1(1,ia) = zz1/xlo
  300. xperp1(2,ia) = zz2/xlo
  301. xperp1(3,ia) = zz3/xlo
  302. endif
  303. if (iextr(ib).eq.0) then
  304. xperp1(1,ib) = xab
  305. xperp1(2,ib) = yab
  306. xperp1(3,ib) = zab
  307. else
  308. xxx = 1.d0
  309. if (iextr(ia).eq.2) xxx = -1.d0
  310. zz1 = xab*xxx + xperp1(1,ib)
  311. zz2 = yab*xxx + xperp1(2,ib)
  312. zz3 = zab*xxx + xperp1(3,ib)
  313. xlo = sqrt(zz1*zz1 + zz2*zz2 + zz3*zz3)
  314. if (xlo.le.0.d0) then
  315. call erreur(945)
  316. return
  317. endif
  318. xperp1(1,ib) = zz1/xlo
  319. xperp1(2,ib) = zz2/xlo
  320. xperp1(3,ib) = zz3/xlo
  321. endif
  322. iextr(ia) = iextr(ia) + 1
  323. iextr(ib) = iextr(ib) + 2
  324. enddo
  325. if (nbsz2.ne.0) segdes,meleme
  326. enddo
  327. do 3659 j = 1, nbpti
  328. if (iextr(j).eq.0) goto 3659
  329. if (abs(xperp1(1,j)).gt.2.D-2 .or.
  330. & abs(xperp1(2,j)).gt.2.D-2) then
  331. xperp2(1,j) = -xperp1(2,j)
  332. xperp2(2,j) = xperp1(1,j)
  333. else
  334. xperp2(1,j) = -xperp1(3,j)
  335. xperp2(3,j) = xperp1(1,j)
  336. endif
  337. if (idim.eq.3) then
  338. xlo = sqrt ( xperp2(1,j)**2 + xperp2(2,j)**2
  339. $ + xperp2(3,j)**2 )
  340. xperp2(1,j) = xperp2(1,j) / xlo
  341. xperp2(2,j) = xperp2(2,j) / xlo
  342. xperp2(3,j) = xperp2(3,j) / xlo
  343. xaa = xperp1(3,j)*xperp2(2,j)-xperp1(2,j)*xperp2(3,j)
  344. xbb = xperp1(1,j)*xperp2(3,j)-xperp1(3,j)*xperp2(1,j)
  345. xcc = xperp1(2,j)*xperp2(1,j)-xperp1(1,j)*xperp2(2,j)
  346. xperp1(1,j) = xaa
  347. xperp1(2,j) = xbb
  348. xperp1(3,j) = xcc
  349. else
  350. xperp1(1,j) = xperp2(1,j)
  351. xperp1(2,j) = xperp2(2,j)
  352. endif
  353. 3659 continue
  354. ENDIF
  355. *---- fin du traitement specifique accrochage glissant ----
  356. C----------------------------------------------------------
  357.  
  358. C ON VA CHANGER LE MAILLAGE IPP2 A ACCROCHER EN ELEMENTS POI1
  359. C S'IL NE L EST PAS DEJA
  360. IPT2 = IPP2
  361. c* SEGACT,IPT2
  362. IF (IPT2.ITYPEL.NE.1) THEN
  363. CALL CHANGE(IPP2,1)
  364. SEGDES,IPT2
  365. IPT2 = IPP2
  366. SEGACT,IPT2
  367. ENDIF
  368. c On a maintenant : IPT2 = IPP2 qui est maillage de POI1 !
  369. NODES = IPT2.NUM(/2)
  370.  
  371. SEGINI,iexclu
  372. IEXX = iexclu
  373.  
  374. c On a toujours : ipt1 = MELEMI
  375. ipt1 = MELEMI
  376.  
  377. C On veut savoir s'il n'y a pas des noeuds de IPP2 qui ne sont pas
  378. C deja des noeuds d'un element de MELEMI, car pour RELA ACCRO, on n'a pas
  379. C besoin de les traiter (et de construire une matrice de rigidite).
  380. C Pour cela, pour tous les noeuds (ip) de MELEMI, idejvu(ip) = 2.
  381. meleme = ipt1
  382. DO iob = 1, NBSZ
  383. IF (NBSZ1.NE.0) THEN
  384. meleme = ipt1.LISOUS(iob)
  385. c* SEGACT,meleme
  386. ENDIF
  387. NBNN1 = meleme.num(/1)
  388. NBEL1 = meleme.num(/2)
  389. DO iem = 1, NBEL1
  390. DO ipt = 1, NBNN1
  391. ip = meleme.NUM(ipt,iem)
  392. iexclu.IDEJVU(ip) = 2
  393. ENDDO
  394. ENDDO
  395. c* IF (NBSZ1.NE.0) SEGDES,meleme
  396. ENDDO
  397.  
  398. SEGINI,icount
  399. ICC = icount
  400. icount.NBPTOT = 0
  401.  
  402. SEGINI,mtrav
  403.  
  404. C idejvu(ip) = 1 si le point ip dans un element de MELEMI
  405. C---- BOUCLE SUR LES SOUS ZONES ----
  406. meleme = ipt1
  407. ITR=0
  408. DO iob = 1, NBSZ
  409. IF (NBSZ1.NE.0) THEN
  410. meleme = ipt1.LISOUS(iob)
  411. c* SEGACT,meleme
  412. ENDIF
  413. C ICI LE ZONAGE (Les segments sont actifs en E/S)
  414. CALL ZONAG(iob,meleme,IPP2,ICC,IEXX,ITR)
  415. c* IF (NBSZ1.NE.0) SEGDES,meleme
  416. ENDDO
  417. C---- FIN DE LA BOUCLE SUR LES SOUS ZONES ----
  418.  
  419. IACCRO = icount.NBPTOT
  420.  
  421. C---- TABLEAU des POINT PAS VUS ----
  422. C-----------------------------------------------------------
  423. NBOUB = NODES - IACCRO
  424. * WRITE(ioimp,*) 'nombre de points exterieurs ' ,nboub
  425. SEGINI,ioubli
  426. iuo = 0
  427. DO j = 1, NODES
  428. ip = IPT2.NUM(1,j)
  429. IF (iexclu.IDEJVU(ip).EQ.0) THEN
  430. iuo = iuo + 1
  431. ioubli.IPTOU(iuo) = ip
  432. ENDIF
  433. ENDDO
  434. * WRITE(ioimp,*) 'nombre de points exterieurs ' ,nboub
  435. IF (iuo.EQ.0) GOTO 9000
  436.  
  437. C-----------------------------------------------------------
  438. C---- RATTRAPAGE EVENTUEL DE POINTS SUR LA LIMITE ----
  439. IPT1 = MELEMI
  440. c* SEGACT,IPT1
  441. krat = 0
  442. C---- boucle sur les sous zones ----
  443. meleme = IPT1
  444. DO iob = 1, NBSZ
  445. IF (NBSZ1.NE.0) THEN
  446. meleme = IPT1.LISOUS(iob)
  447. c* SEGACT,meleme
  448. ENDIF
  449. NBNN1 = meleme.num(/1)
  450. NBEL1 = meleme.num(/2)
  451. IELE = meleme.itypel
  452. kd = kdegre(IELE)
  453. C element quadratiques
  454. if (kd.eq.3) then
  455. inomi1 = 0
  456. nso = nbsom(IELE)
  457. iad = nspos(IELE)-1
  458. do 762 i = 1, NBNN1
  459. do j = 1, nso
  460. iso = ibsom(iad+j)
  461. if (i.eq.iso) goto 762
  462. enddo
  463. inomi1 = inomi1 + 1
  464. mtrav.inomil(inomi1) = i
  465. 762 continue
  466. C elements lineaires
  467. else if (kd.eq.2) then
  468. inomi1 = NBNN1
  469. do i = 1, NBNN1
  470. mtrav.inomil(i) = i
  471. enddo
  472. endif
  473.  
  474. C...... boucle sur les elements de la sz .......
  475. DO 1121 iem = 1, NBEL1
  476. CALL DOXE(xcoor,idim,NBNN1,meleme.num,iem,mtrav.xel)
  477.  
  478. C.......boucle sur les points candidats au rattrapage ......
  479. DO 1120 IPT = 1, IUO
  480.  
  481. ip = ioubli.IPTOU(ipt)
  482. IF (ip.LT.0) GOTO 1120
  483. c***c* Cas deja vu car correspond idejvu(ip) = 2 !
  484. c*** do jl = 1, NBNN1
  485. c*** if (ip.eq.meleme.num(jl,iem)) goto 1120
  486. c*** enddo
  487. IREFP = (IP-1)*IDIMP1
  488. xpo(1) = XCOOR(IREFP+1)
  489. xpo(2) = XCOOR(IREFP+2)
  490. xpo(3) = 0.D0
  491. IF (IDIM.EQ.3) xpo(3) = XCOOR(IREFP+3)
  492. C coordonnes intrinseques et fonctions de forme : qsijs
  493. C WRITE(6,5777 ) ((mtrav.xel(ia,ib),ia=1,idim),ib=1,nbnn1)
  494. C 5777 format(6e12.5)
  495. CALL QSIJS(mtrav.xel,IELE,NBNN1,idim,xpo(1),
  496. & mtrav.SHPP,mtrav.qsi,iretou)
  497. IF (iretou.NE.0) goto 1120
  498. DO i = 1, inomi1
  499. ilp = inomil(i)
  500. IF (mtrav.shpp(1,ilp).lt.CCRIT) goto 1120
  501. ENDDO
  502. C ce point exterieur est rattrapable
  503. krat = krat + 1
  504. ioubli.iptou(ipt) = -ip
  505. NNP = icount.NBPTOT + 1
  506. icount.IEINT(1,NNP) = ip
  507. icount.IEINT(2,NNP) = iem
  508. icount.IEINT(3,NNP) = iob
  509. DO i = 1, idim
  510. icount.QQQ(i,NNP) = mtrav.qsi(i)
  511. ENDDO
  512. icount.MAXPZ(iob) = icount.MAXPZ(iob) + 1
  513. icount.NBPTOT = NNP
  514. C WRITE(6,2375)ip,iem,(xpo(i1),i1=1,idim),(qsi(i2),i2=1,idim)
  515. 2375 format(2i4,6e12.5)
  516. 1120 continue
  517. C.......fin de la boucle sur les points candidats au rattrapage ......
  518. 1121 CONTINUE
  519. C...... fin de la boucle sur les elements enveloppe de la sz .......
  520. c* IF (NSZ1.NE.0) SEGDES,meleme
  521. ENDDO
  522. C---- fin de la boucle sur les zones ----
  523.  
  524. C comptage du nombre de points rattrapes
  525. k = 0
  526. do i = 1, IUO
  527. IF (ioubli.IPTOU(I).LT.0) k = k+1
  528. ENDDO
  529. if (k.ne.krat) then
  530. endif
  531.  
  532. * WRITE(6,*) ' nombre de points rattrapes ' ,k
  533. C on donne le nombre de points non projetés
  534. C if (k.ne.iuo) then
  535. C INTERR(1)=iuo-k
  536. C CALL ERREUR(-322)
  537. C endif
  538. IACCRO = IACCRO + k
  539.  
  540. C---- FIN DE LA SEQUENCE DE RATTRAPAGE
  541. C-----------------------------------------------------------
  542. 9000 CONTINUE
  543.  
  544. c MESSAGE : Nombre de points accrochés %i1 sur %i2 proposés
  545. INTERR(1) = IACCRO
  546. INTERR(2) = NODES
  547. CALL ERREUR(-319)
  548. C
  549. IPT1 = MELEMI
  550. C NOMBRE DE SOUS ZONES TOUCHÉES ET NOMBRE TOTAL DE NOEUDS CONCERNÉS
  551. NBSZR = 0
  552. ktot = 0
  553. DO iob = 1, NBSZ
  554. k = icount.MAXPZ(iob)
  555. ktot = ktot + k
  556. if (k.GT.0) NBSZR = NBSZR + 1
  557. ENDDO
  558. if (ktot.NE.icount.NBPTOT) then
  559. write(ioimp,*) 'pb ktot nbptot',ktot,icount.NBPTOT
  560. CALL ERREUR(5)
  561. endif
  562.  
  563. SEGACT,MCOORD*MOD
  564. C dans le cas IFROCA ne 0 les multiplicateurs existent deja
  565. IF (IFROCA.EQ.0) THEN
  566. IF (igliss.EQ.1) NINC = NINC - 1
  567. NBPTS = NBPTI + (icount.NBPTOT * NINC)
  568. SEGADJ,MCOORD
  569. ENDIF
  570.  
  571. C INITIALISATION DE MRIGID
  572. NRIGEL = NBSZR * NINC
  573. NRIGE = 8
  574. SEGINI,MRIGID
  575. mrigid.IFORIG = IFOMOD
  576. mrigid.MTYMAT ='RIGIDITE'
  577.  
  578. ipt1 = MELEMI
  579. meleme = ipt1
  580. ICZ = 0
  581.  
  582. C________________________________________________________________
  583. C____ BOUCLE SUR LES SOUS ZONES __________________________
  584. DO 400 IOB = 1, NBSZ
  585.  
  586. NBELEM = icount.MAXPZ(iob)
  587. * WRITE(ioimp,*) 'SOUS ZONE ',IOB ,'NOMBRE DE POINTS ',NBELEM
  588. IF (NBELEM.EQ.0) GOTO 395
  589.  
  590. IF (NBSZ1.NE.0) meleme = ipt1.lisous(iob)
  591. ICZ = ICZ + 1
  592.  
  593. KEL = meleme.ITYPEL
  594. NBMA1 = meleme.NUM(/1)
  595.  
  596. NHD = NBMA1+1
  597.  
  598. C Dimensions des maillage (IPT3), xmatri, DESCR
  599. NBNN = NBMA1 + 2
  600. NBSOUS = 0
  601. NBREF = 0
  602.  
  603. NELRIG = NBELEM
  604. NLIGRD = NBNN
  605. IF (igliss.EQ.1) NLIGRD = 1 + IDIM*(NBMA1+1)
  606. NLIGRP = NLIGRD
  607.  
  608. * WRITE(ioimp,*) ' ninc' , ninc
  609. DO iop = 1, NINC
  610. SEGINI,IPT3
  611. IPT3.ITYPEL = 22
  612. SEGINI,XMATRI
  613. XMATRI.SYMRE = 0
  614. IMEL(iop) = IPT3
  615. IMTT(iop) = XMATRI
  616. ENDDO
  617.  
  618. C.... boucle sur les elements .......
  619. xpo(3)= 0.D0
  620. ino = 0
  621. DO 300 K = 1, NODES
  622. if (icount.ieint(3,k).NE.iob) goto 300
  623. ino = ino + 1
  624. IP = icount.IEINT(1,K)
  625. IEL = icount.IEINT(2,K)
  626. irefp = (ip-1) * idimp1
  627. xpo(1) = xcoor(irefp+1)
  628. xpo(2) = xcoor(irefp+2)
  629. IF (IDIM.EQ.3) xpo(3) = xcoor(irefp+3)
  630. qsi(1) = icount.QQQ(1,k)
  631. qsi(2) = icount.QQQ(2,k)
  632. qsi(3) = icount.QQQ(3,k)
  633.  
  634. CALL SHAPE(QSI(1),QSI(2),QSI(3),KEL,SHPP,iret)
  635.  
  636. CCOE(1) = 1.D0
  637. NHD = NBMA1 + 1
  638. DO KM = 1, NBMA1
  639. CCOE(1+KM) = -SHPP(1,KM)
  640. ENDDO
  641.  
  642. DO 59 IK = 1, NINC
  643.  
  644. IF (IFROCA.EQ.0) THEN
  645. IDE = NBPTI * IDIMP1
  646. XCOOR(IDE+1) = xpo(1)
  647. XCOOR(IDE+2) = xpo(2)
  648. XCOOR(IDE+3) = xpo(3)
  649. NBPTI = NBPTI + 1
  650. ELSE
  651. ibp = icpr(ip)
  652. NBPTI = imul(ik,ibp)
  653. ENDIF
  654. C REMPLISSAGE DU MELEME IPT3
  655. ipt3 = imel(ik)
  656. ipt3.NUM(1,ino) = NBPTI
  657. ipt3.NUM(2,ino) = IP
  658. DO i = 1, NBMA1
  659. ipt3.NUM(2+i,ino) = meleme.NUM(i,IEL)
  660. ENDDO
  661.  
  662. C REMPLISSAGE DU XMATRI
  663. XMATRI = imtt(ik)
  664.  
  665. c Cas glissant
  666. IF (igliss.eq.1) THEN
  667. xaa = xperp1(1,ip)
  668. xbb = xperp1(2,ip)
  669. if (idim.eq.3) xcc = xperp1(3,ip)
  670. if (ik.eq.2) then
  671. xaa = xperp2(1,ip)
  672. xbb = xperp2(2,ip)
  673. if (idim.eq.3) xcc = xperp2(3,ip)
  674. endif
  675. IF (IFROCA.NE.0.AND.ik.EQ.ninc) then
  676. C petit calcul du vecteur tangent
  677. if( idim.eq.2) then
  678. XAA = xperp1(2,ip)
  679. XBB =-xperp1(1,ip)
  680. else
  681. XAA = xperp1(2,ip)*xperp2(3,ip)-xperp1(3,ip)*xperp2(2,ip)
  682. XBB = xperp1(3,ip)*xperp2(1,ip)-xperp1(1,ip)*xperp2(3,ip)
  683. XCC = xperp1(1,ip)*xperp2(2,ip)-xperp1(2,ip)*xperp2(1,ip)
  684. endif
  685. endif
  686. ikm = 1
  687. DO 337 km=1,nhd
  688. CCC = -CCOE(KM)
  689. ikm = ikm + 1
  690.  
  691. RE(1,ikm,INO) = -CCC*XAA
  692. RE(ikm,1,INO) = -CCC*XAA
  693. ikm = ikm + 1
  694. RE(1,ikm,INO) = -CCC*XBB
  695. RE(ikm,1,INO) = -CCC*XBB
  696. IF (IDIM.EQ.3) THEN
  697. ikm = ikm + 1
  698. RE(1,ikm,INO) = -CCC*XCC
  699. RE(ikm,1,INO) = -CCC*XCC
  700. ENDIF
  701. 337 CONTINUE
  702.  
  703. c Cas non glissant
  704. ELSE
  705. DO KM = 1, NHD
  706. CCC = CCOE(KM)
  707. RE(1,KM+1,INO) = CCC
  708. RE(1+KM,1,INO) = CCC
  709. ENDDO
  710. ENDIF
  711. 59 continue
  712. c - fin de boucle sur les xmatri -
  713. 300 CONTINUE
  714. C.... fin de boucle sur les elements .......
  715. if (ino.ne.nbelem) then
  716. write(ioimp,*) 'Pb element zone ',iob,ino,nbelem
  717. call erreur(5)
  718. endif
  719.  
  720. INRI = NINC*(ICZ-1) + 1
  721.  
  722. DO iop = 1, NINC
  723.  
  724. SEGINI DESCR
  725. descr.LISINC(1) = 'LX'
  726. descr.LISDUA(1) = 'FLX'
  727. descr.NOELEP(1) = 1
  728. descr.NOELED(1) = 1
  729. c - Cas non glissant -
  730. IF (igliss.eq.0) then
  731. DO k = 2, NBNN
  732. descr.LISINC(k) = NOMDD(ICOR(iop))
  733. descr.LISDUA(k) = NOMDU(ICOR(iop))
  734. descr.NOELEP(k) = k
  735. descr.NOELED(k) = k
  736. ENDDO
  737. c - Cas glissant -
  738. ELSE
  739. DO km1 = 1, nhd
  740. j = 1+(km1-1)*idim
  741. descr.lisinc(j+1)= NOMDD(ICOR(1))
  742. descr.lisinc(j+2)= NOMDD(ICOR(2))
  743. descr.LISDUA(j+1)= NOMDU(ICOR(1))
  744. descr.LISDUA(j+2)= NOMDU(ICOR(2))
  745. descr.NOELEP(j+1)= KM1+1
  746. descr.NOELEP(j+2)= KM1+1
  747. descr.NOELED(j+1)= KM1+1
  748. descr.NOELED(j+2)= KM1+1
  749. if (idim.eq.3) then
  750. descr.lisinc(j +3)= NOMDD(ICOR(3))
  751. descr.LISDUA(j +3)= NOMDU(ICOR(3))
  752. descr.NOELEP(j +3)= KM1+1
  753. descr.NOELED(j +3)= KM1+1
  754. endif
  755. ENDDO
  756. endif
  757. SEGDES,DESCR
  758. IPT3 = IMEL(iop)
  759. SEGDES,IPT3
  760. XMATRI = IMTT(iop)
  761. SEGDES,xMATRI
  762.  
  763. jop = INRI + iop - 1
  764. IRIGEL(1,jop) = IPT3
  765. c IRIGEL(2,jop) = 0
  766. IRIGEL(3,jop) = DESCR
  767. IRIGEL(4,jop) = XMATRI
  768. IRIGEL(5,jop) = NIFOUR
  769. c IRIGEL(6,jop) = 0
  770. IF (IFROCA.NE.0.AND.iop.EQ.ninc) IRIGEL(6,jop) = 2
  771. c IRIGEL(7,jop) = 0
  772. c IRIGEL(8,jop) = 0
  773. COERIG(jop) = 1.D0
  774. ENDDO
  775.  
  776. 395 CONTINUE
  777. IF (NBSZ1.NE.0) SEGDES,meleme
  778.  
  779. 400 CONTINUE
  780. C____ FIN DE BOUCLE SUR LES SOUS ZONES __________________________
  781.  
  782. if (icz.ne.nbszr) then
  783. write(ioimp,*) 'Pb Rigi ',ICZ,NBSZR
  784. call erreur(5)
  785. endif
  786.  
  787. C**** MENAGE ***********************************************************
  788. IF (IFROCA.EQ.0) THEN
  789. SEGACT,MCOORD*NOMOD
  790. ENDIF
  791. SEGSUP,icount,iexclu,mtrav
  792. SEGDES,IPT1,IPT2
  793. SEGDES,MRIGID
  794.  
  795. C**** ECRITURE DE LA RIGIDTE RESULTAT **********************************
  796. CALL ECROBJ('RIGIDITE',MRIGID)
  797.  
  798. RETURN
  799. END
  800.  
  801.  
  802.  
  803.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales