Télécharger rigix.eso

Retour à la liste

Numérotation des lignes :

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

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