Télécharger accro.eso

Retour à la liste

Numérotation des lignes :

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

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