Télécharger rigix.eso

Retour à la liste

Numérotation des lignes :

  1. C RIGIX SOURCE BP208322 18/04/12 21:15:10 9802
  2. subroutine RIGIX (ivamat,ivacar,NMATT,CMATE,
  3. & IMODEL,IPINF,LOCIRI,NBGMAT,NELMAT,IMAT)
  4. C*********************************************************
  5. c
  6. C PROCEDURE UTILISEE DANS LE CAS D'ELEMENTS XFEM
  7. c POUR LE CALCUL DE RIGIDITE ELEMENTAIRES
  8. C
  9. C*********************************************************
  10. C Liste des modifications :
  11. C -as 2009/09/03 : ajout de IMAT en entrée de RIGIX pour le traitement
  12. C d'une matrice de Hooke en entrée (cas IMAT=2)
  13. C -YT 2010/02/11 : Developpement XFEM 3D
  14. C -BP 2011/01 : modification des boucles pour eviter trop de segadj
  15. C -BP 2013/10 : ajout de DIM3 pour les contraintes planes
  16. C -GG 2017/09 : Modification de mise à 0 les DDL non enrichis
  17. C
  18. C*********************************************************
  19. C PARTIE DECLARATIVE
  20. C*********************************************************
  21.  
  22. C
  23. IMPLICIT INTEGER(I-N)
  24. IMPLICIT REAL*8(A-H,O-Z)
  25. C
  26. CHARACTER*8 CMATE
  27. CYT PARAMETER (NDDLMAX=20,NBNIMAX=10)
  28. PARAMETER (NDDLMAX=30,NBNIMAX=10)
  29. CHARACTER*4 MOTINC(NDDLMAX),MOTDUA(NDDLMAX)
  30. DATA MOTINC/'UX ','UY ','UZ ','AX ','AY ','AZ ',
  31. >'B1X ','B1Y ','B1Z ','C1X ','C1Y ','C1Z ','D1X ','D1Y ',
  32. >'D1Z ','E1X ','E1Y ','E1Z ','B2X ','B2Y ','B2Z ','C2X ',
  33. >'C2Y ','C2Z ','D2X ','D2Y ','D2Z ','E2X ','E2Y ','E2Z '/
  34. DATA MOTDUA/'FX ','FY ', 'FZ ','FAX ','FAY ','FAZ ',
  35. >'FB1X','FB1Y','FB1Z','FC1X','FC1Y','FC1Z','FD1X','FD1Y',
  36. >'FD1Z','FE1X','FE1Y','FE1Z','FB2X','FB2Y','FB2Z','FC2X',
  37. >'FC2Y','FC2Z','FD2X','FD2Y','FD2Z','FE2X','FE2Y','FE2Z'/
  38. C
  39. CYT DATA MOTINC/'UX ','UY ','AX ','AY ',
  40. CYT >'B1X ','B1Y ','C1X ','C1Y ','D1X ','D1Y ','E1X ','E1Y ',
  41. CYT >'B2X ','B2Y ','C2X ','C2Y ','D2X ','D2Y ','E2X ','E2Y '/
  42. CYT DATA MOTDUA/ 'FX ','FY ','FAX ','FAY ',
  43. CYT >'FB1X','FB1Y','FC1X','FC1Y','FD1X','FD1Y','FE1X','FE1Y',
  44. CYT >'FB2X','FB2Y','FC2X','FC2Y','FD2X','FD2Y','FE2X','FE2Y'/
  45.  
  46. PARAMETER (NBENRMAX=30)
  47. C NBENRMAX deja defini dans rigixr.eso
  48. INTEGER LOCIRI(10,(1+NBENRMAX))
  49. DIMENSION MLRE(NBENRMAX+1)
  50. C
  51. -INC CCOPTIO
  52. -INC SMCOORD
  53. -INC SMMODEL
  54. -INC SMCHAML
  55. -INC SMINTE
  56. -INC SMRIGID
  57. -INC SMELEME
  58. -INC SMLREEL
  59. -INC CCHAMP
  60. C
  61. POINTEUR MCHEX1.MCHELM,MINT1.MINTE,MINT2.MINTE
  62. C
  63. C Segment (type LISTENTI) contenant les informations sur un element
  64. SEGMENT INFO
  65. INTEGER INFELL(JG)
  66. ENDSEGMENT
  67.  
  68. SEGMENT WRK1
  69. REAL*8 DDHOOK(LHOOK,LHOOK),XE(3,NBBB)
  70. REAL*8 REL(LRE,LRE),RINT(LRE,LRE)
  71. ENDSEGMENT
  72. C
  73. SEGMENT WRK2
  74. REAL*8 SHPWRK(6,NBNO),BGENE(LHOOK,LRE)
  75. c REAL*8 LV7WRK(NBENRMA2,2,6,NBNO)
  76. REAL*8 LV7WRK(NBENRMA2,2,6,NBBB)
  77. ENDSEGMENT
  78. C
  79. SEGMENT MVELCH
  80. REAL*8 VALMAT(NV1)
  81. ENDSEGMENT
  82. C
  83. SEGMENT MPTVAL
  84. INTEGER IPOS(NS),NSOF(NS)
  85. INTEGER IVAL(NCOSOU)
  86. CHARACTER*16 TYVAL(NCOSOU)
  87. ENDSEGMENT
  88. C
  89. SEGMENT MRACC
  90. INTEGER TLREEL(NBENRMA2,NBI)
  91. INTEGER MELRIG(NBELEM)
  92. ENDSEGMENT
  93.  
  94. C
  95. C
  96. C*********************************************************
  97. c
  98. C RECUPERATION + ACTIVATIONS + VALEURS UTILES
  99. c
  100. C*********************************************************
  101. C sont deja actifs dans rigi1.eso :
  102. C SEGACT MMODEL,IMODEL,MELVAL
  103. c IVAMAT
  104. MPTVAL=IVAMAT
  105.  
  106.  
  107. C++++ Recup + Activation de la geometrie ++++++++++++++++
  108. MELEME= IMAMOD
  109. SEGACT MELEME
  110.  
  111. C++++ RECUP DES INFOS EF ++++++++++++++++++++++++++++++++
  112. c + OBTENUES PAR ELQUOI DANS RIGI1 PENDANT PHASE 1
  113. C segment INFO deja actif dans RIGI1
  114. c + rigi1 n appelle pas elquoi, c'est modeli qui l'a fait
  115. c mais du coup, on na pas de segment minte
  116. c (car depend de si pt de g pour rigi, pour sigma....)
  117. c c'est + simple de rappeler elquoi ici
  118. MELE = NEFMOD
  119. call elquoi(MELE,0,3,IPINF,IMODEL)
  120. INFO = IPINF
  121. c MELE = INFELL(1)
  122. c NBPGA2= INFELL(2)
  123. c NBPGAU= INFELL(3)
  124. c NBPGAU= INFELL(4)
  125. c ICARA = INFELL(5)
  126. NGAU1 = INFELL(6)
  127. c LW = INFELL(7)
  128. LRE = INFELL(9)
  129. LHOOK = INFELL(10)
  130. MINT1 = INFELL(11)
  131. segact,MINT1
  132. MINT2 = INFELL(12)
  133. c if(MINT2.ne.0) segact,MINT2
  134. MFR = INFELL(13)
  135. IELE = INFELL(14)
  136. NDDL = INFELL(15)
  137. NSTRS = INFELL(16)
  138. c write(ioimp,*)'-> RIGIX infell',(infell(iou),iou=1,16)
  139.  
  140. c + AUTRES INFOS
  141. C nbre de noeuds par element
  142. NBNN1 = NUM(/1)
  143. C nbre d elements
  144. NBEL1 = NUM(/2)
  145.  
  146. c REM: pour se passer du dimensionnement du nbre d'enrichissement dans
  147. c elquoi et le realiser localement , on pourrait ecrire:
  148. c LRE = NDDLMAX*NBNN1
  149. c NDDL= NDDLMAX
  150.  
  151. C sous decoupage et points de Gauss de l element geometrique de base
  152. CYT if(MELE.eq.263) then
  153. cbp if(MELE.eq.263.or.MELE.eq.264) then
  154. if (MINT2.ne.0) then
  155. segact,MINT2
  156. NGAU2 = MINT2.POIGAU(/1)
  157. else
  158. NGAU2=0
  159. endif
  160. c write(*,*) 'dim de MINT2=6,',(MINT2.SHPTOT(/2)),(MINT2.SHPTOT(/3))
  161. c write(*,*) 'MINT2',(MINT2.QSIGAU(iou),iou=1,NGAU2)
  162. c segdes,MINT2
  163.  
  164.  
  165. c nbre maxi de fonction de forme par noeud (fonction std comprise)
  166. c NBNI = NDDL/IDIM inutile!
  167.  
  168.  
  169. C++++ Recup des infos d enrichissement +++++++++++++++++++
  170. c recup du MCHEX1 d'enrichissement
  171. NBENR1=0
  172. MCHAM1=0
  173. NOBMO1=IVAMOD(/1)
  174. if (NOBMO1.ne.0) then
  175. do iobmo1=1,NOBMO1
  176. if ((TYMODE(iobmo1)).eq.'MCHAML') then
  177. MCHEX1 = IVAMOD(iobmo1)
  178. segact,MCHEX1
  179. if ((MCHEX1.TITCHE).eq.'ENRICHIS') then
  180. MCHAM1 = MCHEX1.ICHAML(1)
  181. segact,MCHAM1
  182. segdes,MCHEX1
  183. goto 1000
  184. endif
  185. segdes,MCHEX1
  186. endif
  187. enddo
  188. write(ioimp,*) 'Le modele est vide (absence d enrichissement)'
  189. * return
  190. else
  191. write(ioimp,*) 'Aucun MCHAML enrichissement dans le Modele'
  192. * return
  193. endif
  194.  
  195. 1000 continue
  196. c niveau d enrichissement(s) du modele (ddls std u exclus)
  197. c NBENR1= 0 si std, 1 si H seul, 2 si H+F1, 3 si H+F1+F2, etc...
  198. if (MCHAM1.ne.0) NBENR1=MCHAM1.IELVAL(/1)
  199. c* write(ioimp,*) 'niveau d enrichissement(s) du modele',NBENR1
  200.  
  201. C*********************************************************
  202. c
  203. c PREPARATION DES OBJETS RESULTATS
  204. c
  205. C*********************************************************
  206. C
  207. C++++ RECHERCHE DES NOMS D'INCONNUES ET DES DUAUX
  208. c MOTINC et MODUA en dur pour l instant
  209.  
  210. C++++ REMPLISSAGE DE DESCR de taille maxi
  211. c (on ajustera en enlevant si besoin)
  212. NLIGRP = LRE
  213. NLIGRD = LRE
  214. SEGINI DESCR
  215. IPDSCR = DESCR
  216. C
  217. KINC = 0
  218. C----->> BOUCLE SUR LES fonctions de Forme
  219. DO IENR=1,NBNIMAX
  220. C------->> BOUCLE SUR LES NOEUDS
  221. DO I=1,NBNN1
  222. C--------->> BOUCLE SUR LA DIMENSION
  223. DO KDIM=1,IDIM
  224. KINC = KINC + 1
  225. CYT LISINC(KINC) = MOTINC((IDIM*(IENR-1))+KDIM)
  226. CYT LISDUA(KINC) = MOTDUA((IDIM*(IENR-1))+KDIM)
  227. LISINC(KINC) = MOTINC((3*(IENR-1))+KDIM)
  228. LISDUA(KINC) = MOTDUA((3*(IENR-1))+KDIM)
  229. NOELEP(KINC) = I
  230. NOELED(KINC) = I
  231. ENDDO
  232. ENDDO
  233. ENDDO
  234.  
  235. IF (KINC.NE.LRE) THEN
  236. WRITE(ioimp,*) 'IL Y A UNE ERREUR DANS DESCR'
  237. WRITE(ioimp,*) 'KINC , LRE = ',KINC,LRE
  238. RETURN
  239. ENDIF
  240. C
  241. SEGDES DESCR
  242. C
  243. C
  244. C++++ INITIALISATIONS...
  245. C
  246. c ... des tables LOCIRI et MLRE
  247. c MLRE contient le nombre d'inconnues de chaque sous-zone
  248. c determinee depuis le nombre de fonctions de forme
  249. c ienr= 1: U+H(1+1=2), 2: U+H+F1(2+4=6), 3: U+H+F1+F2(6+4=10)
  250. if (NBENR1.ne.0) then
  251. do ienr=1,NBENR1
  252. nbniJ = 2 + ((ienr-1)*4)
  253. MLRE(1+ienr) = IDIM*NBNN1*nbniJ
  254. c -LOCIRI
  255. LOCIRI(5,1+ienr)= NIFOUR
  256. enddo
  257. endif
  258. C Tables + longues car 1er indice correspond au fontion de forme std
  259. MLRE(1) = IDIM*NBNN1*1
  260. LOCIRI(5,1)= NIFOUR
  261. c on complete avec des 0
  262. if (NBENR1.lt.(NBENRMAX+1)) then
  263. do ienr=(NBENR1+1),(NBENRMAX)
  264. MLRE(1+ienr) = 0
  265. enddo
  266. endif
  267. c
  268. c ...DU SEGMENT WRK1
  269. NBENRMA2 = NBENRMAX
  270. NBBB = NBNN1
  271. SEGINI,WRK1
  272.  
  273. c ...DU SEGMENT WRK2
  274. c NBNO = NBNI
  275. NBNO = LRE/IDIM
  276. SEGINI,WRK2
  277.  
  278. c ...DU SEGMENT MVELCH
  279. NV1 = NMATT
  280. SEGINI,MVELCH
  281. C
  282. c ...DU SEGMENT MRACC
  283. NBENRMA2 = NBENRMAX
  284. NBI = NBNN1
  285. NBELEM = NUM(/2)
  286. segini,MRACC
  287. C
  288. C
  289.  
  290. C*********************************************************
  291. C
  292. C>>>>>>>>>>>>>>>>>>>>>>>>>>> 1ere BOUCLE SUR LES ELEMENTS
  293. C
  294. NBENR = 0
  295. DO 2000 J=1,NBEL1
  296. c write(ioimp,*) '===boucle1=== EF',J,' /NBEL1 ================'
  297. C
  298. C
  299. C*********************************************************
  300. C POUR CHAQUE ELEMENT, ...
  301. C*********************************************************
  302. C
  303. C++++ NBENRJ = niveau d enrichissement de l element ++++
  304. C =0 si EF std =1 si U+H =2 si U+H+F1 =3 si U+H+F1+F2
  305. NBENRJ=0
  306. if (NBENR1.ne.0) then
  307. do IENR=1,NBENR1
  308. MELVA1 = MCHAM1.IELVAL(IENR)
  309. segact,MELVA1
  310. do I=1,NBNN1
  311. mlree1 = MELVA1.IELCHE(I,J)
  312. if(mlree1.ne.0) NBENRJ=max(NBENRJ,IENR)
  313. enddo
  314. segdes,MELVA1
  315. enddo
  316. endif
  317. NBENR=max(NBENRJ,NBENR)
  318. c
  319. c++++ Selon le niveau d enrichissement NBENRJ...
  320. C
  321. c++++ ...on ajuste la 3eme dimension de XMATRI.RE
  322. c (=nbre d element de ce niveau d enrichissement)
  323. C on la stocke temporairement dans LOCIRI(4, ) meme si c est pas propre
  324. LOCIRI(4,1+NBENRJ) = LOCIRI(4,1+NBENRJ) + 1
  325. c
  326. C++++ ...on copie cet element a la geo correspondante
  327. IPT1 = LOCIRI(1,1+NBENRJ)
  328. c si elle n existe pas, il faut la creer
  329. if (IPT1 .eq. 0) then
  330. NBNN = NBNN1
  331. NBELEM = 1
  332. NBSOUS = 0
  333. NBREF = 0
  334. segini,IPT1
  335. IPT1.ITYPEL = ITYPEL
  336. LOCIRI(1,1+NBENRJ)=IPT1
  337. c si elle existe, il faut l agrandir
  338. else
  339. NBNN = NBNN1
  340. NBELEM = (IPT1.NUM(/2)) + 1
  341. NBSOUS = 0
  342. NBREF = 0
  343. segadj,IPT1
  344. endif
  345.  
  346. c Ajout de la connectivite de l'element courant
  347. do I=1,NBNN1
  348. IPT1.NUM(I,NBELEM) = NUM(I,J)
  349. enddo
  350.  
  351. c on retient la position NBELEM de l element J
  352. MELRIG(J)= NBELEM
  353.  
  354. C++++ ...on crée le descr qui va bien si il n existe pas
  355. DES1 = LOCIRI(3,1+NBENRJ)
  356. if (DES1.eq.0) then
  357. segini,DES1=DESCR
  358. NLIGRP = MLRE(1+NBENRJ)
  359. NLIGRD = MLRE(1+NBENRJ)
  360. segadj,DES1
  361. LOCIRI(3,1+NBENRJ) = DES1
  362. segdes,DES1
  363. endif
  364.  
  365. 2000 CONTINUE
  366. C FIN DE 1ere BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<<
  367. C
  368. C*********************************************************
  369.  
  370.  
  371. C*********************************************************
  372. C POUR CHAQUE NIVEAU D ENRICHISSEMENT, ...
  373. C*********************************************************
  374. do IENR=1,(NBENR1+1)
  375. c attention : IENR = IENRvrai + 1
  376. NELRIG = LOCIRI(4,IENR)
  377. C++++ ...ON CRÉE XMATRI DE LA BONNE TAILLE TOUT DE SUITE
  378. if (NELRIG.ne.0) then
  379. NLIGRP = MLRE(IENR)
  380. NLIGRD = MLRE(IENR)
  381. segini,XMATRI
  382. LOCIRI(4,IENR)=XMATRI
  383. C++++ ...ON DESACTIVE IPT1
  384. IPT1 = LOCIRI(1,IENR)
  385. segdes,IPT1
  386. endif
  387. enddo
  388.  
  389.  
  390.  
  391. C*********************************************************
  392. C
  393. C>>>>>>>>>>>>>>>>>>>>>>>>>>> 2eme BOUCLE SUR LES ELEMENTS
  394. C
  395. DO 2001 J=1,NBEL1
  396. c write(*,*) '===boucle 2=== EF',J,' /NBEL1 ================'
  397. C
  398. C*********************************************************
  399. C POUR CHAQUE ELEMENT, ...
  400. C*********************************************************
  401. C
  402. C++++ ON RECUPERE LES COORDONNEES DES NOEUDS DE L ELEMENT
  403. CALL DOXE(XCOOR,IDIM,NBNN1,NUM,J,XE)
  404.  
  405. C++++ NBENRJ = niveau d enrichissement de l element ++++
  406. C =0 si EF std =1 si U+H =2 si U+H+F1 =3 si U+H+F1+F2
  407. NBENRJ=0
  408. if (NBENR1.ne.0) then
  409. do IENR=1,(NBENR1)
  410. MELVA1 = MCHAM1.IELVAL(IENR)
  411. segact,MELVA1
  412. do I=1,NBNN1
  413. mlree1 = MELVA1.IELCHE(I,J)
  414. C on en profite pour remplir MRACC table de raccourcis pour cet element
  415. TLREEL(IENR,I) = mlree1
  416. if (mlree1.ne.0) then
  417. NBENRJ = max(NBENRJ,IENR)
  418. C et on active la listreel
  419. segact,mlree1
  420. endif
  421. enddo
  422. segdes,MELVA1
  423. enddo
  424. endif
  425. if (NBENRMA2.gt.NBENR1) then
  426. do IENR2=(NBENR1+1),NBENRMA2
  427. do I=1,NBNN1
  428. TLREEL(IENR2,I) = 0
  429. enddo
  430. enddo
  431. endif
  432. C
  433. c++++ on recupere XMATRI
  434. XMATRI = LOCIRI(4,1+NBENRJ)
  435.  
  436. c++++ quelques indices pour dimensionner au plus juste
  437. c nbre total de ddl de l'élément considéré
  438. NLIGRD = MLRE(1+NBENRJ)
  439. NLIGRP = MLRE(1+NBENRJ)
  440. c nbre fonction de forme=((Ni_std+Ni_enrichi)*nbnoeud)=Ninconnues/idim
  441. NBNI = (MLRE(1+NBENRJ)) / IDIM
  442. c position de cet element pour cet enrichissement
  443. IELRIG = MELRIG(J)
  444. c write(*,*) 'EF',J,' NBENRJ=',NBENRJ,
  445. c &'donc',NLIGRD,' ddls et ',NBNI,' fonctions de forme'
  446. c
  447.  
  448. C++++ ON MET A ZERO LA SOMME QUE L ON VA CALCULER
  449. c CALL ZERO(RINT,LRE,NLIGRP)
  450. * CALL ZERO(RINT,NLIGRD,NLIGRP)
  451. c |-> impossible car ZERO considere qu il sagit de la dimension de SHPWRK
  452. CALL ZEROX(RINT,NLIGRD,NLIGRP,LRE)
  453. c
  454. c REM : il est interessant au niveau du tems cpu d'utiliser moins de pts
  455. c de Gauss lorsque l element n est pas enrichi car cela n'apporte
  456. c rien de plus.
  457. c On utilise MINT2 = infell(12) qui contient le segment d'integration
  458. C de l'EF std (QUA4 par ex. pour XQ4R)
  459. if ((NBENRJ.eq.0).and.(MINT2.ne.0)) then
  460. MINTE = MINT2
  461. NBPGAU= NGAU2
  462. else
  463. MINTE = MINT1
  464. NBPGAU= NGAU1
  465. endif
  466.  
  467. C*********************************************************
  468. C
  469. C>>>>>>>>>> BOUCLE SUR LES POINTS DE GAUSS
  470. C
  471. DO 2500 KGAU=1,NBPGAU
  472. c
  473. c write(ioimp,*) '-pt de G ',KGAU,' / ',NBPGAU,'----'
  474. C
  475. C*********************************************************
  476. C Initialisation à 0
  477. C*********************************************************
  478.  
  479. c CALL ZERO(SHPWRK,6,NBNO)
  480. CALL ZERO(SHPWRK,6,NBNI)
  481. C
  482. i6zz = 3
  483. IF (IDIM.EQ.3) i6zz = 4
  484. C
  485. c do ienr7=1,NBENRMAX
  486. do ienr7=1,NBENRJ
  487. do inod7=1,NBNN1
  488. c do i6=1,6
  489. CYT do i6=1,3
  490. do i6=1,i6zz
  491. LV7WRK(ienr7,1,i6,inod7) = 0.D0
  492. LV7WRK(ienr7,2,i6,inod7) = 0.D0
  493. enddo
  494. enddo
  495. enddo
  496.  
  497.  
  498. c*********************************************************
  499. c Calcul des fonction de forme std dans repere local
  500. c*********************************************************
  501.  
  502. ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccccc
  503. c (et donc sur les Ni std)
  504. DO 2510 I=1,NBNN1
  505.  
  506. C++++ Calcul des Ni std
  507. c (rappel: 1:Ni 2:Ni,qsi 3:Ni,eta avec i=1,4)
  508. SHPWRK(1,I) = SHPTOT(1,I,KGAU)
  509. SHPWRK(2,I) = SHPTOT(2,I,KGAU)
  510. SHPWRK(3,I) = SHPTOT(3,I,KGAU)
  511. IF (IDIM.EQ.3) SHPWRK(4,I) = SHPTOT(4,I,KGAU)
  512.  
  513. 2510 CONTINUE
  514. ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  515.  
  516.  
  517.  
  518. c*********************************************************
  519. c Passage des fonctions de forme std dans repere global
  520. c*********************************************************
  521.  
  522. C++++ CALCUL DES Ni,x Ni,y (i=1,4) + CALCUL DE det(J)
  523. CALL JACOBI(XE,SHPWRK,IDIM,NBNN1,DJAC)
  524. c if(J.eq.1.and.KGAU.eq.1)
  525. c &write(*,*) 'Ni(i=1,4)=',(SHPWRK(1,iou),iou=1,NBNN1)
  526.  
  527. c*********************************************************
  528. c Si on est pas du tout enrichi on peut sauter une grosse partie
  529. c*********************************************************
  530. if(NBENRJ.eq.0) goto 2999
  531.  
  532. c*********************************************************
  533. c Calcul des level set + leurs derivees dans repere global
  534. c*********************************************************
  535.  
  536. ccccc BOUCLE SUR LES enrichissements ccccccccccccccccccc
  537. do 2520 ienr=1,NBENRJ
  538.  
  539. c MELVA1=MCHAM1.IELVAL(IENR)
  540. c segact,MELVA1
  541.  
  542. ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccc
  543. DO 2521 I=1,NBNN1
  544.  
  545. C++++ Le I eme noeud est-il ienr-enrichi?
  546. mlree1= TLREEL(IENR,I)
  547. if(mlree1.eq.0) goto 2521
  548.  
  549.  
  550. c Calcul du repere local de fissure(=PSI,PHI)
  551. c (rappel: 1,1:psi 1,2:phi 2,1 psi,x ...etc...)
  552. do 2522 inode=1,NBNN1
  553. c pour le H-enrichissement, on n a pas gardé PSI (inutile)
  554. if (ienr.ne.1) then
  555. c valeur de PSI au inode^ieme noeud
  556. xpsi1 = mlree1.PROG(inode)
  557. c qu on multiplie par la valeur de Ni^std au pt de G considéré
  558. LV7WRK(ienr,1,1,I)= LV7WRK(ienr,1,1,I)
  559. & + (SHPWRK(1,inode)*xpsi1)
  560. LV7WRK(ienr,1,2,I)= LV7WRK(ienr,1,2,I)
  561. & + (SHPWRK(2,inode)*xpsi1)
  562. LV7WRK(ienr,1,3,I)= LV7WRK(ienr,1,3,I)
  563. & + (SHPWRK(3,inode)*xpsi1)
  564. IF (IDIM.EQ.3) LV7WRK(ienr,1,4,I)= LV7WRK(ienr,1,4,I)
  565. & + (SHPWRK(4,inode)*xpsi1)
  566. c valeur de PHI au inode^ieme noeud
  567. xphi1 = mlree1.PROG(NBNN1+inode)
  568. else
  569. xphi1 = mlree1.PROG(inode)
  570. endif
  571. LV7WRK(ienr,2,1,I)= LV7WRK(ienr,2,1,I)
  572. & + (SHPWRK(1,inode)*xphi1)
  573. LV7WRK(ienr,2,2,I)= LV7WRK(ienr,2,2,I)
  574. & + (SHPWRK(2,inode)*xphi1)
  575. LV7WRK(ienr,2,3,I)= LV7WRK(ienr,2,3,I)
  576. & + (SHPWRK(3,inode)*xphi1)
  577. IF (IDIM.EQ.3) LV7WRK(ienr,2,4,I)= LV7WRK(ienr,2,4,I)
  578. & + (SHPWRK(4,inode)*xphi1)
  579. 2522 continue
  580.  
  581. 2521 continue
  582. ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  583.  
  584.  
  585. 2520 CONTINUE
  586. ccccc fin de BOUCLE SUR LES enrichissements cccccccccccccccc
  587.  
  588. c on a construit
  589. C LV7WRK(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD)
  590.  
  591.  
  592. c*********************************************************
  593. c Ajout des fonctions de forme d enrichissement
  594. c + leurs derivees dans repere global
  595. c*********************************************************
  596. CALL SHAPX(MELE,LV7WRK,NBENRMAX,NBENRJ,TLREEL,SHPWRK,IRET)
  597.  
  598. c retour a la partie commune aux EF enrichi et non enrichi
  599. 2999 continue
  600.  
  601. C*********************************************************
  602. c write(*,*) 'C CALCUL DE B'
  603. C*********************************************************
  604. call ZERO(BGENE,LHOOK,LRE)
  605. KB=1
  606. C boucle sur tous les Ni
  607. DO 3001 II=1,NBNI
  608.  
  609. BGENE(1,KB) = SHPWRK(2,II)
  610. BGENE(2,KB+1) = SHPWRK(3,II)
  611. BGENE(4,KB) = SHPWRK(3,II)
  612. BGENE(4,KB+1) = SHPWRK(2,II)
  613.  
  614. IF (IDIM.EQ.3) THEN
  615. BGENE(3,KB+2)=SHPWRK(4,II)
  616. BGENE(5,KB) =SHPWRK(4,II)
  617. BGENE(5,KB+2)=SHPWRK(2,II)
  618. BGENE(6,KB+1)=SHPWRK(4,II)
  619. BGENE(6,KB+2)=SHPWRK(3,II)
  620. ENDIF
  621. KB = KB + IDIM
  622. 3001 CONTINUE
  623.  
  624.  
  625. C*********************************************************
  626. c write(*,*) 'C CALCUL DE D (=MATRICE DE HOOKE)'
  627. C*********************************************************
  628. c
  629. c+++++cas DOHOOO
  630. IF (IMAT.EQ.2) THEN
  631. MELVAL=IVAL(1)
  632. IBMN=MIN(J ,IELCHE(/2))
  633. IGMN=MIN(KGAU,IELCHE(/1))
  634. MLREEL=IELCHE(IGMN,IBMN)
  635. SEGACT MLREEL
  636. IF (KGAU.LE.NBGMAT.AND.(J.LE.NELMAT.OR.NBGMAT.GT.1))
  637. 1 CALL DOHOOO(PROG,LHOOK,DDHOOK)
  638. SEGDES MLREEL
  639.  
  640. c+++++cas DOHMAS
  641. ELSE IF (IMAT.EQ.1) THEN
  642. DO 4001 IM=1,NMATT
  643. IF (IVAL(IM).NE.0) THEN
  644. MELVAL=IVAL(IM)
  645. IBMN=MIN(J ,VELCHE(/2))
  646. IGMN=MIN(KGAU,VELCHE(/1))
  647. if (ibmn.gt.0.and.igmn.gt.0) then
  648. VALMAT(IM)=VELCHE(IGMN,IBMN)
  649. else
  650. VALMAT(IM)=0.D0
  651. endif
  652. ELSE
  653. VALMAT(IM)=0.D0
  654. ENDIF
  655. 4001 CONTINUE
  656. C
  657. CALL DOHMAS(VALMAT,CMATE,IFOUR,LHOOK,1,DDHOOK,IRET)
  658. IF (IRET.EQ.0) THEN
  659. MOTERR(1:8)=CMATE
  660. MOTERR(9:16)=NOMFR(MFR/2+1)
  661. INTERR(1)=IFOUR
  662. CALL ERREUR(81)
  663. RETURN
  664. ENDIF
  665. ENDIF
  666. C
  667. C*********************************************************
  668. c write(*,*) 'C CALCUL DE Bt.D.B'
  669. C*********************************************************
  670.  
  671. C PRISE EN COMPTE DU POIDS D'INTEGRATION
  672. DJAC = ABS(DJAC) * POIGAU(KGAU)
  673.  
  674. c CAS DES CONTRAINTES PLANES
  675. C recuperation de l'epaisseur: DIM3 donnée facultative du materiau
  676. IF (IFOUR.EQ.-2) THEN
  677. MPTVAL=IVACAR
  678. IF (IVACAR.NE.0) THEN
  679. MELVAL=IVAL(1)
  680. IF (MELVAL.NE.0) THEN
  681. IGMN=MIN(KGAU,VELCHE(/1))
  682. IBMN=MIN(J,VELCHE(/2))
  683. DIM3=VELCHE(IGMN,IBMN)
  684. ELSE
  685. DIM3=1.D0
  686. ENDIF
  687. ENDIF
  688. DJAC=DJAC*DIM3
  689. MPTVAL=IVAMAT
  690. ENDIF
  691.  
  692. CALL ZEROX(REL,NLIGRD,NLIGRP,LRE)
  693.  
  694. CALL BDBST(BGENE,DJAC,DDHOOK,LRE,LHOOK,REL)
  695. * CALL BDBST(BGENE,DJAC,DDHOOK,NLIGRP,LHOOK,REL)
  696. c |-> impossible car BDBST considere qu il sagit de la dimension de bgene
  697. c pour gagner du temps cpu il faudrait qqchose du type:
  698. c CALL BDBSTX(BGENE,DJAC,DDHOOK,NLIGRP,LHOOK,LRE,REL)
  699.  
  700. C*********************************************************
  701. C CUMUL DANS RINT(II,JJ)
  702. C*********************************************************
  703. DO 4201 II=1,NLIGRD
  704. DO 4201 JJ=1,NLIGRP
  705. RINT(II,JJ) = RINT(II,JJ) + REL(II,JJ)
  706. 4201 CONTINUE
  707. c
  708. c
  709. 2500 CONTINUE
  710. C FIN DE BOUCLE SUR LES POINTS DE GAUSS <<<<<<<<<<<<<<
  711. C
  712. C*********************************************************
  713. C write(*,*) 'C STOCKAGE DANS XMATRI de RE(II,JJ)'
  714. c et de XMATRI dans IMATRI
  715. C*********************************************************
  716. C
  717. DO 6001 II=1,NLIGRD
  718.  
  719. DO 6002 JJ=1,NLIGRP
  720.  
  721. c XMATRI.RE(II,JJ,nelrig) = RINT(II,JJ)
  722. c RE(II,JJ,NELRIG) = RINT(II,JJ)
  723. RE(II,JJ,IELRIG) = RINT(II,JJ)
  724.  
  725. c si NON enrichi, on met tout a 0 sauf la diago
  726. c if(mlree1.eq.0) RE(II,JJ) = 0. inutile
  727.  
  728. 6002 CONTINUE
  729.  
  730.  
  731. 6001 CONTINUE
  732. c
  733. c quelques desactivations
  734. do IENR=1,NBENR1
  735. do I=1,NBNN1
  736. mlree1 = TLREEL(IENR,I)
  737. if(mlree1.ne.0) segdes,mlree1
  738. enddo
  739. enddo
  740. c
  741. c
  742. 2001 CONTINUE
  743. C FIN DE BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<<
  744.  
  745.  
  746. C
  747. SEGSUP,WRK1,WRK2,MVELCH
  748. SEGSUP,MRACC
  749.  
  750. segdes,MELEME
  751. segdes,MINT1
  752. if(MINT2.ne.0) segdes,MINT2
  753. c
  754. c desactivation des maillages et segment imatri
  755. do ienr=1,(1+NBENR1)
  756. c IPT1 = LOCIRI(1,ienr)
  757. c if(IPT1.ne.0) segdes,IPT1
  758. XMATRI = LOCIRI(4,ienr)
  759. if(XMATRI.ne.0) segdes,XMATRI
  760. enddo
  761.  
  762. c SEGDES dans RIGI1 = IMODEL
  763. C
  764. RETURN
  765. END
  766.  
  767.  
  768.  
  769.  
  770.  
  771.  

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