Télécharger accro.eso

Retour à la liste

Numérotation des lignes :

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

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