Télécharger accro3.eso

Retour à la liste

Numérotation des lignes :

accro3
  1. C ACCRO3 SOURCE PV090527 26/04/30 21:15:02 12529
  2. SUBROUTINE ACCRO3
  3. C========================================================================
  4. C Cree la matrice de liaison entre le champ u et w
  5. C avec une formulation forte => On utilise les fonctions de forme
  6. C et on identifie le deplacement aux noeuds de la fissure
  7. c
  8. c On prend en compte les enrichissements XFEM
  9. c
  10. C Creation : BP, decembre 2012
  11. C Modifications : ...
  12. C
  13. C========================================================================
  14.  
  15. IMPLICIT INTEGER(I-N)
  16. IMPLICIT REAL*8 (A-H,O-Z)
  17.  
  18. -INC PPARAM
  19. -INC CCOPTIO
  20. -INC CCREEL
  21. -INC CCHAMP
  22. -INC CCGEOME
  23.  
  24. -INC SMCOORD
  25. -INC SMELEME
  26. -INC SMRIGID
  27. -INC SMLMOTS
  28. -INC SMMODEL
  29. -INC SMINTE
  30. -INC SMCHAML
  31. -INC SMLREEL
  32.  
  33. POINTEUR MCHEX1.MCHELM
  34. external shape
  35. DATA EPSI/1.D-9/
  36. DIMENSION ICOR(6),IMEL(6),imtt(10)
  37. DIMENSION qsi(3),XPO(3)
  38.  
  39. C INITIALISATION DES INCONNUES obligatoires et facultatives
  40. PARAMETER (NOBL=3,NFAC=9)
  41. CHARACTER*4 DDLOBL(NOBL),DDLFAC(NFAC),MODDL,MODDL2
  42. CHARACTER*4 DUAOBL(NOBL),DUAFAC(NFAC)
  43.  
  44. DATA DDLOBL/'UX ','UY ','UZ '/
  45. DATA DDLFAC/'AX ','AY ','AZ ',
  46. >'B1X ','B1Y ','B1Z ',
  47. >'B2X ','B2Y ','B2Z '/
  48. DATA DUAOBL/'FX ','FY ','FZ '/
  49. DATA DUAFAC/'FAX ','FAY ','FAZ ',
  50. >'FB1X','FB1Y','FB1Z',
  51. >'FB2X','FB2Y','FB2Z'/
  52.  
  53.  
  54. c-----------------------------------------------------
  55. c Segment de travail
  56. c-----------------------------------------------------
  57. SEGMENT MTRAV
  58. REAL*8 XE2(3,NBNN2)
  59. REAL*8 SHPP1(6,NBNN1)
  60. REAL*8 XE1(3,NBNN1)
  61. INTEGER IDEJVU(NBPTS)
  62. ENDSEGMENT
  63.  
  64. c tableaux comptant le nbre d EF de chaque ddl
  65. PARAMETER(NDDLMAX=6)
  66. INTEGER NELDDL(NDDLMAX)
  67.  
  68. if(iimpi.ge.2) then
  69. write(ioimp,*) '-----------------------------------------------'
  70. write(ioimp,*) ' ENTREE dans ACCRO3'
  71. write(ioimp,*) '-----------------------------------------------'
  72. endif
  73.  
  74. c Preliminaires
  75. idim1 = idim + 1
  76. SEGACT,MCOORD*mod
  77.  
  78. c-----------------------------------------------------
  79. c RECUPERATION DU MAILLAGE MASSIF
  80. c-----------------------------------------------------
  81. CALL LIROBJ('MMODEL ',IPMODL,1,IRETOU)
  82. CALL ACTOBJ('MMODEL ',IPMODL,1)
  83. IF(IERR.NE.0) RETURN
  84. MMODE1=IPMODL
  85. segact,MMODE1
  86. c Récupération du nombre de zones du modèle
  87. N1 = MMODE1.KMODEL(/1)
  88. if(N1.gt.1) write(ioimp,*) 'attention 1 seule zone a ce jour!'
  89. IMODE1 = MMODE1.KMODEL(1)
  90. segdes,MMODE1
  91. C Récupération du maillage et du numéro d'élément du modèle
  92. segact,IMODE1
  93. nele1 = IMODE1.NEFMOD
  94. IPT1 = IMODE1.IMAMOD
  95. SEGACT IPT1
  96. C Récupération du numéro d'élément du maillage, du nombre de noeuds et d'éléments
  97. iele1 = IPT1.itypel
  98. nbnn1 = IPT1.num(/1)
  99. nbel1 = IPT1.num(/2)
  100. mele1 = IMODE1.INFELE(1)
  101.  
  102. c-----------------------------------------------------
  103. c RECUPERATION DU MAILLAGE INTERFACE
  104. c-----------------------------------------------------
  105. IPMAI2 = 0
  106. CALL LIROBJ('MMODEL ',IPMODL,0,IRETOU)
  107. IF(IERR.NE.0) RETURN
  108. if (IRETOU.EQ.1)then
  109. CALL ACTOBJ('MMODEL ',IPMODL,1)
  110. MMODE2=IPMODL
  111. segact,MMODE2
  112. N2 = MMODE2.KMODEL(/1)
  113. if(N2.gt.1) write(ioimp,*) 'attention 1 seule zone a ce jour!'
  114. IMODE2 = MMODE2.KMODEL(1)
  115. segdes,MMODE2
  116. segact,IMODE2
  117. nele2 = IMODE2.NEFMOD
  118. IPT2 = IMODE2.IMAMOD
  119. SEGACT IPT2
  120. c pour l'instant on dit que nele = iele (marche pour iele entre 2 et 26, voir bdata.eso)
  121. iele2 = IPT2.itypel
  122. nbnn2 = IPT2.num(/1)
  123. nbel2 = IPT2.num(/2)
  124. ngau2 = IMODE2.INFELE(6)
  125. segdes,IMODE2
  126. else
  127. C Dans le cas où on a un maillage en entrée
  128. CALL LIROBJ('MAILLAGE',IPMAI2,1,IRETOU)
  129. IF(IERR.NE.0) RETURN
  130. IPT2 = IPMAI2
  131. SEGACT IPT2
  132. c pour l'instant on dit que nele = iele (marche pour iele entre 2 et 26, voir bdata.eso)
  133. iele2 = IPT2.itypel
  134. nele2 = iele2
  135. if (nele2.lt.2.or.nele2.gt.26) then
  136. write(ioimp,*)'element geometrique different de l element fini'
  137. call erreur(16)
  138. endif
  139. nbnn2 = IPT2.num(/1)
  140. nbel2 = IPT2.num(/2)
  141. c SEG2
  142. if (nele2.EQ.2) ngau2 = 2
  143. c SEG3
  144. if (nele2.EQ.3) ngau2 = 3
  145. c TRI3
  146. if (nele2.EQ.4) ngau2 = 1
  147. c TRI6
  148. if (nele2.EQ.6) ngau2 = 4
  149. c QUA4
  150. if (nele2.EQ.8) ngau2 = 4
  151. c QUA8
  152. if (nele2.EQ.10) ngau2 = 9
  153. endif
  154.  
  155. call RESHPT(ngau2,nbnn2,iele2,nele2,0,IPTR2,IRET)
  156. MINTE2 = IPTR2
  157. segact,MINTE2
  158.  
  159.  
  160. c-----------------------------------------------------
  161. c RECHERCHE DU MCHAML ISSU MCHEX1 D ENRICHISSEMENT
  162. c-----------------------------------------------------
  163. MCHAM1=0
  164. NBENR2=0
  165. segact,IMODE1
  166. NOBMOD = IMODE1.IVAMOD(/1)
  167. IF (NOBMOD.NE.0) THEN
  168. DO 1002 iobmo1=1,NOBMOD
  169. if((IMODE1.TYMODE(iobmo1)).ne.'MCHAML') goto 1002
  170. MCHEX1 = IMODE1.IVAMOD(iobmo1)
  171. segact,MCHEX1
  172. if((MCHEX1.TITCHE).ne.'ENRICHIS') goto 1003
  173. MCHAM1 = MCHEX1.ICHAML(1)
  174. segact,MCHAM1
  175. NBENR2 = MCHAM1.IELVAL(/1)
  176. do ienr2=1,NBENR2
  177. MELVA1=MCHAM1.IELVAL(IENR2)
  178. if(MELVA1.ne.0) segact,MELVA1
  179. enddo
  180. 1003 continue
  181. segdes,MCHEX1
  182. 1002 CONTINUE
  183. ENDIF
  184.  
  185. c-------------------------------------
  186. c INITIALISATION DES OBJETS DE TRAVAIL : MTRAV et ITYMAT
  187. c-------------------------------------
  188. segini,MTRAV
  189. *bp : NTYMAT = (U ou H ou HB1 ou HB1B2)
  190. * nbre de types de matrices = NRIGEL = NTYMAT * idim
  191. c NTYMAT = 4
  192. NTYMAT = 1+NBENR2
  193.  
  194. c-------------------------------------
  195. c INITIALISATION DU MCOORD
  196. c-------------------------------------
  197. NBPTS0 = NBPTS
  198. NBPTS = NBPTS0 + ((nbel2 + 1)*idim*3)
  199. SEGADJ,MCOORD
  200. NBPTS = NBPTS0
  201. * on compte le vrai nombre de LX ajouté via NBPTS
  202.  
  203. c--------------------------------------------------------------------
  204. c INITIALISATION DU SEGMENT MRIGID
  205. c--------------------------------------------------------------------
  206. NRIGEL = NTYMAT*idim
  207. segini,MRIGID
  208. IFORIG = IFOUR
  209. MTYMAT ='RIGIDITE'
  210. c -on prepare le meleme
  211. NBSOUS = 0
  212. NBREF = 0
  213.  
  214. ityty=0
  215. c on initialise la taille matrice en fonction du type de matrice
  216. do ity=1,NTYMAT
  217. do iidim=1,idim
  218.  
  219. ityty=ityty+1
  220. COERIG(ityty) = 1.D0
  221.  
  222. * dim de la matrice RE elementaire
  223. * = 1 (LX) + 1 (UX) + nbnoeud (UX)
  224. * = 1 (LX) + 1 (AX) + nbnoeud (AX)
  225. * = 1 (LX) + 1 (AX) + 2*nbnoeud (AX+B1X)
  226. * = 1 (LX) + 1 (AX) + 3*nbnoeud (AX+B1X+B2X)
  227. nbno1 = nbnn1
  228. nenr1 = ity
  229. if(nenr1.le.2) then
  230. NLIGRP = 1 + 1 + nbno1
  231. else
  232. NLIGRP = 1 + 1 + ((nenr1-1)*nbno1)
  233. endif
  234. NLIGRD = NLIGRP
  235.  
  236. c -creation du MELEME
  237. NBNN = 2 + nbno1
  238. NBELEM=0
  239. SEGINI,MELEME
  240. c ITYPEL=28
  241. ITYPEL=22
  242. IRIGEL(1,ityty) = MELEME
  243.  
  244. c -remplissage du DESCR
  245. SEGINI,DESCR
  246. IRIGEL(3,ityty) = DESCR
  247. iddl=0
  248. c remplissage des ddl LX de la fissure
  249. iddl=iddl+1
  250. LISINC(iddl)='LX'
  251. LISDUA(iddl)='FLX'
  252. NOELEP(iddl)=1
  253. NOELED(iddl)=1
  254. c remplissage des ddl UX de la fissure (ici appelé WX) []*WX=TX
  255. iddl=iddl+1
  256. if (nenr1.eq.1) then
  257. LISINC(iddl)=DDLOBL(iidim)
  258. LISDUA(iddl)=DUAOBL(iidim)
  259. else
  260. LISINC(iddl)=DDLFAC(iidim)
  261. LISDUA(iddl)=DUAFAC(iidim)
  262. endif
  263. NOELEP(iddl)=2
  264. NOELED(iddl)=2
  265. c remplissage des ddl de la structure
  266. if (nenr1.eq.1) then
  267. do ino1=1,nbno1
  268. iddl=iddl+1
  269. LISINC(iddl)=DDLOBL(iidim)
  270. LISDUA(iddl)=DUAOBL(iidim)
  271. NOELEP(iddl)=2+ino1
  272. NOELED(iddl)=2+ino1
  273. enddo
  274. else
  275. do ini1=1,(nenr1-1)
  276. do ino1=1,nbno1
  277. iddl=iddl+1
  278. LISINC(iddl)=DDLFAC(iidim+(3*(ini1-1)))
  279. LISDUA(iddl)=DUAFAC(iidim+(3*(ini1-1)))
  280. NOELEP(iddl)=2+ino1
  281. NOELED(iddl)=2+ino1
  282. enddo
  283. enddo
  284. endif
  285. if(iimpi.ge.3) write(ioimp,*) ityty,(LISINC(iou),iou=1,NLIGRP)
  286. if(iimpi.ge.3) write(ioimp,*) ityty,(NOELEP(iou),iou=1,NLIGRP)
  287. SEGDES,DESCR
  288.  
  289. c -initialisation du XMATRI
  290. NELRIG=0
  291. rigrel=0
  292. SEGINI,XMATRI
  293. IRIGEL(4,ityty) = XMATRI
  294. IRIGEL(5,ityty) = NIFOUR
  295. IRIGEL(6,ityty) = 0
  296. IRIGEL(7,ityty) = 0
  297. IRIGEL(8,ityty) = 0
  298.  
  299. enddo
  300. enddo
  301.  
  302.  
  303. c----------------------------------------------------------------------
  304. c 1. RECHERCHE DES ELEMENTS DE STRUCTURE CONTENANT DES POINTS DE GAUSS
  305. c DES ELEMENTS DE LA FISSURE
  306. c 2. REMPLISSAGE DU MRIGID (XMATRI et MELEME)
  307. c----------------------------------------------------------------------
  308. iaccro=0
  309. NODES=0
  310. C
  311. c==== Boucle sur les elements de fissure ==============================
  312. DO 1100 iem2=1,nbel2
  313.  
  314. call doxe(xcoor,idim,nbnn2,ipt2.num,iem2,xe2)
  315. nbenrj = 0
  316.  
  317. c======= Boucle sur les noeuds de fissure ============================
  318. DO 1132 ino2=1,nbnn2
  319.  
  320. c on n'attache qu'une seule fois chaque noeud
  321. inode2 = IPT2.NUM(ino2,iem2)
  322. if(IDEJVU(inode2).ne.0) goto 1132
  323. IDEJVU(inode2)=1
  324. NODES=NODES+1
  325.  
  326. c récupération des coordonnees du point de gauss dans le repère global
  327. XPO(1) = xe2(1,ino2)
  328. XPO(2) = xe2(2,ino2)
  329. XPO(3) = xe2(3,ino2)
  330.  
  331. c---------- Boucle sur les elements de structure ----------------------
  332. DO 1131 iem1=1,nbel1
  333.  
  334. c si pas d'enrichissement, on travaille sur tous les elements
  335. if(MCHAM1.eq.0) goto 1133
  336. c on saute les elements non enrichi car a priori ne contiennent pas la fissure
  337. do ienr2=1,NBENR2
  338. MELVA1=MCHAM1.IELVAL(IENR2)
  339. if(MELVA1.ne.0) then
  340. do inode1=1,nbnn1
  341. if(MELVA1.IELCHE(inode1,iem1).ne.0) goto 1133
  342. enddo
  343. endif
  344. enddo
  345. goto 1131
  346. 1133 continue
  347.  
  348. c recuperation des coordonnées des noeuds de IPT1 : xe1 (dans le repère x,y,z)
  349. call doxe(xcoor,idim,nbnn1,ipt1.num,iem1,xe1)
  350.  
  351. qsi(1) = 0.D0
  352. qsi(2) = 0.D0
  353. qsi(3) = 0.D0
  354. c calcul des fonctions de formes de IPT1 au pt de Gauss de IPT2
  355. call QSIJS(xe1,iele1,nbnn1,idim,XPO,SHPP1,qsi,IRET)
  356.  
  357. c test pour savoir si PG est dans EF de IPT1
  358. DO ino1=1,NBNN1
  359. if (SHPP1(1,ino1).LT.-1.01D-7) then
  360. go to 1131
  361. endif
  362. ENDDO
  363. c ON a trouvé : l'iem1 élément de structure contient ce noeud de fissure
  364. IDEJVU(inode2)=10
  365. iaccro=iaccro+1
  366.  
  367. c DETECTION DU TYPE D'ENRICHISSEMENT MAX DE CET ELEMENT = nbenrj
  368. DO 3001 IENR2=1,NBENR2
  369. MELVA1=MCHAM1.IELVAL(IENR2)
  370. IF(MELVA1.eq.0) GOTO 3001
  371. DO 3002 ino1=1,nbnn1
  372. MLREEL = MELVA1.IELCHE(ino1,iem1)
  373. c Test pour savoir si le noeud est enrichi
  374. IF(MLREEL.eq.0) GOTO 3002
  375. nbenrj=max(nbenrj,IENR2)
  376. 3002 continue
  377. 3001 continue
  378.  
  379. if(iimpi.ge.3) write(ioimp,*) 'EF fissure ',iem2,
  380. & ' ptdeG ',ino2,' -> EF MASSIF ',iem1,' nbenrj=',nbenrj
  381.  
  382.  
  383. c Remplissage du MRIGID
  384.  
  385. c ---Boucle sur la Partie standard et enrichie ---
  386. do 6000 ity=1,min(2,NTYMAT)
  387. c ---Boucle sur la dimension ---
  388. c rem : utile uniquement pour le meleme
  389. c (on pourrait garder le meme xmatri pour les differents iidim)
  390. do 6001 iidim=1,idim
  391.  
  392. if(ity.eq.1) then
  393. ityty = iidim
  394. else
  395. ityty = (nbenrj*idim) + iidim
  396. endif
  397. MELEME = IRIGEL(1,ityty)
  398. XMATRI = IRIGEL(4,ityty)
  399. NBELEM = NUM(/2)+1
  400. NLIGRD = RE(/1)
  401. NLIGRP = RE(/2)
  402. NELRIG = RE(/3)+1
  403. segadj,MELEME
  404. rigrel=0
  405. segadj,XMATRI
  406.  
  407. c Remplissage du MELEME
  408. c traitement du LX
  409. NBPTS = NBPTS + 1
  410. if(NBPTS.gt.(XCOOR(/1)/idim1)) then
  411. NBPTS0 = NBPTS
  412. NBPTS = NBPTS0 + (nbel2 + 1)
  413. SEGADJ,MCOORD
  414. NBPTS = NBPTS0
  415. endif
  416. NUM(1,NBELEM) = NBPTS
  417. XCOOR((NBPTS-1)*idim1 +1) = XCOOR((inode2-1)*idim1 +1)
  418. XCOOR((NBPTS-1)*idim1 +2) = XCOOR((inode2-1)*idim1 +2)
  419. if(idim.eq.3)
  420. & XCOOR((NBPTS-1)*idim1 +3) = XCOOR((inode2-1)*idim1 +3)
  421. c traitement du noeud de la fissure
  422. NUM(2,NBELEM) = inode2
  423. c traitement des noeuds de la structure
  424. inono=2
  425. cbp inutile do j2=1,max(1,nbenrj)
  426. cbp inutile car noelep boucle deja sur les enrichissement
  427. do ino1=1,nbnn1
  428. inono=inono+1
  429. NUM(inono,NBELEM) = IPT1.NUM(ino1,iem1)
  430. enddo
  431. c enddo
  432.  
  433. c Remplissage du XMATRI
  434. c traitement du terme LX - ddl fissure
  435. RE(1,2,NELRIG)=1.d0
  436. RE(2,1,NELRIG)=1.d0
  437. c traitement des terme LX - ddl structure
  438. inono=2
  439. c UX = UX seulement
  440. if(ity.eq.1) then
  441. do ino1=1,nbnn1
  442. inono=inono+1
  443. RE(1,inono,NELRIG)=-1.d0*SHPP1(1,ino1)
  444. RE(inono,1,NELRIG)=-1.d0*SHPP1(1,ino1)
  445. enddo
  446. else
  447. c AX = AX ...
  448. if(nbenrj.ge.1) then
  449. MELVA1 = MCHAM1.IELVAL(1)
  450. do ino1=1,nbnn1
  451. MLREEL = MELVA1.IELCHE(ino1,iem1)
  452. inono=inono+1
  453. if(MLREEL.ne.0) then
  454. RE(1,inono,NELRIG)=-1.d0*SHPP1(1,ino1)
  455. RE(inono,1,NELRIG)=-1.d0*SHPP1(1,ino1)
  456. endif
  457. enddo
  458. endif
  459. c ... + B1X + B2X
  460. if(nbenrj.ge.2) then
  461. do jenrj=2,nbenrj
  462. c on ecrit B1X pour les nbnn1 noeuds, puis B2X ...
  463. MELVA1 = MCHAM1.IELVAL(jenrj)
  464. do ino1=1,nbnn1
  465. MLREEL = MELVA1.IELCHE(ino1,iem1)
  466. if(MLREEL.eq.0) then
  467. c pas d'enrichissement => on met 0
  468. RX05 = 0.d0
  469. else
  470. SEGACT,MLREEL
  471. PSIX = 0.d0
  472. do iii0 = 1,nbnn1
  473. PSIX = PSIX + (SHPP1(1,iii0) * PROG(iii0))
  474. enddo
  475. SEGDES,MLREEL
  476. RX05= -1.d0*SQRT(ABS(PSIX))
  477. endif
  478. inono=inono+1
  479. RE(1,inono,NELRIG)=RX05*SHPP1(1,ino1)
  480. RE(inono,1,NELRIG)=RX05*SHPP1(1,ino1)
  481. enddo
  482. enddo
  483. endif
  484. endif
  485.  
  486. 6001 continue
  487. 6000 continue
  488. c ---fin de Boucle sur les XMATRI et MELEME
  489. c des Partie standard et enrichie et sur les dimensions ---
  490.  
  491. c on a trouvé le iem1 élément et on a fait le travail :
  492. c on passe au point de gauss suivant
  493. goto 1132
  494.  
  495. 1131 CONTINUE
  496. c---------- fin de la Boucle sur les elements de structure -------------
  497. if(iimpi.ge.2) write(ioimp,*)
  498. 1 'Fin de la boucle sur les elements de structure'
  499.  
  500. if(IDEJVU(inode2).ne.10) then
  501. write(ioimp,*) 'Attention le noeud ',inode2,' est hors support'
  502. endif
  503.  
  504. 1132 CONTINUE
  505. c======= fin de la Boucle sur les noeuds de fissure ===================
  506.  
  507. 1100 CONTINUE
  508. c==== fin de la Boucle sur les elements de fissure =====================
  509.  
  510. c MESSAGE : Nombre de points accrochés %i1 sur %i2 proposés
  511. INTERR(1)=iaccro
  512. INTERR(2)=NODES
  513. CALL ERREUR(-319)
  514.  
  515. c------------------------------------
  516. c AJUSTEMENT AVANT DE QUITTER
  517. c------------------------------------
  518.  
  519. if(iimpi.ge.2) write(ioimp,*) 'AJUSTEMENT AVANT DE QUITTER'
  520. c NBPTS est deja le vrai nomber de noeuds
  521. SEGADJ,MCOORD
  522.  
  523. C BOUCLE SUR LES SOUS RIGIDITÉS
  524. ityok = 0
  525. DO 2000 ityty=1,(idim*NTYMAT)
  526. MELEME = IRIGEL(1,ityty)
  527. DESCR = IRIGEL(3,ityty)
  528. XMATRI = IRIGEL(4,ityty)
  529. NBELEM = NUM(/2)
  530. if(NBELEM.ne.0) then
  531. ityok = ityok + 1
  532. IRIGEL(1,ityok) = MELEME
  533. IRIGEL(3,ityok) = DESCR
  534. IRIGEL(4,ityok) = XMATRI
  535. segdes,MELEME,XMATRI
  536. else
  537. segsup,MELEME,DESCR,XMATRI
  538. endif
  539. 2000 continue
  540. NRIGEL = ityok
  541. segadj,MRIGID
  542.  
  543. c-------------------------------
  544. c MENAGE AVANT DE QUITTER
  545. c-------------------------------
  546. segsup,MTRAV
  547. segdes IPT1,IPT2
  548. if(MCHAM1.ne.0) segdes,MCHAM1
  549. SEGDES,MRIGID
  550.  
  551. if(iimpi.ge.3) write(ioimp,*) 'ecriture du MRIGID',MRIGID
  552. CALL ECROBJ('RIGIDITE',MRIGID)
  553.  
  554. END
  555.  
  556.  
  557.  
  558.  
  559.  
  560.  

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