Télécharger accro2.eso

Retour à la liste

Numérotation des lignes :

accro2
  1. C ACCRO2 SOURCE CB215821 19/08/20 21:15:00 10287
  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.  
  25. -INC PPARAM
  26. -INC CCOPTIO
  27. -INC SMCOORD
  28. -INC SMELEME
  29. -INC SMRIGID
  30. -INC SMLMOTS
  31. -INC CCHAMP
  32. -INC CCGEOME
  33. -INC SMMODEL
  34. -INC SMINTE
  35. -INC SMCHAML
  36. -INC SMLREEL
  37.  
  38.  
  39. C Segment contenant les informations sur un element
  40. SEGMENT INFO
  41. INTEGER INFELL(JG)
  42. ENDSEGMENT
  43.  
  44. POINTEUR INFO1.INFO,INFO2.INFO
  45. POINTEUR MCHEX1.MCHELM
  46. external shape
  47. DATA EPSI/1.D-9/
  48. DIMENSION ICOR(6),IMEL(6),imtt(10)
  49. DIMENSION QSI(3),XPO(3)
  50.  
  51. C INITIALISATION DES INCONNUES obligatoires et facultatives
  52. PARAMETER (NOBL=3,NFAC=9,NMOCLE=3)
  53. CHARACTER*4 DDLOBL(NOBL),DDLFAC(NFAC),MODDL,MODDL2
  54. CHARACTER*4 DUAOBL(NOBL),DUAFAC(NFAC),MOCLE(NMOCLE)
  55. CHARACTER*4 DDLOB2(NOBL),DDLFA2(NOBL)
  56. CHARACTER*4 DUAOB2(NOBL),DUAFA2(NOBL)
  57.  
  58. DATA DDLOBL/'UX ','UY ','UZ '/
  59. DATA DDLFAC/'AX ','AY ','AZ ',
  60. >'B1X ','B1Y ','B1Z ',
  61. >'B2X ','B2Y ','B2Z '/
  62. DATA DUAOBL/'FX ','FY ','FZ '/
  63. DATA DUAFAC/'FAX ','FAY ','FAZ ',
  64. >'FB1X','FB1Y','FB1Z',
  65. >'FB2X','FB2Y','FB2Z'/
  66.  
  67. c nom de composantes speciales
  68. DATA DDLOB2/'UX ','UY ','UZ '/
  69. DATA DDLFA2/'AX ','AY ','AZ '/
  70. DATA DUAOB2/'FX ','FY ','FZ '/
  71. DATA DUAFA2/'FAX ','FAY ','FAZ '/
  72.  
  73. DATA MOCLE /'PENA','STAB','NOMU'/
  74.  
  75. c-----------------------------------------------------
  76. c Segment de travail
  77. c-----------------------------------------------------
  78. SEGMENT MTRAV
  79. REAL*8 SHPP2(6,NBNN2)
  80. REAL*8 XE2(3,NBNN2),XEL2(3,NBNN2),NGENE2(1,NBNN2)
  81. REAL*8 XE3(3,NBNN2)
  82. REAL*8 SHPP1(6,NBNN1)
  83. REAL*8 XE1(3,NBNN1),XEL1(3,NBNN1)
  84. REAL*8 NGENE1(1,NBNN1)
  85. REAL*8 REL(nbnn2,nbnn2),REL1(nbnn2,nbnn1)
  86. REAL*8 RINT(nbnn2,nbnn2)
  87. REAL*8 RINT0(nbnn2,ngau2*nbnn1),RINT1(nbnn2,ngau2*nbnn1)
  88. REAL*8 RINT2(nbnn2,ngau2*nbnn1),RINT3(nbnn2,ngau2*nbnn1)
  89. INTEGER ITYMAT(2,4*nbnn1*ngau2)
  90. INTEGER ITAMYT(ngau2*NBNN1,4)
  91. INTEGER PRELEM(3*nbnn1*ngau2)
  92. INTEGER ElemPG(ngau2)
  93. ENDSEGMENT
  94.  
  95. c tableaux comptant le nbre d EF de chaque ddl
  96. PARAMETER(NDDLMAX=6)
  97. INTEGER NELDDL(NDDLMAX)
  98.  
  99. if(iimpi.ge.2) then
  100. write(ioimp,*) '-----------------------------------------------'
  101. write(ioimp,*) ' ENTREE dans ACCRO2'
  102. write(ioimp,*) '-----------------------------------------------'
  103. endif
  104.  
  105. SEGACT MCOORD
  106.  
  107. c-----------------------------------------------------
  108. c LECTURE DES MOT CLES PENAlite et STABilistation
  109. c-----------------------------------------------------
  110. RLATIN=0.d0
  111. XISTAB=0.d0
  112. INOMU=0
  113. 700 ICLE=0
  114. CALL LIRMOT(MOCLE,NMOCLE,ICLE,0)
  115. IF (ICLE.EQ.1) GOTO 701
  116. IF (ICLE.EQ.2) GOTO 702
  117. IF (ICLE.EQ.3) THEN
  118. INOMU=1
  119. GOTO 700
  120. ENDIF
  121. GOTO 799
  122.  
  123. c LECTURE DE LA DIRECTION DE RECHERCHE (FLOTTANT)
  124. 701 CALL LIRREE(RLATIN,1,IRETOU)
  125. IF (IERR.NE.0) RETURN
  126. if(iimpi.ge.2) write(ioimp,*) 'On a lu la raideur :',RLATIN
  127. GOTO 700
  128.  
  129. c LECTURE DE LA STABILISATION (FLOTTANT)
  130. 702 CALL LIRREE(XISTAB,1,IRETOU)
  131. IF (IERR.NE.0) RETURN
  132. if(iimpi.ge.2) write(ioimp,*) 'On a lu la STABILISATION :',XISTAB
  133. GOTO 700
  134.  
  135. 799 CONTINUE
  136.  
  137.  
  138. c-----------------------------------------------------
  139. c RECUPERATION DU MAILLAGE MASSIF
  140. c-----------------------------------------------------
  141. CALL LIROBJ('MMODEL ',IPMODL,1,IRETOU)
  142. CALL ACTOBJ('MMODEL ',IPMODL,1)
  143. IF(IERR.NE.0) RETURN
  144. MMODE1=IPMODL
  145. segact,MMODE1
  146. c Récupération du nombre de zones du modèle
  147. N1 = MMODE1.KMODEL(/1)
  148. if(N1.gt.1) write(ioimp,*) 'attention 1 seule zone a ce jour!'
  149. IMODE1 = MMODE1.KMODEL(1)
  150. C Récupération du maillage et du numéro d'élément du modèle
  151. segact,IMODE1
  152. nele1 = IMODE1.NEFMOD
  153. IPT1 = IMODE1.IMAMOD
  154. SEGACT IPT1
  155. C Récupération du numéro d'élément du maillage, du nombre de noeuds et d'éléments
  156. iele1 = IPT1.itypel
  157. nbnn1 = IPT1.num(/1)
  158. nbel1 = IPT1.num(/2)
  159. c récupération des caractéristique EF IPT1
  160. cbp : on essaie d'abord avec INFELE
  161. mele1 = IMODE1.IMAMOD
  162. MINTE1 = IMODE1.INFELE(11)
  163. cbp : si echec, on retente avec elquoi
  164. if(mele1.eq.0.or.MINTE1.eq.0) then
  165. call elquoi(nele1,0,3,IPTR1,IMODE1)
  166. INFO = IPTR1
  167. segact,INFO
  168. mele1 = INFELL(1)
  169. MINTE1 = INFELL(11)
  170. segdes,INFO
  171. endif
  172. cbp : si toujours pb, erreur...
  173. if(mele1.eq.0.or.MINTE1.eq.0) then
  174. write(ioimp,*) 'donnees relatives a l element fini inconnues'
  175. call erreur(591)
  176. endif
  177.  
  178.  
  179. c-----------------------------------------------------
  180. c RECUPERATION DU MAILLAGE INTERFACE
  181. c-----------------------------------------------------
  182. IPMAI2 = 0
  183. CALL LIROBJ('MMODEL ',IPMODL,0,IRETOU)
  184. IF(IERR.NE.0) RETURN
  185.  
  186. C -DANS LE CAS OU ON A UN MODELE EN ENTREE
  187. IF (IRETOU.EQ.1)THEN
  188. CALL ACTOBJ('MMODEL ',IPMODL,1)
  189. MMODE2=IPMODL
  190. segact,MMODE2
  191. N2 = MMODE2.KMODEL(/1)
  192. if(N2.gt.1) write(ioimp,*) 'attention 1 seule zone a ce jour!'
  193. IMODE2 = MMODE2.KMODEL(1)
  194. segact,IMODE2
  195. nele2 = IMODE2.NEFMOD
  196. IPT2 = IMODE2.IMAMOD
  197. SEGACT IPT2
  198. c pour l'instant on dit que nele = iele (marche pour iele entre 2 et 26, voir bdata.eso)
  199. iele2 = IPT2.itypel
  200. nbnn2 = IPT2.num(/1)
  201. nbel2 = IPT2.num(/2)
  202.  
  203. c recuperation des caracteritiques de l'element
  204. cbp : on essaie d'abord avec INFELE
  205. MINTE2 = IMODE2.INFELE(11)
  206. ngau2 = IMODE2.INFELE(6)
  207. cbp : si echec, on retente avec elquoi
  208. if(MINTE2.eq.0.or.ngau2.eq.0) then
  209. c write(ioimp,*) 'appel elquoi',iele2,nele2
  210. call elquoi(nele2,0,3,IPTR2,IMODE2)
  211. INFO = IPTR2
  212. segact,INFO
  213. MINTE2 = INFELL(11)
  214. ngau2 = INFELL(6)
  215. segdes,INFO
  216. endif
  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. CALL ACTOBJ('MAILLAGE',IPMAI2,1)
  228. IF(IERR.NE.0) RETURN
  229. IPT2 = IPMAI2
  230. SEGACT IPT2
  231. c pour l'instant on dit que nele = iele (marche pour iele entre 2 et 26, voir bdata.eso)
  232. iele2 = IPT2.itypel
  233. nele2 = iele2
  234. if (nele2.lt.2.or.nele2.gt.26) then
  235. write(ioimp,*)'element geometrique different de l element fini'
  236. call erreur(16)
  237. endif
  238. nbnn2 = IPT2.num(/1)
  239. nbel2 = IPT2.num(/2)
  240. ngau2 = -3
  241. c SEG2
  242. if (nele2.EQ.2) ngau2 = 2
  243. c SEG3
  244. if (nele2.EQ.3) ngau2 = 3
  245. c TRI3
  246. if (nele2.EQ.4) ngau2 = 1
  247. c TRI6
  248. if (nele2.EQ.6) ngau2 = 4
  249. c QUA4
  250. if (nele2.EQ.8) ngau2 = 4
  251. c QUA8
  252. if (nele2.EQ.10) ngau2 = 9
  253. if(ngau2.eq.-3) then
  254. write(ioimp,*)'nombre de points d integration inconnu pour',
  255. & 'ce type d element geometrique'
  256. call erreur(16)
  257. endif
  258. ENDIF
  259.  
  260. call RESHPT(ngau2,nbnn2,iele2,nele2,0,IPTR2,IRET)
  261. MINTE2 = IPTR2
  262. segact,MINTE1,MINTE2
  263.  
  264.  
  265. c-----------------------------------------------------
  266. c RECHERCHE DU MCHAML ISSU MCHEX1 D ENRICHISSEMENT
  267. c-----------------------------------------------------
  268. MCHAM1=0
  269. NBENR2=0
  270. segact,IMODE1
  271. NOBMOD = IMODE1.IVAMOD(/1)
  272. IF (NOBMOD.NE.0) THEN
  273. DO 1002 iobmo1=1,NOBMOD
  274. if((IMODE1.TYMODE(iobmo1)).ne.'MCHAML') goto 1002
  275. MCHEX1 = IMODE1.IVAMOD(iobmo1)
  276. segact,MCHEX1
  277. if((MCHEX1.TITCHE).ne.'ENRICHIS') goto 1003
  278. MCHAM1 = MCHEX1.ICHAML(1)
  279. segact,MCHAM1
  280. NBENR2 = MCHAM1.IELVAL(/1)
  281. do ienr2=1,NBENR2
  282. MELVA1=MCHAM1.IELVAL(IENR2)
  283. if(MELVA1.ne.0) segact,MELVA1
  284. enddo
  285. 1003 continue
  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. RX05= SQRT(ABS(PSIX))
  664. c write(ioimp,*) 'RX05 = ', RX05
  665. do II=1,nbnn2
  666. RINT2(II,JJ1) = RINT2(II,JJ1) + RX05*REL1(II,JJ)
  667. enddo
  668. endif
  669. ENDIF
  670.  
  671.  
  672. c B2X
  673. IF(nbenrj.ge.3) THEN
  674. MELVA1 = MCHAM1.IELVAL(3)
  675. MLREEL = MELVA1.IELCHE(JJ,iem1)
  676. if(MLREEL.ne.0) then
  677. PSIX = 0.d0
  678. SEGACT,MLREEL
  679. do iii0 = 1,nbnn1
  680. PSIX = PSIX + (SHPP1(1,iii0) * PROG(iii0))
  681. enddo
  682. RX05= SQRT(ABS(PSIX))
  683. do II=1,nbnn2
  684. RINT3(II,JJ1) = RINT3(II,JJ1) + RX05*REL1(II,JJ)
  685. enddo
  686. endif
  687. ENDIF
  688.  
  689. 4205 CONTINUE
  690.  
  691. if(iimpi.ge.4) then
  692. DO II=1,nbnn2
  693. write(ioimp,*) 'RINT1(',II,',:)=',(RINT1(II,iou),iou=1,nbnno2)
  694. ENDDO
  695. endif
  696.  
  697.  
  698. c on a trouvé le iem1 élément et on a fait le travail :
  699. c on passe au point de gauss suivant
  700. goto 1132
  701.  
  702. 1131 CONTINUE
  703. c---------- fin de la Boucle sur les elements de structure ----------------------
  704. if(iimpi.ge.2) write(ioimp,*)
  705. 1 'Fin de la boucle sur les elements de structure'
  706. if(izok.eq.0) write(ioimp,*) 'Attention le pt de Gauss '
  707. & ,igau2,'de l element ',iem2,' est hors support'
  708.  
  709.  
  710. 1132 CONTINUE
  711. c======= fin de la Boucle sur les pt de G de fissure ============================
  712.  
  713. c --- Remplissage du MELEME et DU XMATRI ---
  714. DO ityp1=1,2
  715.  
  716. c on commence par le MELEME des ddls Obligatoire (U)
  717. izone = 1
  718. c on poursuit avec le MELEME des ddls Sauts (H ou HB1 ou HB1B2)
  719. if(ityp1.eq.2) izone=nbenrj+1
  720. ity = ITAMYT(nbnno2,izone)
  721.  
  722. do iidim=1,idim
  723. ityty = ((ity-1) * idim) + iidim
  724. if(iimpi.ge.4) write(ioimp,*) ity,nbnno2,izone,' -> ',ityty
  725.  
  726. c --- Remplissage du MELEME ---
  727. MELEME = IRIGEL(1,ityty)
  728. NBNN = NUM(/1)
  729. NBELEM = NUM(/2) + 1
  730. if(iimpi.ge.4) then
  731. write(ioimp,*) 'on remplit le type de matrice ',ityty,
  732. 1 ' NBNN,NBELEM=',NBNN,NBELEM
  733. endif
  734. segadj,MELEME
  735. c element de fissure
  736. do ino2=1,nbnn2
  737. NUM(ino2,NBELEM) = IPT2.NUM(ino2,iem2)
  738. c NUM(ino2+nbnn2,NBELEM) = IPT2.NUM(ino2,iem2)
  739. enddo
  740. c element de structure
  741. do iino2=1,nbnno2
  742. c NUM((2*nbnn2)+iino2,NBELEM) = PRELEM(iino2)
  743. NUM(nbnn2+iino2,NBELEM) = PRELEM(iino2)
  744. enddo
  745.  
  746. c --- Remplissage du XMATRI.RE ---
  747. XMATRI = IRIGEL(4,ityty)
  748. NLIGRD = RE(/1)
  749. NLIGRP = RE(/2)
  750. NELRIG = RE(/3) + 1
  751. segadj,XMATRI
  752.  
  753. DO II=1,nbnn2
  754.  
  755. c -Ktt
  756. c on choisit de ne garder que la partie diagonale
  757. JJ=II
  758. XMATRI.RE(II,JJ,NBELEM) = RINT(II,JJ)*XISTAB
  759.  
  760. c -Kwt
  761. DO JJ=1,nbnn2
  762. XMATRI.RE(II,JJ+nbnn2,NBELEM) = RINT(II,JJ)
  763. XMATRI.RE(JJ+nbnn2,II,NBELEM) = RINT(II,JJ)
  764. c write(ioimp,*) 'XMATRI.RE(',II,JJ+nbnn2,NBELEM,') =', RINT(II,JJ)
  765. ENDDO
  766.  
  767. C -Kww
  768. c on choisit de ne garder que la partie diagonale
  769. JJ=II
  770. XMATRI.RE(nbnn2+II,nbnn2+JJ,NBELEM) = RINT(II,JJ)*RLATIN
  771.  
  772. C -Kut (on a deja fait le travail de positionnement avant)
  773. if(ityp1.eq.1) then
  774. c UX
  775. JJdec = 2*nbnn2
  776. DO JJ=1,nbnno2
  777. XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT0(II,JJ)
  778. XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT0(II,JJ)
  779. c write (6,* ) 'XMATRI.RE(',II,JJdec+JJ,NBELEM,') =', XMA
  780. c 1TRI.RE(II,JJdec+JJ,NBELEM)
  781. ENDDO
  782. else
  783. c AX
  784. JJdec = 2*nbnn2
  785. DO JJ=1,nbnno2
  786. XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT1(II,JJ)
  787. XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT1(II,JJ)
  788. c write (6,* ) 'XMATRI.RE(',II,JJdec+JJ,NBELEM,') =', XMA
  789. c 1TRI.RE(II,JJdec+JJ,NBELEM)
  790. ENDDO
  791.  
  792. c B1X
  793. IF(nbenrj.ge.2) THEN
  794. JJdec=JJdec+nbnno2
  795. DO JJ=1,nbnno2
  796. XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT2(II,JJ)
  797. XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT2(II,JJ)
  798.  
  799. c Cas ou les inconnues aux noeuds ne sont pas présentes
  800. c mis à 0 de ces inconnues en mettant un 1 sur la diagonale
  801.  
  802. c Début de la boucle pour tester les noeuds de structure d'ancien point de gauss
  803. cbp DO 5002 iigau2 = 1,igau2
  804. DO 5002 iigau2 = 1,ngau2
  805. iiem1 = elemPG(iigau2)
  806. DO JJ1=1,nbnn1
  807. iino1=IPT1.NUM(JJ1,iiem1)
  808. if(PRELEM(JJ).eq.iino1) then
  809. MELVA1 = MCHAM1.IELVAL(2)
  810. MLREEL = MELVA1.IELCHE(JJ1,iiem1)
  811. if(MLREEL.ne.0) goto 5002
  812. XMATRI.RE(JJdec+JJ,JJdec+JJ,NBELEM) = 1
  813. c write(ioimp,*) 'XMATRI.RE(',JJd ec+JJ,JJdec+JJ,NBEL
  814. c 1EM,') =', XMATRI.RE(JJdec+JJ,JJdec+JJ,NBELEM)
  815. endif
  816. ENDDO
  817. 5002 CONTINUE
  818. ENDDO
  819.  
  820.  
  821.  
  822. ENDIF
  823. c B2X
  824. IF(nbenrj.ge.3) THEN
  825. JJdec=JJdec+nbnno2
  826. DO JJ=1,nbnno2
  827. XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT3(II,JJ)
  828. XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT3(II,JJ)
  829. ENDDO
  830. ENDIF
  831. endif
  832.  
  833. ENDDO
  834.  
  835. enddo
  836.  
  837.  
  838. ENDDO
  839.  
  840.  
  841.  
  842. 1100 CONTINUE
  843. c==== fin de la Boucle sur les elements de fissure ==============================
  844. c----------------------------------------------------------------------
  845.  
  846.  
  847. c------------------------------------
  848. c AJUSTEMENT AVANT DE QUITTER
  849. c------------------------------------
  850.  
  851. if(iimpi.ge.2) write(ioimp,*) 'AJUSTEMENT AVANT DE QUITTER'
  852.  
  853.  
  854. C BOUCLE SUR LES SOUS RIGIDITÉS
  855. ityok = 0
  856. DO 2000 ityty=1,(idim*NTYMAT)
  857.  
  858. MELEME = IRIGEL(1,ityty)
  859. DESCR = IRIGEL(3,ityty)
  860. XMATRI = IRIGEL(4,ityty)
  861. NBELEM = NUM(/2)
  862. if(NBELEM.ne.0) then
  863. ityok = ityok + 1
  864. IRIGEL(1,ityok) = MELEME
  865. IRIGEL(3,ityok) = DESCR
  866. IRIGEL(4,ityok) = XMATRI
  867. segdes,XMATRI
  868. else
  869. segsup,MELEME,DESCR,XMATRI
  870. endif
  871. 2000 continue
  872. NRIGEL = ityok
  873. segadj,MRIGID
  874.  
  875. c-------------------------------
  876. c MENAGE AVANT DE QUITTER
  877. c-------------------------------
  878.  
  879. segsup,MTRAV
  880.  
  881. c-------------------------------
  882. c ON REMPLACE LES FX PAR DES LX
  883. c-------------------------------
  884. IF(INOMU.eq.0) THEN
  885. JGN=4
  886. JGM=2*idim
  887. segini,MLMOT1
  888. do iidim=1,idim
  889. MLMOT1.MOTS(iidim) = DUAOBL(iidim)
  890. MLMOT1.MOTS(iidim+idim) = DUAFAC(iidim)
  891. enddo
  892. ILMOT1 = MLMOT1
  893. IPRIG1 = MRIGID
  894. if(iimpi.ge.2) write(ioimp,*)' APPEL A FX2LX '
  895. CALL FX2LX(IPRIG1,ILMOT1,IPRIG2)
  896. MRIGID = IPRIG2
  897. SEGSUP,MLMOT1
  898. ENDIF
  899.  
  900. c write(ioimp,*) 'desactivation du MRIGID',MRIGID
  901. SEGDES,MRIGID
  902.  
  903. if(iimpi.ge.3) write(ioimp,*) 'ecriture du MRIGID',MRIGID
  904. CALL ECROBJ('RIGIDITE',MRIGID)
  905.  
  906. END
  907.  
  908.  
  909.  

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