Télécharger accro2.eso

Retour à la liste

Numérotation des lignes :

  1. C ACCRO2 SOURCE BP208322 16/11/18 21:15:01 9177
  2. SUBROUTINE ACCRO2
  3. C========================================================================
  4. C Cree la matrice de liaison entre le champ u et w
  5. C avec une formulation faible a 3 champs : u w et t
  6. C FAIBLE => On utilise les fonctions de forme
  7. c
  8. C Option PENA et STAB pour créer les matrices de penalite Kww
  9. C et de stabilisation Ktt (utile si LATIN utilisée)
  10. c
  11. C En sortie, on a Kfaibl = [Ktt Ktu Ktw ; Kut 0 0 ; Kwt 0 Kww]
  12. C [0 Ktu Ktw ; sym. ] = \int t * (u-w) ds
  13. c
  14. c ajout des enrichissements XFEM
  15. c
  16. C Creation : Benoit Trolle (SNCF-Lamcos), avril 2012
  17. C Integration dans Cast3m : BP, decembre 2012
  18. C========================================================================
  19.  
  20. IMPLICIT INTEGER(I-N)
  21. IMPLICIT REAL*8 (A-H,O-Z)
  22.  
  23. -INC CCREEL
  24. -INC CCOPTIO
  25. -INC SMCOORD
  26. -INC SMELEME
  27. -INC SMRIGID
  28. -INC SMLMOTS
  29. -INC CCHAMP
  30. -INC CCGEOME
  31. -INC SMMODEL
  32. -INC SMINTE
  33. -INC SMCHAML
  34. -INC SMLREEL
  35.  
  36.  
  37. C Segment contenant les informations sur un element
  38. SEGMENT INFO
  39. INTEGER INFELL(JG)
  40. ENDSEGMENT
  41.  
  42. POINTEUR INFO1.INFO,INFO2.INFO
  43. POINTEUR MCHEX1.MCHELM
  44. external shape
  45. DATA EPSI/1.D-9/
  46. DIMENSION ICOR(6),IMEL(6),imtt(10)
  47. DIMENSION QSI(3),XPO(3)
  48.  
  49. C INITIALISATION DES INCONNUES obligatoires et facultatives
  50. PARAMETER (NOBL=3,NFAC=9,NMOCLE=3)
  51. CHARACTER*4 DDLOBL(NOBL),DDLFAC(NFAC),MODDL,MODDL2
  52. CHARACTER*4 DUAOBL(NOBL),DUAFAC(NFAC),MOCLE(NMOCLE)
  53. CHARACTER*4 DDLOB2(NOBL),DDLFA2(NOBL)
  54. CHARACTER*4 DUAOB2(NOBL),DUAFA2(NOBL)
  55.  
  56. DATA DDLOBL/'UX ','UY ','UZ '/
  57. DATA DDLFAC/'AX ','AY ','AZ ',
  58. >'B1X ','B1Y ','B1Z ',
  59. >'B2X ','B2Y ','B2Z '/
  60. DATA DUAOBL/'FX ','FY ','FZ '/
  61. DATA DUAFAC/'FAX ','FAY ','FAZ ',
  62. >'FB1X','FB1Y','FB1Z',
  63. >'FB2X','FB2Y','FB2Z'/
  64.  
  65. c nom de composantes speciales
  66. DATA DDLOB2/'UX ','UY ','UZ '/
  67. DATA DDLFA2/'AX ','AY ','AZ '/
  68. DATA DUAOB2/'FX ','FY ','FZ '/
  69. DATA DUAFA2/'FAX ','FAY ','FAZ '/
  70.  
  71. DATA MOCLE /'PENA','STAB','NOMU'/
  72.  
  73. c-----------------------------------------------------
  74. c Segment de travail
  75. c-----------------------------------------------------
  76. SEGMENT MTRAV
  77. REAL*8 SHPP2(6,NBNN2)
  78. REAL*8 XE2(3,NBNN2),XEL2(3,NBNN2),NGENE2(1,NBNN2)
  79. REAL*8 XE3(3,NBNN2)
  80. REAL*8 SHPP1(6,NBNN1)
  81. REAL*8 XE1(3,NBNN1),XEL1(3,NBNN1)
  82. REAL*8 NGENE1(1,NBNN1)
  83. REAL*8 REL(nbnn2,nbnn2),REL1(nbnn2,nbnn1)
  84. REAL*8 RINT(nbnn2,nbnn2)
  85. REAL*8 RINT0(nbnn2,ngau2*nbnn1),RINT1(nbnn2,ngau2*nbnn1)
  86. REAL*8 RINT2(nbnn2,ngau2*nbnn1),RINT3(nbnn2,ngau2*nbnn1)
  87. INTEGER ITYMAT(2,4*nbnn1*ngau2)
  88. INTEGER ITAMYT(ngau2*NBNN1,4)
  89. INTEGER PRELEM(3*nbnn1*ngau2)
  90. INTEGER ElemPG(ngau2)
  91. ENDSEGMENT
  92.  
  93. c tableaux comptant le nbre d EF de chaque ddl
  94. PARAMETER(NDDLMAX=6)
  95. INTEGER NELDDL(NDDLMAX)
  96.  
  97. if(iimpi.ge.2) then
  98. write(ioimp,*) '-----------------------------------------------'
  99. write(ioimp,*) ' ENTREE dans ACCRO2'
  100. write(ioimp,*) '-----------------------------------------------'
  101. endif
  102.  
  103. SEGACT MCOORD
  104.  
  105. c-----------------------------------------------------
  106. c LECTURE DES MOT CLES PENAlite et STABilistation
  107. c-----------------------------------------------------
  108. RLATIN=0.d0
  109. XISTAB=0.d0
  110. INOMU=0
  111. 700 ICLE=0
  112. CALL LIRMOT(MOCLE,NMOCLE,ICLE,0)
  113. IF (ICLE.EQ.1) GOTO 701
  114. IF (ICLE.EQ.2) GOTO 702
  115. IF (ICLE.EQ.3) THEN
  116. INOMU=1
  117. GOTO 700
  118. ENDIF
  119. GOTO 799
  120.  
  121. c LECTURE DE LA DIRECTION DE RECHERCHE (FLOTTANT)
  122. 701 CALL LIRREE(RLATIN,1,IRETOU)
  123. IF (IERR.NE.0) RETURN
  124. if(iimpi.ge.2) write(ioimp,*) 'On a lu la raideur :',RLATIN
  125. GOTO 700
  126.  
  127. c LECTURE DE LA STABILISATION (FLOTTANT)
  128. 702 CALL LIRREE(XISTAB,1,IRETOU)
  129. IF (IERR.NE.0) RETURN
  130. if(iimpi.ge.2) write(ioimp,*) 'On a lu la STABILISATION :',XISTAB
  131. GOTO 700
  132.  
  133. 799 CONTINUE
  134.  
  135.  
  136. c-----------------------------------------------------
  137. c RECUPERATION DU MAILLAGE MASSIF
  138. c-----------------------------------------------------
  139. CALL LIROBJ ('MMODEL',IPMODL,1,IRETOU)
  140. IF(IERR.NE.0) RETURN
  141. MMODE1=IPMODL
  142. segact,MMODE1
  143. c Récupération du nombre de zones du modèle
  144. N1 = MMODE1.KMODEL(/1)
  145. if(N1.gt.1) write(ioimp,*) 'attention 1 seule zone a ce jour!'
  146. IMODE1 = MMODE1.KMODEL(1)
  147. segdes,MMODE1
  148. C Récupération du maillage et du numéro d'élément du modèle
  149. segact,IMODE1
  150. nele1 = IMODE1.NEFMOD
  151. IPT1 = IMODE1.IMAMOD
  152. SEGACT IPT1
  153. C Récupération du numéro d'élément du maillage, du nombre de noeuds et d'éléments
  154. iele1 = IPT1.itypel
  155. nbnn1 = IPT1.num(/1)
  156. nbel1 = IPT1.num(/2)
  157. c récupération des caractéristique EF IPT1
  158. cbp : on essaie d'abord avec INFELE
  159. mele1 = IMODE1.IMAMOD
  160. MINTE1 = IMODE1.INFELE(11)
  161. cbp : si echec, on retente avec elquoi
  162. if(mele1.eq.0.or.MINTE1.eq.0) then
  163. call elquoi(nele1,0,3,IPTR1,IMODE1)
  164. INFO = IPTR1
  165. segact,INFO
  166. mele1 = INFELL(1)
  167. MINTE1 = INFELL(11)
  168. segdes,INFO
  169. endif
  170. segdes,IMODE1
  171. cbp : si toujours pb, erreur...
  172. if(mele1.eq.0.or.MINTE1.eq.0) then
  173. write(ioimp,*) 'donnees relatives a l element fini inconnues'
  174. call erreur(591)
  175. endif
  176.  
  177.  
  178. c-----------------------------------------------------
  179. c RECUPERATION DU MAILLAGE INTERFACE
  180. c-----------------------------------------------------
  181. IPMAI2 = 0
  182. CALL LIROBJ ('MMODEL',IPMODL,0,IRETOU)
  183. IF(IERR.NE.0) RETURN
  184.  
  185. C -DANS LE CAS OU ON A UN MODELE EN ENTREE
  186. IF (IRETOU.EQ.1)THEN
  187. MMODE2=IPMODL
  188. segact,MMODE2
  189. N2 = MMODE2.KMODEL(/1)
  190. if(N2.gt.1) write(ioimp,*) 'attention 1 seule zone a ce jour!'
  191. IMODE2 = MMODE2.KMODEL(1)
  192. segdes,MMODE2
  193. segact,IMODE2
  194. nele2 = IMODE2.NEFMOD
  195. IPT2 = IMODE2.IMAMOD
  196. SEGACT IPT2
  197. c pour l'instant on dit que nele = iele (marche pour iele entre 2 et 26, voir bdata.eso)
  198. iele2 = IPT2.itypel
  199. nbnn2 = IPT2.num(/1)
  200. nbel2 = IPT2.num(/2)
  201.  
  202. c recuperation des caracteritiques de l'element
  203. cbp : on essaie d'abord avec INFELE
  204. MINTE2 = IMODE2.INFELE(11)
  205. ngau2 = IMODE2.INFELE(6)
  206. cbp : si echec, on retente avec elquoi
  207. if(MINTE2.eq.0.or.ngau2.eq.0) then
  208. c write(ioimp,*) 'appel elquoi',iele2,nele2
  209. call elquoi(nele2,0,3,IPTR2,IMODE2)
  210. INFO = IPTR2
  211. segact,INFO
  212. MINTE2 = INFELL(11)
  213. ngau2 = INFELL(6)
  214. segdes,INFO
  215. endif
  216. segdes,IMODE2
  217. cbp : si toujours pb, erreur...
  218. if(mele1.eq.0.or.MINTE1.eq.0) then
  219. write(ioimp,*) 'donnees relatives a l element fini INTERFACE'
  220. & ,'inconnues'
  221. call erreur(591)
  222. endif
  223.  
  224. C -DANS LE CAS OU ON A UN MAILLAGE EN ENTREE
  225. ELSE
  226. CALL LIROBJ('MAILLAGE',IPMAI2,1,IRETOU)
  227. IF(IERR.NE.0) RETURN
  228. IPT2 = IPMAI2
  229. SEGACT IPT2
  230. c pour l'instant on dit que nele = iele (marche pour iele entre 2 et 26, voir bdata.eso)
  231. iele2 = IPT2.itypel
  232. nele2 = iele2
  233. if (nele2.lt.2.or.nele2.gt.26) then
  234. write(ioimp,*)'element geometrique different de l element fini'
  235. call erreur(16)
  236. endif
  237. nbnn2 = IPT2.num(/1)
  238. nbel2 = IPT2.num(/2)
  239. ngau2 = -3
  240. c SEG2
  241. if (nele2.EQ.2) ngau2 = 2
  242. c SEG3
  243. if (nele2.EQ.3) ngau2 = 3
  244. c TRI3
  245. if (nele2.EQ.4) ngau2 = 1
  246. c TRI6
  247. if (nele2.EQ.6) ngau2 = 4
  248. c QUA4
  249. if (nele2.EQ.8) ngau2 = 4
  250. c QUA8
  251. if (nele2.EQ.10) ngau2 = 9
  252. if(ngau2.eq.-3) then
  253. write(ioimp,*)'nombre de points d integration inconnu pour',
  254. & 'ce type d element geometrique'
  255. call erreur(16)
  256. endif
  257. ENDIF
  258.  
  259. call RESHPT(ngau2,nbnn2,iele2,nele2,0,IPTR2,IRET)
  260. MINTE2 = IPTR2
  261. segact,MINTE1,MINTE2
  262.  
  263.  
  264. c-----------------------------------------------------
  265. c RECHERCHE DU MCHAML ISSU MCHEX1 D ENRICHISSEMENT
  266. c-----------------------------------------------------
  267. MCHAM1=0
  268. NBENR2=0
  269. segact,IMODE1
  270. NOBMOD = IMODE1.IVAMOD(/1)
  271. IF (NOBMOD.NE.0) THEN
  272. DO 1002 iobmo1=1,NOBMOD
  273. if((IMODE1.TYMODE(iobmo1)).ne.'MCHAML') goto 1002
  274. MCHEX1 = IMODE1.IVAMOD(iobmo1)
  275. segact,MCHEX1
  276. if((MCHEX1.TITCHE).ne.'ENRICHIS') goto 1003
  277. MCHAM1 = MCHEX1.ICHAML(1)
  278. segact,MCHAM1
  279. NBENR2 = MCHAM1.IELVAL(/1)
  280. do ienr2=1,NBENR2
  281. MELVA1=MCHAM1.IELVAL(IENR2)
  282. if(MELVA1.ne.0) segact,MELVA1
  283. enddo
  284. 1003 continue
  285. segdes,MCHEX1
  286. 1002 CONTINUE
  287. ENDIF
  288.  
  289. c-------------------------------------
  290. c INITIALISATION DES OBJETS DE TRAVAIL : MTRAV et ITYMAT
  291. c-------------------------------------
  292.  
  293. segini,MTRAV
  294.  
  295. *bp : NTYMAT = (U ou H ou HB1 ou HB1B2) * dim * nbre de configs possibles
  296. * NTYMAT = nbre de types de matrices
  297. * nbre de configs varie entre :
  298. * -on peut avoir tous les pt de Gauss dans 1 seul elements de structure
  299. * -jusqu'à chaque pt de Gauss dans un element distinct
  300. ity=0
  301. c boucle sur les enrichissements possibles U, H , HB1, HB1B2
  302. if(iimpi.ge.3) write(ioimp,*) 'ity , ino1 , ienr1'
  303. do ienr1=1,4
  304. c boucle sur les nombre de noeuds "structure"
  305. do ino1=(nbnn1),(ngau2*nbnn1)
  306. ity=ity+1
  307. ITYMAT(1,ity)=ino1
  308. ITYMAT(2,ity)=ienr1
  309. ITAMYT(ino1,ienr1)=ity
  310. if(iimpi.ge.3) write(ioimp,*) ity,ino1,ienr1
  311. enddo
  312. enddo
  313. NTYMAT = 4*(nbnn1*(ngau2-1)+1)
  314. if(iimpi.ge.3) write(ioimp,*) 'ity*idim=',ity,NTYMAT
  315.  
  316.  
  317. c--------------------------------------------------------------------
  318. c INITIALISATION DU SEGMENT MRIGID
  319. c--------------------------------------------------------------------
  320. NRIGEL = NTYMAT*idim
  321. segini,MRIGID
  322. IFORIG = IFOUR
  323. MTYMAT ='RIGIDITE'
  324. c -on prepare le meleme
  325. NBSOUS = 0
  326. NBREF = 0
  327.  
  328. ityty=0
  329. c on initialise la taille matrice en fonction du type de matrice
  330. do ity=1,NTYMAT
  331. do iidim=1,idim
  332.  
  333. ityty=ityty+1
  334. COERIG(ityty) = 1.D0
  335.  
  336. * dim de la matrice RE elementaire = nbnoeud*TX + nbnoeud*AX + nbnoeud*AX
  337. nbno1 = ITYMAT(1,ity)
  338. nenr1 = ITYMAT(2,ity)
  339. if(nenr1.le.2) then
  340. NLIGRP = nbnn2 + nbnn2 + nbno1
  341. else
  342. NLIGRP = nbnn2 + nbnn2 + ((nenr1-1)*nbno1)
  343. endif
  344. NLIGRD = NLIGRP
  345.  
  346. c -creation du MELEME
  347. NBNN = nbnn2 + nbno1
  348. NBELEM=0
  349. SEGINI,MELEME
  350. ITYPEL=28
  351. IRIGEL(1,ityty) = MELEME
  352.  
  353. c -remplissage du DESCR
  354. SEGINI,DESCR
  355. IRIGEL(3,ityty) = DESCR
  356. iddl=0
  357. c remplissage des ddl LX de la fissure (ici appelé FX)
  358. c on iverse donc inconnue et duale []*FX=UX
  359. do ino2=1,nbnn2
  360. iddl=iddl+1
  361. if (nenr1.eq.1) then
  362. LISINC(iddl)=DUAOB2(iidim)
  363. LISDUA(iddl)=DDLOB2(iidim)
  364. else
  365. LISINC(iddl)=DUAFA2(iidim)
  366. LISDUA(iddl)=DDLFA2(iidim)
  367. endif
  368. NOELEP(iddl)=ino2
  369. NOELED(iddl)=ino2
  370. enddo
  371. c remplissage des ddl UX de la fissure (ici appelé WX) []*WX=TX
  372. do ino2=1,nbnn2
  373. iddl=iddl+1
  374. if (nenr1.eq.1) then
  375. LISINC(iddl)=DDLOB2(iidim)
  376. LISDUA(iddl)=DUAOB2(iidim)
  377. else
  378. LISINC(iddl)=DDLFA2(iidim)
  379. LISDUA(iddl)=DUAFA2(iidim)
  380. endif
  381. NOELEP(iddl)=ino2
  382. NOELED(iddl)=ino2
  383. enddo
  384. c remplissage des ddl de la structure
  385. if (nenr1.eq.1) then
  386. do ino1=1,nbno1
  387. iddl=iddl+1
  388. LISINC(iddl)=DDLOBL(iidim)
  389. LISDUA(iddl)=DUAOBL(iidim)
  390. NOELEP(iddl)=nbnn2+ino1
  391. NOELED(iddl)=nbnn2+ino1
  392. enddo
  393. else
  394. do ini1=1,(nenr1-1)
  395. do ino1=1,nbno1
  396. iddl=iddl+1
  397. LISINC(iddl)=DDLFAC(iidim+(3*(ini1-1)))
  398. LISDUA(iddl)=DUAFAC(iidim+(3*(ini1-1)))
  399. NOELEP(iddl)=nbnn2+ino1
  400. NOELED(iddl)=nbnn2+ino1
  401. enddo
  402. enddo
  403. endif
  404. if(iimpi.ge.3) write(ioimp,*) ityty,(LISINC(iou),iou=1,NLIGRP)
  405. SEGDES,DESCR
  406.  
  407. c -initialisation du XMATRI
  408. NELRIG=0
  409. SEGINI,XMATRI
  410. IRIGEL(4,ityty) = XMATRI
  411. IRIGEL(5,ityty) = NIFOUR
  412. IRIGEL(6,ityty) = 0
  413. IRIGEL(7,ityty) = 0
  414. IRIGEL(8,ityty) = 0
  415.  
  416. enddo
  417. enddo
  418.  
  419.  
  420. c----------------------------------------------------------------------
  421. c 1. RECHERCHE DES ELEMENTS DE STRUCTURE CONTENANT DES POINTS DE GAUSS
  422. c DES ELEMENTS DE LA FISSURE
  423. c 2. REMPLISSAGE DU MELEME + MRIGID en FONCTION
  424. c----------------------------------------------------------------------
  425.  
  426. C
  427. c==== Boucle sur les elements de fissure ==============================
  428. DO 1100 iem2=1,nbel2
  429.  
  430. call doxe(xcoor,idim,nbnn2,ipt2.num,iem2,xe2)
  431. nbenrj = 0
  432. nbnno2 = 0
  433. call ZERO(RINT,nbnn2,nbnn2)
  434. call ZERO(RINT0,nbnn2,ngau2*nbnn1)
  435. call ZERO(RINT1,nbnn2,ngau2*nbnn1)
  436. call ZERO(RINT2,nbnn2,ngau2*nbnn1)
  437. call ZERO(RINT3,nbnn2,ngau2*nbnn1)
  438.  
  439. c======= Boucle sur les pt de G de fissure ============================
  440. DO 1132 igau2=1,ngau2
  441. izok = 0
  442. c récupération des coordonnees du point de gauss dans le repère global
  443. XPO(1) = 0.D0
  444. XPO(2) = 0.D0
  445. XPO(3) = 0.D0
  446. DO ino=1,nbnn2
  447. XPO(1) = XPO(1) + (MINTE2.SHPTOT(1,ino,IGAU2)*xe2(1,ino))
  448. XPO(2) = XPO(2) + (MINTE2.SHPTOT(1,ino,IGAU2)*xe2(2,ino))
  449. IF(IDIM.EQ.3) then
  450. XPO(3) = XPO(3) + (MINTE2.SHPTOT(1,ino,IGAU2)*xe2(3,ino))
  451. ENDIF
  452. ENDDO
  453.  
  454. c if(iimpi.ge.5) write(ioimp,*) 'Fissure: element',iem2,
  455. c & ' point de Gauss',igau2,' x=',XPO(1),XPO(2)
  456.  
  457. c---------- Boucle sur les elements de structure ----------------------
  458. DO 1131 iem1=1,nbel1
  459.  
  460. c si pas d'enrichissement, on travaille sur tous les elements
  461. if(MCHAM1.eq.0) goto 1133
  462. c on saute les elements non enrichi car a priori ne contiennent pas la fissure
  463. do ienr2=1,NBENR2
  464. MELVA1=MCHAM1.IELVAL(IENR2)
  465. if(MELVA1.ne.0) then
  466. do inode1=1,nbnn1
  467. if(MELVA1.IELCHE(inode1,iem1).ne.0) goto 1133
  468. enddo
  469. endif
  470. enddo
  471. goto 1131
  472. 1133 continue
  473.  
  474. c recuperation des coordonnées des noeuds de IPT1 : xel1 (dans le repère x,y,z)
  475. call doxe(xcoor,idim,nbnn1,ipt1.num,iem1,xel1)
  476.  
  477. c calcul des fonctions de formes de IPT1 au pt de Gauss de IPT2
  478. call QSIJS(xel1,iele1,nbnn1,idim,XPO,SHPP1,qsi,IRET)
  479. c if(iimpi.ge.5) write(ioimp,*) 'Structure: element',iem1,
  480. c & ' noeuds:',(ipt1.num(iou,iem1),iou=1,nbnn1)
  481. c if(iimpi.ge.5) write(ioimp,*)
  482. c & ' N_i=',(SHPP1(1,iou),iou=1,nbnn1)
  483.  
  484.  
  485. c test pour savoir si PG est dans EF de IPT1
  486. DO 1130 ino1=1,NBNN1
  487. if (SHPP1(1,ino1).LT.-1.01D-7) then
  488. go to 1131
  489. endif
  490. 1130 continue
  491. izok = 1
  492. c ON a trouvé l'élément de structure !
  493. c write(ioimp,*) 'ON a trouve l element de structure !'
  494. c DETECTION DU TYPE D'ENRICHISSEMENT MAX DE CET ELEMENT = nbenrj
  495. DO 3001 IENR2=1,NBENR2
  496. MELVA1=MCHAM1.IELVAL(IENR2)
  497. IF(MELVA1.eq.0) GOTO 3001
  498. DO 3002 ino1=1,nbnn1
  499. MLREEL = MELVA1.IELCHE(ino1,iem1)
  500. c Test pour savoir si le noeud est enrichi
  501. IF(MLREEL.eq.0) GOTO 3002
  502. nbenrj=max(nbenrj,IENR2)
  503. 3002 continue
  504. 3001 continue
  505.  
  506. if(iimpi.ge.3) write(ioimp,*) 'EF fissure ',iem2,
  507. & ' ptdeG ',igau2,' -> EF MASSIF ',iem1,' nbenrj=',nbenrj
  508.  
  509. c Préparation du Remplissage du MELEME
  510. do 3003 ino1=1,nbnn1
  511. iino1=IPT1.NUM(ino1,iem1)
  512. if(nbnno2.gt.0) then
  513. do iino2=1,nbnno2
  514. if(PRELEM(iino2).eq.iino1) goto 3003
  515. enddo
  516. endif
  517. c nouveau noeud : on l'ajoute
  518. nbnno2=nbnno2+1
  519. PRELEM(nbnno2)=iino1
  520. 3003 continue
  521. if(iimpi.ge.3) then
  522. write(ioimp,*) ' PRELEM=',( PRELEM(iou),iou=1,nbnno2)
  523. write(ioimp,*) 'EF fissure ',iem2,' ptdeG ',igau2
  524. endif
  525.  
  526. C CALCUL DE REL ET REL1
  527.  
  528. c matrices des fonctions de forme de l'element associé à IPT2
  529.  
  530. do ino2=1,NBNN2
  531. xe3(1,ino2) = sqrt( xe2(1,ino2)*xe2(1,ino2)
  532. 1 + xe2(2,ino2)*xe2(2,ino2) )
  533. NGENE2(1,ino2)=MINTE2.SHPTOT(1,ino2,IGAU2)
  534. do ipp2=1,6
  535. SHPP2(ipp2,ino2)=MINTE2.SHPTOT(ipp2,ino2,IGAU2)
  536. enddo
  537. if(iimpi.ge.4) then
  538. write(ioimp,*) 'SHPP2',(SHPP2(iou,ino2),iou=1,6)
  539. write(ioimp,*) 'xe2',(xe2(iou,ino2),iou=1,2)
  540. endif
  541. enddo
  542. c calcul du jacobien
  543. if(idim.eq.2) then
  544. c calcul du jacobien en 2D pour un seg2
  545. dXdQsi=0.D0
  546. dYdQsi=0.D0
  547. DO i=1,NBNN2
  548. dXdQsi=dXdQsi+SHPP2(2,i)*XE2(1,i)
  549. dYdQsi=dYdQsi+SHPP2(2,i)*XE2(2,i)
  550. ENDDO
  551. c Dans le cas des seg2, poids des pts de gauss = 1
  552. DJAC2=SQRT(dXdQsi*dXdQsi+dYdQsi*dYdQsi)
  553. else
  554. dXdQsi=0.D0
  555. dYdQsi=0.D0
  556. dZdQsi=0.D0
  557. dXdEta=0.D0
  558. dYdEta=0.D0
  559. dZdEta=0.D0
  560. DO i=1,NBNN2
  561. dXdQsi=dXdQsi+SHPP2(2,i)*XE2(1,i)
  562. dXdEta=dXdEta+SHPP2(3,i)*XE2(1,i)
  563. dYdQsi=dYdQsi+SHPP2(2,i)*XE2(2,i)
  564. dYdEta=dYdEta+SHPP2(3,i)*XE2(2,i)
  565. dZdQsi=dZdQsi+SHPP2(2,i)*XE2(3,i)
  566. dZdEta=dZdEta+SHPP2(3,i)*XE2(3,i)
  567. ENDDO
  568. z = (dXdQsi*dYdEta-dXdEta*dYdQsi)
  569. x = (dYdQsi*dZdEta-dYdEta*dZdQsi)
  570. y = (dZdQsi*dXdEta-dZdEta*dXdQsi)
  571. c Dans le cas des tri3, poids du pts de gauss = 0.5
  572. DJAC2 = sqrt(x*x+y*y+z*z)*0.5
  573. endif
  574.  
  575. c matrices des fonctions de forme de l'element associé à IPT1
  576. do ino1=1,NBNN1
  577. NGENE1(1,ino1)=SHPP1(1,ino1)
  578. enddo
  579. if(iimpi.ge.4) then
  580. write(ioimp,*) 'DJAC2=',DJAC2
  581. write(ioimp,*) 'NGENE1=',(NGENE1(1,iou),iou=1,NBNN1)
  582. write(ioimp,*) 'NGENE2=',(NGENE2(1,iou),iou=1,NBNN2)
  583. endif
  584.  
  585. c calcul des matrices elementaires = produits des fonctions de forme
  586. DO II=1,nbnn2
  587. XGENE2 = NGENE2(1,II)*DJAC2
  588. c calcul de K_ww, Kwt et Ktt
  589. DO JJ=1,nbnn2
  590. jjj = JJ
  591. C CC = Jacobien * transposee(N)(II,i) * N(i,JJ) (somme sur i)
  592. CC = XGENE2*NGENE2(1,jjj)
  593. C REL est une matrice symetrique
  594. REL(II,JJ) = CC
  595. REL(JJ,II) = REL(II,JJ)
  596. ENDDO
  597. c calcul de K_tu
  598. DO JJ=1,nbnn1
  599. C CC = Jacobien * transposee(N)(II,i) * N(i,JJ) (somme sur i)
  600. CC=XGENE2*NGENE1(1,JJ)
  601. c write(ioimp,*) ' CC = ', CC
  602. c REL1(II,JJ)=-1.D0*CC
  603. REL1(II,JJ)=CC
  604. c REL(JJdec,II)=REL(II,JJdec)
  605. ENDDO
  606. if(iimpi.ge.4) write(ioimp,*)
  607. 1 'REL1(',II,',:)=',(REL1(II,iou),iou=1,NBNN1)
  608. ENDDO
  609.  
  610. C CUMUL DANS RINT ET RINT1
  611. C CUMUL DANS RINT(II,JJ) pour Kww, Kwt, Ktw et Ktt
  612. c (4 fois la meme matrices = N^{w T} * N^{w} = N^{w T} * N^{t} = ...)
  613. DO II=1,nbnn2
  614. DO JJ=1,nbnn2
  615. RINT(II,JJ) = RINT(II,JJ) + REL(II,JJ)
  616. ENDDO
  617. ENDDO
  618.  
  619. C CUMUL DANS RINT1(II,JJ) pour Kut ( = N^{t T} * N^{u})
  620. C INTEGRATION POSSIBLEMENT INCOMPATIBLE POUR KUT !
  621. c on peut prendre le 1er enrichissement (ou le standard) car on a
  622. c toujours les memes noeuds au debut (sans tenir compte des B1X)
  623.  
  624. c Stockage de l'element iem1 pour les differents PG
  625. elemPG(igau2) = iem1
  626.  
  627. DO 4205 JJ=1,nbnn1
  628. iino1=IPT1.NUM(JJ,iem1)
  629.  
  630. do JJ1 = 1,nbnno2
  631. if(PRELEM(JJ1).eq.iino1) goto 4207
  632. enddo
  633. write(IOIMP,*) 'On n a pas trouvé le noeud ',iino1
  634. call ERREUR(5)
  635. 4207 continue
  636. itoto = 0
  637.  
  638. c UX
  639. DO II=1,nbnn2
  640. RINT0(II,JJ1) = RINT0(II,JJ1) + REL1(II,JJ)
  641. ENDDO
  642. c AX
  643. DO II=1,nbnn2
  644. MELVA1 = MCHAM1.IELVAL(1)
  645. MLREEL = MELVA1.IELCHE(JJ,iem1)
  646. if(MLREEL.ne.0) then
  647. RINT1(II,JJ1) = RINT1(II,JJ1) + REL1(II,JJ)
  648. c write(ioimp,*) ' RINT1(',II,JJ1,') = ', RINT1(II,JJ1)
  649. endif
  650. ENDDO
  651.  
  652.  
  653. c B1X
  654. IF(nbenrj.ge.2) THEN
  655. MELVA1 = MCHAM1.IELVAL(2)
  656. MLREEL = MELVA1.IELCHE(JJ,iem1)
  657. if(MLREEL.ne.0) then
  658. PSIX = 0.d0
  659. SEGACT,MLREEL
  660. do iii0 = 1,nbnn1
  661. PSIX = PSIX + (SHPP1(1,iii0) * PROG(iii0))
  662. enddo
  663. SEGDES,MLREEL
  664. RX05= SQRT(ABS(PSIX))
  665. c write(ioimp,*) 'RX05 = ', RX05
  666. do II=1,nbnn2
  667. RINT2(II,JJ1) = RINT2(II,JJ1) + RX05*REL1(II,JJ)
  668. enddo
  669. endif
  670. ENDIF
  671.  
  672.  
  673. c B2X
  674. IF(nbenrj.ge.3) THEN
  675. MELVA1 = MCHAM1.IELVAL(3)
  676. MLREEL = MELVA1.IELCHE(JJ,iem1)
  677. if(MLREEL.ne.0) then
  678. PSIX = 0.d0
  679. SEGACT,MLREEL
  680. do iii0 = 1,nbnn1
  681. PSIX = PSIX + (SHPP1(1,iii0) * PROG(iii0))
  682. enddo
  683. SEGDES,MLREEL
  684. RX05= SQRT(ABS(PSIX))
  685. do II=1,nbnn2
  686. RINT3(II,JJ1) = RINT3(II,JJ1) + RX05*REL1(II,JJ)
  687. enddo
  688. endif
  689. ENDIF
  690.  
  691. 4205 CONTINUE
  692.  
  693. if(iimpi.ge.4) then
  694. DO II=1,nbnn2
  695. write(ioimp,*) 'RINT1(',II,',:)=',(RINT1(II,iou),iou=1,nbnno2)
  696. ENDDO
  697. endif
  698.  
  699.  
  700. c on a trouvé le iem1 élément et on a fait le travail :
  701. c on passe au point de gauss suivant
  702. goto 1132
  703.  
  704. 1131 CONTINUE
  705. c---------- fin de la Boucle sur les elements de structure ----------------------
  706. if(iimpi.ge.2) write(ioimp,*)
  707. 1 'Fin de la boucle sur les elements de structure'
  708. if(izok.eq.0) write(ioimp,*) 'Attention le pt de Gauss '
  709. & ,igau2,'de l element ',iem2,' est hors support'
  710.  
  711.  
  712. 1132 CONTINUE
  713. c======= fin de la Boucle sur les pt de G de fissure ============================
  714.  
  715. c --- Remplissage du MELEME et DU XMATRI ---
  716. DO ityp1=1,2
  717.  
  718. c on commence par le MELEME des ddls Obligatoire (U)
  719. izone = 1
  720. c on poursuit avec le MELEME des ddls Sauts (H ou HB1 ou HB1B2)
  721. if(ityp1.eq.2) izone=nbenrj+1
  722. ity = ITAMYT(nbnno2,izone)
  723.  
  724. do iidim=1,idim
  725. ityty = ((ity-1) * idim) + iidim
  726. if(iimpi.ge.4) write(ioimp,*) ity,nbnno2,izone,' -> ',ityty
  727.  
  728. c --- Remplissage du MELEME ---
  729. MELEME = IRIGEL(1,ityty)
  730. NBNN = NUM(/1)
  731. NBELEM = NUM(/2) + 1
  732. if(iimpi.ge.4) then
  733. write(ioimp,*) 'on remplit le type de matrice ',ityty,
  734. 1 ' NBNN,NBELEM=',NBNN,NBELEM
  735. endif
  736. segadj,MELEME
  737. c element de fissure
  738. do ino2=1,nbnn2
  739. NUM(ino2,NBELEM) = IPT2.NUM(ino2,iem2)
  740. c NUM(ino2+nbnn2,NBELEM) = IPT2.NUM(ino2,iem2)
  741. enddo
  742. c element de structure
  743. do iino2=1,nbnno2
  744. c NUM((2*nbnn2)+iino2,NBELEM) = PRELEM(iino2)
  745. NUM(nbnn2+iino2,NBELEM) = PRELEM(iino2)
  746. enddo
  747.  
  748. c --- Remplissage du XMATRI.RE ---
  749. XMATRI = IRIGEL(4,ityty)
  750. NLIGRD = RE(/1)
  751. NLIGRP = RE(/2)
  752. NELRIG = RE(/3) + 1
  753. segadj,XMATRI
  754.  
  755. DO II=1,nbnn2
  756.  
  757. c -Ktt
  758. c on choisit de ne garder que la partie diagonale
  759. JJ=II
  760. XMATRI.RE(II,JJ,NBELEM) = RINT(II,JJ)*XISTAB
  761.  
  762. c -Kwt
  763. DO JJ=1,nbnn2
  764. XMATRI.RE(II,JJ+nbnn2,NBELEM) = RINT(II,JJ)
  765. XMATRI.RE(JJ+nbnn2,II,NBELEM) = RINT(II,JJ)
  766. c write(ioimp,*) 'XMATRI.RE(',II,JJ+nbnn2,NBELEM,') =', RINT(II,JJ)
  767. ENDDO
  768.  
  769. C -Kww
  770. c on choisit de ne garder que la partie diagonale
  771. JJ=II
  772. XMATRI.RE(nbnn2+II,nbnn2+JJ,NBELEM) = RINT(II,JJ)*RLATIN
  773.  
  774. C -Kut (on a deja fait le travail de positionnement avant)
  775. if(ityp1.eq.1) then
  776. c UX
  777. JJdec = 2*nbnn2
  778. DO JJ=1,nbnno2
  779. XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT0(II,JJ)
  780. XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT0(II,JJ)
  781. c write (6,* ) 'XMATRI.RE(',II,JJdec+JJ,NBELEM,') =', XMA
  782. c 1TRI.RE(II,JJdec+JJ,NBELEM)
  783. ENDDO
  784. else
  785. c AX
  786. JJdec = 2*nbnn2
  787. DO JJ=1,nbnno2
  788. XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT1(II,JJ)
  789. XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT1(II,JJ)
  790. c write (6,* ) 'XMATRI.RE(',II,JJdec+JJ,NBELEM,') =', XMA
  791. c 1TRI.RE(II,JJdec+JJ,NBELEM)
  792. ENDDO
  793.  
  794. c B1X
  795. IF(nbenrj.ge.2) THEN
  796. JJdec=JJdec+nbnno2
  797. DO JJ=1,nbnno2
  798. XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT2(II,JJ)
  799. XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT2(II,JJ)
  800.  
  801. c Cas ou les inconnues aux noeuds ne sont pas présentes
  802. c mis à 0 de ces inconnues en mettant un 1 sur la diagonale
  803.  
  804. c Début de la boucle pour tester les noeuds de structure d'ancien point de gauss
  805. cbp DO 5002 iigau2 = 1,igau2
  806. DO 5002 iigau2 = 1,ngau2
  807. iiem1 = elemPG(iigau2)
  808. DO JJ1=1,nbnn1
  809. iino1=IPT1.NUM(JJ1,iiem1)
  810. if(PRELEM(JJ).eq.iino1) then
  811. MELVA1 = MCHAM1.IELVAL(2)
  812. MLREEL = MELVA1.IELCHE(JJ1,iiem1)
  813. if(MLREEL.ne.0) goto 5002
  814. XMATRI.RE(JJdec+JJ,JJdec+JJ,NBELEM) = 1
  815. c write(ioimp,*) 'XMATRI.RE(',JJd ec+JJ,JJdec+JJ,NBEL
  816. c 1EM,') =', XMATRI.RE(JJdec+JJ,JJdec+JJ,NBELEM)
  817. endif
  818. ENDDO
  819. 5002 CONTINUE
  820. ENDDO
  821.  
  822.  
  823.  
  824. ENDIF
  825. c B2X
  826. IF(nbenrj.ge.3) THEN
  827. JJdec=JJdec+nbnno2
  828. DO JJ=1,nbnno2
  829. XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT3(II,JJ)
  830. XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT3(II,JJ)
  831. ENDDO
  832. ENDIF
  833. endif
  834.  
  835. ENDDO
  836.  
  837. enddo
  838.  
  839.  
  840. ENDDO
  841.  
  842.  
  843.  
  844. 1100 CONTINUE
  845. c==== fin de la Boucle sur les elements de fissure ==============================
  846. c----------------------------------------------------------------------
  847.  
  848.  
  849. c------------------------------------
  850. c AJUSTEMENT AVANT DE QUITTER
  851. c------------------------------------
  852.  
  853. if(iimpi.ge.2) write(ioimp,*) 'AJUSTEMENT AVANT DE QUITTER'
  854.  
  855.  
  856. C BOUCLE SUR LES SOUS RIGIDITÉS
  857. ityok = 0
  858. DO 2000 ityty=1,(idim*NTYMAT)
  859.  
  860. MELEME = IRIGEL(1,ityty)
  861. DESCR = IRIGEL(3,ityty)
  862. XMATRI = IRIGEL(4,ityty)
  863. NBELEM = NUM(/2)
  864. if(NBELEM.ne.0) then
  865. ityok = ityok + 1
  866. IRIGEL(1,ityok) = MELEME
  867. IRIGEL(3,ityok) = DESCR
  868. IRIGEL(4,ityok) = XMATRI
  869. segdes,MELEME,XMATRI
  870. else
  871. segsup,MELEME,DESCR,XMATRI
  872. endif
  873. 2000 continue
  874. NRIGEL = ityok
  875. segadj,MRIGID
  876.  
  877. c-------------------------------
  878. c MENAGE AVANT DE QUITTER
  879. c-------------------------------
  880.  
  881. segsup,MTRAV
  882. segdes,MINTE1,MINTE2
  883. segdes,IPT1,IPT2
  884. segdes,MCHAM1
  885.  
  886.  
  887. c-------------------------------
  888. c ON REMPLACE LES FX PAR DES LX
  889. c-------------------------------
  890. IF(INOMU.eq.0) THEN
  891. JGN=4
  892. JGM=2*idim
  893. segini,MLMOT1
  894. do iidim=1,idim
  895. MLMOT1.MOTS(iidim) = DUAOBL(iidim)
  896. MLMOT1.MOTS(iidim+idim) = DUAFAC(iidim)
  897. enddo
  898. ILMOT1 = MLMOT1
  899. IPRIG1 = MRIGID
  900. if(iimpi.ge.2) write(ioimp,*)' APPEL A FX2LX '
  901. CALL FX2LX(IPRIG1,ILMOT1,IPRIG2)
  902. MRIGID = IPRIG2
  903. SEGSUP,MLMOT1
  904. ENDIF
  905.  
  906. c write(ioimp,*) 'desactivation du MRIGID',MRIGID
  907. SEGDES,MRIGID
  908.  
  909. if(iimpi.ge.3) write(ioimp,*) 'ecriture du MRIGID',MRIGID
  910. CALL ECROBJ('RIGIDITE',MRIGID)
  911.  
  912. RETURN
  913.  
  914. END
  915.  
  916.  
  917.  
  918.  
  919.  
  920.  
  921.  
  922.  
  923.  
  924.  
  925.  
  926.  

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